source: Sophya/trunk/FrEROS/AnaLC/filtccd.cc@ 3619

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

#include "nbtri.h"
#include "nbmath.h"
#include "matxop.h"
#include "filtccd.h"

static int DEB_FitPan    = 0;
static int DEB_BaseLine  = 0;
static int DEB_Bosse     = 0;
static int DEB_BossePaG  = 0;
static int DEB_SigmaInt  = 0;
static int DEB_FiltMed   = 0;
static int DEB_FiltInt   = 0;

static int Bosse_probnb  = 0;

/*=========================================================================*/
void SetDebFilt(int razer)
{
char *c;
DEB_FitPan = DEB_BaseLine = DEB_Bosse = DEB_BossePaG  
           = DEB_SigmaInt = DEB_FiltMed = DEB_FiltInt = 0;
if( razer < 0 ) return;
if( (c = getenv("DEB_FitPan"))   ) DEB_FitPan   = atoi(c); else DEB_FitPan   = 0;
if( (c = getenv("DEB_BaseLine")) ) DEB_BaseLine = atoi(c); else DEB_BaseLine = 0;
if( (c = getenv("DEB_Bosse"))    ) DEB_Bosse    = atoi(c); else DEB_Bosse    = 0;
if( (c = getenv("DEB_BossePaG")) ) DEB_BossePaG = atoi(c); else DEB_BossePaG = 0;
if( (c = getenv("DEB_SigmaInt")) ) DEB_SigmaInt = atoi(c); else DEB_SigmaInt = 0;
if( (c = getenv("DEB_FiltMed"))  ) DEB_FiltMed  = atoi(c); else DEB_FiltMed  = 0;
if( (c = getenv("DEB_FiltInt"))  ) DEB_FiltInt  = atoi(c); else DEB_FiltInt  = 0;
}

/*=========================================================================*/
double LPacZin(double t,double *Par,double *DgDpar) 
/*                                                  Christophe 11/02/94
fonction de paczinski en -magnitude 2.5*log10(Paczin) au temps t
pour U0,T0,Tau,magfond et derivees de la fonction par rapport
aux parametres [0]=U0,[1]=T0,[2]=Tau,[3]=magfond
*/
{
double Paczin,u2,su2u24,tt0tau,dadu2;

tt0tau = (t-Par[1])/Par[2];
u2 = Par[0]*Par[0] + tt0tau*tt0tau;
su2u24 = sqrt(u2*(u2+4.));
dadu2 = -4./su2u24/su2u24/su2u24;

Paczin = (u2+2.)/su2u24;
DgDpar[0] = DftoDm/Paczin *dadu2* 2.*Par[0];
DgDpar[1] = DftoDm/Paczin *dadu2*(-2.*tt0tau/Par[2]);
DgDpar[2] = DftoDm/Paczin *dadu2*(-2.*tt0tau*tt0tau/Par[2]);
DgDpar[3] = 1.;

Paczin = 2.5*log10(Paczin) + Par[3];
return (Paczin);
}

/*=========================================================================*/
double MPacZin(double t,double *Par,double *DgDpar) 
/*                                                  rcecile 13/06/94
fonction de paczinski a l'envers en -magnitude 2.5*log10(Paczin) au temps t
pour U0,T0,Tau,magfond et derivees de la fonction par rapport
aux parametres [0]=U0,[1]=T0,[2]=Tau,[3]=magfond
*/
{
double Paczin,u2,su2u24,tt0tau,dadu2;

tt0tau = (t-Par[1])/Par[2];
u2 = Par[0]*Par[0] + tt0tau*tt0tau;
su2u24 = sqrt(u2*(u2+4.));
dadu2 = -4./su2u24/su2u24/su2u24;

Paczin = (u2+2.)/su2u24;
DgDpar[0] = DftoDm/Paczin *dadu2* 2.*Par[0];
DgDpar[1] = DftoDm/Paczin *dadu2*(-2.*tt0tau/Par[2]);
DgDpar[2] = DftoDm/Paczin *dadu2*(-2.*tt0tau*tt0tau/Par[2]);
DgDpar[3] = 1.;

Paczin = -2.5*log10(Paczin) + Par[3];
return (Paczin);
}

/*=========================================================================*/
double LPacZin2(double t,double *Par,double *DgDpar) 
/*                                                  Christophe 11/02/94
fonction de paczinski en -magnitude 2.5*log10(Paczin) au temps t
pour 1/U0,T0,1/Tau,magfond et derivees de la fonction par rapport
aux parametres [0]=1/U0,[1]=T0,[2]=1/Tau,[3]=magfond
*/
{
double Paczin,U2,U02,tt0,tt0tau,tt0tau2,dPdU2,U24p1,sU24p1,den;

tt0 = t-Par[1];
tt0tau = tt0*Par[2];
tt0tau2 = tt0tau*tt0tau;
U02 = Par[0]*Par[0];
U2 = U02/(1.+U02*tt0tau2);
U24p1 = 1.+4*U2;  sU24p1 = sqrt(U24p1);
den = 1.+U02*tt0tau2;  den *= den;
dPdU2 = 4.*U2/(U24p1*sU24p1);

Paczin = (1.+2.*U2)/sU24p1;
DgDpar[0] =  DftoDm/Paczin *dPdU2* 2.*Par[0]/den;
DgDpar[1] =  DftoDm/Paczin *dPdU2* 2.*U02*U02*Par[2]*tt0tau/den;
DgDpar[2] = -DftoDm/Paczin *dPdU2* 2.*U02*U02*tt0*tt0tau/den;
DgDpar[3] = 1.;

Paczin = 2.5*log10(Paczin) + Par[3];
return (Paczin);
}

/*=========================================================================*/
double MPacZin2(double t,double *Par,double *DgDpar) 
/*                                                  rcecile 13/06/94
fonction de paczinski a l'envers en -magnitude 2.5*log10(Paczin) au temps t
pour 1/U0,T0,1/Tau,magfond et derivees de la fonction par rapport
aux parametres [0]=1/U0,[1]=T0,[2]=1/Tau,[3]=magfond
*/
{
double Paczin,U2,U02,tt0,tt0tau,tt0tau2,dPdU2,U24p1,sU24p1,den;

tt0 = t-Par[1];
tt0tau = tt0*Par[2];
tt0tau2 = tt0tau*tt0tau;
U02 = Par[0]*Par[0];
U2 = U02/(1.+U02*tt0tau2);
U24p1 = 1.+4*U2;  sU24p1 = sqrt(U24p1);
den = 1.+U02*tt0tau2;  den *= den;
dPdU2 = 4.*U2/(U24p1*sU24p1);

Paczin = (1.+2.*U2)/sU24p1;
DgDpar[0] =  DftoDm/Paczin *dPdU2* 2.*Par[0]/den;
DgDpar[1] =  DftoDm/Paczin *dPdU2* 2.*U02*U02*Par[2]*tt0tau/den;
DgDpar[2] = -DftoDm/Paczin *dPdU2* 2.*U02*U02*tt0*tt0tau/den;
DgDpar[3] = 1.;

Paczin = -2.5*log10(Paczin) + Par[3];
return (Paczin);
}

