!----------------------------------------------------------------------------- ! ! Proper angular diameter distance, returned in units of h^-1 Mpc. ! Note h^-1! ! ! *** Radiation density is ignored in the calculation, so ! this routine cannot be used to calculate the distance to ! the last scattering surface! *** ! ! << sample code >> ! ! USE cosmo ! USE angular_distance ! double precision :: da,redshift ! external da ! ! Specify three cosmological parameters in double precision. ! ! The data type has been defined in MODULE cosmo. ! ode0=0.723d0 ! om0=0.277d0 ! w = -1d0 ! CALL setup_da ! redshift=10d0 ! print*,'omega matter=',om0 ! print*,'omega de=',ode0 ! print*,'w=',w ! print*,'Da(z) at z=',redshift,'is',da(redshift),' h^-1 Mpc' ! end ! ! August 23, 2008 ! E. Komatsu ! !----------------------------------------------------------------------------- MODULE angular_distance INTEGER, parameter :: nw=501 DOUBLE PRECISION, dimension(nw) :: xa,za,z2a contains SUBROUTINE Setup_Da USE cosmo ! tabulate da(x), where x=ln(a) IMPLICIT none INTEGER :: i DOUBLE PRECISION :: x,eta,tol=1d-7,rombint,ok0 EXTERNAL one_over_h ok0=1d0-om0-ode0 do i=1,nw x = (i-nw)*0.01 ! x=ln(a)=-5.00,-4.99,...,0 xa(i) = x eta=rombint(one_over_h,0d0,dexp(-x)-1d0,tol) if(abs(ok0)<1.d-5)then za(i)=2998d0*exp(x)*eta ! flat model else if(ok0>0)then za(i)=2998d0*exp(x)*sinh(sqrt(ok0)*eta)/sqrt(ok0) ! open model else za(i)=2998d0*exp(x)*sin(sqrt(-ok0)*eta)/sqrt(-ok0) ! closed model endif enddo CALL spline(xa,za,nw,1.d30,1.d30,z2a) ! evaluate z2a(i)=d^2[da(i)]/dy^2 return END SUBROUTINE Setup_Da END MODULE angular_distance !----------------------------------------------------------------------------- DOUBLE PRECISION FUNCTION da(z) ! proper distance in units of h^-1 Mpc. USE angular_distance ! Returns da(z) by interpolating the table generated by Setup_Da. IMPLICIT none DOUBLE PRECISION, intent(IN) :: z DOUBLE PRECISION :: a,b,hh,x INTEGER :: jlo x=-dlog(1.+z) ! x=ln(a)=-ln(1+z) CALL hunt(xa,nw,x,jlo) hh=xa(jlo+1)-xa(jlo) a=(xa(jlo+1)-x)/hh b=(x-xa(jlo))/hh da=(a*za(jlo)+b*za(jlo+1)+((a**3-a)*z2a(jlo)+(b**3-b)*z2a(jlo+1))*(hh**2)/6.) return END FUNCTION Da !----------------------------------------------------------------------------- DOUBLE PRECISION FUNCTION one_over_h(z) USE cosmo ! one_over_h(z) = H_0/H(z) (dimensionless) ! x=a IMPLICIT none DOUBLE PRECISION, intent(IN) :: z double precision :: x,ok0,or0=0d0 ! radiation density is ignored! ok0=1d0-om0-ode0 x=1d0/(1d0+z) one_over_h = 1d0/sqrt(om0/x**3d0+ok0/x**2d0+or0/x**4d0+ode0/x**(3d0+3d0*w)) return END FUNCTION one_over_h