!-----------------------------------------------------------------------------
!
! 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