Changeset 3285 in Sophya for trunk/Cosmo/SimLSS

Timestamp:

Jul 29, 2007, 9:01:31 PM (18 years ago)

Author:

cmv

Message:

update logique Init_Spline de CosmoCalc cmv 29/07/07

Location:

trunk/Cosmo/SimLSS

Files:

• cmvobserv3d.cc (modified) (2 diffs)
• cmvtuniv.cc (modified) (5 diffs)
• cosmocalc.cc (modified) (6 diffs)
• cosmocalc.h (modified) (1 diff)

trunk/Cosmo/SimLSS/cmvobserv3d.cc

@@ -60 +60 @@
  double h100=0.71, om0=0.267804, or0=7.9e-05, ol0=0.73,w0=-1.;
  double zref = 0.5;
- double perc=0.01,dzinc=-1.,dzmax=5.; unsigned short glorder=4;
+ double perc=0.01,dzinc=-1.,dzmax=-1.; unsigned short glorder=4;
 
  // *** Spectrum and variance definition
@@ -193 +193 @@
  univ.SetInteg(perc,dzinc,dzmax,glorder);
  univ.SetDynParam(h100,om0,or0,ol0,w0);
+ univ.PrtInteg();
  univ.Print();
  double loscomref = univ.Dloscom(zref);

trunk/Cosmo/SimLSS/cmvtuniv.cc

@@ -42 +42 @@
  cout<<"z1="<<z1<<"  z2="<<z2<<"  dz="<<dz<<endl;
 
- double perc=0.01,dzinc=-1.,dzmax=5.; unsigned short glorder=4;
+ double perc=0.01,dzinc=-1.,dzmax=-1.; unsigned short glorder=4;
  if(narg>2) sscanf(arg[2],"%lf,%lf,%lf,%hu",&perc,&dzinc,&dzmax,&glorder);
- if(dzinc<=0.) dzinc = dz/2.;
  cout<<" perc="<<perc<<" dzinc="<<dzinc<<" dzmax="<<dzmax<<" glorder="<<glorder<<endl;
 
@@ -50 +49 @@
  univ.SetInteg(perc,dzinc,dzmax,glorder);
  univ.SetDynParam(h100,om0,or0,ol0,w0);
+ univ.PrtInteg();
  univ.Print();
 
@@ -55 +55 @@
  univ1.SetInteg(perc,dzinc,dzmax,glorder);
  univ1.SetDynParam(h100,om0,or0,ol0,w0);
+ univ1.PrtInteg();
  univ1.Print();
 
@@ -202 +203 @@
  vector<double> Dinterp;
  double zmin = Z[0], zmax = Z[0];
- for(int i=0;i<Z.size();i+=ninc) {
+ for(unsigned int i=0;i<Z.size();i+=ninc) {
    Dinterp.push_back(D[i]);
    zmax = Z[i];
@@ -214 +215 @@
  double xnt[n];
    
- for(int i=0;i<Z.size();i++) {
+ for(unsigned int i=0;i<Z.size();i++) {
    if(Z[i]>zmax) break;
    xnt[0] = Z[i];

trunk/Cosmo/SimLSS/cosmocalc.cc

@@ -106 +106 @@
 void CosmoCalc::SetInteg(double dperc,double dzinc,double dzmax,unsigned short order)
 {
-  _dperc = fabs(dperc);
-
-  _dzinc = fabs(dzinc);
-  if(_dzinc>_zmax/5.) _dzinc = _zmax/5.; // Protection against big dzinc
-
-  _dzmax = fabs(dzmax);
+  _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;
 
@@ -119 +121 @@
 
   _computespl=true;
+}
+
+
+void CosmoCalc::PrtInteg(void)
+{
+ printf("CosmoCalc::PrtInteg: dperc=%g dzinc=%g dzmax=%g glorder=%hu\n",_dperc,_dzinc,_dzmax,_glorder);
 }
 
@@ -369 +377 @@
   if(_usespline) {
     if( (!_computespl) && (z==_zold) ) return _integval;
-    if(z>_zmax) {_zmax = z; _computespl = true;}
+    if(z>_zmax) {_zmax = z+_dzmax; _computespl = true;}
     if(_computespl) Init_Spline();
     _integval = _spline->CSplineInt(z);
@@ -392 +400 @@
   //-- Look for intervalles and integrate
   vector<double> x,y;
+  x.push_back(0.); y.push_back(0.);
   double zbas=0., fbas = Integrand(zbas);
-  x.push_back(zbas); y.push_back(0.);
-  for(double z=_dzinc;z<_zmax+_dzinc && zbas<_zmax;z+=_dzinc) {
+  for(double z=_dzinc;;z+=_dzinc) {
     double f = Integrand(z);
-    bool doit = false;
     //cout<<fbas<<","<<f<<"  :  "<<fabs(f-fbas)<<" >? "<<_dperc*fabs(fbas)<<endl;
-    if( fabs(f-fbas)>_dperc*fabs(fbas) ) doit = true;
-    if( z-zbas>_dzmax ) doit = true;
-    if( z>_zmax ) {doit = true; z=_zmax;}
-    if( ! doit ) continue;
-    x.push_back(z);
-    double sum = 0., dz = z-zbas;
-    for(unsigned short i=0;i<_glorder;i++) sum += _wgausl[i]*Integrand(zbas+_xgausl[i]*dz);
-    y.push_back(sum*dz);
-    //cout<<"...set... "<<zbas<<","<<z<<"  ,  "<<fbas<<","<<f<<endl;
-    zbas = z;
-    fbas = f;
+    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);
+    }
   }
 
@@ -420 +436 @@
     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]);
@@ -426 +454 @@
   _computespl = false;
 
-  cout<<"CosmoCalc::Init_Spline called: (zmax="<<_zmax<<") _nspl="<<_nspl<<endl;
-  for(int i=0;i<(10<_nspl?10:_nspl);i++) cout<<_xspl[i]<<" ";
-  cout<<" ... "<<_xspl[_nspl-1]<<endl;
-
   return _nspl;
 }

trunk/Cosmo/SimLSS/cosmocalc.h

@@ -21 +21 @@
 
   double ZMax(void) { return _zmax;}
-  void SetInteg(double dperc=0.01,double dzinc=0.05,double dzmax=5.,unsigned short order=4);
+  void SetInteg(double dperc=-1.,double dzinc=-1.,double dzmax=-1.,unsigned short order=4);
+  void PrtInteg(void);
 
   void SetDynParam(double h100,double om0,double or0,double ol0,double w0=-1.);