/* ------------- file: -------------------- cosmology.c ------------------- */

/* Evaluates for a given redshift: 

          Hubble expansion factor and its derivative with respect to z
          Hydrogen number density

 * Sara Seager
 * Last modified: May 7, 1999
 */

#include <math.h>

#include "nrutil.h"
#include "constant.h"
#include "input.h"

/* --- Global Variables            ----------------------------------------- */
extern struct Input input;

/* -----------begin -----------------------expansion.c---------------------- */

void H_z(double z, double *Hz, double *DHz)
{
  double b, Fnu, Zeq, TNOW, Const;

  /* -- Calculate H(z)   -----                   --------------------------- */
  TNOW = To;
  /*
  b = 4.0/3.0;
  Fnu = 21.0/8.0*pow((4.0/11.0),b);
  */
  Fnu = 6.813220e-01;
  Zeq = 3.0*SQ(input.H0*cLight)/(8.0*PI*G*aRad*(1.0+Fnu))/(pow(TNOW,4));
  Const = 7.083464e+40;
  Zeq = 3.0*SQ(input.H0)*Const/(1.0+Fnu)/(pow(TNOW,4));
  
  *Hz = input.H0*sqrt(input.OmegaL + (SQ(1.0+z)*(input.OmegaK + (1.0+z)*input.Omega0*(1.0 + (1.0+z)/(1.0+Zeq)))));

  /* --- Calculate d(H(z))/dz     ------------------------------------------ */
   *DHz = input.H0/(2.0*(*Hz))*input.H0*((1.0+z)*(2.0*input.OmegaK + (1.0+z)*input.Omega0*(3.0 + 4.0*(1.0+z)/(1.0+Zeq)))); 

}
/* ----------end ------------------------- expansion.c --------------------- */

/* -----------begin -----------------------density.c------------------------ */

double NHnow(double z) 
{
  double Y = 0.24;
  double mu_H;
  mu_H = 1.0/(1.0 - Y);
  return(3.0*SQ(input.H0)*input.OmegaB/(8.0*PI*G*mHatom*mu_H)*CUBE(1.0+z));

}
/* --------- end ----------------------- density.c --------------- */