import numpy as np

tunit = 6.88027e+14
tyr = tunit/3.15e7
currenttime = 12.756587046328
omega_b = 0.0449
h = 0.71
baryonmass = 1.9e-29*h**2*omega_b*((8*3.08e24)**3*0.2745822715872526)/2e33
mass_solar = 1.998e33
L_ed = 1.3e38 * 6.24e11 # ev
epsilon = 0.2
T_ed = epsilon*mass_solar/1.44e17
L_ed = 1.26e38
e_ph0 = 4.202e3  # 1.2 ev Temperature, median wavelength
skip = 0 # half PopIII stars become binaries
lifetime_max = 10.0 #Myr
mass_min = 40
mass_min1 = 10

f=open('PopIII_haloes_0002_type1.txt','r')
f1=open('PopIII_haloes_0002_type5.txt','r')

f2=open('Overlaped_haloes_list.txt','r')
overlayed=[]
line=f2.readline()
while line !='':
   overlayed.append(int(line))
   line=f2.readline()


fo = open('PopIII_haloes_xray-sources_distribution_forever_40Ms.txt','w')

lifetime_max = lifetime_max*1e6
mass_max = mass_min*np.exp(lifetime_max*3.15e7/T_ed)
mass_max1 = mass_min1*np.exp(lifetime_max*3.15e7/T_ed)
emission_max = epsilon*(mass_max-mass_min)*mass_solar*(3.0e10)**2*6.24e11
emission_max1 = epsilon*(mass_max1-mass_min1)*mass_solar*(3.0e10)**2*6.24e11

source_number = 0

line = f.readline()
line1 = f1.readline()
j=0
totalluminosity=[]
while line !='':
  halomass = float(line.split()[1])
  haloposition = float(line.split()[3]), float(line.split()[4]),float(line.split()[5].replace(":",""))
  energyphoton = []; emission = []
  line = f.readline();line = f.readline()
  while line !="\n":
     skip = skip + 1
     starmass = float(line.split()[1])
     if skip%2 == 0:
      if starmass > mass_min:
        creationtime=float(line.split()[2])
        energyphoton.append(e_ph0*(10.0/(starmass+mass_min)*2)**0.25)
        mass=mass_min*np.exp((currenttime-creationtime)*tunit/T_ed)
        emission.append(L_ed*mass) 
      elif starmass > mass_min1:
        creationtime=float(line.split()[2])
        energyphoton.append(e_ph0*(10.0/(starmass+mass_min1)*2)**0.25)
        mass=mass_min1*np.exp((currenttime-creationtime)*tunit/T_ed)
        emission.append(L_ed*mass)    
     line=f.readline()
  line1 = f1.readline();line1 = f1.readline()
  while line1 !="\n":
     skip = skip + 1
     starmass = float(line1.split()[1])
     if skip%2 == 0:
      if starmass > 1e20: starmass = starmass/1.0e20
      if starmass > mass_min:
        creationtime=float(line1.split()[2])
        energyphoton.append(e_ph0*(10.0/(starmass+mass_min)*2)**0.25)
        mass=mass_min*np.exp((currenttime-creationtime)*tunit/T_ed)
        emission.append(L_ed*mass)
      elif starmass > mass_min1:
        creationtime=float(line1.split()[2])
        energyphoton.append(e_ph0*(10.0/(starmass+mass_min1)*2)**0.25)
        mass=mass_min1*np.exp((currenttime-creationtime)*tunit/T_ed)
        emission.append(L_ed*mass)
     line1=f1.readline() 
  energyphoton = np.array(energyphoton) 
  emission = np.array(emission) 
  if energyphoton.size != 0:
     emission = np.array(emission)
     totalemission = emission.sum()
     totalluminosity.append(totalemission)
     if j not in overlayed:
        fo.write(str(halomass)+" "+str(totalemission)+" "+str((emission>0).sum())+"\n")   
        source_number = source_number + 1
  j = j +1
  line = f.readline()
  if line == "\n": line = f.readline()
  line1 = f1.readline();
  if line1 == "\n": line1 = f1.readline()

totalluminosity=np.array(totalluminosity)
print totalluminosity.sum()
fo.write(str(totalluminosity.sum())+"\n")

f.close()
f1.close()
fo.close()