import numpy as np
from math import exp

XX = 0.76
Omega_B=0.0449
h=0.71
Omega0=0.266

redshiftlist = [15.0,15.4,15.99,16.6,16.99,17.45,17.91,18.50,19.0,19.5,20.0,21.0,22.0,23.0,23.5,24.0,25.0]
redshiftlist1 =  []
i = 0
for z in redshiftlist:
   if i == 0: 
      z1=z
      i+=1 
      continue    
   redshiftlist1.append((z+z1)/2)
   z1=z
   i+=1


#tauHI_coeff  = (XX/0.75d0)*(Omega_B/0.044d0)*(h/0.7d0)*(Omega0/0.27d0)**-0.5d0 &
#         * e0keV**-3.25d0

#tauHeI_coeff = ((1d0-XX)/0.25d0)*(Omega_B/0.044d0)*(h/0.7d0)*(Omega0/0.27d0)**-0.5d0 &
#         * e0keV**-3.22d0

for zobs1 in redshiftlist:
  print zobs1
  for zs1 in redshiftlist1:
    if zs1 <= zobs1: continue
    tauHI_coeff1  = (XX/0.75)*(Omega_B/0.044)*(h/0.7)*(Omega0/0.27)**-0.5*1.0**-3.25
    tauHeI_coeff1 = ((1.0-XX)/0.25)*(Omega_B/0.044)*(h/0.7)*(Omega0/0.27)**-0.5*1.0**-3.22
    tauHI_coeff3  = (XX/0.75)*(Omega_B/0.044)*(h/0.7)*(Omega0/0.27)**-0.5*3.0**-3.25
    tauHeI_coeff3 = ((1.0-XX)/0.25)*(Omega_B/0.044)*(h/0.7)*(Omega0/0.27)**-0.5*3.0**-3.22
    tauHI1  = 4.10863* tauHI_coeff1   * (zs1/(1.0+25))**1.5 * ((zs1/zobs1)**1.75 -1.0)
    tauHeI1 = 10.13573* tauHeI_coeff1 * (zs1/(1.0+25))**1.5 * ((zs1/zobs1)**1.72 -1.0)
    fmod1 = exp(-tauHI1-tauHeI1)
    tauHI3  = 4.10863* tauHI_coeff3   * (zs1/(1.0+25))**1.5 * ((zs1/zobs1)**1.75 -1.0)
    tauHeI3 = 10.13573* tauHeI_coeff3 * (zs1/(1.0+25))**1.5 * ((zs1/zobs1)**1.72 -1.0)
    fmod3 = exp(-tauHI3-tauHeI3)
    print zs1,fmod1, fmod3, fmod3/fmod1