PROGRAM Compute_Profile_Arnaud
! A sample code for computing a GAS pressure profile, Pgas(x) [x=r/r500], and
! the corresponding SZ profile, Tsz(x) [x=theta/theta500], from the generalized 
! NFW model of Arnaud et al., A&A, 517, 92 (2010)
!
! Output file "x_pgas_tsz.txt" contains
! 1st: x
! 2nd: Pgas(x) in units of eV/cm^3, where x=r/r500
! 3rd: Tsz(x) in units of micro Kelvin [Rayleigh-Jeans], where x=theta/theta500
!
! May 15, 2010: E.Komatsu
  IMPLICIT none
  double precision :: xin=0d0,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,rhoc,Ez,da,z
  double precision :: y,pgas3d,pgas2d,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
  print*,'M500=',M500,' h^-1 Msun'
  print*,'r500=',r500,' h^-1 Mpc'
  print*,'theta500=',theta500*10800d0/3.14159d0,' arcmin'
  print*,'pressure is cut at r/r500=',xout
! calculate profiles
  open(1,file='x_pgas_tsz.txt')
  do i=1,1000
     y=xin + 1d-3*dble(i)*(xout-xin) ! y=theta/theta500 or r/r500
     ! See, e.g., Appendix D1 of Komatsu et al., arXiv:1001.4538
     pgas3d=1.65d0*(h0/0.7d0)**2d0*Ez**(8d0/3d0)*(m500/3d14/0.7d0)**(2d0/3d0+0.12d0)   &
           *pgnfw(0d0,y)/0.5176d0
     pgas2d=1.65d0*(h0/0.7d0)**2d0*Ez**(8d0/3d0)*(m500/3d14/0.7d0)**(2d0/3d0+0.12d0)   &
           *rombint(pgnfw,-dsqrt(xout**2d0-y**2d0),dsqrt(xout**2d0-y**2d0),tol,y)      &
           *(r500/h0)/0.5176d0
     Tsz=-2d0*283d0*(pgas2d/50d0) ! uK
     write(1,'(3E18.5)')y,pgas3d,tsz
  enddo
  close(1)
END PROGRAM Compute_Profile_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(l,y) ! x=r/r500
  IMPLICIT none
  double precision :: l,x,y
  double precision :: c=1.177d0,g=0.3081d0,a=1.0510d0,b=5.4905d0,P0=8.403d0
  double precision :: h0=0.702d0
  x=dsqrt(l**2d0+y**2d0)
  pgnfw=P0*(0.7d0/h0)**1.5d0/(c*x)**g/(1d0+(c*x)**a)**((b-g)/a)
  return
END FUNCTION pgnfw