/*=========================================================================*/
double LErrExt(double m,double *Par,double *DgDpar) 
/* fonction des erreurs externes  --   Cecile 23/02/94 */
{
double errext;
errext =exp(Par[1]*m) ;
if ( errext < 1.e-200 ) errext = 0.;
DgDpar[0] = errext;
errext *= Par[0];
DgDpar[1] = m*errext;
DgDpar[2] = 1.;
errext += Par[2];
return (errext);
}

/*=========================================================================*/
void SetBosse_probnb(int iset)
{
if( iset<=0 ) Bosse_probnb=0; else Bosse_probnb=1;
}

/*=========================================================================*/
int_4 Bosse(float *mag,float *err,int_4 ndata
         ,float mean,int_4 nptmin,float scut,int_4 nptscut
         ,float *chi2,int_4 *npoint,float *Pchi2
         ,int_4 *ideb,int_4 *ifin,int_4 *classB,int_4 NbossMX)
/*                                      Christophe 10/11/93 La Silla
**** Recherche de Bosses ****
ENTREE:
 mag: magnitude
 err: erreur sur la magnitude
 ndata: nombre de donnees
 mean: moyenne de reference pour calculer les bosses
 nptmin: nombre minimum de points pour valider une bosse
 scut, nptscut: definition des conditions d'arret de la bosse:
      Soit un point de signe oppose a plus de abs(scut*sigma)
      Soit nptscut points a moins de abs(scut*sigma)
      si scut<=0. la seule interruption est sur le changement de signe
SORTIE:
 chi2: Xi2 des bosses
 npoint: nombre de points dans la bosse
 Pchi2: probabilite de Xi2 des bosses
 ideb: indice du debut de la bosse
 ifin: indice de fin de la bosse (idem)
 classB: tableau des indices de rangement par probabilite de Xi2 croissant
 NbossMX: nombre maximum de bosses permis
 Bosse: nombre de bosse trouvees
*/
{
int_4 deb=DEB_Bosse;
int_4 ipo,ifirst,j,k,npt,i1,i2,sgn,sgnsave,nbosse,np_scut,cont_bosse,pass;
double v,av,ev,c2,sc2,pc2;

ipo = sgnsave = npt = i1 = i2 = ifirst = np_scut = 0;
nbosse = -1;
c2 = sc2= pc2 = 0.;
nptmin = ( nptmin <= 0 ) ? 1 : nptmin ;
nptscut = ( nptscut <= 0 ) ? 1 : nptscut ;
scut = ( scut <= 0. ) ? 0. : scut ;

while ( ipo < ndata ) {
  ev = err[ipo];
  v = 0.;
  if ( ev > 0. ) v = ( mag[ipo] - mean ) / ev;
  av = fabs(v);
  sgn = ( v > 0. ) ? 1 : -1 ;
  if ( deb >= 4 ) printf("ipo=%5d sgn=%2d v=%f ev=%f\n",ipo,sgn,v,ev);
  /* conditions d'increment de la fluctuation*/
  cont_bosse = 1;
  if ( sgn != sgnsave && av > scut ) cont_bosse = 0;
  if ( cont_bosse && ev > 0. ) {
    npt++;
    c2 += v*v;
    pc2 += nberfc(av/Rac2) + Log2;
    i2 = ipo;
    if( av < scut ) { np_scut++;} else { np_scut=0;}
  }
  /* conditions d'arret de la fluctuation */
  if (   ( cont_bosse == 0  && ev>0. ) || 
         ( np_scut>=nptscut && ev>0. ) || 
         ( ipo == ndata-1 )                    ) {
    /* y a t il assez de points ? */
    if ( npt >= nptmin ) {
      nbosse++;
      /* on prend prob(xi2,npt)/2**(npt-1)                            */
      if(Bosse_probnb) 
           pc2 = (double) probnb( (float) c2,npt,&pass) + (npt-1)*Log2;
      if ( nbosse >= NbossMX ) {
        printf("Bosse: Trop de bosses (%d), changez dimension tableaux",nbosse);
        exit(-1);
      } else if ( nbosse == 0 ) {
        classB[0] = 0;
      } else if ( nbosse > 0 ) {
        for ( j=0;j<nbosse;j++) 
          if ( pc2 > fabs( (double) Pchi2[classB[j]] ) ) break;
        for ( k=nbosse-1;k>=j;k--) classB[k+1] = classB[k];
        classB[j] = nbosse;
      }
      Pchi2[nbosse] = sgnsave*pc2;
      sc2 += sgnsave*c2;
      chi2[nbosse] = sgnsave*c2;
      npoint[nbosse] = npt;
      ideb[nbosse] = i1;
      ifin[nbosse] = i2;
      if(deb>=3) printf("Find bosse %5d, npt=%5d lim=%5d %5d classB=%5d pc2=%.2f\n"
                 ,nbosse,npoint[nbosse],ideb[nbosse],ifin[nbosse]
                 ,classB[nbosse],Pchi2[nbosse]);
    }
    i1 = ipo;
    npt = 1;
    sgnsave = ( v > 0. ) ? 1 : -1;
    c2 = v*v;
    pc2 = nberfc(av/Rac2) + Log2;
  }
  ipo++;
}

nbosse++;
if(deb) printf("Bosse: nombre de bosses trouvees=%d, somme des Pchi2=%f\n"
              ,nbosse,sc2);
if ( nbosse > 0 && deb ) {
  k=3;
  if ( deb >= 2 ) k=nbosse;
  for (ipo=0;ipo<k;ipo++) {
    j = classB[ipo];
    printf("Find bosse %4d, npt=%5d lim=%5d %5d pc2=%.2f\n"
          ,j,npoint[j],ideb[j],ifin[j],Pchi2[j]);
  }
}
return(nbosse);
}

/*=========================================================================*/
int BosseN(float *mag,float *err,int_4 ndata
         ,float mean,float smin,int_4 nptsmin,float scut,int_4 nptscut
         ,float *chi2,int_4 *npoint,float *Pchi2
         ,int_4 *ideb,int_4 *ifin,int_4 *classB,int_4 NbossMX)
