void Initialize()
{
  /* Cosmological parameters: */
  omega_matter=OMEGA0=omega;
  omega_baryon=OMEGAB=omegab;
  omega_lambda=OMEGALAMBDA=omegal;
  OMEGAR=0.;
  hubble=hub;
  sigma8=sig8;
  THETA27=2.725/2.7;
  tilt=RUNNINGnS=npower;
  tiltk0=RUNNINGk0=nk0;
  tiltdn=RUNNINGdn=ndndlk;
  omega_curv = 1.0-omega_matter-omega_lambda;

  InitEandH();
  /* power-spectrum normalization according to sigma8 */
  dpp = dp();
  snorm=(dpp/Pi)/sqrt(2.);
}

double wofz(double z)
     /*
       Critical collapse threshold at redshift z
     */
{
  return CDMdeltac(z,omega_matter,omega_lambda);
}

double ageofa(double a)
     /*
        Finds age in years at redshift z=1/a-1.
        Restricted to a flat universe with no radiation.
     */
{
  double t2=a*a*a * ageFac1;
  if (t2 < 1.e-6)
    return ageFac3*(2./3.)*a*sqrt(a/omega_matter)*(1.-t2/6.);
  else
    return ageFac2*log(sqrt(t2)+sqrt(t2+1.));
}

double aofage(double t)
     /*
       Finds scale factor a at age t (in years).
       Restricted to a flat universe with no radiation.
     */
{
  double t2=1.5*t*sqrt(omega_matter)/ageFac3;
  if (SQR(t2)*ageFac1 < 1.e-6)
    return pow(t2,2./3.)*(1.+SQR(t2)*ageFac1/9.);
  else
    {
      t2=exp(2.*t/ageFac2);
      t2=SQR(t2-1.)/(4.*t2);
      return pow(t2/ageFac1,1./3.);
    }
}

double dagedaofa(double a)
/*
  Finds d(age in years)/da at redshift z,
  for any cosmology (but no radiation), where a=1/(1+z).
*/
{
  return (9.7776e9/hubble)*sqrt(a/(omega_matter+a*
				   (omega_lambda*a*a+omega_curv)));
}

void InitSphCol(double z2, double M, double *dbari, double *Rvir, double *dcrit)
{
  double t1,t2,t3;

  SphCz1=1.e4*(1.+z2);
  SphCRi=RofM(M);
  SphCtf=ageofa(1./(1.+z2));
  SphCHi=(1.+SphCz1)/dagedaofa(1./(1.+SphCz1));
  t1=SphCtf/(Pi/SphCHi);
  t2=1./(1.+SphCz1);
  t3=1./(1.+z2);
  *dcrit = del_crit0(fomega(omega,omegal,z2),omegal);
  *dbari = (*dcrit)*ExactGrowth(SphCz1)/ExactGrowth(z2);
  *Rvir = SphCRi/((1.+z2)*pow(18.*SQR(Pi),1./3.));
  SphCk=1.e6;
  t1=1./MeanProf(SphCRi);
  SphCnorm=(*dbari)*t1;
  SphCz2=z2;
}

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)
{
  /* Given t=0 R_i in units of [R_i of the halo] and halo mass M [solar], returns Ri 
     [comoving Mpc], dbar_i, the shell mass Msh [same units as M], R_f [phys Mpc], 
     del-bar_f, delta_i, delta_f, pec V [km/s], total V [km/s]  */
  double t1,t2,t3,ystart[10],a,dbarim,dbarip,Rim,Rip,dbarfm,dbarfp,Rfp,Rfm,
    fracR;
  int nok,nbad,twosided=1;

  *Ri = SphCRi*relRi;
  fracR=0.034;
  Rip=*Ri *(1+fracR);
  Rim=*Ri *(1-fracR);
  if (Rim <= *Ri)
    {
      // One-sided numerical derivative:
      twosided=0;
      fracR=0.0063;
      Rip=*Ri *(1+fracR);
      Rim=*Ri;
    }
  t1=HelpProfile(SphCk,*Ri);
  t2=HelpProfile(SphCk,Rip);
  t3=HelpProfile(SphCk,Rim);
  *dbari=delfac*SphCnorm/t1;
  dbarip=delfac*SphCnorm/t2;
  dbarim=delfac*SphCnorm/t3;
  t1=(CBE(Rip)*MeanProf(Rip)-CBE(Rim)*MeanProf(Rim))/
    (CBE(Rip)-CBE(Rim));
  *deli = SphCnorm*t1;
  *Msh = M*CBE(relRi)*(1.+(*dbari));
  /* Variables: x = (physical R/Ri)/(a/a_i) and dx/d(H0 t) which starts 
     with pure growing mode */
  ystart[0]=1.;
  t3=1./(dagedaofa(1.));
  /*  ystart[1]=0.;  */
  ystart[1]=-(2./9.)*(*dbari)/(ageofa(1./(1.+SphCz1))*t3);
  Pass1=(1.+(*dbari));
  Pass3=t3;
  odeint(ystart,2,ageofa(1./(1.+SphCz1))*t3,ageofa(1./(1.+SphCz2))*t3,1.e-5,
	 (ageofa(1./(1.+SphCz2))-ageofa(1./(1.+SphCz1)))*t3/1.e4,0.,&nok,
	 &nbad,derivsSphCl,rkqs);
  *yf = ystart[0];
  ystart[0] *= (1+SphCz1)/(1+SphCz2);
  ystart[1] = ystart[1]*((1+SphCz1)/(1+SphCz2));
  *Rf = ystart[0] * (*Ri/(1+SphCz1));
  *dbarf = (1 + *dbari) * CBE(*Ri/(*Rf * (1+SphCz2))) - 1.;
  *Vpec = ystart[1] * 100.*hubble*(*Ri)/(1+SphCz1);
  a=1./(1.+SphCz2);
  *V = *Vpec + (*Rf) * 100.*hubble*
    sqrt(omega_matter/CBE(a)+omega_lambda+omega_curv/SQR(a));

  ystart[0]=1.;
  ystart[1]=-(2./9.)*dbarip/(ageofa(1./(1.+SphCz1))*t3);
  Pass1=(1. + dbarip);
  odeint(ystart,2,ageofa(1./(1.+SphCz1))*t3,ageofa(1./(1.+SphCz2))*t3,1.e-5,
	 (ageofa(1./(1.+SphCz2))-ageofa(1./(1.+SphCz1)))*t3/1.e4,0.,&nok,
	 &nbad,derivsSphCl,rkqs);
  ystart[0] *= (1+SphCz1)/(1+SphCz2);
  Rfp = ystart[0] * (Rip/(1+SphCz1));
  dbarfp = (1 + dbarip) * CBE(Rip/(Rfp * (1+SphCz2))) - 1.;
  if (twosided)
    {
      ystart[0]=1.;
      ystart[1]=-(2./9.)*dbarim/(ageofa(1./(1.+SphCz1))*t3);
      Pass1=(1. + dbarim);
      odeint(ystart,2,ageofa(1./(1.+SphCz1))*t3,ageofa(1./(1.+SphCz2))*t3,1.e-5,
	     (ageofa(1./(1.+SphCz2))-ageofa(1./(1.+SphCz1)))*t3/1.e4,0.,&nok,
	     &nbad,derivsSphCl,rkqs);
      ystart[0] *= (1+SphCz1)/(1+SphCz2);
      Rfm = ystart[0] * (Rim/(1+SphCz1));
      dbarfm = (1 + dbarim) * CBE(Rim/(Rfm * (1+SphCz2))) - 1.;
      *delf=(CBE(Rfp)*dbarfp-CBE(Rfm)*dbarfm)/
	(CBE(Rfp)-CBE(Rfm));
    }
  else
    {
      *delf=(CBE(Rfp)*dbarfp-CBE(*Rf)*(*dbarf))/
	(CBE(Rfp)-CBE(*Rf));
    }
}

