source: Sophya/trunk/Cosmo/SimLSS/geneutils.cc@ 3187

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

#include "pexceptions.h"

#include "histos.h"
#include "hisprof.h"
#include "srandgen.h"

#include "geneutils.h"

//-------------------------------------------------------------------
// Classe d'interpolation lineaire:
// Le vecteur y a n elements y_i tels que y_i = f(x_i)
//   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)
{
  if(_xmin>=_xmax || _y.size()<=2) {  // au moins 3 points!
   cout<<"InterpFunc::InterpFunc : bad arguments values"<<endl;
   throw ParmError("InterpFunc::InterpFunc : bad arguments values");
  }
  _nm1   = _y.size()-1;
  _dx    = (_xmax-_xmin)/(double)_nm1;
}

double InterpFunc::Linear(double x,unsigned short& ok)
{
  ok=0; if(x<_xmin) ok=1; else if(x>_xmax) ok=2;
  x -= _xmin;
  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);
}

double InterpFunc::Parab(double x,unsigned short& ok)
{
  ok=0; if(x<_xmin) ok=1; else if(x>_xmax) ok=2;
  x -= _xmin;
  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;
}

//-------------------------------------------------------------------
int FuncToHisto(GenericFunc& func,Histo& h,bool logaxex)
// Remplit l'histo 1D "h" avec la fonction "func"
// INPUT:
// logaxex = false : remplissage lineaire
//        les abscisses "x" des bins sont remplis avec f(x)
// logaxex = true : remplissage logarithmique (base 10)
//        les abscisses "x" des bins sont remplis avec f(10^x)
// RETURN:
//       0 = OK
//       1 = error
{
  if(h.NBins()<=0) return 1;

  h.Zero();

  for(int_4 i=0;i<h.NBins();i++) {
    double x = h.BinCenter(i);
    if(logaxex) x = pow(10.,x);
    h.SetBin(i,func(x));
  }

  return 0;
}

int FuncToVec(GenericFunc& func,TVector<r_8>& v,double xmin,double xmax,bool logaxex)
// Remplit le TVector avec la fonction "func"
// INPUT:
// logaxex = false : remplissage lineaire
//        les abscisses "x" des bins sont remplis avec f(x)
// logaxex = true : remplissage logarithmique (base 10)
//        les abscisses "x" des bins sont remplis avec f(10^x)
// RETURN:
//       0 = OK
//       1 = error
// Remarque:
//  v(i) = f(xmin+i*dx) avec dx = (xmax-xmin)/v.NElts()
{
  if(v.NElts()<=0 || xmax<=xmin) return 1;

  v = 0.;
  double dx = (xmax-xmin)/v.NElts();
   
  for(int_4 i=0;i<v.NElts();i++) {
    double x = xmin + i * dx;;
    if(logaxex) x = pow(10.,x);
    v(i) = func(x);
  }

  return 0;
}

//-------------------------------------------------------------------
double AngSol(double dtheta,double dphi,double theta0)
// Retourne l'angle solide d'un "rectangle" et coordonnees spheriques
// de DEMI-COTE "dtheta" x "dphi" et centre en "theta0"
// Attention: Le "theta0" de l'equateur est Pi/2  (et non pas zero)
//            Les unites des angles sont en radians
//               theta0 in [0,Pi]
//               dtheta in [0,Pi]
//               dphi   in [0,2Pi]
// Return: l'angle solide en steradian
{
  double theta1 = theta0-dtheta, theta2 = theta0+dtheta;
  if(theta1<0.) theta1=0.;
  if(theta2>M_PI) theta2=M_PI;

  return   2.*dphi * (cos(theta1)-cos(theta2));
}

double AngSol(double dtheta)
// Retourne l'angle solide d'une calotte spherique de demi-ouverture "dtheta"
// Attention: Les unites des angles sont en radians
//               dtheta in [0,Pi]
// Return: l'angle solide en steradian
// Approx pour theta petit: PI theta^2
{
  return   2.*M_PI * (1.-cos(dtheta));
}

//-------------------------------------------------------------------
unsigned long PoissRandLimit(double mu,double mumax)
{
  double pp,ppi,x;
  unsigned long n;

  if(mu>=mumax) {
    pp = sqrt(mu);
    while( (x=pp*NorRand()) < -mu );
    return (unsigned long)(mu+x+0.5);
  }

  ppi = pp = exp(-mu);
  x = drand01();
  n = 0;
  while (x > ppi) {
    n++;
    pp = mu*pp/(double)n;
    ppi += pp;
  }
  return n;
}