/*                                      Christophe 27/12/95 Saclay
**** Recherche de Bosses sur courbe de lumiere ****
ENTREE:
 mag: tableau des magnitudes
 err: tableau des erreurs sur la magnitude (<=0 point non considere)
 ndata: nombre de donnees dans les tableaux (mag,err)
 mean: ligne de base par rapport a laquelle les bosses sont determinees
 smin, nptsmin: condition de validation d une bosse:
       nombre minimum de points au dessus (dessous)
       de smin*err (-smin*err) pour valider une bosse
 scut : condition de demarrage d une bosse:
       1 point > scut (ou < -scut)
 scut, nptscut: definition des conditions d'arret de la bosse:
      nptscut points a moins de  scut**err si bosse positive
      nptscut points a plus  de -scut**err si bosse negative
      Conseil: scut = 1.
 NbossMX: nombre maximum permis de bosses
SORTIE:
 chi2: tableau des Xi2 des bosses
 npoint: tableau des nombres de points dans les bosses (du meme signe que la bosse)
 Pchi2: tableau des probabilites de Xi2 des bosses
 ideb: tableau des indices de debut des bosses
 ifin: tableau des indices de fin des bosses (idem)
 classB: tableau des indices de rangement par probabilite de Xi2 croissant
        ex: pchi2 de la 3ieme bosse = Pchi2[classB[2]]
 return: Bosse: nombre de bosses trouvees
REMARQUE:
 npoints != ifin-ideb+1 car
           1-/ il y a des points a erreur <=0
           2-/ les points de signes opposes
               au signe de la bosse ne sont pas comptabilises
         c'est le nombre de points de meme signe que la bosse dans la bosse
         le nombre total de points dans la bossse est ifin-ideb+1
*/
{
int deb=DEB_Bosse;               /* 0,1,2,3 */
int j,k,ipo,npt=0,npt_cut=0,npt_c2_cut=0,i1=-1,i2=-1,sgn,sgnsave,nbosse,npt_valid,pass;
double v,av,ev,c2,pc2,c2_cut,pc2_cut;
char s;

if(smin<0.) smin = 1.;
if(nptsmin<=0) nptsmin = 1;
if(scut<0.) scut = 1.;
if(nptscut<=0) nptscut = 1;

if(deb>=1) printf("***** BosseN ndata=%d smin=%d, %f scut=%d, %f\n"
                 ,ndata,nptsmin,smin,nptscut,scut);

if( ndata <= 1 ) return(-1);

nbosse = -1;
ipo = sgnsave = 0;

for (;;) {
  /* caracteristiques du point ipo */
  ev = err[ipo];
  if ( ev > 0. ) v = ( mag[ipo] - mean ) / ev; else v = 0.;
  if ( v > 0. ) { av = v; sgn = 1; } else { av = -v; sgn = -1; }
  if (deb>=3)
      printf("     ipo=%5d sgn=%2d v=%f ev=%f c2=%f\n",ipo,sgn,v,ev,v*v);
  if ( ev>0. && sgnsave==0 && av>=scut ) {
    /* essai d'une nouvelle bosse (sgnsave==0), on demarre si > scut */
    i1 = i2 = ipo;
    sgnsave = sgn;
    c2 = v*v;
    pc2 = nberfc(av/Rac2) + Log2;
    npt = 1;
    if( av >= smin ) npt_valid = 1; else npt_valid = 0;
    c2_cut = pc2_cut = 0.;
    npt_cut = npt_c2_cut = 0;
    if(deb>=3)
      printf("* Dep bosse [%d] i1=%d n,c2,pc2= %d, %.3f, %.3f\n"
            ,sgnsave,i1,npt,c2,pc2);
  } else if( sgnsave!=0 &&
             ( (ev>0. && sgnsave*v<scut) || ipo==ndata-1 ) ) {
    /* c'est une possibilite de fin de bosse             */
    /* c a d soit un point normal en position de coupure */
    /*       soit le dernier point des tableaux          */
    if ( ev>0. && sgnsave*v<scut ) {
      /* point normal ou dernier point dans des conditions d arret de bosse */
      s = '-';
      npt_cut++;
      if( sgn == sgnsave ) {
        npt_c2_cut++;
        c2_cut += v*v;
        pc2_cut += nberfc(av/Rac2) + Log2;
      }
    } else if ( ev>0. && sgnsave*v>=scut ) {
      /* c-a-d le dernier point (ipo==ndata-1)             */
      /* avec une amplitude (grande) continuant la bosse   */
      /* --> dans ce cas la bosse doit inclure ce point    */
      /* et eventuellement les precedents de cette coupure */
      s = '+';
      npt += 1 + npt_c2_cut;
      c2 += v*v + c2_cut;
      pc2 += nberfc(av/Rac2) + Log2 + pc2_cut;
      if( av >= smin ) npt_valid++;
      i2 = ipo;
    } else {
      /* cas  ipo==ndata-1 et ev<=0. : point ignore */
      s = '0';
    }
    if(deb>=3)
       printf("- Fin%c ncut=%d n,c2,pc2= %d, %.3f, %.3f cut= %d, %.3f, %.3f\n"
             ,s,npt_cut,npt,c2,pc2,npt_c2_cut,c2_cut,pc2_cut);
    if ( npt_cut>=nptscut || ipo==ndata-1 ) {
      /* c'est une fin de bosse */
      if ( npt_valid >= nptsmin ) {
        /* validation de la bosse */
        nbosse++;
        /* on prend prob(xi2,npt)/2**(npt-1)                            */
      if(Bosse_probnb) 
           pc2 = (double) probnb((float) c2,npt,&pass) + (npt-1)*Log2;
        if ( nbosse >= NbossMX ) {
          printf("Bosse: Trop de bosses (%d), changez dimension tableaux",nbosse);
          return(-NbossMX);
        } else if ( nbosse == 0 ) {
          classB[0] = 0;
        } else if ( nbosse > 0 ) {
          for ( j=0;j<nbosse;j++) 
            if ( pc2 > fabs( (double) Pchi2[classB[j]] ) ) break;
          for ( k=nbosse-1;k>=j;k--) classB[k+1] = classB[k];
          classB[j] = nbosse;
        }
        Pchi2[nbosse] = sgnsave*pc2;
        chi2[nbosse] = sgnsave*c2;
        npoint[nbosse] = npt;
        ideb[nbosse] = i1;
        ifin[nbosse] = i2;
        if(deb>=3)
           printf("* Bosse %d, lim=%d %d classB=%d n,c2,pc2=%d %.3f %.3f\n"
                 ,nbosse,ideb[nbosse],ifin[nbosse],classB[nbosse]
                 ,npoint[nbosse],chi2[nbosse],Pchi2[nbosse]);
      } else {
        if(deb>=3) printf("* Bosse rejetee npt=%d /%d\n",npt,nptsmin);
      }
      sgnsave = 0;
      ipo = i2;
    }
  } else if ( ev>0. && sgnsave!=0 ) {
    /* c'est un point normal dans la bosse */
    npt += 1 + npt_c2_cut;
    c2 += v*v + c2_cut;
    pc2 += nberfc(av/Rac2) + Log2 + pc2_cut;
    if( av >= smin ) npt_valid++;
    i2 = ipo;
    if(deb>=3)
      printf("- add point i2=%d n,c2,pc2= %d, %.3f, %.3f cut= %d, %.3f, %.3f\n"
            ,i2,npt,c2,pc2,npt_c2_cut,c2_cut,pc2_cut);
    npt_cut = npt_c2_cut = 0;
    c2_cut = pc2_cut = 0.;
  } else {
    /* point a erreur negative ou points a moins de scut sans bosse demarree */
    if(deb>=3)
       printf("- pas de bosse demarree (%d) ou erreur negative (%f)\n"
             ,sgnsave,ev);
  }
  ipo++;
  if(ipo >= ndata ) break;
}
nbosse++;

if(deb>=1) {
  printf("Bosse: nombre de bosses trouvees=%d\n",nbosse);
  if ( nbosse>0 ) {
    printf("********************************************************\n");
    k=3;
    if ( deb>=2 ) k=nbosse;
    for (ipo=0;ipo<k;ipo++) {
      j = classB[ipo];
      printf("* bosse %4d, lim=%5d %5d n,c2,pc2=%d %.3f %.3f\n"
            ,j,ideb[j],ifin[j],npoint[j],chi2[j],Pchi2[j]);
    }
    printf("********************************************************\n");
  }
}

return(nbosse);
}

