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

#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)
// Les x_i sont regulierement espaces
//   et x_0=xmin et x_{n-1}=xmax
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;
  _xlarg = _xmax-_xmin;
  _dx    = _xlarg/(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/_xlarg*_nm1);  // 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/_xlarg*_nm1+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);
}

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