*** Comoving Volume in units of h^-3 Mpc^3 ***
August 23, 2008: E.Komatsu

Here we provide a simple subroutine, volume.f90, for computing 
the comoving volume in units of h^-3 Mpc^3 per area.

The subroutine may be called as

CALL volume(zmin,zmax,area,vol)

zmin: input, the minimum redshift, double precision
zmax: input, the maximum redshift, double precision
area: input, area on the sky in deg^2, double precision
vol: output, volume in units of h^-3 Mpc^3

IMPORTANT: the radiation density is ignored in the calculation, and
therefore this routine cannot be used to compute the comoving
volume out 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 angular_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_da
(5) CALL volume(zmin,zmax,area,vol) where all of the arguments are 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_volume"

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

 USE cosmo
 USE angular_distance
 double precision :: vol,redshift1,redshift2,area ! in deg^2
 ! 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
 redshift1=1.8d0
 redshift2=3.5d0
 area=426d0
 print*,'omega matter=',om0
 print*,'omega de=',ode0
 print*,'w=',w
 CALL volume(redshift1,redshift2,area,vol)
 print*,'Volume(z) per',area,' deg^2 between z=' &
       ,redshift1,' and',redshift2,' is',vol,' h^-3 Mpc^3'
 end

Another program, "compute_volume.f90," for generating V(z) over many
redshifts is

PROGRAM Compute_volume
  USE cosmo
  USE angular_distance
  ! A sample program for computing the comoving volume per square 
  ! degrees out to a given z, integrated from the present epoch, 
  ! in units of h^-3 Mpc^3 deg^-2 (note h^-3!)
  ! August 23, 2008: E.Komatsu
  IMPLICIT none
  integer :: j
  double precision :: vol,z,area,zmin
! specify three cosmological parameters
  ode0=0.723d0
  om0=0.277d0
  w=-1d0
! tabulate da(z) by calling "setup_da"
  call setup_da
! now output volume(z) per deg^2 [h^-3 Mpc^3] as a function of z...
  open(1,file='redshift_volume.txt')
  zmin=0d0
  area=1d0 ! deg^2
  z=0d0
  do j=1,10000
     CALL volume(zmin,z,area,vol)
     write(1,'(1F18.8,1E18.8)')z,vol
     z = z+0.01d0
  enddo
  close(1)
  stop
END PROGRAM Compute_volume