Changeset 3157 in Sophya for trunk/Cosmo/SimLSS/geneutils.cc

Timestamp:

Jan 25, 2007, 6:04:48 PM (19 years ago)

Author:

cmv

Message:

intro du facteur de croissance dans la simul cmv 25/01/2007

File:

• trunk/Cosmo/SimLSS/geneutils.cc (modified) (4 diffs)

trunk/Cosmo/SimLSS/geneutils.cc

@@ -19 +19 @@
 // Classe d'interpolation lineaire:
 // Le vecteur y a n elements y_i tels que y_i = f(x_i)
-// Les x_i sont regulierement espaces
-//   et x_0=xmin et x_{n-1}=xmax
+//   ou les x_i sont regulierement espaces
+//   et x_0=xmin et x_{n-1}=xmax   (xmax inclus!)
 InterpFunc::InterpFunc(double xmin,double xmax,vector<double>& y)
   : _xmin(xmin), _xmax(xmax), _y(y)
@@ -29 +29 @@
   }
   _nm1   = _y.size()-1;
-  _xlarg = _xmax-_xmin;
-  _dx    = _xlarg/(double)_nm1;
+  _dx    = (_xmax-_xmin)/(double)_nm1;
 }
 
@@ -37 +36 @@
   ok=0; if(x<_xmin) ok=1; else if(x>_xmax) ok=2;
   x -= _xmin;
-  long i = long(x/_xlarg*_nm1);  // On prend le "i" juste en dessous
+  long i = long(x/_dx);  // On prend le "i" juste en dessous
   if(i<0) i=0; else if(i>=_nm1) i=_nm1-1;
   return _y[i] + (_y[i+1]-_y[i])/_dx*(x-i*_dx);
@@ -46 +45 @@
   ok=0; if(x<_xmin) ok=1; else if(x>_xmax) ok=2;
   x -= _xmin;
-  long i = long(x/_xlarg*_nm1+0.5);  // On prend le "i" le + proche
+  long i = long(x/_dx+0.5);  // On prend le "i" le + proche
   if(i<1) i=1; else if(i>=_nm1-1) i=_nm1-2;
   double a = (_y[i+1]-2.*_y[i]+_y[i-1])/(2.*_dx*_dx);
   double b = (_y[i+1]-_y[i-1])/(2.*_dx);
   return _y[i] + (x-i*_dx)*(a*(x-i*_dx)+b);
+}
+
+//-------------------------------------------------------------------
+// Classe d'inversion d'une fonction STRICTEMENT MONOTONE CROISSANTE
+// - Le vecteur y a "Nin" elements y_i tels que "y_i = f(x_i)"
+// - On a x(i) < x(i+1) et y(i) < y(i+1)
+// - La classe renvoie ymin=y(0) , ymax=y(Nin -1)
+//   et le vecteur x = f^-1(y) de "Nout" elements
+//   Les y_i sont regulierement espaces et ymin et ymax
+//   La re-interpolation inverse est faite par lineaire
+InverseFunc::InverseFunc(vector<double>& x,vector<double>& y)
+  : _ymin(0.) , _ymax(0.) , _x(x) , _y(y)
+{
+  int_4 ns = _x.size();
+  if(ns<3 || _y.size()<=0 || ns!=_y.size())
+    throw ParmError("InverseFunc::InverseFunc_Error: bad array size");
+
+  // Check "x" strictement monotone croissant
+  for(int_4 i=0;i<ns-1;i++)
+    if(_x[i+1]<=_x[i]) {
+      cout<<"InverseFunc::InverseFunc_Error: _x array not stricly growing"<<endl;
+      throw ParmError("InverseFunc::InverseFunc_Error: _x array not stricly growing");
+    }
+
+  // Check "y" monotone croissant
+  for(int_4 i=0;i<ns-1;i++)
+    if(_y[i+1]<_y[i]) {
+      cout<<"InverseFunc::InverseFunc_Error: _y array not growing"<<endl;
+      throw ParmError("InverseFunc::InverseFunc_Error: _y array not growing");
+    }
+
+  // define limits
+  _ymin = _y[0];
+  _ymax = _y[ns-1];
+
+}
+
+InverseFunc::~InverseFunc(void)
+{
+}
+
+int InverseFunc::ComputeLinear(long n,vector<double>& xfcty)
+{
+  if(n<3) return -1;
+
+  xfcty.resize(n);
+
+  long i1,i2;
+  double x;
+  for(int_4 i=0;i<n;i++) {
+    double y = _ymin + i*(_ymax-_ymin)/(n-1.);
+    find_in_y(y,i1,i2);
+    double dy = _y[i2]-_y[i1];
+    if(dy==0.) {
+      x = (_x[i2]+_x[i1])/2.; // la fct a inverser est plate!
+    } else {
+      x = _x[i1] + (_x[i2]-_x[i1])/dy * (y-_y[i1]);
+    }
+    xfcty[i] = x;
+  }
+
+  return 0;
+}
+
+int InverseFunc::ComputeParab(long n,vector<double>& xfcty)
+{
+  if(n<3) return -1;
+
+  xfcty.resize(n);
+
+  long i1,i2,i3;
+  double x;
+  for(int_4 i=0;i<n;i++) {
+    double y = _ymin + i*(_ymax-_ymin)/(n-1.);
+    find_in_y(y,i1,i2);
+    // On cherche le 3ieme point selon la position de y / au 2 premiers
+    double my = (_y[i1]+_y[i2])/2.;
+    if(y<my) {i3=i2; i2=i1; i1--;} else {i3=i2+1;}
+    // Protection
+    if(i1<0) {i1++; i2++; i3++;}
+    if(i3==_y.size()) {i1--; i2--; i3--;}
+    // Interpolation parabolique
+    double dy = _y[i3]-_y[i1];
+    if(dy==0.) {
+      x = (_x[i3]+_x[i1])/2.; // la fct a inverser est plate!
+    } else {
+      double X1=_x[i1]-_x[i2], X3=_x[i3]-_x[i2];
+      double Y1=_y[i1]-_y[i2], Y3=_y[i3]-_y[i2];
+      double den = Y1*Y3*dy;
+      double a = (X3*Y1-X1*Y3)/den;
+      double b = (X1*Y3*Y3-X3*Y1*Y1)/den;
+      y -= _y[i2];
+      x = (a*y+b)*y + _x[i2];
+    }
+    xfcty[i] = x;
+  }
+
+  return 0;
 }