** Temperature and Polarization Radiation Transfer Function **
September 21, 2008: E.Komatsu
Ref: Komatsu & Spergel, PRD 63, 063002 (2001)

THIS PROGRAM IS BASED UPON "CMBFAST" CODE (VERSION 4.0) DEVELOPED BY 
UROS SELJAK AND MATIAS ZALDARRIAGA. FOR COPYRIGHT INFORMATION, 
SEE THE FILE CALLED "LICENSE" IN THIS DIRECTORY.

Here we provide a program for computing the temperature, gT(l,k), and
polarization, gP(l,k), radiation transfer functions of Cosmic Microwave
Background Radiation anisotropies. Here, l is the multipole, while k
is the wavenumber. For definition of gT(l,k), see Eq.(10) of 
Komatsu & Spergel, PRD 63, 063002 (2001). The polarization transfer 
function, gP(l,k), is defined similarly. 

A few remarks for your information:
- "cmbflat.f" in the original CMBFAST package has been replaced by "gTflat.f". 
- Only flat universes and scalar perturbations are supported. 
  (I.e., no tensor modes.)

To understand what these functions are, it is useful to relate them to
the angular power spectra of temperature, C_l^TT, E-polarization, 
C_l^EE, and their cross-correlation, C_l^TE. The relations are:

C_l^TT = (2/pi) int k^2 dk P(k) [gT(l,k)]^2
C_l^EE = [(l-1)l(l+1)(l+2)] (2/pi) int k^2 dk P(k) [gP(l,k)]^2
C_l^TE = sqrt[(l-1)l(l+1)(l+2)] (2/pi) int k^2 dk P(k) gT(l,k)gP(l,k)

where P(k) is the power spectrum of Bardeen's curvature perturbations, 
Phi, during the matter dominated era, which is related to the primordial
curvature perturbation, zeta, as Phi = (3/5)zeta, at the linear order.

For comparison, we also provide the data for C_l's, 
"wmap5baosn_max_likelihood_scalCls.dat," which was
generated using CAMB code for the maximum likelihood parameters
given in Table I of Komatsu et al.(2008) [WMAP 5-year interpretation paper] 
with "WMAP5+BAO+SN". The input file for CAMB is also provided 
(wmap5baosn_max_likelihood_params.ini). 

NOTE: this program is based upon CMBFAST version 4.0, which is older than
version 4.5.1 for which the accuracy has been improved significantly.
(CMBFAST version 4.5.1 and CAMB yield consistent results; 
see Seljak et al., PRD 68, 083507 (2003).)
Therefore, C_l's computed from gT and gP of this program are slightly 
different from those given in "wmap5baosn_max_likelihood_scalCls.dat". 
Nevertheless they are close enough for the purpose of computing non-Gaussian 
statistics such as the bispectrum, Minkowski functions, wavelets, etc.

- To compile and use the  program, edit Makefile
and simply "./make"
- It will generate executable called "jlgen", "compute_gT", and "compute_cl".
- First, run "./jlgen". (You may also use "./sample_jlgen.cmd"), which will
generate a data file containing the spherical Bessel function with the filename
of your choice.
- Second, run "./compute_gT". (You may also use "./sample_gT.cmd"), which will
generate a data file containing l, k, gT, and gP, with the filename of your
choice. 
- Optionally you may run "./compute_cl" to compute the angular power spectra,
C_l^TT, C_l^EE, and C_l^TE, from the radiation transfer functions.

=========================================================================
A simple program to generate the angular power spectra, C_l^TT, C_l^EE, and
C_l^TE, from the radiation transfer functions, gT and gP, generated by
"compute_gT", is given as "compute_cl.f90":

PROGRAM Compute_Cl
  ! A sample program for computing the angular power spectrum of temperature,
  ! C_l^TT, E polarization, C_l^EE, as well as their cross-correlation, C_l^TE,
  ! from the radiation transfer functions. 
  ! September 21, 2008: E.Komatsu
  IMPLICIT none
  integer :: nk,lmax
  double precision :: dl,ak,dak,gT,gP,cl(3)
  double precision :: akpivot,deltaR2,ns,norm
  double precision :: dummy
  character(len=70) :: filename,header
  integer :: i,ik,l
  ! parameters regarding the power spectrum of curvature perturbations
  akpivot=0.002d0 ! pivot wavenumber for the primordial power spectrum
  deltaR2=2.46d-9 ! <Delta_R^2> at the pivot wavenumber
  ns=0.962d0      ! primordial tilt
  ! calculate the proper normalization factor
  norm=(2.726d6)**2d0 & ! conversion factor to uK^2
       *(3d0/5d0)**2d0& ! conversion from R^2 to (Phi during matter era)^2
       *(deltaR2*2d0*3.1415926535d0**2d0/akpivot**3d0)
                        ! conversion from Delta_R^2(k_pivot) tp P_R(k_pivot)
  ! read in temperature (gT) and polarization (gP) radiation transfer functions
  filename='wmap5baosn_max_likelihood_radiation_transfer_function.txt'
  nk=6200 ! # of k-sampling points
  lmax=1500
  open(1,file=filename,status='old')
  ! skip first 5 rows
  do i=1,5
     read(1,'(1A70)')header
     print*,header
  enddo
  ! integrate over wavenumbers and output the angular power spectra of
  ! cl(1)=TT, cl(2)=EE, cl(3)=TE
  open(2,file='cl.txt',status='unknown')
  do l=2,lmax
     cl(1:3)=0d0
     do ik=1,nk-1
        read(1,*)dl,ak,dak,gT,gP
        ! TT
        cl(1)=cl(1)+(2d0/3.1415926535d0)*gT**2d0   &
             *ak**2d0*dak*(ak/akpivot)**(ns-4d0)*norm
        ! EE
        cl(2)=cl(2)+(2d0/3.1415926535d0)*gP**2d0   &
             *(dl-1d0)*dl*(dl+1d0)*(dl+2d0)        & !  (l+2)!/(l-2)!
             *ak**2d0*dak*(ak/akpivot)**(ns-4d0)*norm 
        ! TE
        cl(3)=cl(3)+(2d0/3.1415926535d0)*gT*gP     &
             *dsqrt((dl-1d0)*dl*(dl+1d0)*(dl+2d0)) & !  sqrt[(l+2)!/(l-2)!]
             *ak**2d0*dak*(ak/akpivot)**(ns-4d0)*norm 
     enddo
     write(2,'(1I5,3E15.5)')l,cl(1:3)*dl*(dl+1d0)/2d0/3.1415926535d0
  enddo
  close(1)
  close(2)
END PROGRAM Compute_Cl