#include "machdefs.h" #include #include #include #include #include #include #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<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="<=_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<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)_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 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<? "<<_dperc*fabs(fbas)<_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... "<_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<5)? 5: _nspl; cout<<_nspl<<"..."; for(int i=0;iComputeCSpline(); _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<