Changeset 3285 in Sophya for trunk/Cosmo/SimLSS/cosmocalc.cc

Timestamp:

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

Author:

cmv

Message:

update logique Init_Spline de CosmoCalc cmv 29/07/07

File:

• trunk/Cosmo/SimLSS/cosmocalc.cc (modified) (6 diffs)

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;
 }