/*=========================================================================*/
int_4 BossePaG(float *temps,float *mag,float *err,int_4 ndata,float mean
            ,int_4 I1,int_4 I2,float *U0,float *T0,float *Tau)
/*                                      Christophe 11/11/93 La Silla
**** Recherche des parametres U0,T0,Tau d'une bosse ****
          bosse + ou - en magnitude ou en flux
ENTREE:
 temps: temps
 mag: magnitude
 err: erreur sur la magnitude
 ndata: nombre de donnees (si >0 mag=magnitudes, si <0 mag=flux)
 mean: valeur de la ligne de base
 I1: indice de depart de la bosse
 I2: indice de fin de la bosse
 U0: si >0 la bosse est positive, si <0 la bosse est negative
SORTIE:
 U0: estimation du parametre d'impact U0
 T0: estimation du temps de minimum d'approche T0
 Tau: estimation de la largeur Tau
 BossePaG: retourne 0 si calcul ok, <0 sinon
Methode:
 Les parametres T0 et Tau sont estimes a partir de la
 distribution en temps des points dans la bosse,
 U0 par l'estimation de la hauteur de la gaussienne de sigma
 calcule pour la valeur de Tau et de l'aire de la bosse.
*/
{
int_4 deb=DEB_BossePaG;
int_4 i,flu=0,bossp=1;
float tprev,magprev;
double v,m,s,norm,sint,hgauss,a,aa,u12,hmax;
int ifirst;

if( ndata<0 ) { flu=1; ndata *= -1;}
if( *U0<0 ) bossp = -1;
*U0 = *T0 = *Tau = 0.;
if ( I2 <= I1 ) {printf("BossePaG: I1 <= I2 %d %d\n",I1,I2); return(-1);}
if ( mean <= 0. && flu ) return(-11);
if(deb) printf("BossePaG: de I1=%d a I2=%d n=%d mean=%g flu=%d\n"
              ,I1,I2,ndata,mean,flu);

/* Calcul de T0 et Tau */
hmax= -1.e30;
tprev = temps[I1]; magprev = mean; /* $CHECK$ EA  pas cool l'init de magprev */
/* magprev peut etre une amplitude, mais pas mean !! */
/* $CHECK$ EA si on tue le point I1, sint delire. Je rajoute ifirst */
m = s = norm = sint = 0.;
ifirst = 1;
for (i=I1;i<=I2;i++) {
  if ( err[i] <= 0. ) continue;
  if(flu && bossp<0 && mag[i]<=0.) continue;
  if(flu && bossp>=0) a = mag[i]/mean-1.;
    else if(flu && bossp<0) a = 1.-mean/mag[i];
      else a = mag[i]-mean;
  aa = bossp*a;
  if( aa<=0 ) continue;
  if( aa > hmax ) hmax=aa;
  v = temps[i];
  m += v*aa;
  s += v*v*aa;
  norm += aa;
  if ( !ifirst )
    sint += (temps[i]-tprev)*(a+magprev)/2.;
  else
    ifirst = 0;
  tprev = temps[i];
  magprev = a;
}
if(norm <= 0. ) return(-2);
sint = fabs(sint);
m = m/norm;
s = s/norm - m*m;
if(s <= 0.) return(-3);
s = sqrt(s);
hgauss = sint / S2Pi / s;
if (deb>=2) printf("m=%f s=%f int=%f hgauss=%f hmax=%f\n",m,s,sint,hgauss,hmax);

/* T0 = moyenne */
*T0 = m;

/* U0 = hauteur de la gaussienne de sigma s */
if(flu) a = hgauss+1.; else a = pow(10.,0.4*hgauss);
if(a<=1. || a>= sqrt(GRAND2)) return(-4);
*U0 = 2.*(a/sqrt(a*a-1.)-1.);
if(deb>=2) printf("A0=%f u0^2=%f",a,*U0);

/* Tau = en considerant que 2.36*s est la largeur 
   a mi-hauteur (magnitude) de la courbe de pazcinski */
if(flu) a = hgauss/2.+1.; else a = pow(10.,0.2*hgauss);
if(a<1.) return(-5);
u12 = 2*(a/sqrt(a*a-1.)-1.);
if (deb>=2) printf(" a1/2=%f u1/2^2=%f\n",a,u12);
v = 2.36*s/2.;
if ( u12 <= *U0 ) return(-6);
v = v*v/(u12-*U0);
*Tau = sqrt(v);
*U0 = sqrt(*U0);
if (deb) printf("BossePaG: U0=%f, T0=%f, Tau=%f\n",*U0,*T0,*Tau);

return(0);
}

/*=========================================================================*/
float SigmaInt(float *temps,float *mag,float *err,int_4 *ndata,float nsigma)
/*                                 Christophe 12/11/93 La Silla
**** Calcul du sigma interne d'une courbe de lumiere ****
ENTREE:
 temps: temps
 mag: magnitude
 err: erreur sur la magnitude
 ndata: nombre de donnees
 nsigma: nombre de sigmas pour la coupure de la 2sd passe
         si <=0 une seule passe
SORTIE:
 SigmaInt: retourne le sigma Interne ou <0 si il y a eu un probleme
- Methode: On calcul le sigma avec la dispersion d'un point Y(t) et de la
    prediction lineaire a partir de ses 2 voisins Y0(t) et Y1(t):
    M = Y - Ypred = Y - (Y1 +(Y2-Y1)/(t2-t1)*(t-t1))
    le calcul de propagation des erreurs donne:
    SM = 2*(R^2-R+1)*SY avec SY=SY1=SY2 et R=(t-t1)/(t2-t1)  (<1.)
    donc pour avoir SY on calcul le sigma de la quantite:
    M et on divise le resultat par sqrt( <2*(R^2-R+1)> )
    ce qui dans le cas ou les points sont equidistants 
    donne comme R=1/2: SM/sqrt(3/2)
*/
{
int_4 deb=DEB_SigmaInt;
int_4 pass,passmx,i,n=0;
float sigma;
double m,s,fpred,R,sR;

passmx = ( nsigma > 0. ) ? 2 : 1 ;
if ( *ndata < 3 ) return(-1.);
sigma = GRAND;
if( nsigma <= 0. ) nsigma = 1.;
for (pass=1;pass<=passmx;pass++) {
  n=0;
  m = s = sR = 0.;
  for (i=1;i<*ndata-1;i++) {
    if ( err[i] > 0. && err[i-1] > 0. && err[i+1] > 0. ) {
      R = ( temps[i]  - temps[i-1] ) / ( temps[i+1]- temps[i-1] );
      fpred = mag[i]  -  ( mag[i-1] + (mag[i+1]- mag[i-1] )*R );
      if ( fabs(fpred) < nsigma*sigma ) {
        m += fpred;
        s += fpred*fpred;
        sR += 2.*(R*R-R+1.);
        n++;
  } } }
  if ( n <= 0 ) return(-2.);
  m = m/n;
  sR = sR/n;
  sigma = sqrt(s/n) / sqrt(sR);
  if ( deb ) printf("SigmaInt: pass=%d sigma=%f m=%f sR=%f npoints=%d\n"
                   ,pass,sigma,m,sR,n);
}
*ndata = n;
return(sigma);
}

