PROGRAM Compute_Transform_Arnaud
! A sample code for computing a 2d Fourier transform of the Sunyaev-Zel'dovich effect 
! profile, Tsz(x) [x=theta/theta500], from the generalized NFW model of Arnaud et al., 
! A&A, 517, 92 (2010). The trasnformation is done via
!
! trasnform=(fourpi*r500)/l500^2 
!          *int_{xin}^{xout}[x^2*Tsz(x)*sin(l*x/l500)/(l*x/l500)]dx
!
! Output file "x_transform.txt" contains
! 1st: x
! 2nd: 2d Fourier transform of Tsz(x) in units of micro Kelvin [Rayleigh-Jeans], 
!      where x=l/l500
!
! May 15, 2010: E.Komatsu
  IMPLICIT none
  double precision :: xin=1d-6,xout=6d0 ! pressure profile is cut at r=xout*r500
  double precision :: om0=0.277d0,h0=0.702d0,tol=1d-6
  double precision :: m500,r500,theta500,l500,rhoc,Ez,da,z
  double precision :: y,transform2d,tsz
  double precision :: chi,pgnfw,rombint
  integer :: i
  external chi,pgnfw,rombint
! Sample parameters: Coma cluster
  z=0.0231d0
  m500=6.61d14 ! h^-1 Msun, computed from M500-Tx relation of Vikhlinin et al with Tx=8.4keV
  da=rombint(chi,1d0/(1d0+z),1d0,tol,om0)
  da=da/(1d0+z) ! proper distance, h^-1 Mpc
  print*,'redshift=',z
  print*,'Da(physical)=',da,' h^-1 Mpc'
  Ez=dsqrt(om0*(1d0+z)**3d0+1d0-om0)
  rhoc=2.775d11*Ez**2d0 ! critical density in units of h^2 Msun Mpc^-3
  r500=(3d0*m500/4d0/3.14159d0/500d0/rhoc)**(1d0/3d0) ! h^-1 Mpc
  theta500=r500/da ! radians
  l500=da/r500
  print*,'M500=',M500,' h^-1 Msun'
  print*,'r500=',r500,' h^-1 Mpc'
  print*,'theta500=',theta500*10800d0/3.14159d0,' arcmin'
  print*,'l500=',l500
  print*,'pressure is cut at r/r500=',xout
! calculate profiles
  open(1,file='x_transform.txt')
  do i=1,1000
     y=xin + 1d-3*dble(i)*(xout-xin) ! y=l/l500
     ! See, e.g., Appendix D1 of Komatsu et al., arXiv:1001.4538
     !            Eq.(2) of Komatsu&Seljak (2002)
     transform2d=1.65d0*(h0/0.7d0)**2d0*Ez**(8d0/3d0)*(m500/3d14/0.7d0)**(2d0/3d0+0.12d0)&
           *rombint(pgnfw,xin,xout,tol,y)*(r500/h0)/0.5176d0*(4d0*3.14159d0)/l500**2d0
     Tsz=-2d0*283d0*(transform2d/50d0) ! uK
     write(1,'(2E18.5)')y,tsz
  enddo
  close(1)
END PROGRAM Compute_Transform_Arnaud
!-----------------------------------------------------------------------------
DOUBLE PRECISION FUNCTION chi(x,om0)
! chi(x) = 1/[H(x)*x^2] in units of h^-1 Mpc
  IMPLICIT none
  DOUBLE PRECISION, intent(IN) :: x,om0
  chi = 2998d0/dsqrt(om0*x+(1d0-om0)*x**4d0)
  return
END FUNCTION chi
!-----------------------------------------------------------------------------
DOUBLE PRECISION FUNCTION pgnfw(x,ell) ! x=r/r500 & ell=l/l500
  IMPLICIT none
  double precision :: x,ell
  double precision :: c=1.177d0,g=0.3081d0,a=1.0510d0,b=5.4905d0,P0=8.403d0
  double precision :: h0=0.702d0
  pgnfw=P0*(0.7d0/h0)**1.5d0/(c*x)**g/(1d0+(c*x)**a)**((b-g)/a) &
       *x**2d0*dsin(ell*x)/(ell*x)
  return
END FUNCTION pgnfw