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