/*=========================================================================*/
int_4 FiltMed(float *mag,float *err
             ,float *magflt, float* errflt, int_4 ndata)
/*                                   Christophe 13/11/93 La Silla
**** Filtre Statistique Median a 3 points ****
ENTREE:
 mag: magnitude
 err: erreur sur la magnitude
 ndata: nombre de donnees
SORTIE:
 magflt: magnitude filtree
 errflt: erreur sur la magnitude filtree
 FiltMed: retourne 0 si ok, <0 sinon.
*/
{
int deb=DEB_FiltMed;
int j,k0,k1,k2,imin=-1,imax=-1;

if ( ndata<3 ) {
  for (k0=0;k0<ndata;k0++) 
    { magflt[k0] = mag[k0]; errflt[k0] = err[k0];}
  return(-1);
}

/* filtre median sur tous les points sauf le premier et le dernier */
for (k0=1;k0<ndata-1;k0++) {
  magflt[k0] = mag[k0]; errflt[k0] = err[k0];
  if ( err[k0] <= 0. ) continue;
  /*recherche du 1er point inferieur utilisable */
  k1 = k0-1; while ( err[k1] <= 0. && k1 > 0       ) k1--;
   /*  recherche du 1er point superieur utilisable */
  k2 = k0+1; while ( err[k2] <= 0. && k2 < ndata-1 ) k2++;
  if ( err[k1]<=0. || err[k2]<=0. ) continue;
  if     ((mag[k1]-mag[k0])*(mag[k2]-mag[k0])<0.) j=k0;
  else if((mag[k0]-mag[k1])*(mag[k2]-mag[k1])<0.) j=k1;
  else if((mag[k1]-mag[k2])*(mag[k0]-mag[k2])<0.) j=k2;
  else {
    if     (mag[k0]==mag[k1]) {if(err[k0]<=err[k1]) j=k0; else j=k1;}
    else if(mag[k1]==mag[k2]) {if(err[k1]<=err[k2]) j=k1; else j=k2;}
    else                      {if(err[k2]< err[k0]) j=k2; else j=k0;}
  }
  magflt[k0] = mag[j]; errflt[k0] = err[j];
  if( imin == -1 ) imin = k0; imax = k0;
}

/* cas du premier point */
k0 = 0; magflt[k0] = mag[k0]; errflt[k0] = err[k0];
if( err[k0]>0. && imin>=0 )
  { magflt[k0] = magflt[imin]; errflt[k0] = errflt[imin];}

/* cas du dernier point */
k0 = ndata-1; magflt[k0] = mag[k0]; errflt[k0] = err[k0];
if( err[k0]>0. && imax>=0 )
  { magflt[k0] = magflt[imax]; errflt[k0] = errflt[imax];}

/* debug? */
if (deb>=1) {
  printf("FiltMed:");
  for(k0=0;k0<ndata+10;k0+=10) {
     for(j=k0; j<k0+10 && j<ndata; j++)
        { if(k0==j) printf("\nF "); printf(" %10.3f",mag[j]);}
     for(j=k0; j<k0+10 && j<ndata; j++)
        { if(k0==j) printf("\nE "); printf(" %10.3f",err[j]);}
     for(j=k0; j<k0+10 && j<ndata; j++)
        { if(k0==j) printf("\nFf"); printf(" %10.3f",magflt[j]);}
     for(j=k0; j<k0+10 && j<ndata; j++)
        { if(k0==j) printf("\nEf"); printf(" %10.3f",errflt[j]);}
  }
  printf("\n");
}

return(0);
}

/*=========================================================================*/
int_4 FiltInt(float *temps,float *mag,float *err,int_4 ndata,int_4 I1,int_4 I2
             ,float lbosse,int_4 npomin,int_4 *ideb,int_4 *ifin)
/*                                   Christophe 15/11/93 La Silla
**** Calcul de l'intervalle pour le fit ****
ENTREE:
 temps: temps
 mag: magnitude
 err: erreur sur la magnitude
 ndata: nombre de donnees
 I1: 1er indice de la bosse
 I2: dernier indice de la bosse
 lbosse: largeur (en unite de largeur de bosse)
         a rajouter de part et d'autre de la bosse
 npomin: nombre de points minimum avant/apres la bosse
SORTIE:
 ideb: indice du debut de l'intervalle
 ifin: indice de fin de l'intervalle
 FiltInt: retourne nombre de points dans l'intervalle et hors bosse,
          <0 si probleme.
*/
{
int_4 deb=DEB_FiltInt;
int_4 i,n,ntot;
float t1,t2;
double dt;

if ( I2 <= I1 ) {printf("FiltInt: I2 <= I1 %d %d\n",I1,I2); return(-1);}

ntot = 0;
t1 = temps[I1];
t2 = temps[I2];
lbosse = (lbosse<=0.) ? 1. : lbosse ;
dt = (t2-t1)*lbosse;
if(deb>1) printf("FiltInt: limites bosse %f %f (%f)\n",t1,t2,dt);
t1 -= dt;
t2 += dt;
if(deb>1) printf("FiltInt: limites intervalle de fit %f %f\n",t1,t2);

/* calcul de la borne inferieure de fit */
n=0;
for (i=I1;i>=0;i--) {
  if(err[i] > 0. ) {
    n++;
    if( temps[i] < t1 && n > npomin ) break;
} }
ntot += n;
*ideb= ( i<0 ) ? 0 : i;

/* calcul de la borne superieure de fit */
n=0;
for (i=I2;i<ndata;i++) {
  if(err[i] > 0. ) {
    n++;
    if( temps[i] > t2 && n > npomin ) break;
} }
ntot += n;
*ifin= ( i>= ndata ) ? ndata-1 : i;

if(deb) printf("FiltInt: limites intervalle de fit %d %d\n",*ideb,*ifin);

return(ntot);
}

/*======================================================================================*/
double chi2_exterieur(float *mag,float *err,int_4 ndata,float mean
                     ,int_4 I1,int_4 I2,int_4 *npt)
