source: Sophya/trunk/Cosmo/SimLSS/pkspectrum.cc@ 3315

#include "sopnamsp.h"
#include "machdefs.h"
#include <iostream>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <unistd.h>

#include "pexceptions.h"

#include "constcosmo.h"
#include "geneutils.h"
#include "pkspectrum.h"


///////////////////////////////////////////////////////////
//******************** InitialSpectrum ******************//
///////////////////////////////////////////////////////////

InitialSpectrum::InitialSpectrum(double n,double a)
  : n_(n), A_(a)
{
}

InitialSpectrum::InitialSpectrum(InitialSpectrum& pkinf)
  : n_(pkinf.n_), A_(pkinf.A_)
{
}

InitialSpectrum::~InitialSpectrum(void)
{
}

void InitialSpectrum::SetNorm(double a)
{
  A_ = a;
}

void InitialSpectrum::SetSlope(double n)
{
  n_ = n;
}


///////////////////////////////////////////////////////////
//****************** TransfertEisenstein ****************//
///////////////////////////////////////////////////////////

// From Eisenstein & Hu ApJ 496:605-614 1998 April 1 (ou astro-ph/9709112)
TransfertEisenstein::TransfertEisenstein(double h100,double OmegaCDM0,double OmegaBaryon0,double tcmb,bool nobaryon,int lp)
  : lp_(lp)
  , Oc_(OmegaCDM0) , Ob_(OmegaBaryon0) , h_(h100) , tcmb_(tcmb)
  , nobaryon_(nobaryon) , nooscenv_(0), retpart_(0)
{
 zero_();
 Init_();
}

TransfertEisenstein::TransfertEisenstein(TransfertEisenstein& tf)
  : lp_(tf.lp_)
  ,Oc_(tf.Oc_) , Ob_(tf.Ob_) , h_(tf.h_) , tcmb_(tf.tcmb_)
  , nobaryon_(tf.nobaryon_) , nooscenv_(tf.nooscenv_), retpart_(tf.retpart_)
{
 zero_();
 Init_();
}

