/* Code to calculate sigma8 for COBE-normalized CDM and MDM models, 
using the approximation of Eisenstein & Hu (1997). */

/* This code will calculate sigma_8 and other quantities for CDM and
MDM models.  Because the fitting form for the transfer function does not
include the baryon oscillations, the answer will not be accurate for
large baryon fractions.  The COBE normalization is from the papers
on Bunn & White (1997) and Hu & White (1997). */

/* Accuracy is generally better than 5%.  Recall that the COBE normalization
is only 7% accurate (1-sigma).  This program gives low answers compared to
CMBfast for lower values of Omega_0*h^2, and vice versa for higher values. 
Values of Omega_0 near 0.2 can produce answers that are 5% low. */

/* The integrator used here is a bit crude, but suffices for rapid
10^-3 accuracy.  You could get better performance by substituting an
integrator from a numerical package.  See the end of the file for details. */

/* This program uses the routines TFmdm_set_cosm() and TFmdm_onek_hmpc(),
available from http://www.sns.ias.edu/~whu/power/power.html in the file
power.c.  To compile this program, you'll need that file as well, e.g.
	cc -O sigma8.c power.c -o sigma8 -lm 		*/

/* Lists of versions and bugs:
v1.0 -- 10/23/97 -- Original release version.
v1.1 -- 11/13/97 -- We correct a missing factor of omega_m in the
		    conversion between mass and length scales [lines
		    32 and 34 of main()].  This only affects the results 
		    if you explicitly used the -mass option.
v1.2 -- 04/08/98 -- We correct a typo in the second-order tilt scaling
		    of the COBE normalization for flat models with tensors.
v1.3 -- 09/15/98 -- We correct a typo in the tilt scaling of the COBE
		    normalization for open models with tensors.
v1.4 -- 12/02/98 -- We correct the COBE normalization for tilted open models,
		    undoing the change of v1.3 and fixing what was
		    wrong there to begin with.
*/

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

int TFmdm_set_cosm(double omega_matter, double omega_baryon, double omega_hdm,
	int degen_hdm, double omega_lambda, double hubble, double redshift);
double TFmdm_onek_hmpc(double kk);

double sig2tophat(double x);
double sig2gauss(double x);
double sig2sharpk(double x);
double sig2tophatnu(double x);
double sig2gaussnu(double x);
double sig2sharpknu(double x);

#define WEIRD  -2538	/* A strange number, unlikely to be entered */
#define RHOCRIT 2.78e11	/* The critical density, for h=1, in Msun/Mpc^3 */

/* Convenience from Numerical Recipes in C, 2nd edition
------------ COMPLETEMENT DEBILE (cmv) ????????? ------------
static double sqrarg;
#define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg)
*/
#define SQR(a) (a*a)

/* Global variable */
double   omega_m, omega_b, omega_n, 
	numdegen, lambda, redshift, 
	hubble, scale, mass, tilt;
double   deltaH;		/* COBE-normalization */
double   sig2_scale;	/* For passing a length scale to the integrator */
int 	qtensors, quiet, qredshift;

/* ========================== User interaction routines ============== */
void usage()
/* Print out information about the allowed options and then quit. */
{
    fprintf(stderr,"Command line options:\n");
    fprintf(stderr,"-om or -omega <omega_matter, including CDM, baryons, massive neutrinos>\n");
    fprintf(stderr,"-ob <omega_baryon> -fb <baryon_fraction>\n");
    fprintf(stderr,"-on <omega_neutrino> -fn <neutrino_fraction>\n");
    fprintf(stderr,"-numdegen or -Nnu <number of degenerate neutrino species>, default 1\n");
    fprintf(stderr,"-lambda <omega_lambda> -flat (choose lambda appropriately)\n");
    fprintf(stderr,"-h or -hubble <H_0, units of 100 km/s/Mpc>\n");
    fprintf(stderr,"-H or -Hubble <H_0, units of 1 km/s/Mpc>\n");
    fprintf(stderr,"-tilt or -n <scalar spectral tilt>, no tensors.\n");
    fprintf(stderr,"-Tilt or -N <scalar spectral tilt>, power-law tensors.\n");
    fprintf(stderr,"-z <redshift>, if none given then appropriate value assumed.\n");
    fprintf(stderr,"The following request the rms fluctuations inside a tophat window function:\n");
    fprintf(stderr,"-scale <comoving length scale, h^-1 Mpc>\n");
    fprintf(stderr,"-mass <mass scale, h^-1 Msun>\n");
    fprintf(stderr,"If neither -scale or -mass is given, sigma_8 and other quantities are computed.\n");
    fprintf(stderr,"-quiet (print only answer to stdout)\n");
    exit(1);
    return;
}