/*                         Cecile
calcul du Chi2 par degre de liberte a l'exterieur de la bosse delimitee par i1 et I2
ENTREE :
 mag : magnitude
 err : erreur sur la magnitude
 ndata : nombre de donnees
 mean : moyenne des mesures
 I1, I2 : limites de la bosse
SORTIE:
 chi2_exterieur : valeur du Chi2 reduit a l'exterieur de la bosse
*/
{
int_4 i;
double chi2,d;

if ( I2 <= I1 ) {printf("chi2_ext. : I2 <= I1 %d %d\n",I1,I2); return(-1.);}

*npt=0;
chi2=0.;

if(I1>0) {
  for(i=0;i<I1;i++) {
    if( err[i] > 0.) {
      d = (mag[i] - mean) / err[i];
      chi2+= d*d;
      (*npt)++;
} } }

if(I2<ndata) {
  for(i=I2+1;i<ndata;i++) {
    if( err[i] > 0.) {
      d = (mag[i] - mean) / err[i];
      chi2+= d*d;
      (*npt)++;
} } }

return(chi2);
}

/*===================================================================*/
double Dt_PacZin(double u0,double tau,double A)
/* pour calculer le DT ou la paczinsky a une ampli de A */
{
double x;
if( A <= 1. ) return(-1.);
x = 2.*( A / sqrt(A*A-1.) - 1. ) - u0*u0;
if( x <= 0. ) return(-2.);
x = tau * sqrt(x);
return(x);
}

/*==================================================================*/
int_4 Most_3D2D(float *Y,float *EY,int_4 *ndata,float Ymin,float Ymax
               ,float *Res_2D,float *Res_3D)
/*
Pour calculer la valeur la + probable par estimation des parametres
  d'une distribution parabolique / polynome de degre 3
INPUT:
 Y: tableau des valeurs
 EY: tableau des erreurs (<=0 point non traite)
 ndata: nombre de valeurs
 Ymin,Ymax: limites des valeurs a considerees pour Y
OUTPUT:
 ndata: nombre de valeurs utilisees pour le calcul
 Res_2D: resultat du fit parabolique
 Res_3D: resultat du fit polynome d'ordre 3
 rc: validation des resultats
       =0 : probleme
       =1 : fit parabolique ok
       =2 : fit polynome_ordre_3 ok
       =3 : fit parabolique + polynome_ordre_3 ok
*/
{
int i,n,Flag,lp=0;
double xmoy,x2moy,x3moy;
double v1,v2,VN,det;
double M3[3][3],M2[2][2],B[3],AX[4];

if(lp) printf("Most_3D2D: ndata=%d Ymin=%g Ymax=%g\n",*ndata,Ymin,Ymax);

*Res_3D = 0.;
*Res_2D = 0.;
Flag = 0;

if( *ndata<=0 || Ymax <= Ymin ) goto END;

/* on recentre tout entre 0 et 1 donc Ymin=0 YMax=1 */
VN = Ymax-Ymin;
n = 0;
xmoy = x2moy = x3moy = 0.;

for(i=0;i<*ndata;i++)
  {
  if( EY[i]>0. && Y[i]>Ymin && Y[i]<Ymax )
    {
    v1 = (Y[i]-Ymin)/VN;
    xmoy += v1;
    x2moy += v1*v1;
    x3moy += v1*v1*v1;
    n++;
    }
  }
if(n<=0) goto END;
xmoy /= n;
x2moy /= n;
x3moy /= n;
if(lp) printf("    <x>=%g <x2>=%g <x3>=%g\n",xmoy,x2moy,x3moy);

/* On remplit la matrice 2x2 pour le fit parabolique */
M2[0][0] = 1./5. -1./9.;
M2[0][1] = M2[1][0] = M2[1][1] = 1./12.; /* = 1./4. -1./6. = 1./3. -1./4. */
B[0] = x2moy - 1./3.;
B[1] = xmoy - 0.5;
det = SolveDLinSyst(&(M2[0][0]), B, AX, 2);
if(lp) printf("Solve_2D: det=%g sol=%g %g\n",det,AX[0],AX[1]);
if ( det != 0. ) 
  {
  if ( AX[0] != 0. )
    {
    Flag += 1;
    *Res_2D = -AX[1]/(2.*AX[0]);
    if(lp) printf("       s1=%g\n",*Res_2D);
    *Res_2D = Ymin + *Res_2D*VN;
    }
  }

/* On remplit la matrice 3x3 pour le fit polynome d'ordre 3 */
M3[0][0] = M3[2][2] = 1./5. -1./8.;
/* 1./12. = 1./4. -1./6. = 1./3. -1./4. = 1./6. -1./12 */
M3[0][1] = M3[1][2] = M3[0][2] = M3[1][0] = M3[2][1] = 1./12.;
M3[1][1] = 1./5. -1./9.;
M3[2][0] = 1./7. -1./16.;
B[0] = xmoy - 0.5;
B[1] = x2moy - 1./3.;
B[2] = x3moy - 0.25;
det = SolveDLinSyst(&(M3[0][0]), B, AX, 3);
if(lp) printf("Solve_3D: det=%g sol=%g %g %g\n",det,AX[0],AX[1],AX[2]);
if ( det != 0. ) 
  {
  /* on cherche les 2 racines de l'annulation de la derivee */
  det = AX[1]*AX[1] - 3.*AX[0]*AX[2];
  if(lp) printf("       det=%g\n",det);
  if ( det >= 0. && AX[0] != 0. )
    {
    xmoy =  (-AX[1] + sqrt(det))/(3.*AX[0]);
    v1 = 6.*AX[0]*xmoy+2.*AX[1];
    x2moy = (-AX[1] - sqrt(det))/(3.*AX[0]);
    v2 = 6.*AX[0]*x2moy+2.*AX[1];
    if(lp) printf("       s1=%g (%g) s2=%g (%g)\n",xmoy,v1,x2moy,v2);
    Flag += 2;
    *Res_3D = xmoy;
    if(v2<0.) *Res_3D = x2moy;
    *Res_3D = Ymin + *Res_3D*VN;
    }
  }


END:
if(lp) printf("Most_3D2D: *Res_2D=%g *Res_3D=%g Flag=%d\n"
             ,*Res_2D,*Res_3D,Flag);
return(Flag);
}

/*==================================================================*/
double correlation(int_4 npts,double *x,double *y)
{
/*                        Cecile    Numerical Recipies p. 503
Calcule le coefficient de correlation de x et y sur npts points
ENTREE :
 npts nombre de points
 x, y tableau des valeurs
SORTIE:
 correlation : coefficient de correlation
               -2 si nb de points OK < 3
               -3 si produit des (x[i]-ax)(y[i]-ay)<=0
*/
int_4 i,n;
double ax,ay,xt,yt,sxx,syy,sxy;
n=0; ay = ax= 0.;
for(i=0;i<npts;i++) {ax += x[i]; ay += y[i]; n++;}
if( n<3 ) {
  printf("correlation: moins de 3 points pour le calcul de la correlation (%d)\n",n);
  return(-2.);
}
ax /= n; ay /= n;
xt = yt = syy = sxy = sxx = 0.;
for(i=0;i<npts;i++) {
  xt = x[i]-ax;  yt = y[i]-ay;
  sxx += xt*xt;  syy += yt*yt;  sxy += xt*yt;
} 
return( ((sxx*syy>0.) ? sxy/sqrt(sxx*syy) : -3.) );
}