void TransfertEisenstein::Init_(void)
{

  O0_ = Oc_ + Ob_;
  if(nobaryon_) {O0_ = Oc_; Ob_ = 0.;}
  if(lp_) cout<<"h100="<<h_<<" Omatter="<<O0_<<" Ocdm="<<Oc_<<" Ob="<<Ob_<<endl;

  double H0 = 100. * h_, h2 = h_*h_;

  if(tcmb_<0.) tcmb_ = T_CMB_Par;
  th2p7_ = tcmb_/2.7;
  double th2p7P4 = th2p7_*th2p7_*th2p7_*th2p7_;
  if(lp_) cout<<"tcmb = "<<tcmb_<<" K = "<<th2p7_<<" *2.7K "<<endl;

  // Formule 2 p 606
  zeq_ = 2.50e4 * O0_ * h2 / th2p7P4;
  if(lp_) cout<<"zeq = "<<zeq_<<" (redshift of matter-radiation equality)"<<endl;

  // Formule 3 p 607
  // (attention ici C=1 : H0 -> H0/C si on utilise la premiere formule)
  //  keq_ = sqrt(2.*O0_*H0*H0*zeq_) / SpeedOfLight_Cst;
  keq_ = 7.46e-2 * O0_ * h2 / (th2p7_*th2p7_);
  if(lp_) cout<<"keq = "<<keq_<<" Mpc^-1 (scale of equality)"<<endl;

  // On s'arrete ici si pas de baryons
  if(nobaryon_) return;

  // Formule 4 p 607
  double b1_eq4 = 0.313*pow(O0_*h2,-0.419)*(1. + 0.607*pow(O0_*h2,0.674));
  double b2_eq4 = 0.238*pow(O0_*h2,0.223);
  zd_ = 1291. * pow(O0_*h2,0.251) / (1.+0.659* pow(O0_*h2,0.828))
              * (1. + b1_eq4*pow(Ob_*h2,b2_eq4));
  if(lp_) cout<<"zd = "<<zd_<<" (Redshift of drag epoch)"<<endl;

  // Formule 5 page 607    (R = 3*rho_baryon/4*rho_gamma)
  Req_ = 31.5*Ob_*h2 / th2p7P4 * (1.e3/zeq_);
  //WARNING: W.Hu code (tf_fit.c) en des-accord avec l'article: zd -> (1+zd)
  Rd_  = 31.5*Ob_*h2 / th2p7P4 * (1.e3/zd_);
  //in tf_fit.c: Rd_  = 31.5*Ob_*h2 / th2p7P4 * (1.e3/(1.+zd_));
  if(lp_) {
    cout<<"Req = "<<Req_<<" Rd = "<<Rd_
        <<" (Photon-baryon ratio at equality/drag epoch)"<<endl;
    cout<<"Sound speed at equality "<<1./sqrt(3.*(1.+Req_))
        <<", at drag "<<1./sqrt(3.*(1.+Rd_))<<" in unit of C"<<endl;
  }

  // Formule 6 p 607
  s_ = 2./(3.*keq_) * sqrt(6./Req_)
      * log( (sqrt(1.+Rd_) + sqrt(Rd_+Req_)) / (1.+sqrt(Req_)) );
  if(lp_) cout<<"s = "<<s_<<" Mpc (sound horizon at drag epoch)"<<endl;

  // Formule 7 page 607
  ksilk_ = 1.6*pow(Ob_*h2,0.52)*pow(O0_*h2,0.73) * (1. + pow(10.4*O0_*h2,-0.95));
  if(lp_) cout<<"ksilk = "<<ksilk_<<" Mpc^-1 (silk damping scale)"<<endl;

  // Formules 10 page 608
  double a1 = pow(46.9*O0_*h2,0.670) * (1. + pow(32.1*O0_*h2,-0.532));
  double a2 = pow(12.0*O0_*h2,0.424) * (1. + pow(45.0*O0_*h2,-0.582));
  alphac_ = pow(a1,-Ob_/O0_) * pow(a2,-pow(Ob_/O0_,3.));
  double b1 = 0.944 / (1. + pow(458.*O0_*h2,-0.708));
  double b2 = pow(0.395*O0_*h2,-0.0266);
  betac_ = 1 / ( 1. + b1*(pow(Oc_/O0_,b2) - 1.) );
  if(lp_) cout<<"alphac = "<<alphac_<<" betac = "<<betac_
               <<" (CDM suppression/log shift)"<<endl;

  // Formule 23 page 610
  bnode_ = 8.41 * pow(O0_*h2,0.435);
  if(lp_) cout<<"bnode = "<<bnode_<<" (sound horizon shift)"<<endl;

  // Formule 14 page 608
  //WARNING: W.Hu code (tf_fit.c) en des-accord avec l'article: (1+zeq) -> zeq
  double y = (1.+zeq_)/(1.+zd_);
  //in tf_fit.c: double y = zeq_/(1.+zd_);
  double s1py = sqrt(1.+y);
  double Gy = y*( -6.*s1py + (2.+3.*y)*log((s1py+1.)/(s1py-1.)) );
  alphab_ = 2.07*keq_*s_*pow(1.+Rd_,-3./4.)*Gy;

  // Formule 24 page 610
  betab_ = 0.5 + Ob_/O0_
         + (3.-2.*Ob_/O0_) * sqrt(pow(17.2*O0_*h2,2.) + 1.);
  if(lp_) cout<<"alphab = "<<alphab_<<" betab = "<<betab_
               <<" (Baryon suppression/envelope shift)"<<endl;

  // Formule 31 page 612
  alphag_ = 1.
          - 0.328*log(431.*O0_*h2)*Ob_/O0_
          + 0.38*log(22.3*O0_*h2)*pow(Ob_/O0_,2.);
  if(lp_) cout<<"alphag = "<<alphag_<<" (gamma suppression in approximate TF)"<<endl;

 // The approximate value of the sound horizon, formule 26 page 611
 sfit_ = 44.5*log(9.83/(O0_*h2)) / sqrt(1.+10.*pow(Ob_*h2,3./4.));  // Mpc
 if(lp_) cout<<"sfit="<<sfit_<<" Mpc (fit to sound horizon)"<<endl;

 // La positoin du premier pic acoustique, formule 25 page 611
 kpeak_ = 5*M_PI/(2.*sfit_) * (1.+0.217*O0_*h2);  // 1/Mpc
 if(lp_) cout<<"kpeak="<<kpeak_<<" Mpc^-1 (fit to wavenumber of first peak)"<<endl;

  return;
}

