/* ------ File: -------------------- main.c ------------------------------ */

/* This is the calling routine for the program RECFAST.
   RECFAST calcualtes H and He I recombination, following the 
   paper Seager, Scott, Sasselov. If this code is used, refer to this
   paper. 

   RECFAST uses a 3-level effective H atom: ground state, 1st excited
    state, and continuum. and He atom, along with 
   the matter temperature equation.
   The 3-level atom:
   the ground state, first excited state and continuum.  

 * INPUT: file in.dat:
     starting redshift
     ending redshift
     stepsize
     Omega total
     Omega baryon
     Omega Lambda
     Omega curvature
     h
     fudge factor for H recombination
 * OUTPUT: file out.dat: 
     redshiftz
     x_e
 *
 * Sara Seager
 * Last modified: May 5, 1999.
 */
 
#include <stdio.h>

#include "integrate.h"
#include "input.h"
#include "nrutil.h"

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

main()
{
  long int i;
  double ystart[4], T;
  double fHe, hl;
  double *zout, **yout, zprint;
  int nout, istart;
  FILE *f1;

  f1 = fopen("out.dat","w");
  zout = dvector(1,20000);
  yout = dmatrix(1,3,1,20000);
  
  /* --- Read in starting values             ------------------------------- */
  printf("Using Hummer's case B recombination rates\n");
  printf("Using HeII singlet only recombination rates\n");
  printf("Using T_0 = 2.728\n");
  printf("Using Y_P = 0.24\n\n");
   
  printf("Enter starting redshift\n");
  scanf("%le", &input.zstart);

  printf("Enter final redshift\n");
  scanf("%le", &input.zend);

  printf("Redshift step (10 maximum)\n");
  scanf("%le", &input.hstep);

  printf("Omega total\n");
  scanf("%le", &input.Omega0);

  printf("Omega in baryons\n");
  scanf("%le", &input.OmegaB);

  printf("Omega in the cosmological constant\n");
  scanf("%le", &input.OmegaL);

  printf("Omega from curvature\n");
  scanf("%le", &input.OmegaK);

  printf("Hubble constant in units of 100kms-1Mpc-1\n");
  scanf("%le", &hl);

  input.H0 = hl*100.0/(1.0e+06*3.0856775807e+13);

  printf("Fudge factor for H recombination? (1.14)\n");
  scanf("%le", &input.F);

  fHe = YP/(4.0*(1.0-YP));

  /* --- Finished entering starting data       ----------------------------- */

  /* --- Only need to start integrating before He I recombination.
         Anything before that time is just a constant ionization 
         fraction. xe = 1.0 + fHe                                        --- */

  /* --- To save time here, try a lower value of z. For example, for cdm
     use z = 3800                                                  --------- */
  if (input.zstart > 5000) {
    istart=input.zstart;
    for (i=istart; i>=5010; i-=10) {
      zprint = i;
      fprintf(f1, "%e %e\n", zprint, 1.0+fHe);
    }
    input.zstart = 5000.0;
  }

  /* --- Set starting solution. Everything is ionized, and the 
     matter temperature = the radiation temperature           
     1 = xHe = nHe+/nHTot
     2 = xp = np/nHTot
     3 = TM (matter temperature)                           ------------  --- */
  ystart[1] = fHe;
  ystart[2] = 1.0;
  ystart[3] = To*(1.0+input.zstart);

printf("ystart %e %e %e\n", ystart[1], ystart[2],ystart[3]);  
ystart[3] = 13632.727;
  /* --- Can substitute any integrator for stiff equations at this point --- */

  integrate(input.zstart, ystart, 3, &nout, zout, yout);  
 
  /* --- Print out           ----------------------------------------------  */

  for (i=1; i<=nout; i++) {
    fprintf(f1,"%e %e\n",zout[i], yout[2][i]+yout[1][i]);
  }

  /* --- Clean up              --------------------------------------------- */
  fclose(f1);
  free_dvector(zout,1,20000); free_dmatrix(yout,1,20000,1,20000);
  return 0;

}