/*==================================================================*/
double CorrelatioN(int_4 nptX,float *X,float *Y,int_4* IndXvY)
/*                                 cmv 17/4/99
Calcule le coefficient de correlation de X[] et Y[] sur nptX points
ENTREE :
  nptX nombre de points dans le tableau X[]
  X[], Y[] : tableau des valeurs a correler
  IndXvY[] : tableau des indices de Y[] a correler a X[]
SORTIE:  coefficient de correlation
         -2 si nb de points OK (>3)
         -3 si somme des (X[i]-ax)^2 ou (Y[i]-ay)^2 <=0
EXPLICATION:
  X[1,nptX-1], Y[....], IndXvY[1,nptX-1]
  - le point X[i] correspond au point Y[IndXvY[i]] si IndXvY[i]>=0
  - le point X[i] n'est pas pris en compte si IndXvY[i]<0
*/
{
int_4 i,n;
double ax,ay,xt,yt,sxx,syy,sxy;
n=0; ay=ax=0.;
for(i=0;i<nptX;i++) if(IndXvY[i]>=0) {ax+=X[i]; ay+=Y[IndXvY[i]]; n++;}
if(n<3)
  {printf("CorrelatioN: %d<3 points pour le calcul\n",n); return(-2.);}
ax /= n; ay /= n;
xt = yt = syy = sxy = sxx = 0.;
for(i=0;i<nptX;i++) if(IndXvY[i]>=0)
  { xt = X[i]-ax; yt = Y[IndXvY[i]]-ay;
    sxx += xt*xt; syy += yt*yt;  sxy += xt*yt; }
return( ((sxx*syy>0.) ? sxy/sqrt(sxx*syy) : -3.) );
}

/*=========================================================================*/
int_4 MeanLine(float *mag,float *err,int_4 *ndata
              ,int passmx,float nsigma,float *mean,float *sigma)
/*
Calcul de la moyenne et du sigma en plusieurs passes avec nettoyage
ENTREE:
 mag: magnitude
 err: erreur sur la magnitude (Si NULL pas d'erreur)
 ndata: nombre de donnees
 nsigma: nombre de sigma dans lequel on calcule le mean
 passmx: nombre de passes
SORTIE:
 ndata: nombre de data utilises pour calculer le most
 mean: moyenne des donnees
 sigma: sigma des donnees
 MeanLine: 0 si calcul ok, <0 sinon
*/
{
int_4 pass,n,i;
double m,s2,scut,v;
if( passmx<=0 ) passmx = 1;

*mean = *sigma = 0.;
for (pass=1;pass<=passmx;pass++) {
  m = s2 = 0.; n=0;
  if( pass>=2 ) scut=*sigma*nsigma; else scut=GRAND;
  for (i=0;i<*ndata;i++) {
    v = mag[i];
    if(err!=NULL) if( err[i] <= 0. ) continue;
    if( fabs(v-*mean) >= scut ) continue;
    n++;
    m += v;
    s2 += v * v;
  }
  if ( n>0 ) {
    m /= (double) n;
    v = s2/n - m*m;
    *mean = m;
    *sigma = (v>0) ? sqrt(v) : 0.;
  } else {
    *mean = *sigma = 0.;
    *ndata=n;
    return(-1);
  }
  /* printf("MeanLine: pass=%d mean=%f sigma=%f n=%d\n"
           ,pass,*mean,*sigma,n); */
}
return(0);
}

/*=========================================================================*/
int_4 BaseLine(float *mag,float *err,int_4 *ndata,float nsigma
            ,float binini,float *mean,float *sigma,float *most)
/*                  Christophe 10/11/93 La Silla + modif 28/7/94 Saclay
**** calcul de la ligne de base d'une etoile par maxi d'histogramme ****
Methode: 1-/ moyenne, 2-/ moyenne tronquee,
         3-/ maximum d'histogramme,
         4-/ puis deuxieme passe sur histo plus lache
ENTREE:
 mag: magnitude
 err: erreur sur la magnitude
 ndata: nombre de donnees
 nsigma: nombre de sigma dans lequel on calcule le maximum
         si <=0, une seule passe pour la moyenne et dynamique histo = 5.
 binini: bin minimum de la methode d'histogramme
SORTIE:
 ndata: nombre de data utilises pour calculer le most
 mean: moyenne des donnees
 sigma: sigma des donnees
 most: maximum de probabilite
RETURN: nombre d entrees dans le bin du maximum, <0 si echec
*/
{
int_4 deb=DEB_BaseLine;
int_4 pass,passmx,i,ib0,ib1,n=0,ibin,nbin=0,tab0mx,nbinmx=100,tab0[101];
float mini,maxi,binw;
double m,s2,v,vv,scut;
double SX4,SX3,SX2,SX,S1,SYX2,SYX,SY;
double MX33[3][3],BB[3],AX[3],det;

if(deb) printf("BaseLine: n=%d nsigma=%f binini=%f\n",*ndata,nsigma,binini);

/*********** D'abord une moyenne et sigma ************/
*mean = *sigma = 0.;
passmx = ( nsigma <= 0. ) ? 1 : 2 ;
for (pass=1;pass<=passmx;pass++) {
  m = s2 = 0.; n=0;
  scut = GRAND;
  if( pass == 2 ) scut = nsigma* *sigma;
  for (i=0;i<*ndata;i++) {
    v = mag[i];
    if( err[i] <= 0. || fabs(v-*mean) >= scut ) continue;
    n++;  m += v;  s2 += v*v;
  }
  if ( n >= 2 ) {
    *mean = m / n;
    v = s2 / n - m/n * m/n;
    *sigma = ( v>0. ) ? sqrt(v) : 0.;
  } else {
    *mean = *sigma = 0.;
    *ndata=n;
    return(-1);
  }
  if(deb) printf("BaseLine: pass=%d mean=%f sigma=%f n=%d\n"
                ,pass,*mean,*sigma,n);
  if( *sigma == 0. ) return(-1);
}

/***** Puis une methode de maximum d'histogramme *****/
if ( nsigma <= 0. ) nsigma = 5.;
mini = *mean - nsigma * *sigma;
maxi = *mean + nsigma * *sigma;

for (pass=1;pass<=2;pass++) {

  if( binini > 0. ) { 
    if ( pass == 1 ) nbin = (maxi-mini)/(2.*binini) + 0.5; /* 1ere pass a 2*binini */
    if ( pass == 2 ) nbin = (maxi-mini)/binini + 0.5;      /* 2sd  pass a   binini */
    if ( nbin >= nbinmx ) nbin = nbinmx;
    if ( nbin < 3 ) nbin = 3;
  } else nbin = nbinmx;
  binw = (maxi-mini)/nbin;
  if(deb) printf("mini=%f maxi=%f binw=%f nbin=%d\n",mini,maxi,binw,nbin);

  /*fill histogramme tab0*/
  for (i=0;i<nbin;i++) tab0[i] = 0;
  for (i=0;i<*ndata;i++) {
    if( err[i] <= 0. ) continue;
    ibin = (mag[i]-mini)/binw;
    if( ibin >= 0 && ibin < nbin ) tab0[ibin] += 1;
  }

  /*maximum de l'histogramme tab0*/
  ibin=0;
  tab0mx = 0;
  for (i=0;i<nbin;i++) if ( tab0[i] > tab0mx ) {tab0mx = tab0[i]; ibin=i;}
  ib0 = ibin-1; ib1 = ibin+1;
  if(ib0<0) { ib0 = 0; ib1++;}
  if(ib1>=nbin) { ib1 = nbin-1; ib0--;}
  if (deb>1) {
    printf("Maximum tab0=%d (%d [%d,%d])\n",tab0mx,ibin,ib0,ib1);
    if (deb>2) { printf("Tableau final tab0:\n");
               for (i=0;i<nbin;i++) printf(" %d=%d",i,tab0[i]);
               printf("\n"); }
    }

  /* calcul de la valeur maximale par fit parabole +/- 1 bin autour du maximum */
  S1 = SY = SYX = SYX2 = SX4 = SX3 = SX2 = SX = 0.;
  for (i=ib0;i<=ib1;i++) {
    v = mini + (i+0.5)*binw;
    vv = v*v;
    m = tab0[i];
    s2=1.;
    S1 += s2;
    SX += s2*v;
    SY += s2*m;
    SYX += s2*v*m;
    SYX2 += s2*m*vv;
    SX2 += s2*vv;
    SX3 += s2*vv*v;
    SX4 += s2*vv*vv;
  }
  MX33[0][0] = S1;
  MX33[1][1] = SX2;
  MX33[2][2] = SX4;
  MX33[0][1] = MX33[1][0] = SX;
  MX33[0][2] = MX33[2][0] = SX2;
  MX33[1][2] = MX33[2][1] = SX3;
  BB[0] = SY;
  BB[1] = SYX;
  BB[2] = SYX2;
  det = SolveDLinSyst(&(MX33[0][0]), BB, AX, 3);
  if( det == 0. ) return(-2);
  if( AX[2] == 0. ) return(-2);
  *most = -AX[1]/(2.*AX[2]);

  if(deb)
    printf("BaseLine: pass=%d most_probable=%f ibin=%d\n",pass,*most,ibin);

  if( *most<=mini || *most >= maxi ) return(-3);

  /*nouveaux parametres de l'histogramme tab0 */
  mini   = *most-2.5*binw;
  maxi   = *most+2.5*binw;

}

*ndata=n;
return(tab0mx);
}

