import numpy as np
import glob
import re
from math import log10
import commands


def pop2_lum(mass, nion_mass):
    return mass * nion_mass

def pop3_lum(mass):
    _mass = max(min(mass, 500), 5)
    x = log10(_mass)
    x2 = x*x
    if _mass > 9 and _mass <= 500:
        Q = 10.0**(43.61 + 4.9*x - 0.83*x2)
    elif _mass >= 5 and _mass <= 9:
        Q = 10.0**(39.29 + 8.55*x)
    return  Q

year = 3.1557e7
timeunit = 6.88027e+14

fo=open("Ionising_photons.txt",'w')

dlist = glob.glob('Redshift*')
dlist.sort()
dlist.reverse()
for folder in dlist: 
  filename = glob.glob(folder+"/RedshiftOutput????")
  n = int(re.sub(r'\D',"",filename[0].split('/')[1]))
  #if n <=0: continue
  a=commands.getoutput('grep InitialTime '+folder+'/RedshiftOutput00??')
  currenttime=float(a.split()[2])
  a=commands.getoutput('grep StarClusterIonizingLuminosity '+folder+'/RedshiftOutput00??')
  nion_mass = float(a.split()[2])
  a=commands.getoutput('grep CosmologyCurrentRedshift '+folder+'/RedshiftOutput00??')
  redshift = float(a.split()[2])

  f=open(folder+"/PopIII_list_%04i_type5.txt"%n,"r")
  lines_PopIII5=f.readlines()
  lines_PopIII5.pop(0)
  f.close()

  f=open(folder+"/PopIII_list_%04i_type1.txt"%n,"r")
  lines_PopIII1=f.readlines()
  lines_PopIII1.pop(0)
  f.close()

  f=open(folder+"/PopII_list_%04i.txt"%n,"r")
  lines_PopII=f.readlines()
  lines_PopII.pop(0)
  f.close()

  Nphoton = 0

  for line in lines_PopIII5:
    mass,creationtime=line.split()[1:3]
    mass=float(mass);creationtime=float(creationtime)
    if mass > 1e20: mass=mass/1e20
    logm = np.log10(mass)
    lifetime = 10**(9.785 - 3.759*logm + 1.413*logm*logm - \
                              0.186*logm*logm*logm) 
    if lifetime > (currenttime - creationtime)*timeunit/year:
       lifetime = (currenttime - creationtime)*timeunit/year    
    Nphoton = Nphoton+pop3_lum(mass)*lifetime*year

  for line in lines_PopIII1:
    mass,creationtime=line.split()[1:3]
    mass=float(mass);creationtime=float(creationtime)
    logm = np.log10(mass)
    lifetime = 10**(9.785 - 3.759*logm + 1.413*logm*logm - \
                              0.186*logm*logm*logm)         
    if lifetime > (currenttime - creationtime)*timeunit/year:
       lifetime = (currenttime - creationtime)*timeunit/year
    Nphoton = Nphoton+pop3_lum(mass)*lifetime*year

  for line in lines_PopII:
    mass,creationtime=line.split()[1:3]
    mass=float(mass);creationtime=float(creationtime)
    lifetime=2e7
    if lifetime > (currenttime - creationtime)*timeunit/year:
       lifetime = (currenttime - creationtime)*timeunit/year
    Nphoton = Nphoton+pop2_lum(mass,nion_mass)*lifetime*year

  fo.write("%8.5f %5.2f %7.5e\n"%(currenttime,redshift,Nphoton))


fo.close()