PROGRAM szfast
  USE cosmo
  USE linearpk
  USE sigma
  USE growth
  USE angular_distance
  USE multipole
! This code calculates the power spectrum of the Sunyaev-Zel'dovich effect.
! The methodology is described in Komatsu & Seljak, 336, 1256 (2002).
! The output, "multipole_szpower.txt", contains:
! 1st: multipole
! 2nd: l(l+1)Cl/twopi [uK^2] in the Rayleigh-Jeans limit [i.e., dT/T=-2y]
! December 14, 2010: E.Komatsu
  IMPLICIT none
  double precision :: z1=1d-5,z2=5d0 ! redshift limits
  character(len=128) :: filename
  integer :: npk,i,j,n=1001 ! 201 is good enough for l>100
  double precision :: cl,dcldlnx,lnx,dlnx,lnx1,lnx2
  external dcldlnx
! Specify three cosmological parameters
! The data type has been defined in MODULE cosmo.
  obh2=0.02262d0
  om0=0.277d0
  ode0=0.723d0
  h0=0.702d0
  w=-1d0
! read in and tabulate P(k)
  filename='wmap5baosn_max_likelihood_matterpower_at_z=30.dat'
  npk=896 ! # of lines in the file
  CALL open_linearpk(filename,npk)
! fit sigma^2(R) to Chebyshev polynomials
  CALL compute_sigma2
  CALL close_linearpk
! compute the growth factor
  CALL setup_growth
! compute and tabulate da(z) 
  CALL setup_da
! compute the SZ power spectrum, l(l+1)Cl/twopi [uK^2] 
! in the RJ limit, dT/T=-2y
  open(1,file='multipole_szpower.txt')
  lnx1=dlog(1d0+z1)
  lnx2=dlog(1d0+z2)
  dlnx=(lnx2-lnx1)/dble(n-1)
  do j=1,16
     ell=10d0**(1d0+dble(j-1)*0.2)
     ! integrate by the trapezoidal rule
     cl=0d0
     cl=dcldlnx(lnx1)*(0.5d0*dlnx)
     cl=cl+dcldlnx(lnx2)*(0.5d0*dlnx)
     do i=2,n-1
        lnx=lnx1+dble(i-1)*dlnx
        cl=cl+dcldlnx(lnx)*dlnx
     enddo
     print*,ell,cl*ell*(ell+1d0)/2d0/3.14159d0
     write(1,'(2E15.5)')ell,cl*ell*(ell+1d0)/2d0/3.14159d0
  enddo
  close(1)
END PROGRAM szfast
!--------------------------------------------------------------
double precision function dcldlnx(lnx)
  ! dcldlnx = (1+z) (dV/dz) int dlnM dn/dlnM Tsz^2
  ! lnx=ln(1+z)
  IMPLICIT none
  double precision, intent(IN) :: lnx
  double precision :: M1=5d11,M2=5d15 ! mass limits, h^-1 Msun
  double precision :: integrand,lnM1,lnM2,lnM,z
  integer :: i
  external integrand
  z=exp(lnx)-1d0
  lnM1=dlog(M1) ! minimum mass, h^-1 Msun
  lnM2=dlog(M2) ! maximum mass, h^-1 Msun
  CALL qgaus1(integrand,lnM1,lnM2,dcldlnx,z) ! Gaussian quadrature
  dcldlnx=dcldlnx*(1d0+z)
  return
end function dcldlnx