import numpy as np
import math

h=0.71
omegam=0.266
c=3.0e5

currentredshift = 16.0

f=open("LW_intensity1.txt","r")
line = f.readline()
line = f.readline()
line = f.readline()

redshift = []
LW_source = []
while line != "":
    redshift.append(float(line.split()[0]))
    LW_source.append(float(line.split()[4]))
    line = f.readline()   

def r_os(sourceredshift):
 r_os=2*c*(100.0*h)**-1*omegam**(-0.5)*(1+sourceredshift)**-0.5*(((1+currentredshift)/(1+sourceredshift))**-0.5-1)
 return r_os

def f_m(sourceredshift):
   alpha = (h/0.7)**-1*(omegam/0.27)**-0.5*((1+sourceredshift)/21.0)**-0.5
   r=r_os(sourceredshift)
   if r > 97.39*alpha:
     return 0.0
   else:
     return 1.7*math.exp(-(r/(116.29*alpha))**0.68)-0.7

for currentredshift in redshift[1:]:
 J_LW_now = 0
 z1 = redshift[0]; J_LW1 = LW_source[0]
 for z2,J_LW2 in zip(redshift[1:],LW_source[1:]):
    if z2 < currentredshift: break
    J_LW_now += ((1+currentredshift)/(1+(z2+z1)/2))**3 * (r_os(z1)-r_os(z2)) / (1+(z2+z1)/2.0) * J_LW1 * f_m((z2+z1)/2.0)
    #print (z2+z1)/2,J_LW_now,J_LW1,z1,z2
    z1=z2;J_LW1=J_LW2
 print currentredshift,J_LW_now*1.4212581669612953e-10