/* The following gets rid of longs, which are not needed (ints are 32 bit). */
#define long int
#define IntN3 151
#define Pi 3.1415926535897932385
#define PI 3.1415926535897932385

#define SAFETY 0.9
#define PGROW -0.2
#define PSHRNK -0.25
#define ERRCON 1.89e-4
#define MAXSTP 300000
#define TINY 1.0e-30
#define ITMAX 600
#define EPS 3.0e-8
#define NR_END 1
#define FREE_ARG char*

/* Conveniences: */
static double sqrarg;
#define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg)
static double cubarg;
#define CBE(a) ((cubarg=(a)) == 0.0 ? 0.0 : cubarg*cubarg*cubarg)
#define pi 3.1415926535897932385

static double dmaxarg1,dmaxarg2;
#define DMAX(a,b) (dmaxarg1=(a),dmaxarg2=(b),(dmaxarg1) > (dmaxarg2) ?\
        (dmaxarg1) : (dmaxarg2))

static double dminarg1,dminarg2;
#define DMIN(a,b) (dminarg1=(a),dminarg2=(b),(dminarg1) < (dminarg2) ?\
        (dminarg1) : (dminarg2))

/* 
   Function prototypes: 
*/
double CDMdeltac(double , double , double);
double del_crit0(double , double);
double fomega(double , double , double);
double rombint(double (*func)(double), double , double , double);
double integrate(double (*func)(double), double , double , double);
int dround(double);
double dmin2(double, double);
double sigmatop(double);
double powerspec(double, double *);
double sigmatop2(double);
double dsigmatop(double);
double dsigmatop2(double);
double sigmaOfM(double);
double corr(double , double);
double DELc(double , double , double);
void crash(char *);
double brent(double ax, double bx, double cx, double (*f)(double), double tol,
	     double *xmin);
double gammln(double);
double **matrix(long nrl, long nrh, long ncl, long nch);
void free_matrix(double **m, long nrl, long nrh, long ncl, long nch);
void mnbrak(double *ax, double *bx, double *cx, double *fa, double *fb, double *fc,
	    double (*func)(double));
double f1dim(double x);
void linmin(double p[], double xi[], int n, double *fret, double (*func)(double []));
int *ivector(long nl, long nh);
void free_ivector(int *v, long nl, long nh);
void Initialize();
double MofS(double);
double ageofa(double);
double aofage(double);
double dagedaofa(double);
double HoverH0(double a);
void derivs(double , double [], double []);
double MofR(double);
double RofM(double);
double ExactGrowth(double z);
double dEGrowth(double a);
double FindMvir(double z2);
double delbOfRi(double relRi);
double del_v(double om0, double omLam);
void crash(char *s);
double integrate(double (*func)(double),double a,double b,double tol);
double rombint(double (*func)(double), double a, double b, double tol);
double xirIntegrand(double lkp);
void nrerror(char error_text[]);
double dp();
int min2(int a, int b);
double window(double x);
double windowk(double x);
double windowprime(double x);
double Tcdm2EandH(double k);
void InitEandH();
double PkRIntegrand(double lr);
double PkgIntegrand(double lr);
double sig8Integrand(double lkp);
double dsign(double x);
double wofz(double z);
double ageofa(double a);
double aofage(double t);
double dagedaofa(double a);
double HoverH0(double a);
void InitSphCol(double z2, double M, double *dbari, double *Rvir, double *dcrit);
void SphColODE(double relRi, double M, double *Ri, double *dbari, double *Msh, 
	       double *Rf, double *dbarf, double *deli, double *delf, double *Vpec, 
	       double *V, double *yf, double delfac);
void derivsSphCl(double H0t, double q[], double dq[]);
double XiOfR1R2(double R1,double R2);
double xir(double kl);
double xir2(double kl);
double GenMeanProf(double R, double x1, double *x2, double z);
double MeanProf(double R);
double NewInfall(double S0oS, double nu2oS0);
double NewInfallm1sig(double alpha, double beta);
double NewInfallp1sig(double alpha, double beta);
double FInfall(double gamma);
double HelpProfile(double k, double R);
double fomega(double om0,double omLam,double z);
double del_crit0(double om0, double omLam);
double CDMdeltac(double z, double omega, double omegal);
void rkqs(double y[], double dydx[], int n, double *x, double htry, double eps,
	  double yscal[], double *hdid, double *hnext,
	  void (*derivs)(double, double [], double []));
void rkck(double y[], double dydx[], int n, double x, double h, double yout[],
	  double yerr[], void (*derivs)(double, double [], double []));
void odeint(double ystart[], int nvar, double x1, double x2, double eps, double h1,
	    double hmin, int *nok, int *nbad,
	    void (*derivs)(double, double [], double []),
	    void (*rkqs)(double [], double [], int, double *, double, double, double [],
			 double *, double *, void (*)(double, double [], double [])));
double *vector(long nl, long nh);
void free_vector(double *v, long nl, long nh);
void spline(double x[], double y[], int n, double yp1, double ypn, double y2[]);
void splint(double xa[], double ya[], double y2a[], int n, double x, double *y);
double zbrent(double (*func)(double), double x1, double x2, double tol,
	      double value);
void linint(double *xa,double *ya,int n,double x,double *y, double *ra, 
	    double *rb, int *kLow);