void derivsSphCl(double H0t, double q[], double dq[])
{
  /*
    Solve for spherical collapse, including vacuum energy
    (flat universe). These are the derivatives of the
    two variables, which are: 
    x = (physical R/Ri)/(a/a_i) and dx/d(H0 t)

    Note that this is like Fortran, with q[1] and dq[1], 
    although input to odeint is ystart[0].                    */
  double a=aofage(H0t/Pass3),H;

  if (q[1] < 1.e-8) 
    dq[1]=dq[2]=0.;
  else
    {
      H=1./(dagedaofa(a)*a);
      dq[1]=q[2];
      dq[2]=.5*omega_matter*(q[1]-Pass1/SQR(q[1]))/CBE(a)-
	2.*(H/Pass3)*dq[1];
    }
}

double XiOfR1R2(double R1,double R2)
     /*
       Cross-correlation of (linearly extrapolated) overdensity 
       in co-centric spheres of radii R1 and R2 (com Mpc).

       Integration in several pieces ensures an accurate result.
     */
{
  double klimits[7],tol,klim,khi,klhi,kfac,klow,sOM,Rmin,Rmax;
  int i,nk;

  Rmin=DMIN(R1,R2);
  Rmax=DMAX(R1,R2);
  tol=2.e-6;
  klimits[0]=dmin2(1.e-5/Rmax,.2);
  klim=200./Rmax;
  khi=klim;
  klhi=log(khi);
  if (klhi > klimits[0]) klim=dmin2(klim,khi);
  nk=dround(log10(klim/klimits[0]));
  if (nk < 5) nk=5;
  else if (nk > 7) nk=7;
  kfac=log(klim/klimits[0])/(nk-1.);
  klimits[0]=log(klimits[0]);
  for (i=1;i<nk;i++)
    klimits[i]=klimits[i-1]+kfac;
  klow=1.e-9;
  PassXi1=R1;
  PassXi2=R2;
  PassXi4=0.;
  if (khi < klow) sOM=0.;
  else if (klhi < klimits[0]) sOM=integrate(xir2,klow,khi,tol);
  else 
    {
      sOM=integrate(xir2,klow,exp(klimits[0]),tol);
      for (i=1;i<nk;i++)
	sOM += integrate(xir,klimits[i-1],klimits[i],tol);
    }
  sOM=SQR(snorm)*sOM;
  return sOM;
}

