import math
import sys
from yt.mods import *
import numpy as np
from yt.analysis_modules.halo_finding.api import *
import glob
import re
import commands

#massmax = 1.2e9
massmin = 1.0e6
bins = 32  # 300
binwidth = 0.1 #math.log10(massmax/massmin)/bins
massmax = 10**(math.log10(massmin)+bins*binwidth)
eV_erg = 6.24e11

dlist = glob.glob('RD00??')
dlist.sort()
for folder in dlist:
  n = int(re.sub(r'\D',"",folder))
  if glob.glob(folder+'/xray_source_histogram_30Myr_40Ms.txt') != []: continue
  f=open(folder+'/PopIII_haloes_xray-sources_distribution_30Myr_40Ms.txt','r')
  fo = open(folder+'/xray_source_histogram_30Myr_40Ms.txt','w')

  sourcelist = []
  emissionlist = []
  sourcenumber = []

  line = f.readline()
  while line != '':
   halomass = float(line.split()[0])
   emission = float(line.split()[1])
   source = float(line.split()[2])
   sourcelist.append(halomass)
   emissionlist.append(emission/eV_erg)
   sourcenumber.append(source)
   line = f.readline()

  sourcecount = np.zeros(bins+1)
  totalemission = np.zeros(bins+1)
  totalsource = np.zeros(bins+1)
  totalhalo = np.zeros(bins+1)

  for halomass,emission,source in zip(sourcelist,emissionlist,sourcenumber):
   index = int((math.log10(halomass/massmin)+0.5*binwidth)/binwidth)
   if emission <= 0:
     totalhalo[index] = totalhalo[index] + 1
   else:   
     sourcecount[index] = sourcecount[index] + 1
     totalemission[index] = totalemission[index] + emission
     totalsource[index] = totalsource[index] +source 
     totalhalo[index] = totalhalo[index] + 1

  cumulation = 0.0
  bin = math.log10(massmin) + 0.5*binwidth
  for count,emission,source,halo in zip(sourcecount,totalemission,totalsource,totalhalo):
    cumulation = cumulation + emission
    fo.write(str(10**bin)+" "+str(count)+" "+str(emission)+" "+str(cumulation)+" "+str((source/count if count > 0 else 0))+" "+str(halo)+"\n")
    bin =  bin + binwidth

  f.close()
  fo.close()