;------------------------------------------------------------------------------
;+
; NAME:
;  D_c
;
; PURPOSE:
;   Functions to calculate the radial comoving distances for calculating the fkp weight
;   for large-scale structure analysis.
;
; CALLING SEQUENCE:
;   vol = Comov_vol,z1,z2 ; gives you the comoving volume 
;   Dc = Dc_z,z  ; gives you the comoving distances to redshift z
;
; OPTIONAL INPUTS:
;
; OUTPUTS:
;   comoving distances, comoving volume
;                                
; COMMENTS:
;   FIXED cosmology parameters here, one may want to change that, but probably unlikely since this is just for fkp_weight
;
; EXAMPLES:
;
; REQUIREMENTS:
; 
; REVISION HISTORY:
;   20-Apr-2012 Written by Shirley Ho (LBL/CMU)
;-




function Dc_z_int,z

c = 299792.5 ;km/s , so that c/H0 is in units of Mpc
omega_m = 0.274
omega_lambda = 1.0-omega_m
omega_k = 0.0 
H0 = 70.0  ; km/s/Mpc; we can change these numbers if we want
E_z = sqrt(omega_m*(1+z)^3.d + omega_k*(1+z)^2.d + omega_lambda)   ; okay, I am only doing this for LCDM universe
Dc_z_integrand  = c/H0/E_z    

return, Dc_z_integrand
end 

function Dc_z,z   ; return comoving radial distance in Mpc (notice no h!!) 
Dc = qromb('Dc_z_int',0,z,/double)
return,Dc
end 

function Comov_vol,z1,z2  ; in units of Mpc^3 
Dc_z1 = Dc_z(z1)
Vol_z1 = 4./3.*!pi*Dc_z1^3.d
Dc_z2 = Dc_z(z2)
Vol_z2 = 4./3.*!pi*Dc_z2^3.d

Diff_vol = abs(Vol_z2-Vol_z1)
return, Diff_vol
end