double xir(double kl)
{
  return xirIntegrand(kl);
}

double xir2(double k)
{
  return xir(log(k))/k;
}

double GenMeanProf(double R, double x1, double *x2, double z)
{
  /* Generate average profile around a region of density
     sufficient to form a virialized halo at redshift z.
     R is the initial comoving radius [Mpc], x1 and x2
     give the log-spaced range of radii.                  
     Quantity used is not density, but mean enclosed
     density.                                         
     If InfSig != 0, use +/- 1-sigma profile.
     If profile becomes non-monotonic, stop and set *x2  */
  int i,imono=1;
  double y,t7,M,t4,t3,t6,nu,t1=0.,M0,t8;

  IntProf=IntN3;
  MProfR=R;
  M=MofR(R);
  M0=M;
  t7=sigmaOfM(M);
  t1=sqrt(XiOfR1R2(R,R));
  nu=wofz(z);
  for (i=0; (i<IntProf) && imono; i++)
    {
      y=R*(1.+exp(log(x1)+log((*x2)/x1)*(i/(IntProf-1.))));
      M=MofR(y);
      MProfx[i]=log(y/R-1.);
      t6=XiOfR1R2(R,y);
      t3=sigmaOfM(M);
      if (iInfall)
	{
	  t8=t6;
	  t6 /= SQR(t7);
	  t4=SQR(nu)*t6*(1.-t6)/(SQR(t3)-t6*t8);
	  if (InfSig < 0)
	    t1=NewInfallm1sig(t6,t4);
	  else if (InfSig > 0)
	    t1=NewInfallp1sig(t6,t4);
	  else
	    t1=NewInfall(t6,t4);
	  MProfy[i]=t1;
	}
      else
	MProfy[i]=t6/SQR(t7);
      if ((i>0) && (MProfy[i] >= MProfy[i-1]))
	{
	  imono=0;
	  *x2 = exp(MProfx[i-1]);
	  IntProf=i;
	}
    }
  spline(MProfx,MProfy,IntProf,1.e31,1.e31,MProfy2);
  return 1.;
}

double MeanProf(double R)
{
  double t1,R0;

  R0=MProfR*(exp(MProfx[0])+1.);
  if (R <= R0)    
    {
      t1 = MProfy[0];
      t1 = (R <= MProfR) ? 1. :
	t1 + (1.-t1)*(R-R0)/(MProfR-R0);
    }
  else
    {
      splint(MProfx,MProfy,MProfy2,IntProf,log(R/MProfR-1.),&t1);
    }
  return t1;
}

double MofR(double R)
{
  /* M [solar mass] contained in a sphere of radius R [com Mpc] */
  return (4./3.)*Pi*2.775453e11*omega_matter * SQR(hubble)*CBE(R);
}

double RofM(double M)
{
  double t1;
  /* Inverse of MofR */
  t1=pow(3.*M/(4.*Pi*2.775453e11*omega_matter * SQR(hubble)),1./3.);
  return t1;
}

double NewInfall(double S0oS, double nu2oS)
{
  // Mean infall profile of initial density
  // from Barkana (2004), where the inputs are alpha and beta.
  double nu2oS0=nu2oS/S0oS;
  return 1.-(1.-S0oS+1./nu2oS0)*erf(sqrt(nu2oS0*(1.-S0oS)/2.))-
    sqrt(2.*(1.-S0oS)/(Pi*nu2oS0))*exp(-nu2oS0*(1.-S0oS)/2.);
}

double NewInfallm1sig(double alpha, double beta)
{
  Pass1=alpha;
  Pass2=beta;
  return zbrent(FInfall,-10.,1.,1.e-5,.159);
}

double NewInfallp1sig(double alpha, double beta)
{
  Pass1=alpha;
  Pass2=beta;
  return zbrent(FInfall,-10.,1.,1.e-5,.841);
}

double FInfall(double gamma)
{
  double A,B,alp,bet;
  alp=Pass1;
  bet=Pass2;
  A=(gamma-alp)*sqrt(bet/(2.*alp*(1.-alp)));
  B=(gamma+alp-2.)*sqrt(bet/(2.*alp*(1.-alp)));
  return 1.+.5*(erf(A)+erf(B))+(exp(-A*A)-exp(-B*B))*
    sqrt(alp/(2*Pi*bet*(1.-alp)));
}

double HelpProfile(double k, double R)
{
  double t1;
  t1=1./MeanProf(R);
  return t1;
}

double ExactGrowth(double z)
{
  /* Exact growth factor (via integration in LCDM) */
  double t1,t2,a=1./(1.+z);
  if (fabs(1.-a) < 1.e-6)
    return 1.;
  else if ((fabs(omega_matter-1.)+fabs(omega_lambda)) < 1.e-10)
    return a;
  else
    {
      t1=sqrt(omega_lambda*CBE(a)+omega_matter+(1.-omega_lambda-omega_matter)*a)*
	(1./CBE(sqrt(a)))*integrate(dEGrowth,0.,a,1.e-6);
      t2=integrate(dEGrowth,0.,1.,1.e-6);
      return t1/t2;
    }
}

