source: Sophya/trunk/SophyaLib/NTools/fct1dfit.cc@ 307

#include "machdefs.h"
#include <stdio.h>
#include <stdlib.h>
#include <iostream.h>
#include <math.h>
#include "fct1dfit.h"
#include "perrors.h"
#include "nbconst.h"
#include "tabmath.h"

using namespace PlanckDPC;

//define EXPO exp
#define EXPO tabFExp
#define MINEXPM (100.)

//================================================================
// CLASSES DE FONCTIONS 1D AVEC PARAMETRES POUR LE FIT
//================================================================

/////////////////////////////////////////////////////////////////
//++ 
// Module       Classes de fonctions 1D
// Lib  Outils++ 
// include      fct1dfit.h
//--
/////////////////////////////////////////////////////////////////

/////////////////////////////////////////////////////////////////
//++
// Titre        Gauss1DPol
// \index{Gauss1DPol}
//
//|  Gaussienne+polynome:
//|  Si polcenter=true: xc=(x-par[1]), sinon xc=x
//|  f(x) = par[0]*exp[-0.5*( (x-par[1]) / par[2] )**2 ]
//|        +par[3] + par[4]*xc + .... + par[3+NDegPol]*xc**NDegPol
//|  NDegPol = degre du polynome, si <0 pas de polynome
//--
/////////////////////////////////////////////////////////////////

//++
Gauss1DPol::Gauss1DPol(unsigned int ndegpol,bool polcenter)
//
//      Createur.
//--
: GeneralFunction(1,ndegpol+4), NDegPol(ndegpol), PolCenter(polcenter)
{
}

//++
Gauss1DPol::Gauss1DPol(bool polcenter)
//
//      Createur.
//--
: GeneralFunction(1,3), NDegPol(-1), PolCenter(polcenter)
{
}

Gauss1DPol::~Gauss1DPol()
{
}

double Gauss1DPol::Value(double const xp[], double const* Par)
{
double xc = (xp[0]-Par[1])/Par[2];
double e = 0.5*xc*xc;
if( e<MINEXPM ) e = EXPO(-e); else e = 0.;
double f = Par[0]*e;

if(NDegPol>=0) {
  double xcp = (PolCenter) ? Par[1] : 0.;
  double xpow = 1.;
  for(int i=0;i<=NDegPol;i++) {
    f += Par[3+i]*xpow;
    xpow *= xp[0]-xcp;
  }
}
return (f);
}

double Gauss1DPol::Val_Der(double const xp[],double const* Par
                          ,double *DgDpar)
{

 double xc = (xp[0]-Par[1])/Par[2];
 double xc2 = xc*xc;
 double e = 0.5*xc2;
 if( e<MINEXPM ) e = EXPO(-e); else e = 0.;
 double f = Par[0]*e;

 DgDpar[0] = e;
 DgDpar[1] = xc / Par[2] *f;
 DgDpar[2] = xc2/ Par[2] *f;

 if(NDegPol>=0) {
   double xcp = (PolCenter) ? Par[1] : 0.;
   double xpow = 1.;
   for(int i=0;i<=NDegPol;i++) {
     DgDpar[3+i] = xpow;
     f += Par[3+i]*xpow;
     if(PolCenter && i<NDegPol) DgDpar[2] += -(i+1)*xpow*Par[4+i];
     xpow *= xp[0]-xcp;
   }
 }
return f;
}


/////////////////////////////////////////////////////////////////
//++
// Titre        GaussN1DPol
// \index{GaussN1DPol}
//
//| Gaussienne_Normalisee+polynome (par[0]=Volume:
//|  Si polcenter=true: xc=(x-par[1]), sinon xc=x
//|  f(x) = par[0]/(sqrt(2*Pi)*par[2])*exp[-0.5*((x-par[1])/par[2])**2 ]
//|        +par[3] + par[4]*xc + .... + par[3+NDegPol]*xc**NDegPol
//|  NDegPol = degre du polynome, si <0 pas de polynome
//--
/////////////////////////////////////////////////////////////////

//++
GaussN1DPol::GaussN1DPol(unsigned int ndegpol,bool polcenter)
//
//      Createur.
//--
: GeneralFunction(1,ndegpol+4), NDegPol(ndegpol), PolCenter(polcenter)
{
}

//++
GaussN1DPol::GaussN1DPol(bool polcenter)
//
//      Createur.
//--
: GeneralFunction(1,3), NDegPol(-1), PolCenter(polcenter)
{
}

GaussN1DPol::~GaussN1DPol()
{
}