void crash(char *s)
{
    fprintf(stderr,"%s\n", s); exit(1); return;
}

void warn(char *s)
{
    fprintf(stderr,"%s\n", s); return;
}


void parse_command_line(int argc, char *argv[])
/* Parse the command line, and look for illegal entries */
{
    int i;
    if (argc==1) usage();
    i = 1;
    omega_m = 1.0;
    omega_b = 0.05;
    omega_n = 0.0;
    numdegen = 1.0;
    lambda = 0.0;
    redshift = 0.0;
    hubble = 0.5;
    scale = -1.0;
    mass = -1.0;
    qtensors = 0;
    quiet = 0;
    tilt = 1.0;
    qredshift = 0;

    while (i<argc) {
	if (!strcmp(argv[i],"-om")||!strcmp(argv[i],"-omega")) {
	    ++i;
	    if (i>=argc) usage();
	    sscanf(argv[i],"%lf",&omega_m); ++i;
	} else if (!strcmp(argv[i],"-ob")) {
	    ++i;
	    if (i>=argc) usage();
	    sscanf(argv[i],"%lf",&omega_b); ++i;
	} else if (!strcmp(argv[i],"-fb")) {
	    ++i;
	    if (i>=argc) usage();
	    sscanf(argv[i],"%lf",&omega_b); ++i;
	    omega_b *= -1.0;	/* Signal that a fraction was desired */
	} else if (!strcmp(argv[i],"-on")) {
	    ++i;
	    if (i>=argc) usage();
	    sscanf(argv[i],"%lf",&omega_n); ++i;
	} else if (!strcmp(argv[i],"-fn")) {
	    ++i;
	    if (i>=argc) usage();
	    sscanf(argv[i],"%lf",&omega_n); ++i;
	    omega_n *= -1.0;	/* Signal that a fraction was desired */
	} else if (!strcmp(argv[i],"-h")||!strcmp(argv[i],"-hubble")) {
	    ++i;
	    if (i>=argc) usage();
	    sscanf(argv[i],"%lf",&hubble); ++i;
	} else if (!strcmp(argv[i],"-H")||!strcmp(argv[i],"-Hubble")) {
	    ++i;
	    if (i>=argc) usage();
	    sscanf(argv[i],"%lf",&hubble); ++i;
	    hubble /= 100.0;	/* Units were 1 km/s/Mpc */
	} else if (!strcmp(argv[i],"-lambda")) {
	    ++i;
	    if (i>=argc) usage();
	    sscanf(argv[i],"%lf",&lambda); ++i;
	} else if (!strcmp(argv[i],"-flat")) {
	    ++i;
	    lambda = WEIRD;
	} else if (!strcmp(argv[i],"-numdegen")||!strcmp(argv[i],"-Nnu")) {
	    ++i;
	    if (i>=argc) usage();
	    sscanf(argv[i],"%lf",&numdegen); ++i;
	} else if (!strcmp(argv[i],"-n")||!strcmp(argv[i],"-tilt")) {
	    ++i;
	    if (i>=argc) usage();
	    sscanf(argv[i],"%lf",&tilt); ++i;
	    qtensors = 0;
	} else if (!strcmp(argv[i],"-N")||!strcmp(argv[i],"-Tilt")) {
	    ++i;
	    if (i>=argc) usage();
	    sscanf(argv[i],"%lf",&tilt); ++i;
	    qtensors = 1;
	} else if (!strcmp(argv[i],"-z")) {
	    ++i;
	    if (i>=argc) usage();
	    sscanf(argv[i],"%lf",&redshift); ++i;
	    qredshift = 1;
	} else if (!strcmp(argv[i],"-scale")) {
	    ++i;
	    if (i>=argc) usage();
	    sscanf(argv[i],"%lf",&scale); ++i;
	} else if (!strcmp(argv[i],"-mass")) {
	    ++i;
	    if (i>=argc) usage();
	    sscanf(argv[i],"%lf",&mass); ++i;
	} else if (!strcmp(argv[i],"-quiet")) {
	    ++i;
	    quiet = 1;
	} else usage();
    }

    /* Look for illegal entries */
    if (omega_m<=0.0) crash("Omega must be positive.");
    if (hubble<=0.0) crash("Hubble constant must be positive.");
	else if (hubble>10.0) warn("Did you get the Hubble units right?");
    if (numdegen<0) {
	warn("Negative number of degenerate species reset to 1.");
	numdegen = 1.0;
    }
    if (redshift<=-1.0) crash("Redshift must exceed -1.");
    if (scale>0.0 && mass>0.0) {
	warn("Input length OR mass scale.  Defaulting to length.");
	mass = -1.0;
    }

    /* Now look for fractions, if they were input. */
    if (omega_b<0.0) omega_b = omega_m*(-omega_b);	
    if (omega_n<0.0) omega_n = omega_m*(-omega_n);	
    if (lambda== WEIRD) lambda = 1-omega_m;	/* Set to be flat */

    if (omega_b+omega_n>omega_m) crash("Baryons + Neutrinos exceeds Omega_0.");
	else if (omega_b+omega_n>0.6*omega_m) 
	    warn("CDM fraction is below range of TF fitting formula.");
    if (tilt>1.3 || tilt < 0.7) 
	warn("Tilt out of range of Bunn and White fitting formula.");
    if (tilt!=1 && (omega_m<0.2 || omega_m>1.0)) 
	warn("Omega out of range of Bunn and White fitting formula.");
    if (tilt==1 && (omega_m<0.2 || omega_m>1.6 || lambda<0.0 || lambda>0.8))
	warn("Omega or Lambda out of range of Bunn and White fitting formula.");
    if (qtensors!=0&&(tilt>1.0||
	(fabs(lambda)>1e-5&&fabs(1.0-omega_m-lambda)>1e-5))) 
	crash("Tensor formulae only hold for tilts < 1 and simple cosmologies.");
    if (omega_b/omega_m>0.4) warn("Baryon oscillations will be large.  TF fitting form will be inaccurate.");
    if (omega_n/omega_m>0.4) warn("Neutrino fraction outside of range of TF fitting form.");
    if (omega_m*hubble*hubble>0.4) warn("Omega_0 h^2 outside of range of TF fitting form.");
    if (omega_m*hubble*hubble<0.06) warn("Omega_0 h^2 outside of range of TF fitting form.");
    if (redshift>29) warn("Redshift is higher than range of TF fitting form.");

    return;
}

