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