import numpy as np

tunit = 6.88027e+14
tyr = tunit/3.15e7
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
e_ph0 = 4.202e3  # 1.2 ev Temperature, median wavelength
skip = 0 # half PopIII stars become binaries


f=open('PopIII_haloes_0242_type1.txt','r')
f1=open('PopIII_haloes_0242_type5.txt','r')

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

source_number = 0

line = f.readline()
line1 = f1.readline()
while line !='':
  halomass = float(line.split()[1])
  haloposition = float(line.split()[3]), float(line.split()[4]),float(line.split()[5].replace(":",""))
  energyphoton = []; lifetime = []; emission = []; creationtime = []
  line = f.readline();line = f.readline()
  while line !="\n":
     skip = skip + 1
     starmass = float(line.split()[1])
     if skip%2 == 0:
      if starmass > 10:
        creationtime.append(float(line.split()[2]))
        energyphoton.append(e_ph0*(10.0/(starmass+10)*2)**0.25)
        lifetime.append(T_ed*np.log((starmass)/10))
        emission.append(epsilon*(starmass-10.0)*mass_solar*(3.0e10)**2*6.24e11)      
     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 > 10:
        creationtime.append(float(line1.split()[2]))
        energyphoton.append(e_ph0*(10.0/(starmass+10)*2)**0.25)
        lifetime.append(T_ed*np.log((starmass)/10))
        emission.append(epsilon*(starmass-10.0)*mass_solar*(3.0e10)**2*6.24e11)
     line1=f1.readline() 
  energyphoton = np.array(energyphoton); lifetime = np.array(lifetime); 
  emission = np.array(emission) ; creationtime = np.array(creationtime)
  if energyphoton.size != 0:
     totalemission = emission.sum()
     halolifetime = lifetime.max()
     luminosity = totalemission/halolifetime
     halophotonenergy = (energyphoton*emission).sum()/emission.sum()
     luminosity = luminosity/halophotonenergy
     halolifetime = halolifetime/tunit     
     halocreationtime = creationtime.min()
     fo.write(str("PhotonTestSourceType[%d]       = 1"%source_number+"\n"))
     fo.write(str("PhotonTestSourcePosition[%d]   = %f %f %f"%(source_number, \
        haloposition[0],haloposition[1],haloposition[2])+"\n"))
     fo.write(str("PhotonTestSourceLuminosity[%d] = %e"%(source_number,luminosity)+"\n"))
     fo.write(str("PhotonTestSourceLifeTime[%d]   = %f"%(source_number,halolifetime)+"\n"))
     fo.write(str("PhotonTestSourceCreationTime[%d] = %f"%(source_number,halocreationtime)+"\n"))
     fo.write(str("PhotonTestSourceEnergyBins[%d] = 1"%source_number+"\n")) 
     fo.write(str("PhotonTestSourceEnergy[%d] = %f"%(source_number,halophotonenergy)+"\n\n"))     
     source_number = source_number + 1
  line = f.readline()
  if line == "\n": line = f.readline()
  line1 = f1.readline();
  if line1 == "\n": line1 = f1.readline()

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