/* ==================================== MAIN() ========================= */

void main(int argc, char *argv[])
{
    double erfcc(double x);
    double sigma, sigma8;
    extern double growth_to_z0, z_equality, alpha_nu, beta_c;
    double cobenorm();
    double integrate(double (*func)(double));
    parse_command_line(argc, argv);

    /* Set all the background quantities */
    TFmdm_set_cosm(omega_m, omega_b, omega_n, (int) numdegen, lambda, 
			hubble, redshift);
    deltaH = cobenorm();

    if (!quiet) {	/* Echo cosmological information to the screen */
	printf("Omega_0  = %5.3f, Omega_baryon = %5.3f, Omega_neutrino = %5.3f\n",
		omega_m, omega_b, omega_n);
	printf("Hubble   = %5.3f, Lambda       = %5.3f, N_neutrinos    = %d\n",
		hubble, lambda, (int) numdegen);
	printf("Redshift = %5.3f, Growth = %6.4f relative to z=0\n",
		redshift, growth_to_z0);
	printf("Tilt     = %5.3f, COBE normalization: deltaH = %8.3e",
		tilt, deltaH);
	if (tilt!=1.0) 
	    if (qtensors) printf(" with tensors.\n");
	    else printf(" without tensors.\n");
	else printf("\n");
    }
    /* If we've been asked to compute for a particular scale, do that */
    if (scale>0 || mass>0) {
	if (mass<0)
	    mass = 4.0*M_PI/3.0*omega_m*RHOCRIT*scale*scale*scale;
	else if (scale<0) 
	    scale = pow(3.0/4.0/M_PI/RHOCRIT/omega_m*mass, 1.0/3.0);
	/* Integrate tophat for the given scale */
	sig2_scale = scale;
	sigma = sqrt(integrate(sig2tophatnu));

	/* Print out the result */
	if (!quiet) printf("For length scale %6.4f h^-1 Mpc, mass scale %8.3e h^-1 Msun.\nrms tophat fluctuations are ",
		sig2_scale, mass);
	printf("%6.4f\n", sigma);
    } else {	
	/* No scale was requested, so we do the standard tasks */
	/* Calculate sigma_8 */
	sig2_scale = 8.0;	
	sigma8 = sqrt(integrate(sig2tophatnu));
	if (quiet) {
	    printf("%6.4f\n", sigma8);
	    return;	/* This is all we're doing */
	}
	printf("Sigma_8  = %6.4f\n", sigma8);

	/* Calculate sigma_50 */
	sig2_scale = 50.0;
	sigma = sqrt(integrate(sig2tophatnu));
	printf("Sigma_50 = %6.4f   Sigma_50/sigma_8 = %6.4f\n",
		sigma, sigma/sigma8);

	/* We default to z=4 for DLA */
	if (!qredshift) {
	    redshift = 4.0;		
	    TFmdm_set_cosm(omega_m, omega_b, omega_n, (int) numdegen, lambda, 
			hubble, redshift);
	}
	sig2_scale = 50.0/100.0/sqrt(1+redshift)/sqrt(2*3.14159)*
		pow(8.0*3.14159/3.0/sqrt(89.0),1./3.);
	sigma = sqrt(integrate(sig2gaussnu));
	printf("Omega_gas = %8.3e at redshift %5.3f due to \n   sigma(50 km/s) = sigma(%5.3f h^-1 Mpc) = %6.4f (DLA, Klypin et al. 1994)\n", 
	    omega_b*erfcc(1.33/1.414/sigma), redshift, sig2_scale, sigma);

	/* We default to z=3 for Lyman alpha forest */
	if (!qredshift) {
	    redshift = 3.0;		
	    TFmdm_set_cosm(omega_m, omega_b, omega_n, (int) numdegen, lambda, 
			hubble, redshift);
	}
	sig2_scale = 1.0/(34.0*sqrt(omega_m/2.0));	/* Gnedin's k_34 */
	sigma = sqrt(integrate(sig2gauss));
	printf("Sigma_Lya = %6.4f at redshift %5.3f (Ly-alpha forest, Gnedin 1997)\n", 
		sigma, redshift);

	printf("Small-scale: alpha_nu = %6.4f, beta_c = %6.4f\n", 
		alpha_nu, beta_c);
    }
    return;
}