TransfertEisenstein::~TransfertEisenstein(void)
{
}

void TransfertEisenstein::SetNoOscEnv(unsigned short nooscenv)
// To obtain an approximate form of the non-oscillatory part of the transfert function
// nooscenv = 0 : use the baryon oscillatory part of transfert function (full tf)
// nooscenv = 1 : use approx. paragraph 3.3 p610 (middle of right column)
//                Replace  j0(k*stilde)  ->  [1+(k*stilde)^4]^(-1/4)
// nooscenv = 2 : use formulae 29+30+31 page 612
//                The value of an approximate transfer function that captures
//                the non-oscillatory part of a partial baryon transfer function.
//                In other words, the baryon oscillations are left out,
//                but the suppression of power below the sound horizon is included.
{
 if(nooscenv!=1 && nooscenv!=2) nooscenv = 0;
 nooscenv_ = nooscenv;
}

void TransfertEisenstein::SetReturnPart(unsigned short retpart)
// To return only baryon or CDM part part of transfert function
// retpart = 1 : return only CDM part of transfert function
// retpart = 2 : return only Baryon part of transfert function
// retpart = anything else: return only full transfert function
// WARNING: only relevant for nobaryon_=false AND nooscenv!=2
{
 if(retpart!=1 && retpart!=2) retpart = 0;
 retpart_ = retpart;
}

double TransfertEisenstein::T0tild(double k,double alphac,double betac)
{
  // Formule 10 p 608
  //double q = k*th2p7_*th2p7_/(O0_*h_*h_);
  double q = k/(13.41*keq_);
  // Formule 20 p 610
  double C = (14.2/alphac) + 386./(1.+69.9*pow(q,1.08));
  // Formule 19 p 610
  double x = log(M_E+1.8*betac*q);
  return x / (x + C*q*q);
}

double TransfertEisenstein::operator() (double k)
{

  // --- Pour zero baryon
  //  OU Pour function lissee sans oscillation baryon
  if(nobaryon_  || nooscenv_ == 2) {
    double gamma = O0_*h_;
    // Calcul de Gamma_eff, formule 30 page 612 (pour fct lissee)
    if( nobaryon_==false && nooscenv_ == 2 )
      gamma = O0_*h_*(alphag_ + (1.-alphag_)/(1.+pow(0.43*k*sfit_,4.))); // Gamma_eff
    // Formule 28 page 612 : qui est est equivalent a:
    //    q = k / h_ * th2p7_*th2p7_ / gamma;
    // qui est est equivalent a:
    //    q = k / (13.41 * keq)                   pour Ob=0
    //    q = k / (13.41 * keq) * (O0*h/Gamma)    pour le spectre lisse
    // Les resultats sont legerement differents a cause des valeurs approx.
    // des constantes numeriques: on prend comme W.Hu (tf_fit.c)
    //double q = k / h_ * th2p7_*th2p7_ / gamma;  // Mpc^-1
    double q = k/(13.41*keq_) * (O0_*h_/gamma);  // Mpc^-1
    // Formules 29 page 612
    double l0 = log(2.*M_E + 1.8*q);
    double c0 = 14.2 + 731./(1.+62.5*q);
    return l0 / (l0 + c0*q*q);
  }

  // --- Pour CDM + Baryons
  // --- CDM
  double f = 1. / (1. + pow(k*s_/5.4,4.));
  double Tc = f*T0tild(k,1.,betac_) + (1.-f)*T0tild(k,alphac_,betac_);
  if(retpart_ == 1) return Tc; 

  // --- Baryons
  // Formule 22 page 610
  double stilde, ksbnode = k*s_/bnode_;
  if(ksbnode<0.001) stilde =s_ * ksbnode;
     else   stilde = s_ / pow(1. + pow(1./ksbnode,3.), 1./3.);
  // Formule 21 page 610
  double j0kst = 0.;
  if(nooscenv_ == 1) {
    j0kst = pow(1.+pow(k*stilde,4.) , -1./4.); //lissee sans oscillation baryon
  } else {
    double x = k*stilde;
    if(x<0.01) j0kst = 1. - x*x/6.*(1.-x*x/20.);
      else j0kst = sin(x)/x;
    //cout<<"DEBUG: k="<<k<<" stilde="<<stilde<<" x="<<x<<" j0kst="<<j0kst<<endl;
  }
  double Tb = T0tild(k,1.,1.) / (1. + pow(k*s_/5.2,2.));
  Tb += alphab_/(1.+pow(betab_/(k*s_),3.)) * exp(-pow(k/ksilk_,1.4));
  Tb *= j0kst;
  if(retpart_ == 2) return Tb;

  // --- Total
  double T = (Ob_/O0_)*Tb + (Oc_/O0_)*Tc;

  return T;
}

