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
e_ph0 = 4.202e3  # 1.2 ev Temperature, median wavelength
skip = 1 # half PopIII stars become binaries
mass_min = 10
lifetime_max = 3.0e7 # 30 Myr
mass_max = mass_min*np.exp(lifetime_max*3.15e7/T_ed)  # accrete at Eddington limit for max lifetime
emission_max = epsilon*(mass_max-mass_min)*mass_solar*(3.0e10)**2*6.24e11

f=open('PopIII_list_0002_type1.txt','r')
f1=open('PopIII_list_0002_type5.txt','r')

fo = open('PopIII_xray-sources_lifetime.txt','w')

source_number = 0

line = f.readline()
line = f.readline()
while line !='':
  starmass = float(line.split()[1])
  creationtime = float(line.split()[2])
  starposition = float(line.split()[3]), float(line.split()[4]),float(line.split()[5])
  lifetime = T_ed*np.log((starmass)/mass_min)
  if lifetime > lifetime_max*3.15e7: lifetime = lifetime_max*3.15e7
  if starmass > mass_min and (currenttime - creationtime) < lifetime/tunit:
     energyphoton = e_ph0*(10.0/(mass_min+mass_max)*2)**0.25
     Luminosity = epsilon*((starmass if starmass < mass_max else mass_max)-mass_min)*mass_solar*(3.0e10)**2*6.24e11/lifetime/energyphoton
     lifetime = lifetime/tunit 
     fo.write(str("PhotonTestSourceType[%d]       = 1"%source_number+"\n"))
     fo.write(str("PhotonTestSourcePosition[%d]   = %f %f %f"%(source_number, \
        starposition[0],starposition[1],starposition[2])+"\n"))
     fo.write(str("PhotonTestSourceLuminosity[%d] = %e"%(source_number,Luminosity)+"\n"))
     fo.write(str("PhotonTestSourceLifeTime[%d]   = %f"%(source_number,lifetime)+"\n"))
     fo.write(str("PhotonTestSourceCreationTime[%d] = %f"%(source_number,creationtime)+"\n"))
     fo.write(str("PhotonTestSourceEnergyBins[%d] = 1"%source_number+"\n")) 
     fo.write(str("PhotonTestSourceEnergy[%d] = %f"%(source_number,energyphoton)+"\n\n"))     
     source_number = source_number + 1
  for i in xrange(skip):
     line=f.readline()
  line = f.readline()

line = f1.readline()
line = f1.readline()
while line !='':
  starmass = float(line.split()[1])
  creationtime = float(line.split()[2])
  starposition = float(line.split()[3]), float(line.split()[4]),float(line.split()[5])
  if starmass > 1e20: starmass = starmass/1.0e20 
  if starmass > mass_min:
     energyphoton = e_ph0*(10.0/(mass_min+mass_max)*2)**0.25
     lifetime = T_ed*np.log((starmass)/mass_min)
     if lifetime > lifetime_max*3.15e7 : lifetime = lifetime_max*3.15e7
     Luminosity = epsilon*((starmass if starmass < mass_max else mass_max)-mass_min)*mass_solar*(3.0e10)**2*6.24e11/lifetime/energyphoton
     lifetime = lifetime/tunit 
     fo.write(str("PhotonTestSourceType[%d]       = 1"%source_number+"\n"))
     fo.write(str("PhotonTestSourcePosition[%d]   = %f %f %f"%(source_number, \
        starposition[0],starposition[1],starposition[2])+"\n"))
     fo.write(str("PhotonTestSourceLuminosity[%d] = %e"%(source_number,Luminosity)+"\n"))
     fo.write(str("PhotonTestSourceLifeTime[%d]   = %f"%(source_number,lifetime)+"\n"))
     fo.write(str("PhotonTestSourceCreationTime[%d] = %f"%(source_number,creationtime)+"\n"))
     fo.write(str("PhotonTestSourceEnergyBins[%d] = 1"%source_number+"\n"))
     fo.write(str("PhotonTestSourceEnergy[%d] = %f"%(source_number,energyphoton)+"\n\n"))
     source_number = source_number + 1
  for i in xrange(skip):
     line=f1.readline()
  line = f1.readline()

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