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

f=open('PopIII_list_0221_type1.txt','r')
f1=open('PopIII_list_0221_type5.txt','r')

fo = open('PopIII_evolution.txt','w')
starlist = [[],[]]
line = f.readline()
line = f.readline()
while line !='':
  starmass = float(line.split()[1])
  creationtime = float(line.split()[2])
  #if starmass > 40:
  starlist[0].append(starmass) 
  starlist[1].append(creationtime*tyr)  
  line = f.readline()

line = f1.readline()
line = f1.readline()
while line !='':
  starmass = float(line.split()[1])
  creationtime = float(line.split()[2])
  if starmass > 1e20: starmass= starmass/1e20
  starlist[0].append(starmass)
  starlist[1].append(creationtime*tyr)
  line = f1.readline()

timemin = min(starlist[1])
timemax = max(starlist[1])
masscounts = np.zeros((timemax-timemin)/1e6+1)
for starmass,creationtime in zip(starlist[0],starlist[1]):
    masscounts[int((creationtime-timemin)/1e6)] = masscounts[int((creationtime-timemin)/1e6)]+starmass

bin = timemin
cumulativemass = 0.0
for starmass in masscounts:
  cumulativemass = cumulativemass + starmass
  fo.write(str(bin)+" "+str(starmass)+" "+str(cumulativemass)+" "+str(cumulativemass/baryonmass)+"\n")
  bin = bin + 1e6 


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