double TransfertEisenstein::KPeak(void)
// Position du premier pic acoustic
{
 if(nobaryon_) return -1.;
 return kpeak_;
}


///////////////////////////////////////////////////////////
//********************* GrowthFactor ********************//
///////////////////////////////////////////////////////////

// From Eisenstein & Hu ApJ 496:605-614 1998 April 1
// Pour avoir D(z) = 1/(1+z) faire: OmegaMatter0=1 OmegaLambda0=0
GrowthFactor::GrowthFactor(double OmegaMatter0,double OmegaLambda0)
  : O0_(OmegaMatter0) , Ol_(OmegaLambda0) , Ok_(1.-OmegaMatter0-OmegaLambda0)
{
  if(OmegaMatter0==0.) {
    cout<<"GrowthFactor::GrowthFactor: Error bad OmegaMatter0  value : "<<OmegaMatter0<<endl;
    throw ParmError("GrowthFactor::GrowthFactor: Error badOmegaMatter0  value");
  }
  norm_ = 1.; // puisque (*this)(0.) a besoin de norm_
  norm_ = (*this)(0.);
  cout<<"GrowthFactor::GrowthFactor : norm="<<norm_<<endl;
}

GrowthFactor::GrowthFactor(GrowthFactor& d1)
  : O0_(d1.O0_) , Ol_(d1.Ol_) , Ok_(d1.Ok_) , norm_(d1.norm_)
{
}

GrowthFactor::~GrowthFactor(void)
{
}

double GrowthFactor::operator() (double z)
// see Formulae A4 + A5 + A6 page 614
{
 z += 1.;
 double z2 = z*z, z3 = z2*z;
 double den = Ol_ + Ok_*z2 + O0_*z3;
 double o0z = O0_ *z3 / den;
 double olz = Ol_ / den;

 // 4./7. = 0.571429
 double D1z = pow(o0z,0.571429) - olz + (1.+o0z/2.)*(1.+olz/70.);
 D1z = 2.5*o0z / z / D1z;

 return D1z / norm_;
}


///////////////////////////////////////////////////////////
//************** PkSpectrum0 et PkSpectrumZ *************//
///////////////////////////////////////////////////////////

PkSpectrum0::PkSpectrum0(InitialSpectrum& pkinf,TransfertEisenstein& tf)
  : pkinf_(pkinf) , tf_(tf)
{
}

PkSpectrum0::PkSpectrum0(PkSpectrum0& pk0)
  : pkinf_(pk0.pkinf_) , tf_(pk0.tf_)
{
}

PkSpectrum0::~PkSpectrum0(void)
{
}