double GaussN1DPol::Value(double const xp[], double const* Par)
{
 double xc = (xp[0]-Par[1])/Par[2];
 double xc2 = xc*xc;
 double e = 0.5*xc2;
 if( e<MINEXPM ) e = EXPO(-e)/(S2Pi*Par[2]); else e = 0.;
 double f = Par[0]*e;

 if(NDegPol>=0) {
   double xcp = (PolCenter) ? Par[1] : 0.;
   double xpow = 1.;
   for(int i=0;i<=NDegPol;i++) {
     f += Par[3+i]*xpow;
     xpow *= xp[0]-xcp;
   }
 }

return (f);
}

double GaussN1DPol::Val_Der(double const xp[], double const* Par
                           , double *DgDpar)
{
 double xc = (xp[0]-Par[1])/Par[2];
 double xc2 = xc*xc;
 double e = 0.5*xc2;
 if( e<MINEXPM ) e = EXPO(-e)/(S2Pi*Par[2]); else e = 0.;
 double f = Par[0]*e;

 DgDpar[0] = e;
 DgDpar[1] = xc / Par[2] *f;
 DgDpar[2] = (xc2-1.)/ Par[2] *f;

 if(NDegPol>=0) {
   double xcp = (PolCenter) ? Par[1] : 0.;
   double xpow = 1.;
   for(int i=0;i<=NDegPol;i++) {
     DgDpar[3+i] = xpow;
     f += Par[3+i]*xpow;
     if(PolCenter && i<NDegPol) DgDpar[2] += -(i+1)*xpow*Par[4+i];
     xpow *= xp[0]-xcp;
   }
 }

 return f;
}

/////////////////////////////////////////////////////////////////
//++
// Titre        Exp1DPol
// \index{Exp1DPol}
//
//|  Exponentielle+polynome:
//|  xx = x - X_Center
//|  f(x) = exp[par[0]+par[1]*xx]
//|        +par[2] + par[3]*xx + .... + par[2+NDegPol]*xx**NDegPol
//|  NDegPol = degre du polynome, si <0 pas de polynome
//--
/////////////////////////////////////////////////////////////////

//++
Exp1DPol::Exp1DPol(unsigned int ndegpol,double x0)
//
//      Createur.
//--
: GeneralFunction(1,ndegpol+3), NDegPol(ndegpol), X_Center(x0)
{
}

//++
Exp1DPol::Exp1DPol(double x0)
//
//      Createur.
//--
: GeneralFunction(1,2), NDegPol(-1), X_Center(x0)
{
}

Exp1DPol::~Exp1DPol()
{
}

double Exp1DPol::Value(double const xp[], double const* Par)
{
double xc = Par[0]+Par[1]*(xp[0]-X_Center);
double f = ( xc>-MINEXPM ) ? EXPO(xc): 0.;

if(NDegPol>=0) {
  double xpow = 1.;
  for(int i=0;i<=NDegPol;i++) {
    f += Par[2+i]*xpow;
    xpow *= xp[0]-X_Center;
  }
}
return (f);
}

double Exp1DPol::Val_Der(double const xp[],double const* Par
                          ,double *DgDpar)
{
double xc = Par[0]+Par[1]*(xp[0]-X_Center);
double f = ( xc>-MINEXPM ) ? EXPO(xc): 0.;

DgDpar[0] = f;
DgDpar[1] = (xp[0]-X_Center) * f;

if(NDegPol>=0) {
  double xpow = 1.;
  for(int i=0;i<=NDegPol;i++) {
    DgDpar[2+i] = xpow;
    f += Par[2+i]*xpow;
    xpow *= xp[0]-X_Center;
  }
}
return f;
}

/////////////////////////////////////////////////////////////////
//++
// Titre        Polyn1D
// \index{Polyn1D}
//
//| polynome 1D:
//|  xx = x - X_Center
//|  f(x) = par[0] + par[1]*xx + .... + par[NDegPol+1]*xx**NDegPol
//|  NDegPol = degre du polynome
//--
/////////////////////////////////////////////////////////////////

//++
Polyn1D::Polyn1D(unsigned int ndegpol,double x0)
//
//      Createur.
//--
: GeneralFunction(1,ndegpol+1), NDegPol(ndegpol), X_Center(x0)
{
}

Polyn1D::~Polyn1D()
{
}

double Polyn1D::Value(double const xp[], double const* Par)
{
 double xpow = 1.;
 double f = 0.;
 for(int i=0;i<=NDegPol;i++) {
   f += Par[i]*xpow;
   xpow *= xp[0]-X_Center;
 }
 return (f);
}

