!----------------------------------------------------------------------------- ! ! This subroutine converts the virial mass, Mvir, to an over-density mass, ! Mdel, for a given over-density, delrho. The dark matter density profile ! is given by a Navarro-Frenk-White profile, which is parametrized by ! a scale radius (rs) and a concentration parameter (c). ! ! INPUTS: Mvir, rs, c, rhodel [double precision] ! OUTPUT: Mdel [double precision] ! ! For details, see, e.g., Eq (D15) of Komatsu et al. (WMAP 7-year paper), ! arXiv:1001.4538. ! ! December 10, 2010 ! E. Komatsu ! (originally written in February 14, 2002) ! !----------------------------------------------------------------------------- SUBROUTINE mvir2mdel(Mvir,rs,c,delrho,Mdel) IMPLICIT none double precision, intent(IN) :: Mvir,rs,c,delrho double precision, intent(OUT):: Mdel integer :: i,j,iter,imax=50 double precision, allocatable :: mtest(:) double precision :: ltest,m1,m2,ml,mu,lnMdel,var(4) double precision :: zbrent,mv2md external zbrent,mv2md var(1:4)=(/Mvir,rs,c,delrho/) allocate(mtest(0:imax)) ! Find initial mass range for zbrent (non-linear equation solver), ! within which Mdel will be found. mtest(0)=Mvir ! initial guess for the overdensity mass ltest=mv2md(dlog(mtest(0)),var) if(ltest.le.0d0)then do i=1,imax mtest(i)=mtest(i-1)*2.d0 ltest=mv2md(dlog(mtest(i)),var) if (ltest.gt.0d0)then m1=dlog(mtest(i)) m2=dlog(mtest(i-1)) exit endif enddo else do i=1,imax mtest(i)=mtest(i-1)/2.d0 ltest=mv2md(dlog(mtest(i)),var) if(ltest.lt.0d0)then m1=dlog(mtest(i)) m2=dlog(mtest(i-1)) exit endif enddo endif ! Lower & upper limits for M, between which Mdel will be found. ml=min(m1,m2) mu=max(m1,m2) ! Find Mdel using zbrent. lnMdel=zbrent(mv2md,ml,mu,1d-4,iter,var) Mdel=dexp(lnMdel) return END SUBROUTINE mvir2mdel !----------------------------------------------------------------------------- double precision function mv2md(lnM,var) ! zbrent will solve the equation, mv2md(lnMdel)=0, for lnMdel. IMPLICIT none double precision, intent(IN) :: lnM,var(4) double precision :: Mdel,rdel,func,fourpi double precision :: Mvir,rs,c,delrho external func fourpi=4d0*3.141592653589793238d0 Mdel=dexp(lnM) Mvir=var(1) rs=var(2) c=var(3) delrho=var(4) rdel=(Mdel/(fourpi/3d0*delrho))**(1d0/3d0) mv2md=Mdel/Mvir-func(rdel/rs)/func(c) ! See, e.g., Eq.(D15) of Komatsu et al., arXiv:1001.4538. return END function mv2md !----------------------------------------------------------------------------- double precision function func(x) IMPLICIT none double precision, intent(IN) :: x func=dlog(1d0+x)-x/(1d0+x) return END function func