double PkSpectrum0::operator() (double k)
{
  double tf = tf_(k);
  double pkinf = pkinf_(k);
  return pkinf *tf*tf;
}

//------------------------------------
PkSpectrumZ::PkSpectrumZ(PkSpectrum0& pk0,GrowthFactor& d1,double zref)
  : pk0_(pk0) , d1_(d1) , zref_(zref) , scale_(1.) , typspec_(0)
  , zold_(-1.) , d1old_(1.)
{
}

PkSpectrumZ::PkSpectrumZ(PkSpectrumZ& pkz)
  : pk0_(pkz.pk0_) , d1_(pkz.d1_) , zref_(pkz.zref_) , scale_(pkz.scale_) , typspec_(0)
  , zold_(pkz.zold_) , d1old_(pkz.d1old_)
{
}

PkSpectrumZ::~PkSpectrumZ(void) 
{
}

void PkSpectrumZ::SetTypSpec(unsigned short typspec)
  // typsec = 0 : compute Pk(k)
  //        = 1 : compute Delta^2(k) = k^3*Pk(k)/2Pi^2
{
  if(typspec>1) {
    cout<<"PkSpectrumZ::SetTypSpec: Error bad typspec value: "<<typspec<<endl;
    throw ParmError("PkSpectrumZ::SetTypSpec: Error bad typspec value");
  }
  typspec_ = typspec;
}

double PkSpectrumZ::operator() (double k)
{
  return (*this)(k,zref_);
}

double PkSpectrumZ::operator() (double k,double z)
{
  double d1;
  if(z == zold_) d1 = d1old_;
    else {d1 = d1old_ = d1_(z); zold_ = z;}

  double v = pk0_(k) * d1*d1;
  if(typspec_) v *= k*k*k/(2.*M_PI*M_PI);

  return scale_ * v;
}



///////////////////////////////////////////////////////////
//******************* VarianceSpectrum ******************//
///////////////////////////////////////////////////////////

VarianceSpectrum::VarianceSpectrum(GenericFunc& pk,unsigned short typfilter=0)
  : pk_(pk) , R_(0.)
{
  SetFilter(typfilter);
}

VarianceSpectrum::VarianceSpectrum(VarianceSpectrum& vpk)
  : pk_(vpk.pk_) , R_(vpk.R_)
{
  SetFilter(vpk.typfilter_);
}

VarianceSpectrum::~VarianceSpectrum(void)
{
}

//------------------------------------
void VarianceSpectrum::SetFilter(unsigned short typfilter)
// typfilter = 0 : spherical 3D top-hat
//           = 1 : spherical 3D gaussian
//           = 2 : no filter juste integrate spectrum)
{
  if(typfilter>2) {
    cout<<"VarianceSpectrum::SetFilter: Error bad value for type of filter"<<endl;
    throw ParmError("VarianceSpectrum::SetFilter: Error bad value for type of filter");
  }
  typfilter_ = typfilter;
}

void VarianceSpectrum::SetInteg(double dperc,double dlogkinc,double dlogkmax,unsigned short glorder)
// ATTENTION: on n'integre pas f(k)*dk mais k*f(k)*d(log10(k))
// see argument details in function IntegrateFuncLog (geneutils.cc)
{
  dperc_ = dperc;  if(dperc_<=0.) dperc_ = 0.1;
  dlogkinc_ = dlogkinc;
  dlogkmax_ = dlogkmax;
  glorder_ = glorder;
}