double dEGrowth(double a)
{
  // Helper function in ExactGrowth.
  double t1;
  t1=omega_lambda*CBE(a)+omega_matter+(1.-omega_lambda-omega_matter)*a;
  return CBE(sqrt(a/t1));
}

double FindMvir(double z2)
{
  double relRi,dbari;

  dbari=del_v(fomega(omega,omegal,z2),omegal)*
    ExactGrowth(SphCz1)/ExactGrowth(z2);
  /* Find initial comoving radius (in units of initial tophat radius) with 
     the right overdensity needed to correspond to the final virial radius: */
  relRi=zbrent(delbOfRi,1.,2.8,1.e-5,dbari);
  return CBE(relRi)*(1.+dbari);
}

double delbOfRi(double relRi)
{
  return SphCnorm/HelpProfile(SphCk,SphCRi*relRi);
}

double del_v(double om0, double omLam)
     /*
       The initial overdensity that corresponds to the final virial radius.
       Uses formula from Barkana (2004) infall paper
       (valid for om0 >= .05 in LCDM)
     */
{
  double t1,x=log(om0);

  t1=1.5927;
  if (om0 > (1.-1.e-5)) return t1;
  else return t1*(1.+(x/1.e3)*(3.9+x*(-.36-.09*x)));
}

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

double integrate(double (*func)(double),double a,double b,double tol)
     /*
       This calls an integration routine, and integrates the function
       f between a and b, to relative accuracy tol.
     */
{
  if (fabs(a-b) < fabs(a+b)*1.e-12)
    return 0.;
  else if (a > b)
    {
      printf("   %-11.7e  %-11.7e\n",a,b);
      crash("a > b in integration limits. Please fix that. Ending program.");
    }
  return rombint(func,a,b,tol);
}

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

double rombint(double (*func)(double), double a, double b, double tol)
     /* rombint returns the integral from a to b of func using Romberg integration.
	The method converges provided that f(x) is continuous in (a,b).
	tol indicates the desired relative accuracy in the integral. */
{
  int MAXJ=5; 
  int MAXITER=30,nint,i,j,k,jmax;
  double g[MAXJ+1],h,gmax,g0,fourj,g1,error;

  h=0.5*(b-a);
  gmax=h*(func(a)+func(b));
  g[0]=gmax;
  nint=1;
  error=1.e20;
  for (i=0; !((i > MAXITER) || ((i > 5) && (fabs(error) < tol))); i++)
    {
      /*     Calculate next trapezoidal rule approximation to integral. */
      g0=0.;
      for (k=1; k<=nint; k++)
	{
	  g0+=func(a+(k+k-1)*h);
	  /*	  printf("In rombint, %-15.8e\t%-15.8e\n",exp(a+(k+k-1)*h),g0);*/
	}
		   
      g0=0.5*g[0]+h*g0;
      h=0.5*h;
      nint*=2;
      jmax=min2(i,MAXJ);
      fourj=1.;
      for (j=1; j<=jmax; j++)
	/*     Use Richardson extrapolation. */
	{
	  fourj*=4.0;
	  g1=g0+(g0-g[j-1])/(fourj-1.0);
	  g[j-1]=g0;
	  g0=g1;
	}
      /*      if (fabs(g0) > tol) */
      if (fabs(g0) > 1.e-50) 
	error=1.0-gmax/g0;
      else 
	error=gmax; 
      gmax=g0;
      g[jmax]=g0;
    }
  if ((i > MAXITER) && (fabs(error) > tol))
    printf("rombint failed to converge; integral=%g, error=%g\n",
	   g0,error);
  return g0;
}

double xirIntegrand(double lkp)
{
  double k=exp(lkp),t1=k*PassXi4;
  t1 = (t1 < 1.e-2) ? k*(1.+SQR(t1)*(-1./6.+SQR(t1)/120.)) : sin(t1)/PassXi4;
  return Tcdm2EandH(k)*t1*window(k*PassXi1)*window(k*PassXi2)*
    pow(k,RUNNINGnS+1.+.5*RUNNINGdn*log(k/RUNNINGk0));
}

void nrerror(char error_text[])
{
	fprintf(stderr,"Run-time error...\n");
	fprintf(stderr,"%s\n",error_text);
	fprintf(stderr,"...now exiting to system...\n");
	exit(1);
}

int *ivector(long nl, long nh)
/* allocate an int vector with subscript range v[nl..nh] */
{
	int *v;

	v=(int *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(int)));
	if (!v) nrerror("allocation failure in ivector()");
	return v-nl+NR_END;
}

void free_ivector(int *v, long nl, long nh)
/* free an int vector allocated with ivector() */
{
	free((FREE_ARG) (v+nl-NR_END));
}

/* This routine normalizes the power spectrum using sigma 8 */
double dp()
{
  double sig8;       /* sigma 8 as calculated by the integral */

  sig8=sqrt(integrate(sig8Integrand,log(1.e-5),log(1.e4),1.e-4)/
	    (2.*SQR(PI)));
  return(sigma8/sig8);
}
 