/* ========================== Integration Kernals ======================= */

double W2tophat(double x)
/* The square of the tophat window */
{
    double j1onx;
    if (x<0.03) 
	j1onx = (1.0/3.0-x*x*1.0/30.0);
    else j1onx = (cos(x)-sin(x)/x)/x/x;
    return 9.0*j1onx*j1onx;
}

double sig2tophat(double x)
/* This is the integrand for the mass fluctuation integral
sig2 = \int_0^\infty (dk/k) Delta^2(k) W^2(kr)
broken and recombined so that one integrates on (0,1] in a new variable x. */
/* Here, the window is a real-space tophat of radius sig2_scale, the
latter being passed from outside in units of h^-1 Mpc */
{
    double powerspcb(double k);
    return (powerspcb(x/sig2_scale)*W2tophat(x)+
 	    powerspcb(1.0/x/sig2_scale)*W2tophat(1.0/x))/x;
}

double sig2tophatnu(double x)
/* Same as sig2tophat(), but with P_cbnu */
{
    double powerspcbnu(double k);
    return (powerspcbnu(x/sig2_scale)*W2tophat(x)+
 	    powerspcbnu(1.0/x/sig2_scale)*W2tophat(1.0/x))/x;
}

double sig2gauss(double x)
/* This is the integrand for the mass fluctuation integral
sig2 = \int_0^\infty (dk/k) Delta^2(k) W^2(kr)
broken and recombined so that one integrates on (0,1] in a new variable x. */
/* Here, the window is a gaussian of radius sig2_scale, the
latter being passed from outside in units of h^-1 Mpc */
{
    double powerspcb(double k);
    return (powerspcb(x/sig2_scale)*exp(-x*x)+
 	    powerspcb(1.0/x/sig2_scale)*exp(-1.0/x/x))/x;
}