//------------------------------------
double VarianceSpectrum::Filter2(double x)
// ATTENTION: c'est le filtre au carre qui est renvoye
{
  // Just integrate the spectrum without filtering
  if(typfilter_ == 2) return 1.;

  double x2 = x*x;
  // Filtre gaussien G(x) = exp(-x^2/2)
  //        remarque G(x)^2 = exp(-x^2)
  // on prend le DL de G(x)^2 pour x->0 a l'ordre O(x^6)
  //             DL(x) = 1-x^2*(1-x^2/2)
  //             pour x<0.01  |DL(x)-G(X)^2|<2.0e-13
  if(typfilter_ == 1)
    if(x<0.01) return 1.-x2*(1.-x2/2.); else return exp(-x2);

  // Filtre top-hat T(x) = 3*(sin(x)-x*cos(x))/x^3
  // --- Gestion de la pseudo-divergence pour x->0
  // on prend le DL de T(x)^2 pour x->0 a l'ordre O(x^7)
  //             DL(x) = 1-x^2/5*(1-3*x^2/35*(1-4*x^2/81))
  //             pour x<0.1  |DL(x)-T(X)^2|<2.5e-13
  double f2=0.;
  if(x<0.1) {
    f2 = 1.-x2/5.*(1.-3.*x2/35.*(1.-4.*x2/81.));
  } else {
    f2 = 3.*(sin(x)-x*cos(x))/(x2*x);
    f2 *= f2;
  }
  return f2;
  
}

double VarianceSpectrum::Variance(double R,double kmin,double kmax)
// Compute variance of spectrum pk_ by integration
// Input:
//     R = taille du filter top-hat ou gaussien
//     kmin,kmax = bornes en k de l'integrale pour calculer la variance
// Return:
//     valeur de la variance (sigma^2)
// Remarque:
//   la meilleure approximation du filtre top-hat (R) est un filtre gaussien avec (Rg=R/sqrt(5))
//   la variance renvoyee est la variance de la masse
{
  if(R<=0. || kmin<=0 || kmax<=0. || kmin>=kmax) {
    cout<<"VarianceSpectrum::Variance: Error R<=0 or kmin<=0 or kmax<=0 or kmin>=kmax"<<endl;
    throw ParmError("VarianceSpectrum::Variance: Error R<=0 or kmin<=0 or kmax<=0 or kmin>=kmax");
  }

  R_ = R;
  double lkmin = log10(kmin), lkmax = log10(kmax);

  double var = IntegrateFuncLog(*this,lkmin,lkmax,dperc_,dlogkinc_,dlogkmax_,glorder_);

  return var;
}

//------------------------------------
double VarianceSpectrum::FindMaximum(double R,double kmin,double kmax,double eps)
// Retourne le maximum de la fonction a integrer
// La recherche a lieu entre [kmin,kmax] par pas logarithmiques
// Input:
//     R : taille du filter top-hat ou gaussien
//     kmin,kmax : intervalle de recherche
//     eps : precision requise sur les valeurs
// Return:
//     position (en k) du maximum
{
  if(R<=0. || kmin<=0 || kmax<=0. || kmin>=kmax) {
    cout<<"VarianceSpectrum::FindMaximum: Error R<=0 or kmin<=0 or kmax<=0 or kmin>=kmax || eps<=0"<<endl;
    throw ParmError("VarianceSpectrum::FindMaximum: Error R<=0 or kmin<=0 or kmax<=0 or kmin>=kmax || eps<=0");
  }

  R_ = R;

  int n = 10; // toujours >2
  double lkmin = log10(kmin), lkmax = log10(kmax), dlk = (lkmax-lkmin)/n;

  double lkfind=lkmin, pkfind=-1.;
  while(1) {
    for(int i=0; i<=n; i++) {
      double lk = lkmin  + i*dlk;
      double v = (*this)(pow(10.,lk));
      if(v<pkfind) continue;
      pkfind = v; lkfind = lk;
    }
    //cout<<"VarianceSpectrum::FindMaximum: lkfind="<<lkfind<<" pkfind="<<pkfind
    //    <<" lkmin,max="<<lkmin<<","<<lkmax<<" dlk="<<dlk<<endl;
    // --- Convergence si l'encadrement de "kfind" est tel que "dk/kfind<eps"
    // On a dk = 10^(lkfind+dlk) - 10^(lkfind-dlk) = kfind * (10^(dlk) - 10^(-dlk))
    if( pow(10.,dlk)-pow(10.,-dlk) < eps ) break;
    if(lkfind-dlk>lkmin) lkmin = lkfind-dlk;
    if(lkfind+dlk<lkmax) lkmax = lkfind+dlk;
    dlk = (lkmax-lkmin)/n;
  }

  return pow(10.,lkfind);
}