double Polyn1D::Val_Der(double const xp[], double const* Par
                       , double *DgDpar)
{
 double xpow = 1.;
 double f = 0.;
 for(int i=0;i<=NDegPol;i++) {
   DgDpar[i] = xpow;
   f += Par[i]*xpow;
   xpow *= xp[0]-X_Center;
 }
 return f;
}


/////////////////////////////////////////////////////////////////
//++
// Titre        HarmonieNu
// \index{HarmonieNu}
//
//| Analyse harmonique:
//|  f(t) = par(1) + Sum[par(2k)  *cos(2*pi*k*par(0)*(t-t0)]
//|                + Sum[par(2k+1)*sin(2*pi*k*par(0)*(t-t0)]
//|  la somme Sum porte sur l'indice k qui varie de [1,NHarm()+1]
//|  avec la convention   k=1  pour le fondamental
//|                       k>1  pour le (k-1)ieme harmonique
//--
//++
//|  par(0) = inverse de la periode (frequence)
//|  par(1) = terme constant
//|  par(2), par(3) = termes devant le cosinus et le sinus
//|                   du fondamental
//|  par(2(m+1)), par(2(m+1)+1) = termes devant le cosinus
//|                   et le sinus de l'harmonique m (m>=1).
//|  NHarm() = nombre d'harmoniques a fitter.
//|  T0() = centrage des temps, ce n'est pas un parametre du fit.
//      `Conseil:' Avant de faire un fit avec les `NHarm()'
//      harmoniques, il est preferable de faire un pre-fit
//      ou seuls les parametres 1,2 et 3 sont libres et d'injecter
//      le resultat du fit sur ces 3 parametres comme valeurs
//      de depart pour le fit global avec les `NHarm()' harmoniques.
//      De toute facon, le fit ne marchera que si la periode
//      est initialisee de facon tres precise.
//--
/////////////////////////////////////////////////////////////////

//++
HarmonieNu::HarmonieNu(unsigned int nharm,double t0)
//
//      Createur.
//--
: GeneralFunction(1,4+2*nharm), NHarm(nharm), T0(t0)
{
}

HarmonieNu::~HarmonieNu()
{
}

double HarmonieNu::Value(double const xp[], double const* Par)
{
 double a = 2.*Pi*(xp[0]-T0)*Par[0];
 double f = Par[1];
 for(int k=1;k<=NHarm+1;k++) {
   // k=1 fondamental, k=2 a NHarm+1 harmoniques numero 1 a NHarm
   f += Par[2*k]*cos(k*a) + Par[2*k+1]*sin(k*a);
 }
 return f;
}

double HarmonieNu::Val_Der(double const xp[], double const* Par
                        , double *DgDpar)
{
 double cs,sn;
 double dadp0 = 2.*Pi*(xp[0]-T0);
 double a = dadp0*Par[0];
 double f = Par[1];
 DgDpar[0] = 0.;
 DgDpar[1] = 1.;
 for(int k=1;k<=NHarm+1;k++) {
   cs = cos(k*a); sn = sin(k*a);
   f += Par[2*k]*cs + Par[2*k+1]*sn;
   DgDpar[0] += k*dadp0*(-Par[2*k]*sn+Par[2*k+1]*cs);
   DgDpar[2*k]   = cs;
   DgDpar[2*k+1] = sn;
 }
 return f;
}


/////////////////////////////////////////////////////////////////
//++
// Titre        HarmonieT
// \index{HarmonieT}
//
//| Analyse harmonique:
//|  f(t) = par(1) + Sum[par(2k)  *cos(2*pi*k*(t-t0)/par(0)]
//|                + Sum[par(2k+1)*sin(2*pi*k*(t-t0)/par(0)]
//|  la somme Sum porte sur l'indice k qui varie de [1,NHarm()+1]
//|  avec la convention   k=1  pour le fondamental
//|                       k>1  pour le (k-1)ieme harmonique
//--
//++
//|  par(0) = periode
//|  par(1) = terme constant
//|  par(2), par(3) = termes devant le cosinus et le sinus
//|                   du fondamental
//|  par(2(m+1)), par(2(m+1)+1) = termes devant le cosinus
//|                   et le sinus de l'harmonique m (m>=1).
//|  NHarm() = nombre d'harmoniques a fitter.
//|  T0() = centrage des temps, ce n'est pas un parametre du fit.
//--
/////////////////////////////////////////////////////////////////

//++
HarmonieT::HarmonieT(unsigned int nharm,double t0)
//
//      Createur.
//--
: GeneralFunction(1,4+2*nharm), NHarm(nharm), T0(t0)
{
}

HarmonieT::~HarmonieT()
{
}

