*** Proper Luminosity Distance in units of h^-1 Mpc ***
September 27, 2008: E.Komatsu

Here we provide a simple subroutine, dl.f90, for computing 
the proper luminosity distance, Dl(z), in units of h^-1 Mpc.

IMPORTANT: the radiation density is ignored in the calculation, and
therefore this routine cannot be used to compute the luminosity 
distance to the decoupling epoch, z=1090, as radiation
density should not be ignored there. For the other applications
this approximation should be OK.

(1) Put "USE cosmo" at the beginning of the program
(2) Put "USE luminosity_distance" at the beginning of the program
(3) Specify three cosmological parameters:                       
   - matter density parameter, om0
   - dark energy density parameter, ode0
   - dark energy equation of state parameter, w
(4) CALL setup_dl
(5) You are now ready to use a double-precision function "dl(z)" where
   z is also in double precision.

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

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

 USE cosmo
 USE luminosity_distance
 double precision :: dl,redshift
 external dl
 ! 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_dl
 redshift=10d0
 print*,'omega matter=',om0
 print*,'omega de=',ode0
 print*,'w=',w
 print*,'Dl(z) at z=',redshift,'is',dl(redshift),' h^-1 Mpc'
 end

Another program, "compute_dl.f90," for generating Dl(z) over many
redshifts is

PROGRAM Compute_dl
  USE cosmo
  USE luminosity_distance
  ! A sample program for computing the proper luminosity distance
  ! in units of h^-1 Mpc
  ! September 27, 2008: E.Komatsu
  IMPLICIT none
  integer :: j
  double precision :: dl,z
  external dl
! specify four cosmological parameters
  ode0=0.723d0
  om0=0.277d0
  w=-1d0
! tabulate dl(z) by calling "setup_dl"
  call setup_dl
! now output dl(z) as a function of z...
  open(1,file='redshift_dl.txt')
  z=0d0
  do j=1,10000
     write(1,'(2F18.8)')z,dl(z)
     z = z+0.01d0
  enddo
  close(1)
  stop
END PROGRAM Compute_dl