int VarianceSpectrum::FindLimits(double R,double high,double &kmin,double &kmax,double eps)
// Retourne "[kmin,kmax]" tel que la fonction a integrer soit "f(k) <= high"
//   La recherche a lieu entre [kmin,kmax] par pas logarithmiques
// Input:
//     R : taille du filter top-hat ou gaussien
//     kmin,kmax : intervalle de recherche
//     eps : precision requise sur les valeurs kmin et kmax
// Output:
//     kmin,kmax telles que "f(k) <= high"
// Return:
//     rc  = 0 si OK
//     rc |= 1 "f(kmin) >= high"   (bit0 =1)
//     rc |= 2 "f(kmax) >= high"   (bit1 =1)
//     rc |= 4 "f(k) < high pour tout k"   (bit2 =1)
{
  if(R<=0. || kmin<=0 || kmax<=0. || kmin>=kmax  || eps<=0.) {
    cout<<"VarianceSpectrum::FindLimits: Error R<=0 or kmin<=0 or kmax<=0 or kmin>=kmax or eps<=0"<<endl;
    throw ParmError("VarianceSpectrum::FindLimits: Error R<=0 or kmin<=0 or kmax<=0 or kmin>=kmax || eps<=0");
  }

  R_ = R;
  int n = 10; // toujours >2

  int rc = 0;
  double lkmin,lkmax,dlk,lkfind;

  // --- Find kmin
  lkmin=log10(kmin); lkmax=log10(kmax); dlk=(lkmax-lkmin)/n;
  while(1) {
    lkfind = lkmin;
    for(int i=0;i<=n;i++) {
      if( (*this)(pow(10,lkfind)) >= high ) break;
      lkfind = lkmin + i*dlk;
    }
    //cout<<"VarianceSpectrum::FindLimits[kmin]: lkfind="<<lkfind
    //    <<" lkmin,max="<<lkmin<<","<<lkmax<<" dlk="<<dlk<<endl;
    if(fabs(lkfind-lkmax)<dlk/2.) {rc |= 4; return rc;} // protect against f(k)<high for all k
    if( pow(10.,dlk)-pow(10.,-dlk) < eps ) break;
    if(lkfind-dlk>lkmin) lkmin = lkfind-dlk;
    if(lkfind+dlk<lkmax) lkmax = lkfind+dlk;
    dlk = (lkmax-lkmin)/n;
  }
  if(lkfind-lkmin<dlk/2.) rc |= 1;  // f(kmin) >= high
  else kmin = pow(10.,lkmin);
  //cout<<"rc="<<rc<<" lkmin="<<lkmin<<"  pk="<<(*this)(pow(10.,lkmin))<<endl;

  // --- Find kmax
  lkmin=log10(kmin); lkmax=log10(kmax); dlk=(lkmax-lkmin)/n;
  while(1) {
    lkfind=lkmax;
    for(int i=0;i<=n;i++) {
      if( (*this)(pow(10,lkfind)) >= high ) break;
      lkfind -= dlk;
      lkfind = lkmax - i*dlk;
    }
    //cout<<"VarianceSpectrum::FindLimits[kmax]: lkfind="<<lkfind
    //    <<" lkmin,max="<<lkmin<<","<<lkmax<<" dlk="<<dlk<<endl;
    if( pow(10.,dlk)-pow(10.,-dlk) < eps ) break;
    if(lkfind-dlk>lkmin) lkmin = lkfind-dlk;
    if(lkfind+dlk<lkmax) lkmax = lkfind+dlk;
    dlk = (lkmax-lkmin)/n;
  }
  if(lkmax-lkfind<dlk/2.) rc |= 2;  // f(kmax) >= high
  else kmax = pow(10.,lkmax);
  //cout<<"rc="<<rc<<" lkmax="<<lkmax<<"  pk="<<(*this)(pow(10.,lkmax))<<endl;

  return rc;
}