PROGRAM Compute_Profile_Battaglia
  USE battaglia
! 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 Battaglia et al.
! profile [Battaglia et al., ApJ, 758, 75 (2012)]
!
! Output file "xvir_pgas_tsz.txt" contains
! 1st: x
! 2nd: Pgas(x) in units of eV/cm^3, where xvir=r/rvir
! 3rd: Tsz(x) in units of micro Kelvin [Rayleigh-Jeans], where xvir=theta/thetavir
!
! Output file "x500_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
!
! December 11, 2010: E.Komatsu
  IMPLICIT none
  double precision :: xin=0d0,xout=3d0 ! pressure profile is cut at r=xout*rvir
  double precision :: omega,om0=0.277d0,obh2=0.02262d0,fb,tol=1d-6
  double precision :: mvir,rvir,thetavir,rs,delc,c
  double precision :: m500,r500,theta500,rhoc,Ez,da,z
  double precision :: r200,P200
  double precision :: y,pgas3d,pgas2d,tsz
  double precision :: chi,ypgas,rombint,pi=3.1415926535d0
  integer :: i
  external chi,ygas,ypgas,rombint
! Sample parameters: Coma cluster
  h0=0.702d0
  z=0.0231d0
  mvir=1.58d15 ! h^-1 Msun, roughly corresponding to M500=6.61d14 h^-1 Msun 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'
! compute omega's, rhoc, and delc
  Ez=sqrt(om0*(1d0+z)**3d0+1d0-om0) ! E(z)=H(z)/H0
  omega=om0*(1d0+z)**3d0/Ez**2d0 ! Omega(z)
  rhoc=2.775d11*Ez**2d0 ! critical density in units of h^2 M_sun/Mpc^3
  ! Eq.(6) of Bryan & Norman (1998), valid only for a flat universe!
  delc=18d0*pi*pi+82d0*(omega-1d0)-39d0*(omega-1d0)**2d0
  ! One may also use Eq.(C19) of Nakamura & Suto (1997):
  ! delc=omega*18d0*pi*pi*(1d0+0.4093d0*(1d0/omega-1d0)**0.9052d0)
! compute virial and NFW parameters
  rvir=(3d0*mvir/4d0/pi/delc/rhoc)**(1d0/3d0) ! virial radius, h^-1 Mpc
  thetavir=rvir/da ! radians
  print*,'Mvir=',Mvir,' h^-1 Msun'
  print*,'rvir=',rvir,' h^-1 Mpc'
  print*,'thetavir=',thetavir*10800d0/pi,' arcmin'
  print*,'pressure is cut at r/rvir=',xout
  ! Concentration parameter from Duffy et al. (2008)
  c=7.85d0*(mvir/2d12)**(-0.081)/(1d0+z)**0.71 ! Duffy et al. (2008)
  ! One may also use the concentration parameter from Seljak (2000)
  !!$ c=10d0*(mvir/3.42d12)**(-0.2)/(1d0+z)
  rs=rvir/c ! NFW scale radius, h^-1 Mpc
! for convenience, convert the virial mass to M500
  CALL mvir2mdel(Mvir,rs,c,500d0*rhoc,m500)
  r500=(3d0*m500/4d0/pi/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/pi,' arcmin'
! calculate profiles
  CALL mvir2mdel(mvir,rs,c,200d0*rhoc,m200)
  r200=(3d0*m200/4d0/pi/200d0/rhoc)**(1d0/3d0) ! h^-1 Mpc
  fb=obh2/(om0*h0**2d0)
  P200=0.50439d0*200d0*(h0*Ez)**2d0*fb*(m200/1d15)/(2d0*r200) ! eV/cm^3
  zz=z
  open(1,file='xvir_pgas_tsz.txt')
  open(2,file='x500_pgas_tsz.txt')
  do i=1,1000
     y=xin + 1d-3*dble(i)*(xout-xin) ! y=theta/thetavir or r/rvir
     ! See, e.g., Appendix D2 of Komatsu et al., arXiv:1001.4538
     pgas3d=P200*ypgas(0d0,y*rvir/r200)
     pgas2d=P200*(r200/h0) &
     *rombint(ypgas,-dsqrt(xout**2d0-y**2d0)*rvir/r200,dsqrt(xout**2d0-y**2d0)*rvir/r200,tol,y*rvir/r200)
     Tsz=-2d0*283d0*(pgas2d/50d0) ! uK
     write(1,'(3E18.5)')y,pgas3d,tsz
     write(2,'(3E18.5)')y*rvir/r500,pgas3d,tsz
  enddo
  close(1)
  close(2)
END PROGRAM Compute_Profile_Battaglia
!-----------------------------------------------------------------------------
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 ypgas(l,y) ! r/r200
  USE battaglia
  IMPLICIT none
  double precision :: l,x,y
  double precision :: xc,g=-0.3d0,a=1d0,b,P0
  x=dsqrt(l**2d0+y**2d0) ! r/r200
  xc=0.497d0*(m200/h0/1d14)**(-0.00865d0)*(1d0+zz)**(0.731d0)
  P0=18.1d0*(m200/h0/1d14)**(0.154d0)*(1d0+zz)**(-0.758d0)
  b=4.35d0*(m200/h0/1d14)**(0.0393d0)*(1d0+zz)**(0.415d0)
  ypgas=P0*(x/xc)**g/(1d0+(x/xc)**a)**b
  return
END FUNCTION ypgas