/* This is the integrand which is used to calculate sigma 8 */
double sig8Integrand(double lkp)
{
  double k=exp(lkp);
  return Tcdm2EandH(k)*SQR(window(k*8./hubble))*k*SQR(k)*
    pow(k,RUNNINGnS-1.+.5*RUNNINGdn*log(k/RUNNINGk0));
}

int min2(int a, int b)
{
  return (a < b ? a : b);
}


/* This is a spherical top hat window function
   for use by a spherical Abe Lincoln.     */
double window(double x)
{
  double itsit,t;
  if(x>=0.01){
    itsit = 3.*(sin(x)-cos(x)*x);
    itsit = itsit/(x*x*x);
  }
  else{
    t=SQR(x); 
    itsit = 1.+t*(-.1+t/280.);
  }
  return(itsit);

}

/* This is a sharp k-space window function */
double windowk(double x)
{
  double itsit;
  if(x>=1.){
    itsit = 0.;
  }
  else{itsit = 1.;}
  return(itsit);

}

/* This is the derivative of the spherical top hat 
   window function with respect to x */
double windowprime(double x)
{
  double itsit,t=SQR(x);

  if(x>=0.01){
    itsit = 3.*cos(x)*x+sin(x)*(t-3.);
    itsit = itsit*3./SQR(t);
  }
  else{
    itsit = x*(-.2+t/70.);
  }
  return(itsit);
}

/* Run this once to initialize the transfer function 
   from Eisenstein and Hu (assuming no neutrinos)   */
void InitEandH()
{
  double obhh,alphan,b1,b2;
  double x,z_drag,y_drag,k_equality,R_drag,
    R_equality,z_equality;
  double fc,fcb,fb,fnub;  /* Fraction of Omega0 in cold dm, cdm+bary... */
  double pc,pcb;          /* Functions of the f's */

  omhh=OMEGA0*SQR(hubble);
  obhh=OMEGAB*SQR(hubble);
  b1 = 0.313*pow(omhh,-0.419)*(1.+0.607*pow(omhh,0.674));
  b2 = 0.238*pow(omhh,0.223);
  x = 1291.*pow(omhh,0.251)/(1+0.659*pow(omhh,0.828));
  z_drag = x * (1.+b1*pow(obhh,b2));
  
  z_equality = 24986.8*omhh/SQR(SQR(THETA27));     /* Actually 1+z_eq */
  y_drag = z_equality /(1+z_drag);

  fb  = OMEGAB/OMEGA0;
  fcb = 1.;
  fc  = fcb-fb;
  fnub = fb;
  pcb = 0.;
  pc  = 1.25 - sqrt(1.+24.*fc)/4.;
  x  = (fc/fcb)*(5.-2.*(pc+pcb))/(5.-4.*pcb);
  x  *= (1.+fnub*(-0.553+0.126*fnub*fnub))*pow((1.+y_drag),(pcb-pc));
  alphan = x*(1.+((pc-pcb)/2.)*(1.+1./((3.-4.*pc)*(7.-4.*pcb)))/(1.+y_drag));

  betac=1./(1.-0.949*fb);
  alpha_gamma = sqrt(alphan);
  k_equality = 0.0746*omhh/SQR(THETA27);
  R_drag = 31.5*obhh*1000./((1. + z_drag)*SQR(SQR(THETA27)));
  R_equality = 31.5*obhh*1000./((1. + z_equality)*SQR(SQR(THETA27)));
  sound_horizon_fit = (2./(3.*k_equality))*sqrt(6./R_equality)*
    log((sqrt(1.+R_drag)+sqrt(R_drag+R_equality) )
	/(1.+sqrt(R_equality)));
}

/* CDM power spectrum based on the transfer function 
   from Eisenstein and Hu (assuming no neutrinos) */
double Tcdm2EandH(double k)
{
  /* 
     Input k is the wavenumber in units of 1/Mpc   
  */
  double q;   /* The scaled wavenumber */
  double c,l; /* As defined by Eisenstein and Hu */
  double t;   /* The transfer function */
  double t2,gamma_eff;

  gamma_eff =omhh*(alpha_gamma+(1-alpha_gamma)/
		   (1+SQR(SQR(k*sound_horizon_fit*0.43))));
  q = k*SQR(THETA27)/gamma_eff;
  t2=pow(q,1.11);
  c = 14.4+(325./(1.+60.5*t2));
  l = log(exp(1.)+1.84*betac*alpha_gamma*q);
  t = 1./(1.+c*SQR(q)/l);
  return k*SQR(t*PI)/8.; 
}

