*** Converting Virial Mass to Overdensity Mass (in units of h^-1 Msun) ***
December 11, 2010

Here we provide a simple subroutine, mvir2mdel.f90, for converting the 
virial masses, Mvir, to overdensity masses, Mdel, for a given overdensity, delrho.

Here is a call:
   CALL mvir2mdel(Mvir,rs,c,delrho,Mdel)

where rs and c are a scale radius and a concentration parameter of a NFW profile
(Navarro-Frenk-White profile). 

This subroutine solves the following non-linear equation:

Mdel/func(rdel/rs) = Mvir/func(c)

for a given Mvir, rs, and c. Here, the overdensity radis, rdel, is given by

rdel = (Mdel/(fourpi/3d0*delrho))**(1/3).

- To compile and use the sample programs (given below), edit Makefile
and simply "./make"
- It will generate executables called "sample" and "compute_mdel"

=========================================================================
A simple program to use "mvir2mdel.f90" (also included as "sample.f90" in
this directory):

 double precision :: mvir,z,rs,c,rho,del,mdel,rvir
 double precision :: pi=3.1415926535d0,Ez,rhoc,delc
 double precision :: om0=0.277d0,omega
 z=0.1d0
 del=500d0
 mvir=2d14 ! h^-1 Msun
 print*,'z=',z
 print*,'overdensity (del)=',del
 print*,'Mvir=',mvir,' h^-1 Msun'
! 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)
 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)
 print*,'Mdel=',mdel,' h^-1 Msun for Seljak (2000) concentraion'
 ! 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)
 print*,'Mdel=',mdel,' h^-1 Msun for Duffy et al. (2008) concentraion'
 end



Another program, "compute_mdel.f90," for generating Mdel as a function of Mvir is

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