PROGRAM compute_mdel ! A sample program for computing M500 as a function of the virial mass, ! Mvir, in units of h^-1 Msun. Two concentration parameters are used: ! 1. c_seljak = 10d0*(mvir/3.42d12)**(-0.2)/(1d0+z) [Seljak 2000] ! 2. c_duffy = 7.85d0*(mvir/2d12)**(-0.081)/(1d0+z)**0.71 [Duffy et al. 2008] ! January 4, 2010: E.Komatsu IMPLICIT none double precision :: mvir,z,rs,c,rho,del,mdel,rvir double precision :: pi=3.1415926535d0,Ez,rhoc,delc double precision :: om0=0.277d0,omega integer :: i print*,'redshift?' read*,z ! set the overdensity del=500d0 print*,'overdensity=',del ! 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) ! loop over masses open(1,file='mvir_mdel_seljak.txt') open(2,file='mvir_mdel_duffy.txt') do i=1,601 mvir=10d0**(10+dble(i-1)*1d-2) rvir=(3d0*mvir/4d0/pi/delc/rhoc)**(1./3.) ! virial radius, h^-1 Mpc ! Concentration parameter from Seljak (2000) c=10d0*(mvir/3.42d12)**(-0.2)/(1d0+z) rs=rvir/c ! NFW scale radius, h^-1 Mpc CALL mvir2mdel(Mvir,rs,c,del*rhoc,Mdel) write(1,'(2E13.5)') mvir,mdel ! 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 CALL mvir2mdel(Mvir,rs,c,del*rhoc,Mdel) write(2,'(2E13.5)') mvir,mdel enddo close(1) close(2) END PROGRAM compute_mdel