PROGRAM Compute_Transform_KS USE nfw ! A sample code for computing a 2d Fourier transform of the Sunyaev-Zel'dovich effect ! profile, Tsz(x) [x=theta/thetas], from the Komatsu-Seljak profile ! [Komatsu & Seljak, MNRAS, 327, 1353 (2001); MNRAS, 336, 1256 (2002)] ! ! The trasnformation is done via ! ! trasnform=(fourpi*rs)/ls^2 ! *int_{xin}^{xout}[x^2*Tsz(x)*sin(l*x/ls)/(l*x/ls)]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/ls ! ! December 11, 2010: E.Komatsu IMPLICIT none double precision :: xin=1d-6,xout=3d0 ! pressure profile is cut at r=xout*rvir double precision :: omega,om0=0.277d0,h0=0.702d0,obh2=0.02262d0,tol=1d-6 double precision :: mvir,rvir,thetavir,ls,rs,delc double precision :: eta0,rhogas0,kTgas0 double precision :: m500,r500,theta500,l500,rhoc,Ez,da,z double precision :: y,transform2d,tsz double precision :: chi,ygas,ypgas,rombint,pi=3.1415926535d0 integer :: i external chi,ygas,ypgas,rombint ! Sample parameters: Coma cluster 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 Seljak (2000) c=10d0*(mvir/3.42d12)**(-0.2)/(1d0+z) ! One may also use the concentration parameter from Duffy et al. (2008) ! c=7.85d0*(mvir/2d12)**(-0.081)/(1d0+z)**0.71 ! Duffy et al. (2008) rs=rvir/c ! NFW scale radius, h^-1 Mpc ls=da/rs print*,'rs=',rs,' h^-1 Mpc' print*,'ls=',ls ! calculate profiles rhogas0=7.96d13*(obh2/om0)*(mvir/1d15)/rvir**3d0 & *c**2d0/(1d0+c)**2d0/ygas(c)/(dlog(1d0+c)-c/(1d0+c)) eta0=2.235d0+0.202d0*(c-5d0)-1.16d-3*(c-5d0)**2d0 kTgas0=8.8d0*eta0*(mvir/1d15)/rvir open(1,file='x_transform.txt') do i=1,1000 y=xin + 1d-3*dble(i)*(xout-xin) ! y=l/ls ! See, e.g., Appendix D2 of Komatsu et al., arXiv:1001.4538 ! Eq.(2) of Komatsu&Seljak (2002) transform2d=55d0*(rhogas0/1d14)*(kTgas0/8d0)*(rs/h0) & *rombint(ypgas,c*xin,c*xout,tol,y)*(4d0*3.14159d0)/ls**2d0 Tsz=-2d0*283d0*(transform2d/50d0) ! uK write(1,'(2E18.5)')y,tsz enddo close(1) close(2) END PROGRAM Compute_Transform_KS !----------------------------------------------------------------------------- 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 ygas(x) ! r/rs USE nfw IMPLICIT none double precision :: x double precision :: gam,eta0,B gam=1.137d0+8.94d-2*dlog(c/5d0)-3.68d-3*(c-5d0) eta0=2.235d0+0.202d0*(c-5d0)-1.16d-3*(c-5d0)**2d0 B=3d0/eta0*(gam-1d0)/gam/(dlog(1d0+c)/c-1d0/(1d0+c)) ygas=(1d0-B*(1d0-dlog(1d0+x)/x))**(1d0/(gam-1d0)) return END FUNCTION ygas !----------------------------------------------------------------------------- DOUBLE PRECISION FUNCTION ypgas(x,ell) ! x=r/rs, ell=l/ls, ypgas=ygas^gamma USE nfw IMPLICIT none double precision :: x,ell double precision :: gam,eta0,B gam=1.137d0+8.94d-2*dlog(c/5d0)-3.68d-3*(c-5d0) eta0=2.235d0+0.202d0*(c-5d0)-1.16d-3*(c-5d0)**2d0 B=3d0/eta0*(gam-1d0)/gam/(dlog(1d0+c)/c-1d0/(1d0+c)) ypgas=(1d0-B*(1d0-dlog(1d0+x)/x))**(gam/(gam-1d0)) & *x**2d0*dsin(ell*x)/(ell*x) return END FUNCTION ypgas