double sigmaOfM(double M)
/*
  M is in solar mass. Returns the (linearly-extrapolated) density 
  fluctuation sigma(M).
*/
{
  double klimits[7],tol,klim,khi,klhi,kfac,klow,sOM,t1;
  int i,nk;

  /* scale is in comoving Mpc: */
  scale=pow(M*3./(4.*pi*2.775453e11*omhh),1./3.);
  t1=RofM(M);
  tol=2.e-6;
  klimits[0]=dmin2(1.e-5/scale,.2);
  klow=1.e-9;
  klim=200./scale;
  if (klimits[0] <= klow*1.01) klimits[0]=klow*1.01;
  if (klim <= klimits[0]*1.01) klim=klimits[0]*1.01;
  khi=klim;
  klhi=log(khi);
  nk=dround(log10(klim/klimits[0]));
  if (nk < 5) nk=5;
  else if (nk > 7) nk=7;
  kfac=log(klim/klimits[0])/(nk-1.);
  klimits[0]=log(klimits[0]);
  for (i=1;i<nk;i++)
    klimits[i]=klimits[i-1]+kfac;
  if (khi < klow) sOM=0.;
  else if (klhi < klimits[0]) sOM=integrate(sigmatop2,klow,khi,tol);
  else 
    {
      sOM=integrate(sigmatop2,klow,exp(klimits[0]),tol);
      for (i=1;i<nk;i++)
	sOM += integrate(sigmatop,klimits[i-1],klimits[i],tol);
    }
  sOM=snorm * sqrt(sOM);
  return sOM;
}

int dround(double x)
{
  return (int)(x < 0 ? (x-.5) : (x+.5));
}

double dmin2(double a, double b)
{
  return (a < b ? a : b);
}

double sigmatop(double kl)
     /*
       Specialized to z=0, no neutrinos or g-waves.
     */
{
  double k,x;
  double t1;

  k=exp(kl);
  x=scale*k;

  t1=Tcdm2EandH(k);
  return t1*pow(k,(2.+tilt+.5*tiltdn*log(k/tiltk0)))
    * SQR(3.*(x*cos(x) - sin(x))/(x*x*x));
}

double sigmatop2(double k)
{
  return sigmatop(log(k))/k;
}

double fomega(double om0,double omLam,double z)
{
  double f;

  f=omLam/CBE(1.+z);
  return om0/(om0+f+(1.-om0-omLam)/(1.+z));
}

double del_crit0(double om0, double omLam)
     /*
       Critical density for collapse (before linear extrapolation to z=0).
       Uses accurate formula from Barkana (2004)
       (valid in LCDM for om0 >= .05).
     */
{
  double t1,x=log(om0);

  t1=.15*pow(12.*pi,2./3.);
  if (om0 > (1.-1.e-5)) return t1;
  else if (omLam < 1.e-5) return t1*pow(om0,.0185);
  else return t1*(1.+(x/1.e3)*(4.86+x*(-.51-.113*x)));
}

double CDMdeltac(double z, double omega, double omegal)
     /*
       Critical density for collapse delta_c(z) in LCDM
     */
{
  double t1;
  t1=del_crit0(fomega(omega,omegal,z),omegal)/ExactGrowth(z);
  return t1;
}

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 []));
  int i;
  double errmax,h,xnew,*yerr,*ytemp;

  yerr=vector(1,n);
  ytemp=vector(1,n);
  h=htry;
  for (;;) {
    rkck(y,dydx,n,*x,h,ytemp,yerr,derivs);
    errmax=0.0;
    for (i=1;i<=n;i++) errmax=DMAX(errmax,fabs(yerr[i]/yscal[i]));
    errmax /= eps;
    if (errmax > 1.0) {
      h=SAFETY*h*pow(errmax,PSHRNK);
      if (h < 0.1*h) h *= 0.1;
      xnew=(*x)+h;
      if (xnew == *x) nrerror("stepsize underflow in rkqs");
      continue;
    } else {
      if (errmax > ERRCON) *hnext=SAFETY*h*pow(errmax,PGROW);
      else *hnext=5.0*h;
      *x += (*hdid=h);
      for (i=1;i<=n;i++) y[i]=ytemp[i];
      break;
    }
  }
  free_vector(ytemp,1,n);
  free_vector(yerr,1,n);
}

