source: Sophya/trunk/Cosmo/SimLSS/cosmocalc.cc@ 4043

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

#include "pexceptions.h"
#include "cspline.h"

#include "constcosmo.h"
#include "cosmocalc.h"
#include "geneutils.h"

namespace SOPHYA {

///////////////////////////////////////////////////////////
//*********************** CosmoCalc *********************//
///////////////////////////////////////////////////////////

//----------------------------------------------------------
// flat = 0 : do not set Otot=1, compute Otot from others
//      = 1 : change Omatter to get Otot=1
//      = 2 : change Olambda to get Otot=1
//      = 3 : change Orelat to get Otot=1
// usespline = "true" : compute integrale then use spline to extrapolate
//           = "false" : compute integrale everytime
// zmax : maximum redshift(only if usespline=true)
         
CosmoCalc::CosmoCalc(unsigned short flat,bool usespline,double zmax)
  : _flat(flat),
    _zold(-1.) , _integval(0.),
    _usespline(usespline), _computespl(true), _zmax(zmax), _spline(NULL),
    _nspl(0), _xspl(NULL), _yspl(NULL)
{
  if(_usespline && _zmax<=0.) {
    cout<<"CosmoCalc::CosmoCalc: Error bad _zmax value    zmax="<< _zmax<<endl;
    throw ParmError("CosmoCalc::UseSpline: Error bad _zmax value");
  }
  DefaultParam();
  SetInteg();
}

CosmoCalc::CosmoCalc(CosmoCalc& univ)
  : _spline(NULL), _xspl(NULL), _yspl(NULL)
{
  Clone(univ);
  if(_usespline && _zmax<=0.) {
    cout<<"CosmoCalc::CosmoCalc: Error bad _zmax value    zmax="<< _zmax<<endl;
    throw ParmError("CosmoCalc::UseSpline: Error bad _zmax value");
  }
}

CosmoCalc::~CosmoCalc(void)
{
  if(_spline != NULL) delete _spline; _spline = NULL;
  if(_xspl != NULL) delete [] _xspl; _xspl = NULL;
  if(_yspl != NULL) delete [] _yspl; _yspl = NULL;
}

void CosmoCalc::Clone(CosmoCalc& univ)
{
 _flat = univ._flat;
 _h100 = univ._h100;
 _H0 = univ._H0;
 _Olambda0 = univ._Olambda0;
 _W0 = univ._W0;
 _Wa = univ._Wa;
 _Omatter0 = univ._Omatter0;
 _Obaryon0 = univ._Obaryon0;
 _Orelat0 = univ._Orelat0;
 _Otot0 = univ._Otot0;
 _Ocurv0 = univ._Ocurv0;
 _Dhubble = univ._Dhubble;

 _dperc = univ._dperc;
 _dzinc = univ._dzinc;
 _dzmax = univ._dzmax;
 _glorder = univ._glorder;
 if(_xgausl.size()>0) {
   _xgausl.resize(0); _wgausl.resize(0);
   for(int i=0;i<(int)_xgausl.size();i++) {
     _xgausl.push_back(univ._xgausl[i]);
     _wgausl.push_back(univ._wgausl[i]);
   }
 }

 _zold = -1.;
 _integval = -1.;

 _usespline = univ._usespline;
 _computespl = true;
 _zmax = univ._zmax;
 _nspl = 0;
 if(_spline) delete _spline; _spline = NULL;
 if(_xspl) delete [] _xspl; _xspl = NULL;
 if(_yspl) delete [] _yspl; _yspl = NULL;
}

//----------------------------------------------------------
// On change le zmax de la cosmologie
void CosmoCalc::ChangeZmax(double zmax,double dzinc,double dzmax)
{
  _zmax = zmax;
  if(dzinc<=0. ) dzinc = _dzinc;
  if(dzmax<=0.)  dzmax = _dzmax;
  cout<<"CosmoCalc::ChangeZmax: zmax="<<_zmax<<" dzinc="<<dzinc<<" dzmax="<<dzmax<<endl;
  SetInteg(_dperc,dzinc,_dzmax,_glorder);
}

//----------------------------------------------------------
// On va couper l'intervalle entre [0,zmax]:
// On parcourt [0,zmax] par pas de dzinc
// et on cree un intervalle
//     - si la fonction varie de plus de "dperc"
//     - ou si l'increment en z est superieur a "dzmax"
void CosmoCalc::SetInteg(double dperc,double dzinc,double dzmax,unsigned short order)
{
  _dperc = dperc; _dzinc = dzinc; _dzmax = dzmax;

  if(_dperc<=0.) _dperc = 0.01;

  if(_dzinc<=0.) _dzinc = _zmax/100.;
  if(_dzinc>=_zmax/2.) _dzinc = _zmax/2.; // Protection against big dzinc

  if(_dzmax<=0.) _dzmax = _zmax;
  if(_dzmax<_dzinc) _dzmax = _dzinc;

  _glorder = order;
  if(_glorder<=1) _glorder = 4;
  Compute_GaussLeg(_glorder,_xgausl,_wgausl,0.,1.);

  _computespl=true;
}


void CosmoCalc::PrtInteg(void)
{
 printf("CosmoCalc::PrtInteg: dperc=%g dzinc=%g dzmax=%g glorder=%hu\n",_dperc,_dzinc,_dzmax,_glorder);
}

void CosmoCalc::SetObaryon0(double v)
{
  _Obaryon0 = v;
  if(_Obaryon0<0.) {
    cout<<"CosmoCalc::SetObaryon0: Error bad _Obaryon0 value: "<<_Obaryon0<<endl;
    throw ParmError("CosmoCalc::SetObaryon0: Error bad _Obaryon0 value");
  }
}

void CosmoCalc::SetDynParam(double h100,double om0,double or0,double ol0,double w0,double wa)
{
  _computespl=true;

  _h100 = h100;
  _H0 = 100. * _h100 ;
  _Dhubble = SpeedOfLight_Cst / _H0;

  _Omatter0 = om0;
  _Orelat0 = or0;
  _Olambda0 = ol0;
  _W0 = w0;
  _Wa = wa;
  _Otot0 = _Olambda0 + _Omatter0 + _Orelat0;
  _Ocurv0 = 1. - _Otot0;
  if( _h100<0. || _Omatter0<0. || _Orelat0<0. || _Olambda0<0. ) {
    cout<<"CosmoCalc::SetDynParam: Error bad parameter value"<<endl;
    throw ParmError("CosmoCalc::SetDynParam: Error bad parameter value");
  }
  if(_flat==0) return;

  _Otot0 = 1.; _Ocurv0 = 0.;
  if(_flat==1) {
    _Omatter0 = _Otot0 - _Olambda0 - _Orelat0;
    if( _Omatter0<0. ) {
      cout<<"CosmoCalc::SetDynParam: Error bad _Omatter0 value: "<<_Omatter0<<endl;
      throw ParmError("CosmoCalc::SetDynParam: Error bad _Omatter0 value");
    }
    return;
  }

  //--
  if(_flat==2) {
    _Olambda0 = _Otot0 - _Omatter0 - _Orelat0;
    if( _Olambda0<0. ) {
      cout<<"CosmoCalc::SetDynParam: Error bad _Olambda0 value: "<<_Olambda0<<endl;
      throw ParmError("CosmoCalc::SetDynParam: Error bad _Olambda0 value");
    }
    return;
  }

  //--
  if(_flat==3) {
    _Orelat0 = _Otot0 - _Omatter0 - _Olambda0;
    if( _Orelat0<0. ) {
      cout<<"CosmoCalc::SetDynParam: Error bad _Orelat0 value: "<<_Orelat0<<endl;
      throw ParmError("CosmoCalc::SetDynParam: Error bad _Orelat0 value");
    }
    return;
  }

  cout<<"CosmoCalc::SetDynParam: Error bad _flat value: "<<_flat<<endl;
  throw ParmError("CosmoCalc::SetDynParam: Error bad _flat value");

}

void CosmoCalc::DefaultParam(void)
{
  _computespl=true;

  double h100 = 0.71;
  double ol0 = 0.73, w0 = -1., wa = 0.;
  double om0 = 0.135/(h100*h100);
  // Relat = photons (2.725 K) + neutrinos (1.9 K)
  double or0 = 2.47e-5 * (1. + 21./8.*pow(T_NU_Par/T_CMB_Par,4.)) / (h100*h100);
  SetDynParam(h100,om0,or0,ol0,w0,wa);
  _Obaryon0 = 0.0224/(h100*h100);

  return;
}

//----------------------------------------------------------
void CosmoCalc::Print(double z)
{
  if(z<0.) z = 0.;

  printf("CosmoCalc::Print(spl=%d,zmax=%g,flat=%u)  for z=%g (a=%g)\n"
         ,int(_usespline),ZMax(),Flat(),z,1./(1.+z));
  printf("h100=%g H0=%g Dhub=%g  H(z)=%g  Rhoc=%g g/cm^3 =%g Msol/Mpc^3\n"
        ,h100(),H0(),Dhubble(),H(z),Rhoc(z),Rhoc(z)*GCm3toMsolMpc3_Cst);
  printf("Olambda=%g  W=%g   (%g+(%g)*(1-a))\n",Olambda(z),W(z),_W0,_Wa);
  printf("Omatter=%g  Obaryon=%g\n",Omatter(z),Obaryon(z));
  printf("Orelat=%g\n",Orelat(z));
  printf("Otot=%g  Ocurv=%g\n",Otot(z),Ocurv(z));
  if(z <= 0.) return;
  printf("Distance comoving: los=%g Mpc, transv=%g Mpc\n",Dloscom(z),Dtrcom(z));
  printf("Distance: angular=%g Mpc, lum=%g Mpc\n",Dang(z),Dlum(z));
  printf("Volume comoving element: %g Mpc^3/sr/dz=1\n",dVol(z));
  printf("Volume comoving in [0,z] for 4Pi sr: %g Mpc^3\n",Vol4Pi(z));

}
//----------------------------------------------------------
double CosmoCalc::Olambda(double z)
// Equation d'etat de l'energie noire: P=W(z)*Rho*C^2
//   avec W(z) = W0 + Wa*(1-a) = W0 + Wa * z/(1+z) et a=1/(1+z)
{
  double zp1 = 1. + z;
  double v = _Olambda0 * pow(zp1,3.*(1.+_W0+_Wa));
  if(_Wa!=0.) v *= exp(-3.*_Wa*z/zp1);
  return v / E2(z);
}

double CosmoCalc::Omatter(double z)
{
  double zp1 = 1. + z;
  return _Omatter0 * zp1*zp1*zp1 /E2(z);
}

double CosmoCalc::Obaryon(double z)
{
  double zp1 = 1. + z;
  return _Obaryon0 * zp1*zp1*zp1 /E2(z);
}

double CosmoCalc::Orelat(double z)
{
  double z2 = (1.+z)*(1.+z);
  return _Orelat0 * z2*z2 /E2(z);
}

double CosmoCalc::Ocurv(double z)
{
  return _Ocurv0 * (1.+z)*(1.+z) /E2(z);
}

double CosmoCalc::Otot(double z)
{
  return Olambda(z) + Omatter(z) + Orelat(z);
}

double CosmoCalc::Rhoc(double z)
// Densite critique au redshift "z"  en "g/cm^3"
// Attention: c'est une densite cad une masse par volume (non-comobile)
{
  double h2 = H(z) / MpctoMeters_Cst;  h2 *= h2;
  return 3.* h2 / (8.*M_PI*G_Newton_Cst) * 1000.;
}

//----------------------------------------------------------
double CosmoCalc::Dloscom(double z)     /* Mpc comobile */
{
  return _Dhubble * NInteg(z);
}

double CosmoCalc::Dtrcom(double z)     /* Mpc comobile */
{
  double v = NInteg(z);  // Zero curvature

  if(_flat) return _Dhubble * v;

  if(_Ocurv0>0.) {  // hyperbolique
    double s = sqrt(_Ocurv0);
    v = sinh(s*v)/s;
  } else if(_Ocurv0<0.) {  // spherique
    double s = sqrt(-_Ocurv0);
    v = sin(s*v)/s;
  }

  return _Dhubble * v;

}

double CosmoCalc::Dang(double z)     /* Mpc */
{
  return Dtrcom(z) / (1.+z);
}

double CosmoCalc::Dlum(double z)     /* Mpc */
{
  return (1.+z) * Dtrcom(z);
}

double CosmoCalc::dVol(double z)     /* Mpc^3/ sr / unite z */
{
  double d = Dtrcom(z);
  return _Dhubble * d*d / E(z);
}

//----------------------------------------------------------
double CosmoCalc::Vol4Pi(double z)  /* Mpc^3 pour 4Pi sr entre [0,z] */
  // --- on pose x = dm/dh = (1+z)*da/dh, s = sqrt(|Ok|)
  // Ok=0 : V(z)/dh^3 = (4*Pi/3) * x^3
  // Ok>0 : V(z)/dh^3 = (4*Pi/(2*Ok)) * (x*sqrt(1+Ok*x^2) - 1/s * asinh(s*x))
  // Ok<0 : V(z)/dh^3 = (4*Pi/(2*Ok)) * (x*sqrt(1+Ok*x^2) - 1/s * asin(s*x) )
  // --- on pose y = s*x  (y>0)
  // Ok>0 : V(z)/dh^3 = (4*Pi/(2*Ok*s)) * (y*sqrt(1+y^2) - asinh(y))
  // Ok<0 : V(z)/dh^3 = (4*Pi/(2*Ok*s)) * (y*sqrt(1-y^2) - asin(y) )
  // --- a petit "z" on a pour "y -> 0+" en faisant le DL
  // Ok>0 : V(z)/dh^3 =  4*Pi*y^3 / (3*Ok*s) * (1 - 3*y^2/10) + O(y^7)
  // Ok<0 : V(z)/dh^3 = -4*Pi*y^3 / (3*Ok*s) * (1 + 3*y^2/10) + O(y^7)
  //                          (remarque: Ok^(3/2) = s*Ok)
{
  double v,x = Dtrcom(z)/_Dhubble;

  if(_flat) {
    v = 4.*M_PI/3. * x*x*x;
  } else if(_Ocurv0>0.) {
    double s = sqrt(_Ocurv0), y = s*x, y2 = y*y;
    if(y<1e-6) {
      v = 1. - 3.*y2/10.;
      v *= 4.*M_PI*y*y2 / (3.*_Ocurv0*s);
    } else {
      v = y*sqrt(1.+y2) - asinh(y);
      v *= 4.*M_PI/(2.*_Ocurv0*s);
    }
  } else if (_Ocurv0<0.) {
    double s = sqrt(-_Ocurv0), y = s*x, y2 = y*y;
    if(y<1e-6) {
      v = 1. + 3.*y2/10.;
      v *= -4.*M_PI*y*y2 / (3.*_Ocurv0*s);
    } else {
      v = y*sqrt(1.-y2) - asin(y);
      v *= 4.*M_PI/(2.*_Ocurv0*s);
    }
  } else {
    v = 4.*M_PI/3. * x*x*x;
  }

  return v * _Dhubble*_Dhubble*_Dhubble;
}

double CosmoCalc::Vol4Pi(double z1,double z2)  /* Mpc^3 pour 4Pi sr entre [z1,z2] */
{
  return Vol4Pi(z2) - Vol4Pi(z1);
}

double CosmoCalc::E2(double z) const
{
  double zp1 = 1. + z;

  double oldum = _Olambda0;
  if(oldum>0.) {
    oldum *= pow(zp1,3.*(1.+_W0+_Wa));
    if(_Wa!=0.) oldum *= exp(-3.*_Wa*z/zp1);
  }

  return oldum + zp1*zp1*(_Ocurv0 + zp1*(_Omatter0+zp1*_Orelat0));
}

//----------------------------------------------------------
double CosmoCalc::ZFrLos(double loscom /* Mpc com */, int niter)
// Recherche du redshift correspondant a une distance comobile
// le long de la ligne de visee (radiale) egale a "loscom" Mpc
// niter = nomber of iterations for precision measurement
{
  if(niter<3) niter = 6;
  double dz = ZMax()/10.; if(dz<=0.) dz = 0.1;
  double zmin=0., zmax=0.;
  while(Dloscom(zmax)<loscom) zmax += dz;
  if(zmax==0.) return 0.;
  for(int i=0; i<niter; i++) {
    zmin=zmax-dz; if(zmin<0.) zmin=0.;
    dz /= 10.;
    for(double z=zmin; z<zmax+dz; z+=dz) {
      double d = Dloscom(z);
      if(d<loscom) continue;
      zmax = z;
      //cout<<"ZFrLos: z="<<zmax<<" d="<<d<<" / "<<loscom<<endl;
      break;
    }
  }
  return zmax;
}

//----------------------------------------------------------
double CosmoCalc::NInteg(double z)
{
  if(z<0.) {
    cout<<"CosmoCalc::NInteg: Error bad z value  z="<<z<<endl;
    throw ParmError("CosmoCalc::NInteg: Error bad z value");
  }

  // On utilise le spline
  if(_usespline) {
    if( (!_computespl) && (z==_zold) ) return _integval;
    if(z>_zmax) {_zmax = z+_dzmax; _computespl = true;}
    if(_computespl) Init_Spline();
    _integval = _spline->CSplineInt(z);
    _zold = z;
    return _integval;
  }

  // On calcule l'integrale
  if(z != _zold) {
    _integval = IntegrateFunc(*this,0.,z,_dperc,_dzinc,_dzmax,_glorder);
    _zold = z;
  }
  return _integval;

}

int_4 CosmoCalc::Init_Spline(void)
{
  if(_spline!=NULL)
    {delete _spline; delete [] _xspl; delete [] _yspl;}

  //-- Look for intervalles and integrate
  vector<double> x,y;
  x.push_back(0.); y.push_back(0.);
  double zbas=0., fbas = Integrand(zbas);
  for(double z=_dzinc;;z+=_dzinc) {
    double f = Integrand(z);
    //cout<<fbas<<","<<f<<"  :  "<<fabs(f-fbas)<<" >? "<<_dperc*fabs(fbas)<<endl;
    if( z>_zmax || z-zbas>_dzmax || fabs(f-fbas)>_dperc*fabs(fbas) ) {
      double sum=0.,  dz=z-zbas;
      for(unsigned short i=0;i<_glorder;i++) sum += _wgausl[i]*Integrand(zbas+_xgausl[i]*dz);
      x.push_back(z); y.push_back(sum*dz);
      //cout<<"...set... "<<zbas<<","<<z<<"  ,  "<<fbas<<","<<f<<endl;
      zbas = z;  fbas = f;
      if( z>_zmax ) break;
    }
  }

  //-- Protection car il faut au moins 4 points pour un spline cubique
  if(x.size()<4) {
    x.resize(0); y.resize(0);
    x.push_back(0.); y.push_back(0.);
    for(int i=1;i<4;i++) {
      x.push_back(i*_zmax/3.);
      double sum=0., dz=x[i]-x[i-1];
      for(unsigned short i=0;i<_glorder;i++) sum += _wgausl[i]*Integrand(x[i-1]+_xgausl[i]*dz);
      y.push_back(sum*dz);
    }
  }

  //-- Fill spline
  _nspl = x.size();
  _xspl = new double[_nspl];
  _yspl = new double[_nspl];
  for(int i=0;i<_nspl;i++) {
    _xspl[i] = x[i];
    _yspl[i] = y[i];
    if(i!=0) _yspl[i] += _yspl[i-1];
  }

#if 1
  cout<<"CosmoCalc::Init_Spline called: (zmax="<<_zmax<<") _nspl="<<_nspl<<endl;
  int n = (_nspl>5)? 5: _nspl;
  cout<<_nspl<<"...";
  for(int i=0;i<n;i++) cout<<_xspl[i]<<" ("<<_yspl[i]<<") ";
  cout<<"\n... ";
  n = (_nspl<5)? 0: _nspl-5;
  for(int i=n;i<_nspl;i++) cout<<_xspl[i]<<" ("<<_yspl[i]<<") ";
  cout<<endl;
#endif

  double yp1 = Integrand(_xspl[0]);
  double ypn = Integrand(_xspl[_nspl-1]);
  _spline = new CSpline(_nspl,_xspl,_yspl,yp1,ypn,CSpline::NoNatural,false);
  _spline->ComputeCSpline();
  _computespl = false;

  return _nspl;
}

//==========================================================================
//==========================================================================
//==========================================================================
double LargeurDoppler(double v, double nu)
// largeur doppler pour une vitesse v en km/s et une frequence nu
{
 return v / SpeedOfLight_Cst * nu;
}

double DzFrV(double v, double zred)
// largeur en redshift pour une vitesse v en km/s au redshift zred
{
  return v / SpeedOfLight_Cst * (1. + zred);
}

double DNuFrDz(double dzred,double nu_at_0,double zred)
// Largeur DNu pour une largeur en redshift "dzred" au redshift "zred"
//    pour la frequence "nu_at_0" a z=0
// nu =  NuHi(z=0)/(1.+z0)
// dnu = NuHi(z=0)/(1.+z0-dz/2) - NuHi/(1.+z0+dz/2)
//     = NuHi(z=0)*dz/[ (1+z0)^2 - (dz/2)^2 ]
//     = NuHi(z=0)*dz/(1.+z0)^2 / [ 1 - [dz/(1+z0)/2)]^2 ]
//     = NuHi(z=0)*dz/(1.+z0)^2 / [1 - dz/(1+z0)/2] / [1 + dz/(1+z0)/2]
//    ~= NuHi(z=0)*dz/(1.+z0)^2   (approx. pour dz<<z0 a l'ordre (dz/z0)^2)
{
  double zp1 = 1.+zred;
  return nu_at_0*dzred/(zp1*zp1)/(1.-dzred/zp1/2.)/(1.+dzred/zp1/2.);
}

double DzFrDNu(double dnu_at_0,double nu_at_0,double zred)
// Largeur en redshift au redshift "zred" pour une largeur
// en frequence "dnu_at_0" a la frequence "nu_at_0" a z=0
{
  if(dnu_at_0<=0.) return 0.;
  double zp1 = 1.+zred;
  double dnusnu0 = dnu_at_0/nu_at_0;
  return 2./dnusnu0 * (sqrt(1.+(dnusnu0*zp1)*(dnusnu0*zp1)) - 1.);
}
double DzFrDNuApprox(double dnu_at_0,double nu_at_0,double zred)
// idem DzFrDNu mais on utilise l'approximation: dnu=NuHi(z=0)*dz/(1.+z0)^2
{
  double zp1 = 1.+zred;
  return dnu_at_0/nu_at_0 *(zp1*zp1);
}

}  // Fin namespace SOPHYA