double sig2gaussnu(double x)
/* Same as sig2gauss(), but with P_cbnu */
{
    double powerspcbnu(double k);
    return (powerspcbnu(x/sig2_scale)*exp(-x*x)+
 	    powerspcbnu(1.0/x/sig2_scale)*exp(-1.0/x/x))/x;
}

double sig2sharpk(double x)
/* This is the integrand for the mass fluctuation integral
sig2 = \int_0^\infty (dk/k) Delta^2(k) W^2(kr)
broken and recombined so that one integrates on (0,1] in a new variable x. */
/* Here, the window is a k-space tophat of radius 1/sig2_scale, the
latter being passed from outside in units of h^-1 Mpc */
{
    double powerspcb(double k);
    return (powerspcb(x/sig2_scale))/x;
}

double sig2sharpknu(double x)
/* Same as sig2sharpk(), but with P_cbnu */
{
    double powerspcbnu(double k);
    return (powerspcbnu(x/sig2_scale))/x;
}

/* ============================ Power Spectrum and COBE ================ */
double powerspcbnu(double k)
/* Returns Delta^2(k), COBE-normalized, based on the approximations of 
Eisenstein & Hu (1997).  k is in h^-1 Mpc. */
/* TFmdm_set_cosm() and cobenorm() must have been called before this */
/* Returns the density-weighted CDM+Baryon+Neutrino power spectrum */
{
    extern double tf_cbnu, growth_to_z0;
    TFmdm_onek_hmpc(k);  /* This sets the value in tf_cbnu */
    return pow(2997.0*k, tilt+3.0)*SQR(deltaH*tf_cbnu*growth_to_z0);
}

double powerspcb(double k)
/* Returns Delta^2(k), COBE-normalized, based on the approximations of 
Eisenstein & Hu (1997).  k is in h^-1 Mpc. */
/* TFmdm_set_cosm() and cobenorm() must have been called before this */
/* Returns the density-weighted CDM+Baryon power spectrum */
{
    extern double growth_to_z0;
    return pow(2997.0*k, tilt+3.0)*SQR(deltaH*TFmdm_onek_hmpc(k)*growth_to_z0); 
}