void rkck(double y[], double dydx[], int n, double x, double h, double yout[],
  double yerr[], void (*derivs)(double, double [], double []))
{
  int i;
  static double a2=0.2,a3=0.3,a4=0.6,a5=1.0,a6=0.875,b21=0.2,
    b31=3.0/40.0,b32=9.0/40.0,b41=0.3,b42 = -0.9,b43=1.2,
    b51 = -11.0/54.0, b52=2.5,b53 = -70.0/27.0,b54=35.0/27.0,
    b61=1631.0/55296.0,b62=175.0/512.0,b63=575.0/13824.0,
    b64=44275.0/110592.0,b65=253.0/4096.0,c1=37.0/378.0,
    c3=250.0/621.0,c4=125.0/594.0,c6=512.0/1771.0,
    dc5=-277.0/14336.0;
  double dc1=c1-2825.0/27648.0,dc3=c3-18575.0/48384.0,
    dc4=c4-13525.0/55296.0,dc6=c6-0.25;
  double *ak2,*ak3,*ak4,*ak5,*ak6,*ytemp;

  ak2=vector(1,n);
  ak3=vector(1,n);
  ak4=vector(1,n);
  ak5=vector(1,n);
  ak6=vector(1,n);
  ytemp=vector(1,n);
  for (i=1;i<=n;i++)
    ytemp[i]=y[i]+b21*h*dydx[i];
  (*derivs)(x+a2*h,ytemp,ak2);
  for (i=1;i<=n;i++)
    ytemp[i]=y[i]+h*(b31*dydx[i]+b32*ak2[i]);
  (*derivs)(x+a3*h,ytemp,ak3);
  for (i=1;i<=n;i++)
    ytemp[i]=y[i]+h*(b41*dydx[i]+b42*ak2[i]+b43*ak3[i]);
  (*derivs)(x+a4*h,ytemp,ak4);
  for (i=1;i<=n;i++)
    ytemp[i]=y[i]+h*(b51*dydx[i]+b52*ak2[i]+b53*ak3[i]+b54*ak4[i]);
  (*derivs)(x+a5*h,ytemp,ak5);
  for (i=1;i<=n;i++)
    ytemp[i]=y[i]+h*(b61*dydx[i]+b62*ak2[i]+b63*ak3[i]+b64*ak4[i]+b65*ak5[i]);
  (*derivs)(x+a6*h,ytemp,ak6);
  for (i=1;i<=n;i++)
    yout[i]=y[i]+h*(c1*dydx[i]+c3*ak3[i]+c4*ak4[i]+c6*ak6[i]);
  for (i=1;i<=n;i++)
    yerr[i]=h*(dc1*dydx[i]+dc3*ak3[i]+dc4*ak4[i]+dc5*ak5[i]+dc6*ak6[i]);
  free_vector(ytemp,1,n);
  free_vector(ak6,1,n);
  free_vector(ak5,1,n);
  free_vector(ak4,1,n);
  free_vector(ak3,1,n);
  free_vector(ak2,1,n);
}

static int kmax,kount;
static double *xp,**yp,dxsav;	      /* Global variables needed
					 for odeint, but not
					 actually used */

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 [])))
{
  int nstp,i;
  double xsav,x,hnext,hdid,h;
  double *yscal,*y,*dydx;

  yscal=vector(1,nvar);
  y=vector(1,nvar);
  dydx=vector(1,nvar);
  x=x1;
  h=(x2 > x1) ? fabs(h1) : -fabs(h1);
  *nok = (*nbad) = kount = 0;
  for (i=1;i<=nvar;i++) y[i]=ystart[i-1];
  if (kmax > 0) xsav=x-dxsav*2.0;
  for (nstp=1;nstp<=MAXSTP;nstp++) {
    (*derivs)(x,y,dydx);
    for (i=1;i<=nvar;i++)
      yscal[i]=fabs(y[i])+fabs(dydx[i]*h)+TINY;
    if (kmax > 0 && kount < kmax-1 && fabs(x-xsav) > fabs(dxsav)) {
      xp[++kount]=x;
      for (i=1;i<=nvar;i++) yp[i][kount]=y[i];
      xsav=x;
    }
    if ((x+h-x2)*(x+h-x1) > 0.0) h=x2-x;
    (*rkqs)(y,dydx,nvar,&x,h,eps,yscal,&hdid,&hnext,derivs);
    if (hdid == h) ++(*nok); else ++(*nbad);
    if ((x-x2)*(x2-x1) >= 0.0) {
      for (i=1;i<=nvar;i++) ystart[i-1]=y[i];
      if (kmax) {
        xp[++kount]=x;
        for (i=1;i<=nvar;i++) yp[i][kount]=y[i];
      }
      free_vector(dydx,1,nvar);
      free_vector(y,1,nvar);
      free_vector(yscal,1,nvar);
      return;
    }
    if (fabs(hnext) <= hmin) nrerror("Step size too small in odeint");
    h=hnext;
  }
  nrerror("Too many steps in routine odeint");
}

double *vector(long nl, long nh)
/* allocate a double vector with subscript range v[nl..nh] */
{
  double *v;

  v=(double *)malloc((size_t) ((nh-nl+1+NR_END)*sizeof(double)));
  if (!v) nrerror("allocation failure in vector()");
  return v-nl+NR_END;
}

void free_vector(double *v, long nl, long nh)
/* free a double vector allocated with vector() */
{
  free((FREE_ARG) (v+nl-NR_END));
}

void spline(double x[], double y[], int n, double yp1, double ypn, double y2[])
{
  int i,k;
  double p,qn,sig,un,u[1000];

  if (n < 2) return;
  if (yp1 > 0.99e30)
    y2[0]=u[0]=0.0;
  else {
    y2[0] = -0.5;
    u[0]=(3.0/(x[1]-x[0]))*((y[1]-y[0])/(x[1]-x[0])-yp1);
  }
  for (i=1;i<=n-2;i++) {
    sig=(x[i]-x[i-1])/(x[i+1]-x[i-1]);
    p=sig*y2[i-1]+2.0;
    y2[i]=(sig-1.0)/p;
    u[i]=(y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]);
    u[i]=(6.0*u[i]/(x[i+1]-x[i-1])-sig*u[i-1])/p;
  }
  if (ypn > 0.99e30)
    qn=un=0.0;
  else {
    qn=0.5;
    un=(3.0/(x[n-1]-x[n-2]))*(ypn-(y[n-1]-y[n-2])/(x[n-1]-x[n-2]));
  }
  y2[n-1]=(un-qn*u[n-2])/(qn*y2[n-2]+1.0);
  for (k=n-2;k>=0;k--)
    y2[k]=y2[k]*y2[k+1]+u[k];
}