/*=========================================================================*/
int BaseLineS(float *y,float *ey,int ny,int nfiltre,float *mean,float *sigma)
/*    D'apres Alain M. + modifs + code cmv 5/5/97
Determination de la ligne de base comme sequence de nfiltre points
  consecutifs de dispersion minimale:
  y : tableau des flux
  ey : tableau des erreurs sur les flux (pt non util. si <=0)
  ny : nombre de points
  nfiltre : nombre de points consecutifs pour calculer le sigma
            si <0 y est d'abord range par flux croissants
  mean : valeur du fond de ciel retourne
  sigma : valeur du sigma mini trouve
  rc : nombre d essais ( si <=0 echec)
*/
{
int rc=0, isort=0,i,j,npt,n;
double s,m;
float *ys=NULL;

*mean = *sigma = -1.;
if( ny<=1 ) return(-1);

if( nfiltre<0 ) {
  isort = 1;
  nfiltre *= -1;
  npt = 0;
  ys = (float *) malloc(ny*sizeof(float));
  if( ys==NULL ) return(-2);
  for(i=0;i<ny;i++) if(ey[i]>0.) ys[npt++] = y[i];
  qsort(ys,(size_t) npt,(size_t) sizeof(float),qSort_Float);
} else { npt=ny; ys = y;}
if( nfiltre>npt || nfiltre<=1 ) { if(isort) free(ys); return(-3);}

for(i=0;i<npt-nfiltre+1;i++) {
  if( !isort && ey[i]<=0. ) continue;
  s = m = 0.;  n = 0;
  for(j=i;j<npt;j++) {
    if( !isort && ey[j]<=0. ) continue;
    m += ys[j];  s += ys[j]*ys[j];  n++;
    if( n<nfiltre ) continue;
    rc++;
    m /= (double) n;
    s = s/n - m*m;
    s = ( s>0. ) ? sqrt(s) : 0.;
    if( s<*sigma || rc==1 ) { *sigma=s; *mean=m; }
    break;
  }
}

if(isort) free(ys);
return(rc);
}

/*=========================================================================*/
int BaseLineP(float *y,float *ey,int ny,int nfiltre,float *mean)
/*    D'apres Alain M. + modifs + code cmv 5/5/97
Determination de la ligne de base comme endroit de pente minimale
  entre le point [i] et [i+nfiltre-1]. Ceci marche mieux avec
  le tableau trie mais la possibilite de ne pas trier est laissee:
  y : tableau des flux
  ey : tableau des erreurs sur les flux (pt non util. si <=0)
  ny : nombre de points
  nfiltre : nombre de points pour calculer la pente
            si <0 y est d'abord range par flux croissants
  mean : valeur du fond de ciel retourne
  rc : nombre d essais ( si <=0 echec)
*/
{
int rc=0, isort=0,i,j,npt,n;
float p=-1.,pmin=-1.,m;
float *ys=NULL;

*mean = -1.;
if( ny<=1 ) return(-1);

if( nfiltre<0 ) {
  isort = 1;
  nfiltre *= -1;
  npt = 0;
  ys = (float *) malloc(ny*sizeof(float));
  if( ys==NULL ) return(-2);
  for(i=0;i<ny;i++) if(ey[i]>0.) ys[npt++] = y[i];
  qsort(ys,(size_t) npt,(size_t) sizeof(float),qSort_Float);
} else { npt=ny; ys = y;}
if( nfiltre>npt || nfiltre<=1 ) { if(isort) free(ys); return(-3);}

for(i=0;i<npt-nfiltre+1;i++) {
  if( !isort && ey[i]<=0. ) continue;
  m = 0.;  n = 0;
  for(j=i;j<npt;j++) {
    if( !isort && ey[j]<=0. ) continue;
    m += ys[j]; n++;
    p = ys[j];
    if( n<nfiltre ) continue;
    rc++;
    p = fabs(ys[i]-p);
    if( p<pmin || rc==1 ) { pmin = p; *mean = m/n;}
    break;
  }
}

if(isort) free(ys);
return(rc);
}