double cobenorm()
/* Return the Bunn & White (1997) fit for delta_H */
/* Given lambda, omega_m, qtensors, and tilt */
/* Open model with tensors is from Hu & White */
{
    double n;
    n = tilt-1;
    if (fabs(omega_m+lambda-1.0)<1e-5) {	/* Flat universe */
	if (qtensors) 
	    return 1.94e-5*pow(omega_m, -0.785-0.05*log(omega_m))*
		exp(n+1.97*n*n);
	else return 1.94e-5*pow(omega_m, -0.785-0.05*log(omega_m))*
		exp(-0.95*n-0.169*n*n);
	/* Error -- was 0.0169 in earlier version of code */
    } else if (fabs(lambda)<1e-5) {	/* No lambda */
	if (qtensors)
	    return 1.95e-5*pow(omega_m,-0.35-0.19*log(omega_m)-0.15*n)*
		exp(+1.02*n+1.7*n*n);
	 /* Error -- was exp(-1.02n-1.7*n*n) in earlier version (v1-1.3)*/
	else return 1.95e-5*pow(omega_m, -0.35-0.19*log(omega_m)-0.17*n)*
		exp(-n-0.14*n*n);
	 /* Error -- was exp(-n-0.14*n*n) in earlier version (v1-1.2)*/
	 /* Error -- was exp(n+0.14*n*n) in earlier version (v1.3) */
    } else return 1e-5*(2.422-1.166*exp(omega_m)+0.800*exp(lambda)
	+3.780*omega_m-2.267*omega_m*exp(lambda)+0.487*SQR(omega_m)+
	0.561*lambda+3.329*lambda*exp(omega_m)-8.568*omega_m*lambda+
	1.080*SQR(lambda));
}

/* The following routine is from Numerical Recipes in C, 2nd edition,
by Press, Teulkolsky, Vetterling, and Flannery, (Cambridge Univ. Press) */
double erfcc(double x)
{
	double t,z,ans;
	z=fabs(x);
	t=1.0/(1.0+0.5*z);
	ans=t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+
		t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+
		t*(-0.82215223+t*0.17087277)))))))));
	return x >= 0.0 ? ans : 2.0-ans;
}

/* ======================= Integration Engine ======================= */

double midpoint_int(double (*func)(double))
/* Integrate the function func() on the range (0,1] */
/* Use midpoint rule, expanding by factor of 3 each time */
/* Linearly extrapolate to zero stepsize.  Not fancy */
{
    int iter, bins, j;
    double width, xeval, sum[11];

    sum[0] = (*func)(0.5);	/* First step is easy */
    bins = 1; width = 1.0;
    for (iter=1; iter<11; iter++) {	/* Now start increasing the bins */
	bins *= 3;
	width /= 3.0;
	sum[iter]=0.0;
	for (xeval=width/2.0,j=1;j<=bins;xeval+=width,j++) {
	    if (j%3==2) continue;	/* We've already done this one */
	    sum[iter]+=(*func)(xeval);
	}
	/* Add the new evaluations to the old ones. */
	sum[iter]=sum[iter]*width+sum[iter-1]/3.0;
	/* printf("%d iteration, sum = %f\n", iter, sum[iter]); */
	/* Now decide if we're done */
	if (iter<4) continue;	/* We certainly want to do some work */
	/* If error scales as width^2, then linear extrapolation in
	quantity width^2 is sum[iter]+(sum[iter]-sum[iter-1])/8.0 */
	if (fabs( (sum[iter]+(sum[iter]-sum[iter-1])/8.0)/
		  (sum[iter-1]+(sum[iter-1]-sum[iter-2])/8.0)-1.0)<1e-4) 
	    return sum[iter]+(sum[iter]-sum[iter-1])/8.0;
	/* if (fabs(sum[iter-1]/sum[iter]-1)<1e-4) return sum[iter];  */
    }
    crash("Too many steps in the integration.");
}

double integrate(double (*func)(double))
/* Integrate the function func() on the range (0,1] */
{
    double midpoint_int(double (*func)(double));
    return midpoint_int(func);
}

/* The above is a simple integrator.  You could get faster execution times
by using a more accurate integrator.  The following is an example from 
Numerical Recipes, but all that is needed is for the call to integrate()
to return the value of func() integrated on the range (0,1]. */

#ifdef YES_NR
#include "qromo.c"
#include "midpnt.c"
#include "polint.c"
#include "nrutil.c"

float integrate(float (*func)(float))
/* Integrate the function func() on the range (0,1] */
{
    float qromo(float (*func)(float), float a, float b,
	float (*choose)(float(*)(float),float,float,int));
    float midpnt(float (*func)(float), float a, float b, int n);
    return qromo(func,0,1,midpnt);
}
#endif