double HarmonieT::Value(double const xp[], double const* Par)
{
 double a = 2.*Pi*(xp[0]-T0)/Par[0];
 double f = Par[1];
 for(int k=1;k<=NHarm+1;k++) {
   // k=1 fondamental, k=2 a NHarm+1 harmoniques numero 1 a NHarm
   f += Par[2*k]*cos(k*a) + Par[2*k+1]*sin(k*a);
 }
 return f;
}

double HarmonieT::Val_Der(double const xp[], double const* Par
                        , double *DgDpar)
{
 double cs,sn;
 double dadp0 = 2.*Pi*(xp[0]-T0);
 double a = dadp0/Par[0];
 double f = Par[1];
 DgDpar[0] = 0.;
 DgDpar[1] = 1.;
 for(int k=1;k<=NHarm+1;k++) {
   cs = cos(k*a); sn = sin(k*a);
   f += Par[2*k]*cs + Par[2*k+1]*sn;
   DgDpar[0] += k*dadp0*(-Par[2*k]*sn+Par[2*k+1]*cs)/(-Par[0]*Par[0]);
   DgDpar[2*k]   = cs;
   DgDpar[2*k+1] = sn;
 }
 return f;
}

//================================================================
// CLASSES DE FONCTIONS 2D AVEC PARAMETRES POUR LE FIT
//================================================================

/////////////////////////////////////////////////////////////////
//++ 
// Module       Classes de fonctions 2D
// Lib  Outils++ 
// include      fct1dfit.h
//--
/////////////////////////////////////////////////////////////////

/////////////////////////////////////////////////////////////////
//++
// Titre        Polyn2D
// \index{Polyn2D}
//
//| polynome 2D de degre total degre:
//|  NDegPol = degre du polynome (note N dans la suite)
//|  x = x - X_Center,  y = y - Y_Center
//|  f(x,y) = p[0] +sum(k=1,n){ sum(i=0,k){ p[ki]*x^i*y^(k-i) }}
//|  Il y a k+1 termes de degre k  (ex: x^i*y^(k-i))
//|  terme de degre k avec x^i: p[ki] avec  ki = k*(k+1)/2 + i
//| C'est a dire:
//| deg0:   p0
//| deg1: + p1*y + p2*x
//| deg2: + p3*y^2 + p4*x*y + p5*x^2
//| deg3: + p6*y^3 + p7*x*y^2 + p8*x^2*y + p9*x^3
//| deg4: + p10*y^4 + p11*x*y^3 + p12*x^2*y^2 + p13*x^3*y + p14*x^4
//| deg5: + p15*y^5 + ...                       ... + ... + p20*x^5
//| ...
//| degk: + p[k*(k+1)/2]*y^k + ...           ... + p[k*(k+3)/2]*x^k
//| ...
//| degn: + p[n*(n+1)/2]*y^n + ...           ... + p[n*(n+3)/2]*x^n
//--
/////////////////////////////////////////////////////////////////

//++
Polyn2D::Polyn2D(unsigned int ndegpol,double x0,double y0)
//
//      Createur.
//--
: GeneralFunction(2,ndegpol*(ndegpol+3)/2+1), NDegPol(ndegpol), X_Center(x0), Y_Center(y0)
{
}

Polyn2D::~Polyn2D()
{
}

double Polyn2D::Value(double const xp[], double const* Par)
{
 double f = 0.;
 double xpow = 1.;
 for(int i=0;i<=NDegPol;i++) { // On gere la puissance x^i i<ndeg
   int ip;
   double ypow = 1.;
   for(int j=0;j<=NDegPol-i;j++) { // On gere la puissance x^j tq degre=i+j<=ndeg
     ip = (i+j)*(i+j+1)/2+i;
     f += Par[ip]*xpow*ypow;
     ypow *= xp[1]-Y_Center;
   }
   xpow *= xp[0]-X_Center;
 }
 return f;
}

double Polyn2D::Val_Der(double const xp[], double const* Par
                       , double *DgDpar)
{
 double f = 0.;
 double xpow = 1.;
 for(int i=0;i<=NDegPol;i++) { // On gere la puissance x^i i<ndeg
   int ip;
   double ypow = 1.;
   for(int j=0;j<=NDegPol-i;j++) { // On gere la puissance x^j tq degre=i+j<=ndeg
     ip = (i+j)*(i+j+1)/2+i;
     DgDpar[ip] = xpow*ypow;
     f += Par[ip]*DgDpar[ip];
     ypow *= xp[1]-Y_Center;
   }
   xpow *= xp[0]-X_Center;
 }
 return f;
}