source: Sophya/trunk/SophyaPI/PIext/userfitfunex.c@ 4037

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/********************************************************************/
/* fonction de fit "f(x;h,x0,sig,c) = h*exp(-0.5*((x-x0)/sig)^2)+c" */
/* avec nvar=1 npar=4                                               */
/*      x = x[0]                                                    */
/*      p[0]=h, p[1]=x0, p[2]=sig, p[3]=c                           */
/********************************************************************/
double gauss1(double const* x,double const* p)
{
double xc;
xc = (x[0]-p[1])/p[2]; xc *= -xc/2.;
return((xc>=-80.) ? p[0]*exp(xc)+p[3] : p[3]);
}

/* la definition de cette fonction n'est pas obligatoire */
/* dp[0] = df(x;h,x0,sig,c)/dh, dp[1] = etc...           */
double gauss1_der(double const* x,double const* p,double* dp)
{
double xc,xc2,e,f;
xc = (x[0]-p[1])/p[2];
xc2 = xc*xc;
e = -xc2/2.;
if(e>=-80.) e = exp(e); else e = 0.;
f = p[0]*e;
dp[0] = e;
dp[1] = xc  / p[2] *f;
dp[2] = xc2 / p[2] *f;
dp[3] = 1.;
return(f+p[3]);
}

/********************************************************/
/* fonction de fit gaussienne 2D avec correlation       */
/* f(x,y;h,x0,y0,sx,sy,rho,c) =                         */
/*   h * exp{ -0.5*[ (x-x0)^2/sx^2 + (y-y0)^2/s^2       */
/*                   -2*rho*(x-x0)*(y-y0)/(sx*sy) ]}    */
/* avec nvar=2 npar=7                                   */
/*      x = x[0] et y=x[2]                              */
/*      p[0]=h, p[1]=x0, p[2]=y0                        */
/*      p[3]=sx, p[4]=sy, p[5]=rho, p[6]=fond           */
/********************************************************/
double gauss2(double const* x,double const* p)
{
double ax,ay,a;
ax = (x[0]-p[1])/p[3];
ay = (x[1]-p[2])/p[4];
a  = -(ax*ax+ay*ay-2.*p[5]*ax*ay)/2.;
return((a>=-80.) ? p[0]*exp(a)+p[6] : p[6]);
}

/* la definition de cette fonction n'est pas obligatoire */
/* dp[0] = f(x,y;h,x0,y0,sx,sy,rho,c)/dh, dp[1] = etc... */
double gauss2_der(double const* x,double const* p,double* dp)
{
double ax,ay,a,p0ea;
ax   = (x[0]-p[1])/p[3];
ay   = (x[1]-p[2])/p[4];
a    = -(ax*ax+ay*ay-2.*p[5]*ax*ay)/2.;
p0ea = (a>=-80.) ? exp(a) : 0.;
dp[0] = p0ea;
  p0ea *= p[0];
dp[1] = (ax-p[5]*ay)/p[3] *p0ea;
dp[2] = (ay-p[5]*ax)/p[4] *p0ea;
dp[3] = ax*dp[1];
dp[4] = ay*dp[2];
dp[5] = ax*ay *p0ea;
dp[6] = 1.;
return(p0ea+p[6]);
}