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