import numpy as np
import commands

tunit=6.88027e+14/3.156e7/1e6

for n, rs in zip([10,12,14,21,31,41],[20,19,18,17,16,15]):

 fi=open('Redshift%02i/rarepeak_galaxies_properties_%04i'%(rs,n),'r')
 lines1 = fi.readlines()
 fi.close()

 fi=open('Redshift%02i/SED/Luminosity.out'%(rs),'r')
 lines2 = fi.readlines()
 fi.close()

 magnitude = []
 for line in lines2:
     magnitude.append([int(line.split()[0]),float(line.split()[2])])
 magnitude=np.asarray(magnitude)

 fi=open('Redshift%02i/SED/Luminosity_M1600.out'%(rs),'r')
 lines3 = fi.readlines()
 fi.close()

 M1600 = []
 for line in lines3:
     M1600.append([int(line.split()[0]),float(line.split()[2])])
 M1600=np.asarray(M1600)

 fo=open('Redshift%02i/rarepeak_galaxies_properties_luminosity_%04i'%(rs,n),'w')

 a=commands.getoutput('grep CosmologyCurrentRedshift '+'Redshift%02i/RedshiftOutput%04i'%(rs,n))
 redshift=float(a.split()[2])

 a=commands.getoutput('grep InitialTime '+'Redshift%02i/RedshiftOutput%04i'%(rs,n))
 currenttime=float(a.split()[2])*tunit

 fo.write("# Rarepeak, Wise & Abel LW background, z=%4.1f, time=%10.6f Myr \n"%(redshift,currenttime))
 fo.write("%6s %12s %12s %12s %16s %16s %16s %16s %16s %20s %12s %12s %12s %20s %16s %16s %12s %12s \n" % 
             ("#","x", "y", "z", "M_vir [Msun]", "M_gas [Msun]", "M_star [Msun]",
              "Metallicity (VW)", "Metallicity (MW)","Stellar Metallicity", "t_25 [Myr]", "t_50 [Myr]","t_75 [Myr]",  
              'M_hm_stellar [Msun]','M_hm_dm [Msun]','M_hm_gas [M_sun]', 'M_bol','M_1600'))

 
 for line1 in lines1[2:]:
    result = np.zeros([2])
    halo_id = int(line1.split()[0])
    
    if (magnitude[:,0] == halo_id).sum() > 0:
      result[0] = magnitude[magnitude[:,0] == halo_id][0,1]
      result[1] = M1600[M1600[:,0] == halo_id][0,1]
    line1=line1.rstrip('\r\n') 
    fo.write(line1+" %12.8f %12.8f\n" %(result[0],result[1]) )
    
 fo.close()