void splint(double xa[], double ya[], double y2a[], int n, double x, double *y)
{
  // Note: extrapolates only linearly (not with a cubic)
  void nrerror(char error_text[]);
  int klo,khi,k;
  double h,b,a,da,db,t1,t2;

  if (n < 2) 
    {
      *y = ya[0];
      return;
    }
  else if (x <= xa[0])
    {
      linint(xa,ya,2,x,y,&t1,&t2,&k);
      return;
    }
  else if (x >= xa[n-1])
    {
      linint(xa+n-2,ya+n-2,2,x,y,&t1,&t2,&k);
      return;
    }
  klo=0;
  khi=n-1;
  while (khi-klo > 1) {
    k=(khi+klo)/2;
    if (xa[k] > x) khi=k;
    else klo=k;
  }
  h=xa[khi]-xa[klo];
  if (h == 0.) 
    {
      printf("%g %g %g\n",xa[0],xa[n-1],x);
      crash("bad xa input in splint");
    }
  a=(xa[khi]-x)/h;
  b=(x-xa[klo])/h;
  da=-1./h;
  db=1./h;
  *y=a*ya[klo]+b*ya[khi]+((a*a*a-a)*y2a[klo]+(b*b*b-b)*y2a[khi])*(h*h)/6.0;
}

double zbrent(double (*func)(double), double x1, double x2, double tol,
	      double value)
     /*
       Search for a root of the equation func=value
     */
{
  void nrerror(char error_text[]);
  int iter;
  double a=x1,b=x2,c=x2,d,e,min1,min2;
  double fa=(*func)(a)-value,fb=(*func)(b)-value,fc,p,q,r,s,tol1,xm;

  if ((fa > 0.0 && fb > 0.0) || (fa < 0.0 && fb < 0.0))
    {
      printf("a, b, fa, fb, value= %-15.8e\t%-15.8e\t%-15.8e\t%-15.8e\t%-15.8e\n",
	     a,b,fa,fb,value);
      nrerror("Root must be bracketed in zbrent");  
    }
  fc=fb;
  for (iter=1;iter<=ITMAX;iter++) {
    if ((fb > 0.0 && fc > 0.0) || (fb < 0.0 && fc < 0.0)) {
      c=a;
      fc=fa;
      e=d=b-a;
    }
    if (fabs(fc) < fabs(fb)) {
      a=b;
      b=c;
      c=a;
      fa=fb;
      fb=fc;
      fc=fa;
    }
    tol1=DMIN(2.0*EPS,tol)*fabs(b)+0.5*tol;
    xm=0.5*(c-b);
    /* Modified exit condition: */
    if (fabs(xm) <= tol1 || fabs(fb) < EPS) 
      return b;
    if (fabs(e) >= tol1 && fabs(fa) > fabs(fb)) {
      s=fb/fa;
      if (a == c) {
        p=2.0*xm*s;
        q=1.0-s;
      } else {
        q=fa/fc;
        r=fb/fc;
        p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0));
        q=(q-1.0)*(r-1.0)*(s-1.0);
      }
      if (p > 0.0) q = -q;
      p=fabs(p);
      min1=3.0*xm*q-fabs(tol1*q);
      min2=fabs(e*q);
      if (2.0*p < (min1 < min2 ? min1 : min2)) {
        e=d;
        d=p/q;
      } else {
        d=xm;
        e=d;
      }
    } else {
      d=xm;
      e=d;
    }
    a=b;
    fa=fb;
    if (fabs(d) > tol1)
      b += d;
    else
      b += (xm > 0.0 ? fabs(tol1) : -fabs(tol1));
    fb=(*func)(b)-value;
  }
  nrerror("Maximum number of iterations exceeded in zbrent\n");
  return 0.0;
}

void linint(double *xa,double *ya,int n,double x,double *y, double *ra, 
	    double *rb, int *kLow)
     /*
       Do linear interpolation (and extrapolation) instead of cubic spline.
       Can deal with < or > ordering.
     */
{
  int k,khi=n-1,klo=0;
  double a,b,h;
  
  if (xa[khi] >= xa[klo])
    {
      while (khi-klo > 1) 
	{
	  k=(khi+klo)/2;
	  if (xa[k] > x) khi=k;
	  else klo=k;
	}
    }
  else
    {
      while (khi-klo > 1) 
	{
	  k=(khi+klo)/2;
	  if (xa[k] < x) khi=k;
	  else klo=k;
	}
    }
  h=xa[khi]-xa[klo];
  if (h == 0.) 
    {
      printf("%i %g %g %g\n",n,xa[0],xa[n-1],x);
      crash("bad xa input in linint");
    }
  a=(xa[khi]-x)/h;
  b=(x-xa[klo])/h;
  *ra = a;
  *rb = b;
  *y=a*ya[klo]+b*ya[khi];
  *kLow = klo;
}