source: Sophya/trunk/SophyaLib/NTools/nbmath.c@ 2882

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

#include <math.h>
#include "nbmath.h"
#include "matxop.h"
#include "nbinteg.h"
#include "nbtri.h"

#define ITMAX 256
#define EPS 3.0e-7
#define DEB_GausPiv 0
#define DEB_MeanSig 0

/*
++ 
  Module        Fonction mathematiques (C)
  Lib   LibsUtil
  include       nbmath.h

        Fonction mathematiques (C)
--
*/

void FitFun_MrqCof( double (*FunFit) (double ,double *,double *) ,int_4 type
           ,FLOATDOUBLE temps,FLOATDOUBLE mag,FLOATDOUBLE err
           ,int_4 *ndata,int_4 I1,int_4 I2
           ,double *parcur,int_4 npar,int_4 *ind
           ,double *ATGA,double *BETA,double *ci2,int_4 deb);
void nbgcf(double *gammcf,double a,double x,double *gln);
void nbgser(double *gamser,double a,double x,double *gln);


static int FITFUN_DPOL = -1;


/*=========================================================================*/
/*                                      Christophe 8/11/93 La Silla */
/*
++
double FitFun
|        ( double (*FunFit) (double ,double *,double *) ,int_4 type
|        ,FLOATDOUBLE temps,FLOATDOUBLE mag,FLOATDOUBLE err
|        ,int_4 *ndata,int_4 I1,int_4 I2
|        ,double *par,double *epar,double *stepar
|        ,double *minpar,double *maxpar,int_4 npar
|        ,double stochi2,int_4 NstepMX,int_4 deb)

        Pour fitter une fonction parametrique sur des donnees.
        (Obsolete, utilisez GeneralFit en C++).
--
*/
/*
++
| Fonction de fit de la courbe FunFit sur les donnees
| - Methode:   fit des moindres carres dans le cas non lineaire
| - Reference: Statistical and Computational Methods in Data Analysis
|              Siegmund Brandt, North-Holland 1970  p 204-206.
|              Introduction des limites pour la variation des parametres (cmv).
|              Increment des parametres selon la methode de Levenberg-Marquardt
|              (Numerical Receipes in C p 542)
| - Remarques diverses:
|     Les points avec erreurs negatives ou nulles ne sont pas
|                utilises dans le fit.
|     Si le fit est hautement non lineaire, il faut donner de bonnes
|                approximations pour les parametres de depart.
|     Les limites minpar,maxpar sont des limites absolues mais possibles
| - Memo personnel:
|     data = Ti, Fi ERRi  i=0,n-1
|     parametres Pj  j=1,npar
|     fonction FunFit(t,P1,P2,..,Pnpar)
|     matrice A(n,4) = [Aij] = - [dg(T)/dPj]  pour t=Ti
|     matrice G(n,n) = [Gij] = matrice diagonale 1/ERRj**2
|     matrice C(n) = [Ci] = (Fi - FunFit(t,P1,P2,..,Pnpar))
|     Pj(iteration n) = Pj(iteration n-1) + Zj
--
*/
/*
++
| -----------
| | ENTREES |
| -----------
| - FunFit  fonction a fiter
| - type    type des donnees a fiter (union temps,mag,err):
|           float=FITFUN_FLOAT / double=FITFUN_DOUBLE
| - temps:  union FLOATDOUBLE des tableaux des temps:
| - mag:    union FLOATDOUBLE des tableaux des magnitudes:
| - err:    union FLOATDOUBLE des tableaux des erreurs sur la magnitude
| - ndata:  nombre de donnees en entree, donc en c:
| - I1:     indice de depart pour le fit
| - I2:     indice de fin pour le fit
|           (si I2 <= I1 defaut: 1,*ndata)
| - par:    contient les valeurs d'initialisation des parametres
| - stepar: si >0. on fit le parametre, sinon parametre fixe
| - minpar: les valeurs minimum pouvant prendre les parametres
| - maxpar: les valeurs minimum pouvant prendre les parametres
|           (si maxpar[i] <= minpar[i] pas de limitation sur le parametre i)
| - stochi2:la valeur de variation du chi2 d'une iteration l'autre
|           en dessous de laquelle on arrete d'iterer:
|           (si <=0 le default est 0.01).
| - npar:   nombre de parametres du fit *par --> *(par+npar-1)
| - NstepMX:nombre maximum d'iterations
| - deb:    niveau de print (0/1/2/3)
--
*/
/*
++
| ----------
| | SORTIES |
| ----------
| - ndata: nombre de donnees utilisees dans le fit
| - par: contient les valeurs fitees des parametres
| - epar: contient les erreurs sur les valeurs fitees des parametres
| - FitFun: la function elle meme retourne le Xi2 du fit si succes
|       -1. si l'inversion de la matrice des erreurs n'a pas ete possible
|       -2. si les parametres initiaux sont en desaccord avec les limites
|       -3. si le nombre de donnees est inferieur au nombre de parametres
|       -4. si le fit n'a pas converge (nstep>nNstepMX)
--
*/
/*
++
| --------------
| | MEMO PERSO |
| --------------
| la structure FLOATDOUBLE:
| union FloatDouble {
|   float  *fx;
|   double *dx;
| };
| - Si les tableaux temps,mag,emag sont des float alors
|   temps.fx doit contenir l'adresse du tableau des temps
|            (et idem pour mag.fx et err.fx)
| - Si les tableaux temps,mag,emag sont des double alors
|   temps.dx doit contenir l'adresse du tableau des temps
|            (et idem pour mag.dx et err.dx)
--
*/
double FitFun( double (*FunFit) (double ,double *,double *) ,int_4 type
             ,FLOATDOUBLE temps,FLOATDOUBLE mag,FLOATDOUBLE err
             ,int_4 *ndata,int_4 I1,int_4 I2
             ,double *par,double *epar,double *stepar
             ,double *minpar,double *maxpar,int_4 npar
             ,double stochi2,int_4 NstepMX,int_4 deb)
{
int_4 ind[NPAR_FITFUN], i, j, k, nstep, nparsave, nstop, nstopMX;
double LimitC2,lambda,eps;
double ATGA[NPAR_FITFUN][NPAR_FITFUN], COVAR[NPAR_FITFUN][NPAR_FITFUN]
     , ALPHA[NPAR_FITFUN][NPAR_FITFUN];
double DA[NPAR_FITFUN], BETA[NPAR_FITFUN];
double parcur[NPAR_FITFUN],partry[NPAR_FITFUN];
double ci2, ci20;

/* initialisation */
if( type != FITFUN_FLOAT && type != FITFUN_DOUBLE) {
  printf("FitFun: Type de variable non prevue dans union FLOATDOUBLE %d\n",type);
  exit(-1);
}

if(npar>NPAR_FITFUN) {
  printf("FitFun: trop de parametres npar=%d > NPAR_FITFUN=%d\n",npar,NPAR_FITFUN);
  exit(-1);
}
nstopMX = 3;
eps = 1.e-8;
LimitC2=0.01;
if ( stochi2 > 0. ) LimitC2 = stochi2;      /* test des valeurs d'arret */
nstop = 0;

/* fit du centre autour des valeurs estimees */
nparsave = npar;
for (i=0;i<nparsave;i++) parcur[i] = *(par+i);

npar=0;
for (i=0;i<nparsave;i++) {
  ind[i]= -1;
  *(epar+i) = 0.;
  partry[i] = *(parcur+i);
  if ( stepar[i] > 0. ) {
    ind[npar]=i;
    npar++;
} }

if ( I2 <= I1 ) { I1 = 0; I2 = *ndata-1;}
if ( I2-I1+1 < npar ) {
  printf("le nombre de donnees %d est inferieur au nombre de parametres %d\n"
        ,I2-I1+1,npar);
  return(-3.);
}
 
for (j=0;j<npar;j++) {
  i = ind[j];
  if(*(maxpar+i)>*(minpar+i) && (*(parcur+i)<=*(minpar+i) || *(parcur+i)>=*(maxpar+i))) {
    printf("Parametre %d initialise hors limites: %f # [%f,%f]\n"
          ,i,*(parcur+i),*(minpar+i),*(maxpar+i));
    return(-2.);
} }

if ( deb >= 2 ) {
  printf("\n******************* ENTREE DANS FitFun ");
  printf("npar=%d LimitC2=%f\n",npar,LimitC2);
  printf("parametres"); for(i=0;i<nparsave;i++) printf(" %f",*(parcur+i)); printf("\n");
  if ( deb >= 3 ) {
    printf("       ind"); for(i=0;i<nparsave;i++) printf(" %12d",*(ind+i)); printf("\n");
    printf("    stepar"); for(i=0;i<nparsave;i++) printf(" %12.5g",*(stepar+i)); printf("\n");
    printf("    minpar"); for(i=0;i<nparsave;i++) printf(" %12.5g",*(minpar+i)); printf("\n");
    printf("    maxpar"); for(i=0;i<nparsave;i++) printf(" %12.5g",*(maxpar+i)); printf("\n");
} }

/* premiere essai avec les parametres d'initialisation */

FitFun_MrqCof(FunFit,type,temps,mag,err,ndata,I1,I2,parcur,npar,ind,ALPHA[0],BETA,&ci2,deb);

ci20 = ci2;
lambda = 0.001;

if ( *ndata < npar ) {
  printf("le nombre de donnees %d est inferieur au nombre de parametres %d\n",*ndata,npar);
  return(-3.);
}

if ( deb >= 2 ) printf("Passage d'initialisation ci2= %e, lambda=%e\n",ci2,lambda);

/* et maintenant les iterations */

nstep = 0;
/************************ ITERATIONS *******************************************/
ITERATE:
nstep++;

if ( deb >= 2 ) printf("------------------------- pas %d\n",nstep);

if(nstep > NstepMX) {
  printf("FitFun: Le fit n'a pas converge (trop d'iterations %d).\n",nstep);
  if (deb>0 ) {
    printf("== Parametres finals:");
    for(k=0;k<nparsave;k++) printf(" %15.8e",parcur[k]); printf("\n");
    for(i=0;i<npar;i++) { *(epar+ind[i]) = -1.; par[ind[i]] = parcur[ind[i]];}
  }
  return(-4.);
}

for(j=0;j<npar;j++) {
  for(i=0;i<npar;i++) ATGA[j][i] = ALPHA[j][i];
  ATGA[j][j] *= 1.0 + lambda;
}

/* inversion de   ( At * G * A ) */
if( GausPiv(ATGA[0],NPAR_FITFUN,npar,COVAR[0],NPAR_FITFUN,npar,1) == 0. ) {
  printf("La matrice des erreurs n'est pas inversible\n");
  return(-1.);
}
if( deb >= 3 ) {
  printf("matrice 1/( At * G * A )\n");
  for ( i=0; i<npar; i++ ) {
    for ( j=0; j<npar; j++ ) printf(" %15.8e",COVAR[i][j]);
    printf("\n");
} }

/* compute next increment */
for(i=0;i<npar;i++) {
  DA[i] = 0.;
  for(j=0;j<npar;j++) DA[i] += COVAR[i][j] * BETA[j];
}
if ( deb >= 2 ) {
  printf("Correction parametres proposee:");
  for (i=0;i<npar;i++) printf(" %12.5e",DA[i]); printf("\n");
}

/* on s'arrete la, il y a eu convergence */
if ( lambda == 0. ) {
  for(i=0;i<npar;i++) { *(epar+ind[i]) = sqrt(COVAR[i][i]); par[ind[i]] = parcur[ind[i]];}
  if (deb>0 ) {
    printf("===============================================================\n");
    printf("== ci2= %15.8e (%15.8e) ndata= %d nstep= %d\n",ci2,ci2/(*ndata-npar+1),*ndata,nstep);
    printf("== Parametres finals:"); for(k=0;k<nparsave;k++) printf(" %15.8e", par[k]); printf("\n");
    printf("==               +/- "); for(k=0;k<nparsave;k++) printf(" %15.8e",epar[k]); printf("\n");
    printf("===============================================================\n");
  }
  return( (float) ci2 );
}

/* essai des nouveaux parametres */
for(j=0;j<npar;j++) {
  k = ind[j];
  if( nstep == 1 && fabs(DA[j]) > stepar[k] ) {
    partry[k] = parcur[k] + stepar[k] * DA[j] / fabs(DA[j]);
  } else {
    partry[k] = parcur[k] + DA[j];
  }
  if(maxpar[k] > minpar[k] ) {
    if( partry[k] < minpar[k] ) {
      partry[k] = minpar[k];
      if(deb>=2) printf("Parametre %3d limite au minimum %e\n",k,partry[k]);
    }
    if( partry[k] > maxpar[k] ) {
      partry[k] = maxpar[k];
      if(deb>=2) printf("Parametre %3d limite au maximum %e\n",k,partry[k]);
    }
} }
if ( deb>=2)
  {
  printf("essai avec");
  for(k=0;k<nparsave;k++) printf(" %15.8e", partry[k]);
  printf("\n");
  }

/* calcul du nouveau ci2 avec partry */
FitFun_MrqCof(FunFit,type,temps,mag,err,ndata,I1,I2,partry,npar,ind,ATGA[0],DA,&ci2,deb);
if(deb>=2) printf("Ci2: old=%e new=%e %e\n",ci20,ci2,ci2-ci20);

/* test sur l'evolution du ci2 et la strategie a suivre */
if ( ci2 < ci20 ) {
  nstop = 0;
  for (j=0;j<npar;j++) {
    for(i=0;i<npar;i++) ALPHA[j][i] = ATGA[j][i];
    BETA[j] = DA[j];
    parcur[ind[j]] = partry[ind[j]];
  }
  lambda *= 0.1;
  if( ci20 - ci2 < LimitC2 ) lambda = 0.;
  ci20 = ci2;
  if(deb>=2) {
    printf("lambda divided by 10. %e\n",lambda);
    printf("Nouveaux parametres"); for(k=0;k<nparsave;k++) printf(" %12.5e",parcur[k]); printf("\n");
  }
} else {
  /* on arrete si pas assez de variation de chi2 + pas assez de variation de parametres */
  if( ci2 - ci20 < LimitC2 ) {
    k=0;
    for(i=0;i<npar;i++) if( fabs( parcur[ind[i]] - partry[ind[i]] ) < eps ) k++;
    if(k == npar) { nstop++; } else { nstop=0;};
    /* printf("ci2-ci20=%f npar=%d (%d) nstop=%d\n"
           ,ci2-ci20,npar,k,nstop); */
  } else {
    nstop = 0;
  }
  lambda *= 10.;
  if(nstop>=nstopMX) lambda = 0.;
  ci2 = ci20;
  if(deb>=2) printf("Echec essai, lambda multiplied by 10. %e\n",lambda);
}

goto ITERATE;

}

/*=========================================================================*/
void FitFun_MrqCof(double (*FunFit) (double ,double *,double *) ,int_4 type
           ,FLOATDOUBLE temps,FLOATDOUBLE mag,FLOATDOUBLE err
           ,int_4 *ndata,int_4 I1,int_4 I2
           ,double *parcur,int_4 npar,int_4 *ind
           ,double *ATGA,double *BETA,double *ci2,int_4 deb)
{
int_4 i, j, k;
double deriv[NPAR_FITFUN], t, f, e, Gkk, Ck, funfit;

for (i=0; i<npar;i++) { BETA[i] = 0.; for (j=0; j<npar;j++) *(ATGA+i*NPAR_FITFUN+j) = 0.;}
*ci2 = 0.;
*ndata = 0;

/* calcul de  ( At * G * A ) */
for (k=I1; k<=I2; k++) {
  e = (type == FITFUN_FLOAT) ? *(err.fx + k) : *(err.dx + k);
  if ( e > 0. ) {
    (*ndata)++;
    t = (type == FITFUN_FLOAT) ? *(temps.fx + k) : *(temps.dx + k);
    f = (type == FITFUN_FLOAT) ? *(mag.fx   + k) : *(mag.dx + k);
    funfit = FunFit(t,parcur,deriv);
    Gkk = 1./e/e;
    Ck = f - funfit;
    *ci2 += Ck*Ck*Gkk;
    for ( j=0; j<npar; j++ ) {
      for ( i=j; i<npar; i++ ) *(ATGA+i*NPAR_FITFUN+j) += deriv[ind[i]] * Gkk * deriv[ind[j]];
      BETA[j] += deriv[ind[j]] * Gkk * Ck;
} } }

for(j=1;j<npar;j++) for(i=0;i<j;i++) *(ATGA+i*NPAR_FITFUN+j) = *(ATGA+j*NPAR_FITFUN+i);

if( deb >= 3 ) {
  printf("matrice ( At * G * A )\n");
  for ( i=0; i<npar; i++ ) {
    for ( j=0; j<npar; j++ ) printf(" %15.8e",*(ATGA+i*NPAR_FITFUN+j));
    printf("\n");
  }
  printf("BETA:");
  for ( j=0; j<npar; j++ ) printf(" %15.8e",BETA[j]);
  printf("\n");
}

}

/*=========================================================================*/
/*                                          Christophe 8/11/93 La Silla */
/*
++
double GausPiv(double *A,int_4 nca,int_4 n,double *B,int_4 ncb,int_4 m,int_4 Inv)
        Inversion de matrice par la methode du pivot de Gauss.
--
*/
/*
++
| Inv=0: resolution d'un systeme A(n,n) * X(n,m) = B(n,m)  
|        en entree matrices A,B de taille A[?>=n][nca>=n] et B[?>=n][ncb>=m]
|        en sortie la matrice B contient le resultat X
| Inv#0: inversion de la matrice A(n,n)
|        en entree matrice A
|        en sortie la matrice B contient le resultat 1/A 
| GausPiv: retourne le determinant (0. si non-inversible!)
--
*/
/*
++
|Remarque: 
|-matrice A[?>=n][nca>=n],
|       element Aij = "A'[i][j] de ss matrice A'[n][n]" = *(A+i*nca+j)
|       ou 0<=i<n indice de ligne; ou 0<=j<n indice de colonne; 
|-matrice B[?>=n][ncb>=m],
|       element Bij = "B'[i][j] de ss matrice B'[n][m]" = *(B+i*ncb+j)
|       ou 0<=i<n indice de ligne; ou 0<=j<m indice de colonne;
|exemple 1: on veut resoudre le systeme A(4,4) * X = B(4,3)
|           et on definit dans le prog appelant: double A[12][7],B[32][5]
|           on appelle: GausPiv(&A[0][0],7,4,&B[0][0],5,3,0)
|exemple 2: on veut resoudre le systeme A(4,4) * X = B(4)
|           et on definit dans le prog appelant: double A[4][4],B[4]
|           on appelle: GausPiv(A[0],4,4,B,1,1,0)
|exemple 3: on veut resoudre le systeme A(4,4) * X = B(4)
|           et on definit dans le prog appelant: double A[9][4],B[7,3]
|           on appelle: GausPiv(A[0],4,4,B,3,1,0)
|Methode: methode de gauss avec choix du pivot maximum
|Reference: Statistical and Computational Methods in Data Analysis
|           Siegmund Brandt, North-Holland 1970  p 350-361.
|Attention: la matrice A est detruite et la matrice B aussi (evidemment)
--
*/
double GausPiv(double *A,int_4 nca,int_4 n,double *B,int_4 ncb,int_4 m,int_4 Inv)
{
int_4 deb=DEB_GausPiv;
int_4 i,j,k,l,i1,i2,j1,j2,j3,j4,ij,ii,il,lj,kj,k1,nj,kk,ik,nn,kmax,i1max;
double amax,save,det;

if ( n < 1 ) {
  printf("GausPiv: pas de resolution de systeme de dimension %d\n",n);
  return(0.);
}

/* ********************************************************** */
/* cas d'une inversion de matrice: on met B = matrice identite*/
/* ********************************************************** */
if ( Inv != 0 ) {
  if ( m != n ) {
    printf("GausPiv: Vous voulez inverser une matrice non carree %d %d\n",n,m);
    return(0.);
  }
  for (i=0;i<n;i++) {
    for (j=0;j<n;j++) {
      if ( i != j ) { *(B+i*ncb+j) = 0.; } else { *(B+i*ncb+j) = 1.;}
} } }

if ( deb>0 ) {
  printf("** Matrices de depart\n");
  printf("Matrice A %d %d:\n",nca,n);
  for (i=0;i<n;i++) {
    for (j=0;j<n;j++) printf(" %15.8e",*(A+i*nca+j));
    printf("\n");
  }
  printf("Matrice B %d %d:\n",ncb,m);
  for (i=0;i<n;i++) {
    for (j=0;j<m;j++) printf(" %15.8e",*(B+i*ncb+j));
    printf("\n");
} }

/* ***************************************************************** */
/*    reduction de la matrice A a la forme triangulaire superieure   */
/* ***************************************************************** */
det=1.;
kmax = n-1;
for (k=0;k<kmax;k++) {         /* on travaille sur la ligne k */
  if ( deb>0 ) printf("-----> reduction ligne %d\n",k);
  amax=0.;    /*recherche du plus grand coeff pour le pivot*/
  j2=k;       /* dans col k pour les lignes >= k */
  for (j1=k;j1<n;j1++) {
    ik = j1*nca+k;
    if( fabs(amax)-fabs(*(A+ik)) < 0. ) {
      amax = *(A+ik);
      j2 = j1;
  } }
  if ( deb>0 ) printf("pivot sur la ligne %d max= %f\n",j2,amax);
  if ( j2 != k ) {      /* echange des lignes k et j2 si necessaire*/
    for (j=k;j<n;j++) {
      j3=k*nca+j;
      j4=j2*nca+j;
      save = *(A+j3);
      *(A+j3) = *(A+j4);
      *(A+j4) = save;
    }
    for (j=0;j<m;j++) {
      j3=k*ncb+j;
      j4=j2*ncb+j;
      save = *(B+j3);
      *(B+j3) = *(B+j4);
      *(B+j4) = save;
    }
    if ( deb>0 ) {
      printf("Matrice A echange des lignes %d et %d :\n",k,j2);
      for (i=0;i<n;i++) {
        for (j=0;j<n;j++) printf(" %15.8e",*(A+i*nca+j));
        printf("\n");
      }
      printf("Matrice B echange des lignes %d et %d :\n",k,j2);
      for (i=0;i<n;i++) {
        for (j=0;j<m;j++) printf(" %15.8e",*(B+i*ncb+j));
        printf("\n");
  } } }
  k1=k+1;  /* on effectue la substitution-soustraction sur les lignes qui suivent*/
  kk = k*nca+k;      /* c'est le pivot*/
  if ( deb>0 ) printf("valeur du pivot %f\n",*(A+kk));
  if ( *(A+kk) == 0. ) {
    /* printf("Matrice non-inversible: pivot %d nul\n",k); */
    return(0.);
  } else { det *= *(A+kk); }
  for (i=k1;i<n;i++) {
    ik = i*nca+k;  /* c est le premier element de la ligne i */
    for (j=k1;j<n;j++) {
      ij= i*nca+j;
      kj= k*nca+j;
      *(A+ij) += -*(A+kj) * *(A+ik) / *(A+kk);
    }
    for (j=0;j<m;j++) {
      ij= i*ncb+j;
      kj= k*ncb+j;
      *(B+ij) += -*(B+kj) * *(A+ik) / *(A+kk);
  } }
  if ( deb>0 ) {
    printf("Matrice A:\n");
    for (i=0;i<n;i++) {
      for (j=0;j<n;j++) printf(" %15.8e",*(A+i*nca+j));
      printf("\n");
    }
    printf("Matrice B:\n");
    for (i=0;i<n;i++) {
      for (j=0;j<m;j++) printf(" %15.8e",*(B+i*ncb+j));
      printf("\n");
} } }     /* fin du travail sur la ligne k */
det *= *(A+(n-1)*nca+(n-1));
if ( deb>0 ) printf("Determinant de A= %e\n",det);

/* ********************************************************** */
/*       back substitution                                    */
/* ********************************************************** */
nn=(n-1)*nca+(n-1);     /* l'element (n,n) */
for (j=0;j<m;j++) {
  nj=(n-1)*ncb+j;
  if ( *(A+nn) == 0. ) {
    /* printf("Matrice non-inversible pivot (%d,%d) nul\n",n-1,n-1); */
    return(0.);
  }
  *(B+nj) /= *(A+nn);   /* la derniere ligne traitee a part*/
  i1max=n-1;
  if( i1max >= 0 ) {    /* les lignes n-1 a 1 dans cet ordre*/
    for(i1=1;i1<=i1max;i1++) {
      i=n-1-i1;         /* ici i va de n-2 a 0 */
      ij=i*ncb+j;
      ii=i*nca+i;
      i2=i+1;           /* pour la ligne i on somme de i+1 a n */
      for (l=i2;l<n;l++) {
        il=i*nca+l;
        lj=l*ncb+j;
        *(B+ij) -= *(A+il) * *(B+lj);
      }
      if ( *(A+ii) == 0. ) {
        /* printf("Matrice non-inversible pivot (%d,%d) nul\n",i,i); */
        return(0.);
      }
      *(B+ij) /= *(A+ii);
} } }
if ( deb>0 ) {
  printf("-----> Matrice B resultat:\n");
  for (i=0;i<n;i++) {
    for (j=0;j<m;j++) printf(" %15.8e",*(B+i*ncb+j));
    printf("\n");
} }
return(det);
}

/*=========================================================================*/
/*                                       CMV 17/11/92  */
/*
++
int_4 paramga -
  (float sxin, float syin, float rhoin -
  , float *smax, float *axisrat, float *tiltdeg)

        Pour transformer le sx,sy,rho d'une gaussienne en sa,sc/sa,teta.
--
*/
/*
++
| - input : sx ]0,inf] ,sy ]0,inf] ,rho ]-1,+1[
| - output : smax (longeur du grand axe), axisrat (rapport des axes ]0.,1.])
|     tiltdeg (angle entre le grand axe et l'axe ox en degre ]-90,+90])
| - Les fonctions sont:
| exp[-0.5*{[(x/SX)**2-2*RHO/(SX*SY)*x*y+(y/SY)**2]/(1-RHO**2)}]
| exp[-0.5*{ (A/SA)**2                  +(C/SC)**2            }]
| - Remarque: sx,sy representent les sigmas des variables x,y et
| sont differents des sigmas des coupes x=0 et y=0 qui valent:
|    sigmaX(y=0) = sx*sqrt(1-ro^2) different de sx
|    sigmaY(x=0) = sy*sqrt(1-ro^2) different de sy
--
*/
int_4 paramga(float sxin, float syin, float rhoin
             , float *smax, float *axisrat, float *tiltdeg)
/*
Cf Formulae and Methods in Experimental Data Evaluation Vol. 1
rechercher a Bivariate gaussian distribution
*/
{
double sx,sy,rho,a,b,c,alpha,coscar,sincar,cossin,sa,sc;

*smax = 0.; *axisrat = 0.; *tiltdeg = 0.; a = 0.; b = 0.; c = 0.;

/* lecture des arguments */
sx = sxin;  sy = syin;  rho = rhoin;

/* cas des coniques degenerees ou erreurs dans les parametres*/
if ( sx <= 0. || sy <= 0. || rho > 1. || rho < -1. ) return(-1);

/* sortie fit : [(x/SX)**2 - 2*RHO/(SX*SY)*x*y + (y/SY)**2]/(1-RHO**2) = 1
   equation ellipse: A*x**2 + 2*B*x*y + C*y**2  = 1       */
sa = (1.-rho*rho);
a = 1./sx/sx/sa;  b = -rho/sx/sy/sa;  c = 1./sy/sy/sa;

/* cas des coniques degenerees (droites limites) */
if ( rho == 1. || rho == -1. ) return(-2);

/* axes principaux OX OY pour : x = cos(ALPHA)*X - sin(APLHA)*Y
        (ALPHA = (ox,oX))       y = sin(ALPHA)*X + cos(APLHA)*Y
  soit   ( A*cos**2 + C*sin**2 + 2*B*cos*sin )*X**2
     +   ( A*sin**2 + C*cos**2 - 2*B*cos*sin )*Y**2
     + 2*( -A*cos*sin + C*cos*sin + B*(cos**2-sin**2) )*X*Y
                   = (X/SA)**2 + (Y/SC)**2 = 1
  et oX oY axes principaux pour : Tan ( 2*ALPHA ) = 2*B/(A-C) */

if ( a == c && rho == 0. ) {
  /*  cas d'un cercle */
  coscar = 0.5;    sincar = coscar;    cossin = coscar;    alpha = 0.;
} else {
  if ( a == c ) {
    /* cas SX=SY ellipse a +/-45 degres */
    if ( rho > 0. ) {
      /* cas +45 degres */
      alpha = Pi / 4.;
    } else {
      /* cas -45 degres */
      alpha = -Pi / 4.;
    }
  } else {
    /* cas general */
    alpha = atan(2.*b/(a-c)) /2.;
  }
  coscar = cos(alpha);
  sincar = sin(alpha);
  cossin = coscar*sincar;
  coscar *= coscar;
  sincar *= sincar;
  alpha *= 180./Pi;
}

sa = a*coscar+c*sincar+2.*b*cossin;
sc = c*coscar+a*sincar-2.*b*cossin;
sa = ( sa < 0. ) ? 1./sqrt(-sa) : 1./sqrt(sa);
sc = ( sc < 0. ) ? 1./sqrt(-sc) : 1./sqrt(sc);

if ( sa >= sc ) {
  *smax = sa;
  *axisrat = sc/sa;
  *tiltdeg = alpha;
} else {
  *smax = sc;
  *axisrat = sa/sc;
  *tiltdeg = alpha+90.;
  if ( *tiltdeg > 90. ) *tiltdeg = *tiltdeg-180.;
}

return(0);
}

/*========================================================================*/
/*                               CMV 21/01/94 */
/*
++
int_4 gaparam -
  (float smax, float axisrat, float tiltdeg -
  , float *sxin, float *syin, float *rhoin)

        Pour transformer le sa,sc/sa,tiltdeg d'une gaussienne en sx,sy,rho
        Pour calculer les parametres d'une gaussienne
--
*/
/*
++
| - input : sa ]0,inf] ,axisrat ]0,1] ,tiltdeg ]-90,+90]
|         sa (longeur du grand axe), axisrat (rapport des axes )
|         tiltdeg (angle entre le grand axe et l'axe ox en degre)
| - output : sx ]0,inf] , sy ]0,inf] , rho ]-1,1[
| - Les fonctions sont:
| exp[-0.5*{ (A/SA)**2                  +(C/SC)**2            }]
| exp[-0.5*{[(x/SX)**2-2*RHO/(SX*SY)*x*y+(y/SY)**2]/(1-RHO**2)}]
| - Voir aussi la remarque de la fonction paramga.
--
*/
int_4 gaparam(float smax, float axisrat, float tiltdeg
             , float *sxin, float *syin, float *rhoin)
{
double sx,sy,rho,alpha,coscar,sincar,cossin,sa,sc,a,b,c;

*sxin = 0.; *syin = 0.; *rhoin = 0.;

/* cas des coniques degenerees ou erreurs dans les parametres*/
if ( smax <= 0. || axisrat <= 0. || axisrat > 1. ) return(-1);

/* lecture des arguments */
sa = smax;  sc = sa*axisrat;  alpha = tiltdeg*Pi/180.;

coscar = cos(alpha);
sincar = sin(alpha);
cossin = coscar*sincar;
coscar *= coscar;
sincar *= sincar;

a = 1./(coscar/sa/sa + sincar/sc/sc);
  a = sqrt(a);         /* = sqrt(1-ro**2)*sx */
b = 1./(sincar/sa/sa + coscar/sc/sc);
  b = sqrt(b);         /* = sqrt(1-ro**2)*sy */
c = cossin*(-1./sa/sa + 1./sc/sc);
rho = c*a*b;
sx = a / sqrt(1.-rho*rho);
sy = b / sqrt(1.-rho*rho);

*sxin = sx; *syin = sy; *rhoin = rho;

return(0);
}

/*===================================================================*/
/*
++
double nberfc(double x)
        Erreur fonction complementaire: -log10(erfc)
        Pratiquement pour avoir la proba que
        ]-inf,-x] U [x,+inf[ selon une loi gaussiene centree et normee
        il faut appeler nberfc(x/sqrt(2))
        Le job tiens compte de la saturation machine pour x grand.
        La valeur est approximee au developpement limite de erfc
        ( see . Abramowitz p299 7.1.23 )
--
*/
double nberfc(double x)
{
double v;
if(x<0.) return(0.);
if ( x < 10. ) {
  v = erfc(x);
  v = -log10(v);
} else {
  v = 2.*x*x;
  v = 1./v - 3./(v*v) + 15./(v*v*v);
  v += LnPi/2. + log(x) + x*x;
  v /= Ln10;
}
return(v);
}

/*=========================================================================*/
/*
++
float probnb(float sci2,int_4 inddl,int_4 *ipass)
        Retourne -log10(p) ou p est la probabilte que le chi2
        d'une fluctuation inddl points soit superieur a sci2.
        ipass est le type de calcul utilise.
--
*/
float probnb(float sci2,int_4 inddl,int_4 *ipass)
{
int_4 it;
float nddl;
double ci2,a1,a2,p,q,val,terme1,terme2,terme3,terme4,teradd,xnprim;

*ipass = 0;
val = 0.;
nddl = inddl;
ci2 = sci2;
a1 = nddl/2.;
a2 = ci2/2.;

q = nbgammq(a1,a2);
if(q>0. && q<1.) {
  val = -log10(q);   /* moins le log10 de la proba que le chi2 soit > ci2 */
  *ipass = 1;
} else if ( q==1.) {
  p = nbgammp(a1,a2);
  val = p/Ln10;      /* -log10(q) = -log(1-p)/Ln10 = -(p)/Ln10 = p/Ln10 */
  *ipass = 2;
} else {
  /* approximation de mimile */
  xnprim = nddl / 2. - 1.;
  terme1 = nddl / 2. * Log2;        /* 1/2**n */
  terme2 = Hln2pi + (xnprim+0.5)*log(xnprim) - xnprim + 1./(12.*xnprim);
  terme2 /= Ln10;         /* formule de stirling pour Gamma[nddl/2] */
  terme3 = ci2/2. - xnprim*log(ci2);
  terme3 /= Ln10;         /* exp(-ci2/2)*ci2**(n/2-1) */
  terme4 = teradd = 1.;
  for ( it=1; it <= inddl-2; it +=2) {
    teradd *= ( (double) (nddl -1.) - (double) it ) / ci2;
    terme4 += teradd;
    if ( teradd < 1.e-4 ) break;
  }
  terme4 = -log10(2.*terme4);
  val = terme1 + terme2 + terme3 + terme4;
  *ipass = 3;
}
return ((float) val);
}

/*=========================================================================*/
/*
++
double nbgammln(double xx)
        Fonction log de la fonction Gamma
--
*/
double nbgammln(double xx)
{
   double x,tmp,ser;
   static double cof[6]={76.18009173,-86.50532033,24.01409822,
      -1.231739516,0.120858003e-2,-0.536382e-5};
   int_4 j;

   x=xx-1.0;
   tmp=x+5.5;
   tmp -= (x+0.5)*log(tmp);
   ser=1.0;
   for (j=0;j<=5;j++) {
      x += 1.0;
      ser += cof[j]/x;
   }
   return -tmp+log(2.50662827465*ser);
}

/*=========================================================================*/
/*
++
double nbgammq(double a,double x)
        Fonction GammaQ
--
*/
double nbgammq(double a,double x)
{
   double gamser,gammcf,gln;

   if (x < 0.0 || a <= 0.0) { printf("Invalid arguments in routine NBGAMMQ (%f,%f)\n"
                                    ,a,x); exit(-1);}
   if (x < (a+1.0)) {
      nbgser(&gamser,a,x,&gln);
      return 1.0-gamser;
   } else {
      nbgcf(&gammcf,a,x,&gln);
      return gammcf;
   }
}

/*=========================================================================*/
/*
++
double nbgammp(double a,double x)
        Fonction GammaP
--
*/
double nbgammp(double a,double x)
{
   double gamser,gammcf,gln;

   if (x < 0.0 || a <= 0.0) { printf("Invalid arguments in routine NBGAMMP (%f,%f)\n"
                                    ,a,x); exit(-1);}
   if (x < (a+1.0)) {
      nbgser(&gamser,a,x,&gln);
      return gamser;
   } else {
      nbgcf(&gammcf,a,x,&gln);
      return 1.0-gammcf;
   }
}

/*=========================================================================*/
void nbgcf(double *gammcf,double a,double x,double *gln)
{
   int_4 n;
   double gold=0.0,g,fac=1.0,b1=1.0;
   double b0=0.0,anf,ana,an,a1,a0=1.0;

   *gln=nbgammln(a);
   a1=x;
   for (n=1;n<=ITMAX;n++) {
      an=(double) n;
      ana=an-a;
      a0=(a1+a0*ana)*fac;
      b0=(b1+b0*ana)*fac;
      anf=an*fac;
      a1=x*a0+anf*a1;
      b1=x*b0+anf*b1;
      if (a1) {
         fac=1.0/a1;
         g=b1*fac;
         if (fabs((g-gold)/g) < EPS) {
            *gammcf=exp(-x+a*log(x)-(*gln))*g;
            return;
         }
         gold=g;
      }
   }
   printf("a too large, ITMAX too small in routine NBGCF (%f,%f)\n",a,x); exit(-1);
}

/*=========================================================================*/
void nbgser(double *gamser,double a,double x,double *gln)
{
   int_4 n;
   double sum,del,ap;

   *gln=nbgammln(a);
   if (x <= 0.0) {
      if (x < 0.0) { printf("x less than 0 in routine NBGSER (%f,%f)\n",a,x); exit(-1); }
      *gamser=0.0;
      return;
   } else {
      ap=a;
      del=sum=1.0/a;
      for (n=1;n<=ITMAX;n++) {
         ap += 1.0;
         del *= x/ap;
         sum += del;
         if (fabs(del) < fabs(sum)*EPS) {
            *gamser=sum*exp(-x+a*log(x)-(*gln));
            return;
         }
      }
      printf("a too large, ITMAX too small in routine NBGSER (%f,%f)\n",a,x); exit(-1);
      return;
   }
}

/*==========================================================================*/
/*
++
void Set_Ihoq(int degre,int *mIhoqN,double *IhoqNX,double *IhoqNW)
        Set mIhoqN,IhoqNX,IhoqNW pour integration Gauss-Legendre
--
*/
void Set_Ihoq(int degre,int *mIhoqN,double *IhoqNX,double *IhoqNW)
{
switch (degre) {
  case  2:
         *mIhoqN = mIhoq2;
          IhoqNX = Ihoq2X;
          IhoqNW = Ihoq2W;
          break;
  case  3:
         *mIhoqN = mIhoq3;
          IhoqNX = Ihoq3X;
          IhoqNW = Ihoq3W;
          break;
  case  4:
         *mIhoqN = mIhoq4;
          IhoqNX = Ihoq4X;
          IhoqNW = Ihoq4W;
          break;
  case  5:
         *mIhoqN = mIhoq5;
          IhoqNX = Ihoq5X;
          IhoqNW = Ihoq5W;
          break;
  case  6:
         *mIhoqN = mIhoq6;
          IhoqNX = Ihoq6X;
          IhoqNW = Ihoq6W;
          break;
  case  7:
         *mIhoqN = mIhoq7;
          IhoqNX = Ihoq7X;
          IhoqNW = Ihoq7W;
          break;
  case  8:
         *mIhoqN = mIhoq8;
          IhoqNX = Ihoq8X;
          IhoqNW = Ihoq8W;
          break;
  case  9:
         *mIhoqN = mIhoq9;
          IhoqNX = Ihoq9X;
          IhoqNW = Ihoq9W;
          break;
  case 10:
         *mIhoqN = mIhoq10;
          IhoqNX = Ihoq10X;
          IhoqNW = Ihoq10W;
          break;
  case 16:
         *mIhoqN = mIhoq16;
          IhoqNX = Ihoq16X;
          IhoqNW = Ihoq16W;
          break;
  default:
         *mIhoqN = mIhoq8;
          IhoqNX = Ihoq8X;
          IhoqNW = Ihoq8W;
          break;
}
return;
}

/*==========================================================================*/
/*
++
int Get_Ihoq(int degre,double *ihoqnx,double *ihoqnw)
        Ecrit les positions et les poids pour l'integration de Gauss-Legendre
        dans les tableaux *ihoqnx,*ihoqnw.
        Retourne le nombre de points.
        Si *ihoqnx==NULL ou *ihoqnw==NULL
        retourne seulement le nombre de points (cad la taille minimale
        que doivent les tableaux *ihoqnx,*ihoqnw).
--
*/
int Get_Ihoq(int degre,double *ihoqnx,double *ihoqnw)
{
int mIhoqN, i;
double *IhoqNX, *IhoqNW;
switch (degre) {
  case  2:
          mIhoqN = mIhoq2;
          IhoqNX = Ihoq2X;
          IhoqNW = Ihoq2W;
          break;
  case  3:
          mIhoqN = mIhoq3;
          IhoqNX = Ihoq3X;
          IhoqNW = Ihoq3W;
          break;
  case  4:
          mIhoqN = mIhoq4;
          IhoqNX = Ihoq4X;
          IhoqNW = Ihoq4W;
          break;
  case  5:
          mIhoqN = mIhoq5;
          IhoqNX = Ihoq5X;
          IhoqNW = Ihoq5W;
          break;
  case  6:
          mIhoqN = mIhoq6;
          IhoqNX = Ihoq6X;
          IhoqNW = Ihoq6W;
          break;
  case  7:
          mIhoqN = mIhoq7;
          IhoqNX = Ihoq7X;
          IhoqNW = Ihoq7W;
          break;
  case  8:
          mIhoqN = mIhoq8;
          IhoqNX = Ihoq8X;
          IhoqNW = Ihoq8W;
          break;
  case  9:
          mIhoqN = mIhoq9;
          IhoqNX = Ihoq9X;
          IhoqNW = Ihoq9W;
          break;
  case 10:
          mIhoqN = mIhoq10;
          IhoqNX = Ihoq10X;
          IhoqNW = Ihoq10W;
          break;
  case 16:
          mIhoqN = mIhoq16;
          IhoqNX = Ihoq16X;
          IhoqNW = Ihoq16W;
          break;
  default:
          mIhoqN = mIhoq8;
          IhoqNX = Ihoq8X;
          IhoqNW = Ihoq8W;
          break;
}
if( ihoqnx!=NULL && ihoqnw!=NULL ) {
  for(i=0;i<mIhoqN;i++) {
    ihoqnx[i] = IhoqNX[i];
    ihoqnw[i] = IhoqNW[i];
  }
}
return(mIhoqN);
}

/*===================================================================*/
#define EPS_gauleg 3.0e-11
/*
++
void nbgauleg(double x1,double x2,double *x,double *w,int n)
        Integration des fonctions a une dimension y=f(x)
        par la Methode de Gauss-Legendre.
        Calcul des coefficients du developpement
--
*/
/*
++
| INPUT:
|  x1,x2 : bornes de l'intervalle (dans nbinteg.h -> x1=-0.5 x2=0.5)
|  n = degre du developpement
| OUTPUT:
|  x[] = valeur des abscisses ou l'on calcule (dim=n)
|  w[] = valeur des coefficients associes
| REMARQUES:
|  - x et w doivent au moins etre dimensionner a n.
|  - l'integration est rigoureuse si sur l'intervalle d'integration
|    la fonction f(x) peut etre approximee par un polynome
|    de degre 2*m (monome le + haut x**(2*n-1) )
|  - Voir la fonction Integ_Fun pour un calcul d'ordre 8
--
*/
void nbgauleg(double x1,double x2,double *x,double *w,int n)
{
   int m,j,i;
   double z1,z,xm,xl,pp,p3,p2,p1;

   m=(n+1)/2;
   xm=0.5*(x2+x1);
   xl=0.5*(x2-x1);
   for (i=1;i<=m;i++)  {
      z=cos(3.141592654*(i-0.25)/(n+0.5));
      do {
         p1=1.0;
         p2=0.0;
         for (j=1;j<=n;j++) {
            p3=p2;
            p2=p1;
            p1=((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j;
         }
         pp=n*(z*p1-p2)/(z*z-1.0);
         z1=z;
         z=z1-p1/pp;
      } while (fabs(z-z1) > EPS_gauleg);
      x[i-1]=xm-xl*z;
      x[n-i]=xm+xl*z;
      w[i-1]=2.0*xl/((1.0-z*z)*pp*pp);
      w[n-i]=w[i-1];
   }
}
#undef EPS_gauleg

/*==========================================================================*/
/*
++
double Integ_Fun(double xmin,double xmax,double (*fonc)(double),int npas)
        Pour integrer la fonction fonc de xmin a xmax sur npas.
        On emploit une methode Higher-order-gaussienne d'ordre 8, ce qui
        fait un calcul equivalent de N*npas pas.
--
*/
double Integ_Fun(double xmin,double xmax,double (*fonc)(double),int npas)
{
int i,j;
double dlim,sum,xc,xci;

if( xmax <= xmin ) return(0.);
if( npas <= 0 ) npas=1;
sum = 0.;
dlim = (xmax-xmin)/npas;
for(i=0;i<npas;i++) {
  xci = (double) i + 0.5;
  for(j=0;j<mIhoq8;j++) {
    xc = xmin + ( xci + Ihoq8X[j] ) * dlim;
    sum += fonc(xc) * Ihoq8W[j];
} }
return(sum*dlim);
}

/*=========================================================================*/
/*
++
double Integ_Fun_2D -
  (double (*fonc)(double x,double y) -
  ,double xmin,double xmax,double ymin,double ymax -
  ,int npasx,int npasy)

        Integration 2D de fonc(x,y) dans le carre [xmin,xmax] et [ymin,ymax]
--
*/
double Integ_Fun_2D(double (*fonc)(double x,double y)
                ,double xmin,double xmax,double ymin,double ymax
                ,int npasx,int npasy)
{
int i,j;
double x,y,pasx,pasy,sum,sumy;

pasx = (xmax-xmin)/npasx;
pasy = (ymax-ymin)/npasy;

sum = 0.;
for(i=0;i<npasx;i++) {
  x = xmin + ((double) i + 0.5 ) * pasx;
  sumy = 0.;
  for(j=0;j<npasy;j++) {
    y = ymin + ((double) j + 0.5 ) * pasy;
    sumy += fonc(x,y);
  }
  sum += sumy;
}
return( sum*pasx*pasy);
}

/*=========================================================================*/
/*                                    Christophe 01/07/94 */
/*
++
void Set_FitFunDPol(int DPol)
        Pour donner l'ordre du polynome pour le fit.
--
*/
void Set_FitFunDPol(int DPol) 
{
FITFUN_DPOL = DPol;
}

/*==================================================================*/
/*                                             Christophe 01/07/94  */
/*
++
double Gauss1DPolF(double x,double *Par,double *DgDpar) 
        Fonction de fit Gausienne+polynome.
--
*/
/*
++
|  f(x) = par[0]*exp[-0.5*( (x-par[1]) / par[2] )**2 ]
|        +par[3] + par[4]*x + .... + par[3+FITFUN_DPOL]*x**FITFUN_DPOL
|  FITFUN_DPOL peut etre definit avec la routine Set_FitFunDPol
--
*/
double Gauss1DPolF(double x,double *Par,double *DgDpar) 
{
double f,xc,xc2,e,xpow;
int i;

xc = (x-Par[1])/Par[2];
xc2 = xc*xc;
e = exp(-0.5*xc2);
f = Par[0]*e;

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

if(FITFUN_DPOL>=0) {
  xpow = 1.;
  for(i=0;i<=FITFUN_DPOL;i++) {
    DgDpar[3+i] = xpow;
    f += Par[3+i]*xpow;
    xpow *= x;
  }
}

return (f);
}

/*==================================================================*/
/*                                             Christophe 29/08/95  */
/*
++
double GaussI1DPol(double x,double *Par,double *DgDpar) 
        Fonction de fit Gausienne integree+polynome.
--
*/
/*
++
|  f(x) = par[0] / (sqrt(2*Pi)*par[2]) * exp[-0.5*( (x-par[1]) / par[2] )**2 ]
|        +par[3] + par[4]*x + .... + par[3+FITFUN_DPOL]*x**FITFUN_DPOL
|  FITFUN_DPOL peut etre definit avec la routine Set_FitFunDPol
--
*/
double GaussI1DPol(double x,double *Par,double *DgDpar) 
{
double f,xc,xc2,e,xpow;
int i;

xc = (x-Par[1])/Par[2];
xc2 = xc*xc;
e = exp(-0.5*xc2)/(S2Pi*Par[2]);
f = Par[0]*e;

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

if(FITFUN_DPOL>=0) {
  xpow = 1.;
  for(i=0;i<=FITFUN_DPOL;i++) {
    DgDpar[3+i] = xpow;
    f += Par[3+i]*xpow;
    xpow *= x;
  }
}

return (f);
}

/*====================================================================*/
/*                                             Christophe 01/07/94    */
/*
++
double Polyn1D(double x,double *Par,double *DgDpar) 
        Fonction de fit de polynome.
--
*/
/*
++
|  f(x) = par[0] + par[1]*x + .... + par[FITFUN_DPOL]*x**FITFUN_DPOL
|  FITFUN_DPOL peut etre definit avec la routine Set_FitFunDPol
--
*/
double Polyn1D(double x,double *Par,double *DgDpar) 
{
double f,xpow;
int i,DPol;

DPol = (FITFUN_DPOL<0) ? 0 : FITFUN_DPOL ;
xpow = 1.;
f = 0.;
for(i=0;i<=DPol;i++)
  {
  DgDpar[i] = xpow;
  f += Par[i]*xpow;
  xpow *= x;
  }
return (f);
}

/*=========================================================================*/
/*
++
double FitProp(double *x,double *y,double *ey,int *n,double *a1)
        Fit d'une proportion a1: y = a1 * x
--
*/
/*
++
| Input:
|    x = tableau des abscisses
|    y = tableau des ordonnees
|    ey = erreurs sur les y
|    n = nombre de donnees en entree
| Output:
|    n = nombre de points utilises (point non utilise si ey<=0)
|    a1 coefficient
| Return:
|    valeur du xi2 si fit reussi,
|    -1. si pas assez de points
|    -2. si determinant negatif
--
*/
double FitProp(double *x,double *y,double *ey,int *n,double *a1)
{
register int i,np;
register double X2,XY,Y2,w;

np=0;
X2=XY=Y2=*a1=0.;
for(i=0;i<*n;i++) {
  if(ey[i]<=0.) continue;
  np++;
  w = ey[i]*ey[i];
  X2 += x[i]*x[i]/w;
  XY += x[i]*y[i]/w;
  Y2 += y[i]*y[i]/w;
}
*n = np;
if(np<1) return(-1.);
if(X2==0.) return(-2.);
*a1 = XY/X2;
w = Y2 + *a1**a1*X2 -2.**a1*XY;
return(w);
}

/*=========================================================================*/
/*
++
double FitLin(double *x,double *y,double *ey,int *n,double *a0,double *a1)
        Fit d'une droite: y=a0+a1*x
--
*/
/*
++
| Input:
|    x = tableau des abscisses
|    y = tableau des ordonnees
|    ey = erreurs sur les y
|    n = nombre de donnees en entree
| Output:
|    n = nombre de points utilises (point non utilise si ey<=0)
|    a0,a1 coefficients de la droite
| Return:
|    valeur du xi2 si fit reussi,
|    -1. si pas assez de points
|    -2. si determinant negatif
--
*/
double FitLin(double *x,double *y,double *ey,int *n,double *a0,double *a1)
{
register int i,np;
register double I,X,X2,Y,XY,Y2,w;

np=0;
Y2=I=X=X2=Y=XY=*a0=*a1=0.;
for(i=0;i<*n;i++) {
  if(ey[i]<=0.) continue;
  np++;
  w = ey[i]*ey[i];
  I  += 1./w;
  X  += x[i]/w;
  X2 += x[i]*x[i]/w;
  Y  += y[i]/w;
  Y2 += y[i]*y[i]/w;
  XY += x[i]*y[i]/w;
}
*n = np;
if(np<2) return(-1.);
w = X*X-X2*I;
if(w==0.) return(-2.);
*a1 = (Y*X-XY*I)/w;
*a0 = (X*XY-X2*Y)/w;
w = Y2 + *a0**a0*I + *a1**a1*X2
    + 2.*( - Y**a0  - *a1*XY + *a0**a1*X );
return(w);
}

/*=========================================================================*/
/*
++
double FitPar(double *x,double *y,double *ey,int *n,double *a0,double *a1,double *a2)
        Fit d'une parabole: y=a0+a1.x+a2.x^2
--
*/
/*
++
| Input:
|    x = tableau des abscisses
|    y = tableau des ordonnees
|    ey = erreurs sur les y
|    n = nombre de donnees en entree
| Output:
|    n = nombre de points utilises (point non utilise si ey<=0)
|    a0,a1 coefficients de la droite
| Return:
|    valeur du xi2 si fit reussi,
|    -1. si pas assez de points
|    -2. si determinant negatif
--
*/
double FitPar(double *x,double *y,double *ey,int *n,double *a0,double *a1,double *a2)
{
register int i,np;
register double I,X,X2,X3,X4,Y,Y2,XY,X2Y,w,x2;

np=0;
I=X=X2=X3=X4=Y=Y2=XY=X2Y=*a0=*a1=*a2=0.;

for(i=0;i<*n;i++) {
  if(ey[i]<=0.) continue;
  np++;
  x2 = x[i]*x[i];
  w = ey[i]*ey[i];
  I   += 1./w;
  X   += x[i]/w;
  X2  += x2/w;
  X3  += x2*x[i]/w;
  X4  += x2*x2/w;
  Y   += y[i]/w;
  Y2  += y[i]*y[i]/w;
  XY  += x[i]*y[i]/w;
  X2Y += x2*y[i]/w;
}
*n = np;
if(np<3) return(-1.);
w = X2*(X2*X2-X3*X) -X*(X3*X2-X4*X) +I*(X3*X3-X4*X2);
if(w==0.) return(-2.);
*a2 =  (Y*(X2*X2-X3*X)  -XY*(X*X2-X3*I)  +X2Y*(X*X-X2*I)  )/w;
*a1 = -(Y*(X3*X2-X4*X)  -XY*(X2*X2-X4*I) +X2Y*(X2*X-X3*I) )/w;
*a0 =  (Y*(X3*X3-X4*X2) -XY*(X2*X3-X4*X) +X2Y*(X2*X2-X3*X))/w;
w = Y2 + *a0**a0*I + *a1**a1*X2 + *a2**a2*X4
    +2.*( -Y**a0  -*a1*XY -*a2*X2Y
          +*a0**a1*X +*a0**a2*X2+*a1**a2*X3 );
return(w);
}

/*=========================================================================*/
/*
++
double FitParLin -
  (double *xx,double *y,double *ey,int *n,double x0,int Deg_d -
  ,double *a0,double *a1,double *b1,double *b2)

--
*/
/*
++
| Fit d'une parabole: y=a0+b1*(x-x0)+b2*(x-x0)**2
|                                     pour x>x0
|    et d'une droite: y=a0+a1*(x-x0)  pour x<=x0
| Input:
|    x = tableau des abscisses
|    y = tableau des ordonnees
|    ey = erreurs sur les y
|    n = nombre de donnees en entree
|    x0 = abscisse du point pivot (il n'est pas fitte)
|    Deg_d = degre du fit a droite de x0 (1 ou 2)
| Output:
|    n = nombre de points utilises (point non utilise si ey<=0)
|    a0,a1 coefficients de la droite
|    a0,b1,b2 coefficients de la parabole
--
*/
/*
++
| Return:
|    valeur du xi2 si fit reussi,
|    -1. si pas assez de points
|    -2. si determinant negatif
| Remarque:
|    Il faut evidemment au moins 2 points a gauche du pivot x0
|                       et 3 points a droite du pivot x0
| Matrices A * a = B :
|        |   I_   X_g   X_d X2_d |   | a0 |   |    Y_ |
|        |  X_g  X2_g     0    0 | * | a1 | = |  XY_g |
|        |  X_d     0  X2_d X3_d |   | b1 |   |  XY_d |
|        | X2_d     0  X3_d X4_d |   | b2 |   | X2Y_d |
--
*/
double FitParLin(double *xx,double *y,double *ey,int *n,double x0,int Deg_d
                ,double *a0,double *a1,double *b1,double *b2)
{
int i,np,npg,npd,nddl;
double w,x,x2,A[4][4],B[4];
double Y2,I,Y,X_g,X2_g,XY_g,X_d,X2_d,X3_d,X4_d,XY_d,X2Y_d;

if( Deg_d<1 || Deg_d>2 ) Deg_d=2;
nddl = Deg_d + 2;
*a0=*a1=*b1=*b2=0.;
Y2=I=Y=X_g=X2_g=XY_g=X_d=X2_d=X3_d=X4_d=XY_d=X2Y_d=0.;
np=npg=npd=0;
for(i=0;i<*n;i++) {
  if(ey[i]<=0.) continue;
  np++;
  x = xx[i]-x0;
  x2 = x*x;
  w = ey[i]*ey[i];
  I  += 1./w;
  Y  += y[i]/w;
  Y2 += y[i]*y[i]/w;
  if(x<=0.) {
    if(x<0.) npg++;
    X_g  += x/w;
    X2_g += x2/w;
    XY_g += x*y[i]/w;;
  } else {
    npd++;
    X_d   += x/w;
    X2_d  += x2/w;
    X3_d  += x2*x/w;
    X4_d  += x2*x2/w;
    XY_d  += x*y[i]/w;;
    X2Y_d += x2*y[i]/w;;
  }
}
*n = np;
A[0][0] = I;
A[1][1] = X2_g;
A[3][3] = X4_d;
A[0][1] = A[1][0] = X_g;
A[0][2] = A[2][0] = X_d;
A[3][2] = A[2][3] = X3_d;
A[0][3] = A[3][0] = A[2][2] = X2_d;
A[1][2] = A[2][1] = A[1][3] = A[3][1] = 0.;
B[0] = Y;
B[1] = XY_g;
B[2] = XY_d;
B[3] = X2Y_d;
if( np<4 || npg < 1 || npd < 2 ) return(-1.);
w = GausPiv(&A[0][0],4,nddl,B,1,1,0);
if(w==0.) return(-2.);
*a0 = B[0];
*a1 = B[1];
*b1 = B[2];
if(nddl==4) *b2 = B[3];
w = Y2 + *a0**a0*I + *a1**a1*X2_g + *b1**b1*X2_d + *b2**b2*X4_d
  + 2.*( - *a0*Y - *a1*XY_g + *a0**a1*X_g
         - *b1*XY_d - *b2*X2Y_d + *a0**b1*X_d + *a0**b2*X2_d + *b1**b2*X3_d
       );
return(w);
}

/*=========================================================================*/
/*
++
double FitPropClean(double *x,double *y,double *ey,int *n,double per_clean,double *a1)
        Fit de y(i) = a1*x(i) en deux passes:
--
*/
/*
++
|   1-/ fit avec tous les points
|   2-/ fit ou on enleve les per_clean*n points les plus eloignes
| Input:
|    x = tableau des abscisses
|    y = tableau des ordonnees
|    ey = erreurs sur les y
|    n = nombre de donnees en entree
|    per_clean = pourcentage des points a tuer pour la seconde passe
--
*/
/*
++
| Output:
|    n = nombre de points utilises (point non utilise si ey<=0)
|    a1 = coefficients de proportionalite
| Return:
|    valeur du xi2 si fit reussi,
|    -1. si echec fit passe no 1
|    -2. si echec fit passe no 2
|    -3. si probleme malloc
|    -4. si probleme nombre de points a tuer
--
*/
double FitPropClean(double *x,double *y,double *ey,int *n,double per_clean,double *a1)
{
int npt,k,i,nclass;
double *clas,aa1,c2,cut;

*a1 =0.;

/* 1ere passe */
npt = *n;
c2 = FitProp(x,y,ey,&npt,&aa1);
if(c2<0.) {*n = npt; return(-1.);}
*a1 = aa1;
/* printf("pass 1: %g*x c2=%g/%d\n",*a1,c2,npt); */

clas = (double*) malloc( *n * sizeof(double) );
if( clas == NULL ) {*n=npt; return(-3.);}

/* elimination des mauvais points */
nclass = 0;
for(i=0;i<*n;i++) {
  if(ey[i]<=0.) continue;
  c2 = (y[i]-aa1*x[i])/ey[i];
  clas[nclass] = c2*c2;
  nclass++;
}
qsort(clas,(size_t) nclass,(size_t) sizeof(double),qSort_Dble);
k = (int) ( (1. - per_clean ) * (double) nclass );
if(k<0) k=0;
if(k>=nclass) k = nclass-1;
cut = clas[k];
k = 0;
for(i=0;i<*n;i++) {
  clas[i] = ey[i];
  c2 = (y[i]-aa1*x[i])/ey[i];
  c2 *= c2;
  if(ey[i]>0. && c2>cut) {clas[i] = -1.; k++;}
}
/* printf("nombre pt tues %d cut=%g\n",k,cut); */

/* 2sd passe */
npt = *n;
c2 = FitProp(x,y,clas,&npt,&aa1);
if(c2<0.) {*n = npt; free(clas); return(-2.);}
*a1 = aa1;
*n = npt;
/* printf("pass 2: %g*x c2=%g/%d\n",*a1,c2,npt); */

free(clas);
return(c2);
}

/*=========================================================================*/
/*
++
double FitLinClean -
  (double *x,double *y,double *ey,int *n -
  ,double per_clean,double *a0,double *a1)

        Fit de y(i) = a0 + a1*x(i) en deux passes:
--
*/
/*
++
|   1-/ fit avec tous les points
|   2-/ fit ou on enleve les per_clean*n points les plus eloignes
| Input:
|    x = tableau des abscisses
|    y = tableau des ordonnees
|    ey = erreurs sur les y
|    n = nombre de donnees en entree
|    per_clean = pourcentage des points a tuer pour la seconde passe
| Output:
|    n = nombre de points utilises (point non utilise si ey<=0)
|    a0,a1 = coefficients de la droite
--
*/
/*
++
| Return:
|    valeur du xi2 si fit reussi,
|    -1. si echec fit passe no 1
|    -2. si echec fit passe no 2
|    -3. si probleme malloc
|    -4. si probleme nombre de points a tuer
--
*/
double FitLinClean(double *x,double *y,double *ey,int *n
                  ,double per_clean,double *a0,double *a1)
{
int npt,k,i,nclass;
double *clas,aa0,aa1,c2,cut;

*a0 = *a1 =0.;

/* 1ere passe */
npt = *n;
c2 = FitLin(x,y,ey,&npt,&aa0,&aa1);
if(c2<0.) {*n = npt; return(-1.);}
*a0 = aa0;
*a1 = aa1;
/* printf("pass 1: %g + %g*x c2=%g/%d\n",*a0,*a1,c2,npt); */

clas = (double*) malloc( *n * sizeof(double) );
if( clas == NULL ) {*n=npt; return(-3.);}

/* elimination des mauvais points */
nclass = 0;
for(i=0;i<*n;i++) {
  if(ey[i]<=0.) continue;
  c2 = (y[i]-(aa0+aa1*x[i]))/ey[i];
  clas[nclass] = c2*c2;
  nclass++;
}
qsort(clas,(size_t) nclass,(size_t) sizeof(double),qSort_Dble);
k = (int) ( (1. - per_clean ) * (double) nclass );
if(k<0) k=0;
if(k>=nclass) k = nclass-1;
cut = clas[k];
k = 0;
for(i=0;i<*n;i++) {
  clas[i] = ey[i];
  c2 = (y[i]-(aa0+aa1*x[i]))/ey[i];
  c2 *= c2;
  if(ey[i]>0. && c2>cut) {clas[i] = -1.; k++;}
}
/* printf("nombre pt tues %d cut=%g\n",k,cut); */

/* 2sd passe */
npt = *n;
c2 = FitLin(x,y,clas,&npt,&aa0,&aa1);
if(c2<0.) {*n = npt; free(clas); return(-2.);}
*a0 = aa0;
*a1 = aa1;
*n = npt;
/* printf("pass 2: %g + %g*x c2=%g/%d\n",*a0,*a1,c2,npt); */

free(clas);
return(c2);
}

/*===================================================================*/
/*
++
double FitParClean -
  (double *x,double *y,double *ey,int *n -
  ,double per_clean,double *a0,double *a1,double *a2)

        Fit de y(i) = a0 + a1*x(i) + a2*x(i)^2 en deux passes:
--
*/
/*
++
|   1-/ fit avec tous les points
|   2-/ fit ou on enleve les per_clean*n points les plus eloignes
| Input:
|    x = tableau des abscisses
|    y = tableau des ordonnees
|    ey = erreurs sur les y
|    n = nombre de donnees en entree
|    per_clean = pourcentage des points a tuer pour la seconde passe
| Output:
|    n = nombre de points utilises (point non utilise si ey<=0)
|    a0,a1,a2 = coefficients de la parabole
| Return:
|    valeur du xi2 si fit reussi,
|    -1. si echec fit passe no 1
|    -2. si echec fit passe no 2
|    -3. si probleme malloc
|    -4. si probleme nombre de points a tuer
--
*/
double FitParClean(double *x,double *y,double *ey,int *n
                  ,double per_clean,double *a0,double *a1,double *a2)
{
int npt,k,i,nclass;
double *clas,aa0,aa1,aa2,c2,cut;

*a0 = *a1 = *a2 =0.;

/* 1ere passe */
npt = *n;
c2 = FitPar(x,y,ey,&npt,&aa0,&aa1,&aa2);
if(c2<0.) {*n = npt; return(-1.);}
*a0 = aa0;
*a1 = aa1;
*a2 = aa2;
/* printf("pass 1: %g + %g*x + %g*x**2 c2=%g/%d\n",*a0,*a1,*a2,c2,npt); */

clas = (double*) malloc( *n * sizeof(double) );
if( clas == NULL ) {*n=npt; return(-3.);}

/* elimination des mauvais points */
nclass = 0;
for(i=0;i<*n;i++) {
  if(ey[i]<=0.) continue;
  c2 = (y[i]-(aa0+aa1*x[i]+aa2*x[i]*x[i]))/ey[i];
  clas[nclass] = c2*c2;
  nclass++;
}
qsort(clas,(size_t) nclass,(size_t) sizeof(double),qSort_Dble);
k = (int) ( (1. - per_clean ) * (double) nclass );
if(k<0) k=0;
if(k>=nclass) k = nclass-1;
cut = clas[k];
k = 0;
for(i=0;i<*n;i++) {
  clas[i] = ey[i];
  c2 = (y[i]-(aa0+aa1*x[i]+aa2*x[i]*x[i]))/ey[i];
  c2 *= c2;
  if(ey[i]>0. && c2>cut) {clas[i] = -1.; k++;}
}
/* printf("nombre pt tues %d cut=%g\n",k,cut); */

/* 2sd passe */
npt = *n;
c2 = FitPar(x,y,clas,&npt,&aa0,&aa1,&aa2);
if(c2<0.) {*n = npt; free(clas); return(-2.);}
*a0 = aa0;
*a1 = aa1;
*a2 = aa2;
*n = npt;
/* printf("pass 2: %g + %g*x + %g*x**2 c2=%g/%d\n",*a0,*a1,*a2,c2,npt); */

free(clas);
return(c2);
}

/*===========================================================*/
/*
++
float interpole(float m,int n,float *t1,float *t2)
        Interpolation lineaire y=f(x):
--
*/
/*
++
|   t1[n] = tableau des x
|   t2[n] = tableau des y
|   m = valeur ou on interpole
|   RETURN: f(m)
--
*/
float interpole(float m,int n,float *t1,float *t2)
{
int i,n1;
float r;

if ( n == 1 ) return ( t2[0] );

n1=n-1;
if ( m < t1[0] )  m=t1[0];
if ( m > t1[n1] ) m=t1[n1];

for(i=0;i<n;i++) if ( m <= t1[i] ) break;
i = (i==0) ? 1  : i;
i = (i==n) ? n1 : i;

r = t2[i-1] + (t2[i]-t2[i-1]) / (t1[i]-t1[i-1]) * (m-t1[i-1]);

/*
printf("interpole: i=%3d m=%g t1=%g %g t2=%g %g r=%g\n"
      ,i,m,t1[i-1],t1[i],t2[i-1],t2[i],r);
for(i=0;i<n;i++) printf(" (%g,%g)",t1[i],t2[i]); printf("\n");
*/

return(r);
}

/*=========================================================*/
/*                         Christophe 10/11/93 La Silla    */
/*
++
int_4 MeanSig(float *mag,float *err,int_4 *ndata,float nsigma,float *mean,float *sigma)
        Fonction de calcul de la moyenne tronquee et du sigma d'un tableau
--
*/
/*
++
| Methode: 1-/ moyenne-sigma, 2-/ moyenne tronquee et sigma,
| 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
| SORTIE:
|  ndata: nombre de data utilises pour calculer mean et sigma
|  mean: moyenne des donnees
|  sigma: sigma des donnees
|  MeanSig: 0 si calcul ok, <0 sinon
--
*/
int_4 MeanSig(float *mag,float *err,int_4 *ndata,float nsigma,float *mean,float *sigma)
{
int_4 deb=DEB_MeanSig;
int_4 pass,passmx,i,n=0;
double m,s2,v,scut;

*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 ) {
      n++;
      m += v;
      s2 += v * v;
  } }
  if ( n >= 2 ) {
    *mean = m / n;
    *sigma = sqrt( s2 / n - m/n * m/n );
  } else {
    *mean = *sigma = 0.;
    *ndata=n;
    return(-1);
  }
  if ( deb>0 ) printf("MeanSig: pass=%d mean=%f sigma=%f n=%d\n"
                   ,pass,*mean,*sigma,n);
}
*ndata=n;
return(0);
}

/*=========================================================================*/
/*
++
int Most_Probable -
  ( float *val, int nval, float binin, int nmoy -
  , float low, float high, float *most, int deb)

        Pour calculer la valeur la plus probable par maximum d'histogramme.
--
*/
/*
++
| INPUT:
|  - val= pointeur sur un tableau de valeurs (float).
|  - nval= nombre de valeurs a traiter.
|  - binin= taille du bin.
|  - nmoy= la valeur retournee est la moyenne de +/- nmoy bins
|          autour du max de l'histo.
|  - low, high= les valeurs permises mini et maxi des valeurs
|                (si low>high calcul automatique).
|  - deb= niveau de debug.
| OUTPUT:
|  - most = valeur la plus probable.
| RETURN CODE:
|  -1 : nval <=0
|  -2 : bin <=0
|  -3 : low>=high apres cadrage automatique
|  -4 : nombre de valeurs remplies nul au niveau de la 1ere passe
|  -51 : pas trouve de maximum 1ere passe
|  -52 : pas trouve de maximum 2sd passe
--
*/
int Most_Probable( float *val, int nval, float binin, int nmoy
                 , float low, float high, float *most, int deb)
{
int histo[503],i,nbin,n,ib,ibmax,vmax,pass;
float bin,xmax=0.;
double moy,sum;

if(deb>1) printf("Most_Probable: nval=%d binin=%f nmoy=%d low,high=%f %f\n"
                ,nval,binin,nmoy,low,high);
if (nval<=0) return(-1);
if (binin<=0.) return(-2);
nmoy = (nmoy<0) ? 0 : nmoy;
if(low>= high) {
  low = GRAND;
  high = -low;
  for(i=0;i<nval;i++) {if(val[i]>high) high=val[i]; if(val[i]<low) low=val[i];}
}
if(deb>1) printf("nmoy=%d low,high=%f %f\n",nmoy,low,high);
if(low>= high) return(-3);

/* premiere passe "nbin" bins entre low et high, nbin=min(high-low)/bin+1,501) */
/* seconde  passe "nbin" autour du maxi 1ere passe, nbin=min(high-low)/bin,501) */

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

  bin=binin;
  nbin = (int) ( (high-low)/bin ) + 1;
  if(nbin>501) {
    nbin = 501;
    if(pass!=1){
      low  = xmax - 250*binin - binin/2.;
      high = xmax + 250*binin + binin/2.;
  } }
  /* re-scaling du bin pour etre sur */ 
  bin=(high-low)/nbin;
  if(deb>1) printf("pass=%2d low,high=%f %f bin=%f nbin=%d\n"
                  ,pass,low,high,bin,nbin);
  for(i=0;i<nbin;i++) histo[i]=0;
  n=0;
  for(i=0;i<nval;i++) {
    if( val[i]>=high || val[i]<low ) continue;
    ib = (int) ((val[i]-low)/bin);
    if( ib<0 && ib >= nbin ) continue;
    histo[ib]++;
    n++;
  }
  if(n<=0) return (-4);
  vmax = -1; ibmax= -1;
  if(n<=0) return (-4);
  for(i=0;i<nbin;i++) {
    if(histo[i]<vmax) continue;
    vmax = histo[i];
    ibmax = i;
  }
  if(ibmax<0) return (-50-pass);
  xmax = low + ibmax*bin + bin/2.;
  if(deb>1) printf("bin du maximum %5d (%d) abscisse=%f\n",ibmax,vmax,xmax);
  if(deb>3) {
    printf("histogramme:\n");
    for (i= ibmax-25;i<=ibmax+25;i++) {
      if (i<0) continue;
      if (i>=nbin) break;
      printf(" %5d",histo[i]);
      if(i%10==9) printf("\n");
    }
    printf("\n");
  }
}

sum = moy =0.; n=0;
for (i= -nmoy;i<=nmoy;i++) {
  if (ibmax-i>=0 && ibmax+i<nbin) {
    moy += ( xmax + i * bin ) * histo[ibmax+i];
    sum += (double) histo[ibmax+i];
    n++;
} }

*most = moy / sum;
if(deb>0) printf("Most_probable: most=%f (avec %d points)\n",*most,n);

return(0);

}

/*============================================================*/
/*                                      Christophe 27/01/95   */
/*
++
float Ajust_GaFn
|    (int Ns,float *fcfr
|    ,float *haut,float *ehaut,float *mean,float *emean
|    ,float *sigma,float *esigma,float *fond,float *efond
|    ,int fixfond,int NBIN_RESOL
|    ,float perc_net,float frac_sgb,float dyn_sgb,int deb)

        Fonction de fit gaussienne + fond (1D) du tableau fcfr
        avec reglage des parametres d'entree.
--
*/
/*
++
| ATTENTION: le tableau fcfr est modifie !!!              
| ---------                                               
| ENTREES:                                                
|   Ns = nombre d'entrees dans fcfr
|   fcfr = tableau des valeurs a fitter
|   fond = valeur a donner au fond si fixfond = 1
|   fixfond =  1 fond fixe dans le fit a la valeur fond
|           =  2 fond fixe dans le fit a une valeur
|                     calculee automatiquement sur
|                     les bords de l'histogramme 
|           >  2 fond fixe a zero (<=> 1 + fond=0)
|           <= 0 fond libre dans le fit
|   NBIN_RESOL = nombre maxi de bins pour l'histo de fit
|                de la resolution
|   perc_net = pourcentage de donnees a nettoyer pour
|              calculer mean,sigma de la 1ere passe
|   frac_sgb = le bin de fit sera sigma*fracsgb
|              ou sigma est calcule en premiere passe
|   dyn_sgb = limites de l'histo de fit
|        xmin=mean-dyn_sgb*sigma , xmax=mean+dyn_sgb*sigma
|   deb = niveau de debug
--
*/
/*
++
| SORTIES:                                                
|   haut = valeur de la hauteur
|   ehaut = erreur sur la valeur de la hauteur
|   mean = valeur moyenne du decalage
|   emean = erreur sur la valeur moyenne du decalage
|   sigma = valeur de la resolution
|   esigma = erreur sur la valeur de la resolution
|   fond = valeur du fond
|   efond = erreur sur la valeur du fond 
|   return code = chi2 du fit  (<0 si Pb)
--
*/
float Ajust_GaFn(int Ns,float *fcfr
           ,float *haut,float *ehaut,float *mean,float *emean
           ,float *sigma,float *esigma,float *fond,float *efond
           ,int fixfond,int NBIN_RESOL
           ,float perc_net,float frac_sgb,float dyn_sgb,int deb)
{
int_4 i,j,i1,i2,nbin,ibin,entries;
float c2,xmin,xmax,binw,ymax;
float *y,*ey,*x;
double m,s;
FLOATDOUBLE X,Y,EY;
double (*Fun) (double ,double *,double *);
int_4 npar,n;
double par[4], epar[4],minpar[4],maxpar[4],stepar[4],stochi2;

fixfond = ( fixfond <= 0 ) ? 0 : fixfond;
perc_net = ( perc_net >= 1. || perc_net < 0. ) ? 0.1 : perc_net;
NBIN_RESOL = ( NBIN_RESOL <= 0 ) ? 50 : NBIN_RESOL;
frac_sgb = ( frac_sgb <= 0. ) ? 0.5 : frac_sgb;
dyn_sgb  = ( dyn_sgb  <= 0. ) ? 5.0 : dyn_sgb;
x = y = ey = NULL;
x  = (float *) calloc(NBIN_RESOL+2,sizeof(float));
y  = (float *) calloc(NBIN_RESOL+2,sizeof(float));
ey = (float *) calloc(NBIN_RESOL+2,sizeof(float));
if ( x==NULL || y==NULL || ey==NULL ) {c2 = -100.; goto END;}

if(deb>0) printf("Ajust_GaFn[%d] perc_net=%f nbin=%d frac_sgb=%f dyn_sgb=%f\n"
              ,Ns,perc_net,NBIN_RESOL,frac_sgb,dyn_sgb);

/* on vire perc_net % des points a gauche et a droite de la distribution */
qsort(fcfr,(size_t) Ns,(size_t) sizeof(float),qSort_Float);
i = perc_net * Ns;
if(perc_net>0. && i==0) i=1;
i1 = i; i2 = Ns-1-i;
if(deb>1) {
  printf("  df/f=");
  for(j=0;j<10;j++) printf("%8.2f ",fcfr[j]);
  printf("\n");
  printf("selection de %d (%f) a %d (%f)\n",i1,fcfr[i1],i2,fcfr[i2]);
}
if(i2<=i1) { c2= -101.; goto END;}

/* calcul mean et sigma 1ere passe */
m = s = 0.; j = 0;
for(i=i1;i<=i2;i++) {m += fcfr[i]; s += fcfr[i]*fcfr[i]; j++;}
m /= j;
s = sqrt(s/j - m*m);
*mean = m; *sigma = s;
if(deb>1) printf("premiere passe: mean=%f sigma=%f\n",*mean,*sigma);

/* remplissage histo a +/- dyn_sgb * sigma */
xmin = *mean - dyn_sgb* *sigma;
xmax = *mean + dyn_sgb* *sigma;
if(xmin>=xmax) {c2 = -102.; goto END;}
binw = *sigma * frac_sgb;
nbin = (xmax-xmin)/binw;
if(nbin<6) nbin=6; if(nbin>NBIN_RESOL) nbin=NBIN_RESOL;
binw = (xmax-xmin)/nbin;
for(i=0;i<nbin;i++)
    {x[i] = xmin + ((float) i + 0.5)*binw; y[i] = 0.;}
if(deb>1) {
  printf("histo[%d] lim=[%f,%f] bw=%f x=\n"
        ,nbin,xmin,xmax,binw);
  for(j=0;j<nbin;j++)
    {printf("%8.2f ",x[j]); if(j%10==9) printf("\n");}
  if(j%10!=0) printf("\n");
}
entries = 0;
for(i=0;i<Ns;i++) { ibin = (fcfr[i]-xmin)/binw;
         if(ibin>=0 && ibin<nbin) {y[ibin]+=1.; entries++;} }
ymax = 0.;
for(i=0;i<nbin;i++) {ey[i] = (y[i]>1.)? sqrt((double) y[i]) : 1.;
                    if(y[i]>ymax) ymax=y[i];}
if(deb>1) {
  printf("histo (entries=%d) ymax=%f y=\n",entries,ymax);
  for(j=0;j<nbin;j++)
    {printf("%8.0f ",y[j]); if(j%10==9) printf("\n");}
  if(j%10!=0) printf("\n");
  if(deb>2) {
    printf("histo ey=\n");
    for(j=0;j<nbin;j++)
      {printf("%8.0f ",ey[j]); if(j%10==9) printf("\n");}
    if(j%10!=0) printf("\n"); }
}
if(ymax<2.) {c2 = -103.; goto END;}

if(fixfond!=1) {
  *fond = (y[0]+y[1]+y[nbin-1]+y[nbin-2])/4.;
  if(fixfond>2) *fond = 0.;
}
*haut = ymax - *fond;

/* fit gaussienne + fond constant */
Fun = Gauss1DPolF ; Set_FitFunDPol(0);
X.fx  = x;
Y.fx  = y;
EY.fx = ey;
par[0] = *haut;
par[1] = *mean;
par[2] = *sigma;
par[3] = *fond;
stepar[0] = par[0]/10.; stepar[1] = par[2]/5.;
stepar[2] = par[2]/10.; stepar[3] = 0.5;
if(fixfond) { stepar[3] = 0.; epar[3]=0.;}
for(i=0;i<4;i++) {minpar[i] = 1.; maxpar[i] = -1.;}
minpar[2] = 0.00001; maxpar[2] = 9999.;
npar = 4;
stochi2 = 0.0001;
n = nbin;
j = (deb>2)? 1:0;
c2 = FitFun(Fun,FITFUN_FLOAT,X,Y,EY,&n,0,0,par
           ,epar,stepar,minpar,maxpar,npar,stochi2,60,j);

END:

if( c2 > 0. ) {
  *haut  =  par[0]; *mean  =  par[1]; *sigma  =  par[2]; *fond  =  par[3];
  *ehaut = epar[0]; *emean = epar[1]; *esigma = epar[2]; *efond = epar[3];
} else {
  *ehaut = *emean = *esigma = *efond = 0.;
}
if(x !=NULL) free(x);
if(y !=NULL) free(y);
if(ey!=NULL) free(ey);

return(c2);
}

/*==========================================================*/
/*                                    Christophe 27/03/95   */
/*
++
float HFit_GaFn -
  (int Ns,float *fcfr -
  ,float *haut,float *ehaut,float *mean,float *emean -
  ,float *sigma,float *esigma,float *fond,float *efond -
  ,int fixfond,int nbin,float xmin,float xmax,int deb)

        Fonction de fit gaussienne + fond (1D) du tableau fcfr
        dans un histo de "nbin" bins de "xmin a xmax"
--
*/
/*
++
| ENTREES:
|   Ns = nombre d'entrees dans fcfr
|   fcfr = tableau des valeurs a fitter
|   mean = valeur de depart pour la moyenne
|   sigma = valeur de depart pour le sigma
|   fond = valeur de depart pour le fond
|   fixfond > 0  fond fixe dans le fit a la valeur fond
|           <= 0 fond libre dans le fit
|   nbin = nombre de bins pour l'histo de fit
|   xmin = valeur inferieure de l'histo
|   xmax = valeur superieure de l'histo
|   deb = niveau de debug
--
*/
/*
++
| SORTIES:
|   haut = valeur de la hauteur
|   ehaut = erreur sur la valeur de la hauteur
|   mean = valeur moyenne du decalage
|   emean = erreur sur la valeur moyenne du decalage
|   sigma = valeur de la resolution
|   esigma = erreur sur la valeur de la resolution
|   fond = valeur du fond
|   efond = erreur sur la valeur du fond
|   return code = chi2 du fit  (<0 si Pb)
--
*/
float HFit_GaFn(int Ns,float *fcfr
           ,float *haut,float *ehaut,float *mean,float *emean
           ,float *sigma,float *esigma,float *fond,float *efond
           ,int fixfond,int nbin,float xmin,float xmax,int deb)
{
int_4 i,j,ibin,entries;
float c2,binw,ymax;
float *y,*ey,*x;
FLOATDOUBLE X,Y,EY;
double (*Fun) (double ,double *,double *);
int_4 npar,n;
double par[4], epar[4],minpar[4],maxpar[4],stepar[4],stochi2;

if(xmax<=xmin) return(-102);
fixfond = ( fixfond <= 0 ) ? 0 : fixfond;
nbin = ( nbin <= 0 ) ? 50 : nbin;
x = y = ey = NULL;
x  = (float *) calloc(nbin+2,sizeof(float));
y  = (float *) calloc(nbin+2,sizeof(float));
ey = (float *) calloc(nbin+2,sizeof(float));
if ( x==NULL || y==NULL || ey==NULL ) {c2 = -100.; goto END;}

if(deb>0) printf("HFit_GaFn[%d pts] nbin=%d de %g a %g\n"
              ,Ns,nbin,xmin,xmax);

/* remplissage histo */
binw = (xmax-xmin)/nbin;
for(i=0;i<nbin;i++) {x[i]=xmin+((float) i +0.5)*binw; y[i]=0.;}
if(deb>1) {
  printf("histo[%d] lim=[%f,%f] bw=%f x=\n",nbin,xmin,xmax,binw);
  for(j=0;j<nbin;j++){printf("%8.2f ",x[j]); if(j%10==9) printf("\n");}
  if(j%10!=0) printf("\n");
}
entries = 0;
for(i=0;i<Ns;i++) {
  ibin = (fcfr[i]-xmin)/binw;
  if(ibin>=0 && ibin<nbin) {y[ibin]+=1.; entries++;}
}

ymax = 0.;
for(i=0;i<nbin;i++) {ey[i] = (y[i]>1.)? sqrt((double) y[i]) : 1.;
                    if(y[i]>ymax) ymax=y[i];}
if(deb>1) {
  printf("histo (entries=%d) ymax=%f y=\n",entries,ymax);
  for(j=0;j<nbin;j++){printf("%8.0f ",y[j]); if(j%10==9)printf("\n");}
  if(j%10!=0) printf("\n");
  if(deb>2) {
    printf("histo ey=\n");
    for(j=0;j<nbin;j++){printf("%8.0f ",ey[j]); if(j%10==9)printf("\n");}
    if(j%10!=0) printf("\n");
} }
if(ymax<2.) {c2 = -103.; goto END;}

*haut = ymax - *fond;

/* fit gaussienne + fond constant */
Fun = Gauss1DPolF ; Set_FitFunDPol(0);
X.fx  = x;
Y.fx  = y;
EY.fx = ey;
par[0] = *haut;
par[1] = *mean;
par[2] = *sigma;
par[3] = *fond;
stepar[0] = par[0]/10.; stepar[1] = par[2]/5.;
stepar[2] = par[2]/10.; stepar[3] = 0.5;
if(fixfond) { stepar[3] = 0.; epar[3]=0.;}
for(i=0;i<4;i++) {minpar[i] = 1.; maxpar[i] = -1.;}
minpar[2] = 0.00001; maxpar[2] = 9999.;
npar = 4;
stochi2 = 0.0001;
n = nbin;
j = (deb>2)? 1:0;
c2 = FitFun(Fun,FITFUN_FLOAT,X,Y,EY,&n,0,0,par
           ,epar,stepar,minpar,maxpar,npar,stochi2,60,j);

END:

if( c2 > 0. ) {
  *haut  =  par[0]; *mean  =  par[1]; *sigma  =  par[2]; *fond  =  par[3];
  *ehaut = epar[0]; *emean = epar[1]; *esigma = epar[2]; *efond = epar[3];
} else {
  *ehaut = *emean = *esigma = *efond = 0.;
}
if(x !=NULL) free(x);
if(y !=NULL) free(y);
if(ey!=NULL) free(ey);

return(c2);
}

/*==============================================================*/
/*                                    Christophe 15/02/95       */
/*
++
int Ajust_MnSg -
  (int Ns,float *fcfr -
  ,float *mean,float *emean,float *sigma,float *esigma -
  ,float perc_net,float n_sigma,int deb)

        Calcul mean sigma tronques du tableau fcfr.
--
*/
/*
++
| ATTENTION: le tableau fcfr est modifie
| ---------
| Methode:
|  1- calcul de la moyenne Mcent en enlevant perc_net% des pts
|     les + faibles et perc_net% des pts les + forts
|  2- calcul de la valeur cutup correspondant a 2*perc_net%
|     des pts les + hauts en valeur absolue (fcfr-Mcent)
|  3- calcul moyenne et sigma en enlevant les points ayant
|        abs(fcfr-Mcent) > cutup    (1ere passe)
|  4- calcul moyenne et sigma en enlevant les points ayant
|        abs(fcfr-Mcent) > cutup   et
|        abs(fcfr-M(1ere passe)) > n_sigma*sigma(1ere passe)
--
*/
/*
++
| ENTREES:
|   Ns = nombre d'entrees dans fcfr
|   fcfr = tableau des valeurs a fitter
|   perc_net = pourcentage de donnees a nettoyer pour
|              calculer mean,sigma
|   n_sigma = nombre de sigma en passe 2 pour
|              calculer mean,sigma
|   deb = niveau de debug
| SORTIES:
|   mean = valeur moyenne du decalage
|   emean = erreur sur la valeur moyenne du decalage
|   sigma = valeur de la resolution
|   esigma = erreur sur la valeur de la resolution
|   return code = nombre de points utilises (<0 si Pb)
--
*/
int Ajust_MnSg(int Ns,float *fcfr
           ,float *mean,float *emean,float *sigma,float *esigma
           ,float perc_net,float n_sigma,int deb)
{
int i,i1,i2,n,nvire,pass;
double m,s,mcent,cutup,cuts;
float *dum;

/* set up des valeurs par defaut */
perc_net = ( perc_net >= 1. || perc_net < 0. ) ? 0.1 : perc_net;
n_sigma = ( n_sigma <= 0. ) ? 3.0 : n_sigma;
nvire = perc_net * Ns;
if(perc_net>0. && nvire==0) nvire=1;
*emean = *esigma = -1.;
if(deb>0) printf("Ajust_MnSg[%d] perc_net=%f (vire=%d) n_sigma=%f\n"
                ,Ns,perc_net,nvire,n_sigma);

/* on ordonne par ordre croissant de fcfr */
qsort(fcfr,(size_t) Ns,(size_t) sizeof(float),qSort_Float);
i1 = nvire; i2 = Ns-1-nvire;
if(deb>1) printf("    select de %d (%f) a %d (%f)\n",i1,fcfr[i1],i2,fcfr[i2]);
if(i2<=i1) return(-101);

/* on vire perc_net % des points a gauche et a droite de la distribution */
mcent = 0.; n = 0;
for(i=i1;i<=i2;i++) {mcent += fcfr[i]; n++;}
if(n<1) return(-102);
mcent /= (double) n;
if(deb>1) printf("    mean cent=%f sur %d points\n",mcent,n);

/* allocate space and fill absolute value */
if( (dum  = (float *) calloc(Ns,sizeof(float))) == NULL ) return(-103);
for(i=0;i<Ns;i++) dum[i] = fabs( mcent - fcfr[i] );

/* on ordonne par ordre croissant de abs(fcfr-mcent) */
qsort(dum,(size_t) Ns,(size_t) sizeof(float),qSort_Float);
i = Ns-1-2*nvire;
if(i<0) {n= -104; goto END_Ajust_MnSg;}

/* on vire 2.*perc_net % des points les + hauts de la distribution */
cutup = dum[i];
if(deb>1) printf("    cutup=%f sur %d points\n",cutup,i);

/* calcul mean et sigma en 2 passes avec coupure n_sigma*sigma(1ere passe)*/
cuts = GRAND2;
for(pass=1;pass<=2;pass++) {
  m = s = 0.; n = 0;
  for(i=0;i<Ns;i++) {
    if( fabs(mcent-fcfr[i]) > cutup ) continue;
    if( fabs((double) (*mean-fcfr[i])) > cuts ) continue;
    m += fcfr[i];
    s += fcfr[i]*fcfr[i];
    n++;
  }
  if(n<=2) {n= -105; *mean = *sigma = -1.; goto END_Ajust_MnSg;}
  m /= n;
  s = s/n - m*m;
  if(s<=0.) {n= -106; *mean = *sigma = -1.; goto END_Ajust_MnSg;}
  s = sqrt(s);
  *mean = m;
  *sigma = s;
  cuts = n_sigma * s;
  if(deb>1)
    printf("    pass=%d mean=%g sigma=%g (%d pts) cuts=%f\n"
          ,pass,*mean,*sigma,n,cuts);
}

*emean = s/sqrt((double) n);
*esigma = s/sqrt(2. * (double) n);
if(deb>0) printf("mean=%g (+/-%g) sigma=%g (+/-%g) sur %d pts\n"
              ,*mean,*emean,*sigma,*esigma,n);

END_Ajust_MnSg:
if( dum != NULL ) free(dum);
return(n);
}

/*=================================================================*/
/*                                       Nathalie 16/02/95         */
/*
++
float Int3DCube(float(*fonction)(float,float,float),float step,float precision)
        Cette fonction calcule l'integrale 3D: f(x,y,z)dx.dy.dz
        L'integrale est calculee en incrementant des couronnes cubiques
        centrees sur le point de coordonnees (0,0,0)
--
*/
/*
++
| ENTREES:
|   fonction = f(x,y,z)
|   step = pas de l'integration
|   precision = coupure en dV(couronne i)/Vtot(couronnes 0 a i)
| SORTIES:
|   return: l'integrale
| REMARQUES:
|   la fonction f(x,y,z) doit decroitre tres vite en r**2
--
*/
float Int3DCube(float(*fonction)(float,float,float),float step,float precision)
{

  float f;
  double vol,volprec;
  int ecart,i,j,k,l;
  int encore,deb=0;

  vol = volprec = 0.;

  /* contribution du cube central */
  for(l=0;l<nd3d;l++) {
    f = fonction(dx3d[l]*step,dy3d[l]*step,dz3d[l]*step);
    vol += f*w3d[l];
  }
  ecart = 1;

  /* increment en couronnes cubiques de 2*ecart+1 de cote */  
  do {
    volprec = vol;
    for (i= -ecart;i<=ecart;i++) for(l=0;l<nd3d;l++) {
      for(j= -ecart;j<=ecart;j++) {
        f=fonction((i+dx3d[l])*step,(j+dy3d[l])*step,(-ecart+dz3d[l])*step);
        vol += f*w3d[l];
        f=fonction((i+dx3d[l])*step,(j+dy3d[l])*step,( ecart+dz3d[l])*step);
        vol += f*w3d[l];
      }

      for(k= -ecart+1;k<=ecart-1;k++) {
        f=fonction((i+dx3d[l])*step,(-ecart+dy3d[l])*step,(k+dz3d[l])*step);
        vol += f*w3d[l];
        f=fonction((i+dx3d[l])*step,( ecart+dy3d[l])*step,(k+dz3d[l])*step);
        vol += f*w3d[l];
      }
    }
    for (j= -ecart+1;j<=ecart-1;j++) for(k= -ecart+1;k<=ecart-1;k++) 
      for(l=0;l<nd3d;l++) {
        f=fonction((-ecart+dx3d[l])*step,(j+dy3d[l])*step,(k+dz3d[l])*step);
        vol += f*w3d[l];
        f=fonction(( ecart+dx3d[l])*step,(j+dy3d[l])*step,(k+dz3d[l])*step);
        vol += f*w3d[l];
    }

    ecart = ecart + 1;
    encore = ((fabs((vol-volprec)/vol) > precision) || (ecart <= 2) );
    if(deb>1) printf("ec=%d v=%g deltav/v=%g\n"
           ,ecart-1,vol*step*step*step,((vol-volprec)/vol));
  } while(encore);

  vol *= step * step* step;
  if(deb>0) printf("\nVol = %g\n",vol); 
  return(vol);
}


/*=================================================================*/
/*                                             cmv  6/11/97        */
/*
++
int IntFaisDr -
  (int n,float *cs,float *sn,float *p,float *x0,float *y0 -
  ,int npass,float perclean,int lp)

        Pour calculer le point d'intersection moyen
        d'un faisceau de droites.
--
*/
/*
++
| L'equation des droites est sous la forme:
|      cos(alpha)*x + sin(alpha)*y = p
| [cos(alpha),sin(alpha)] coordonnees du vecteur unitaire joignant
| l'origine au point de moindre approche de la droite a l'origine.
| p est la distance de la droite a l'origine.
--
*/
/*
++
| -- Input:
|  n  : nombre de droites
|  cs : tableau des coefficients cos(alpha)
|  sn : tableau des coefficients sin(alpha)
|  p  : tableau des coefficients p
|  npass    : nombre de passes de nettoyage
|  perclean : pourcentage d'etoiles a tuer
|             d'une passe a la suivante
|      ex: perclean=0.2, si il y a 200 etoiles dans une passe
|      la passe suivante en utilisera 160 et la suivante 128.
|  lp : niveau d'impression
--
*/
/*
++
| -- Output:
|  x0 : valeur de l'abscisse du point d'intersection moyen trouve
|  y0 :             ordonnee
| -- Return:
|  nombre de droites utilisees pour le calcul de x0,y0
--
*/
int IntFaisDr(int n,float *cs,float *sn,float *p,float *x0,float *y0
             ,int npass,float perclean,int lp)
{
int i,j,ns, nsuse=-9, ipass;
double sumC2,sumS2,sumCS,sumCP,sumSP,den;
float *d2=NULL, d2cut=-1.;
unsigned short *good=NULL;

*x0 = *y0 = 0.;
if(n<=1) return(-1);
if(lp<0) lp=0;\
if(npass<=0) npass=1;
if(perclean<=0.) { perclean=0.; npass=1;}
if(lp)
  printf("IntFaisDr avec %d points en %d passes avec nettoyage de %g%%\n"
        ,n,npass,perclean*100);

d2 = (float *) malloc(n * sizeof(float));
if(d2==NULL) {
  if(lp)
    printf("allocation de %d float impossible\n",n);
  return(-2);
}
good = (unsigned short *) malloc(n * sizeof(unsigned short));
if(good==NULL) {
  if(lp)
    printf("allocation de %d unsigned short impossible\n",n);
  free(d2);
  return(-3);
}
for(i=0;i<n;i++) good[i]=1;

for(ipass=1;ipass<=npass;ipass++) {

  /* Calcul du point d'intersection pour la passe ipass */
  sumC2=sumS2=sumCS=sumCP=sumSP=0.;
  ns=0;
  for(i=0;i<n;i++) {
    if( !good[i] ) continue;
    sumC2 += cs[i]*cs[i];
    sumS2 += sn[i]*sn[i];
    sumCS += cs[i]*sn[i];
    sumCP += cs[i]*p[i];
    sumSP += sn[i]*p[i];
    ns++;
  }
  den = sumCS*sumCS - sumC2*sumS2;
  if(den==0.) {
    if(lp)
      printf("denominateur nul pass=%d, sumCS=%g sumC2=%g sumS2=%g pour %d dr\n"
            ,ipass,sumCS,sumC2,sumS2,ns);
    free(d2); free(good);
    return(-100-ipass);
  }
  if(ns<2) {
    if(lp) printf("Pas ou plus assez de droites %d\n",ns);
    return(nsuse);
  }
  *x0 = (sumSP*sumCS - sumCP*sumS2)/den;
  *y0 = (sumCP*sumCS - sumSP*sumC2)/den;
  nsuse = ns;
  if(lp) printf("Pass=%d, %d dr x0=%g y0=%g\n",ipass,nsuse,*x0,*y0);
  if(ipass==npass) break;

  /* Calcul de la coupure en distance**2 pour la passe suivante */
  ns = 0;
  for(i=0;i<n;i++) {
    if( !good[i] ) continue;
    d2[ns] = (*x0*cs[i] * *y0*sn[i]-p[i]);
    d2[ns] *= d2[ns];
    ns++;
  }  
  qsort(d2,(size_t) ns,sizeof(float),qSort_Float);
  j = (int) (ns-ns*perclean); if(j<0) j=0;  if(j>=ns) j=ns-1;
  d2cut = d2[j];
  ns=0;
  for(i=0;i<n;i++) {
    if( !good[i] ) continue;
    den = (*x0*cs[i] * *y0*sn[i]-p[i]);
    den *= den;
    if( den>d2cut ) good[i]=0; else ns++;
  }
  if(lp>1)
    printf("  next pass: i=%d  cut=%g -> restera %d dr\n",j,d2cut,ns);
}

free(d2); free(good);
return(nsuse);
}

/*=========================================== cmv 22/11/04 ===*/
/*
++
double log_factoriel(unsigned int n)

        Compute the neperian log of n!
--
*/
double log_factoriel(unsigned int n)
{
  unsigned int i;
  double sum=0.;
  if(n<2) return sum;
  for(i=n;i>1;i--) sum += log((double)i);
  return sum;
}

/*
++
double log_stirling(unsigned int n)

        Compute the neperian log of the Stirling approx of n!
--
*/
double log_stirling(unsigned int n)
/*
Approx: n! ~ sqrt(2Pi n) n^n exp(-n) = sqrt(2Pi) n^(n+1/2) exp(-n)
        log(n!) ~ (n+1/2) log(n) -n + 1/2 log(2Pi)
(ordre 0: log(n!) ~ n (log(n)-1)
*/
{
  if(n<1) return 0.;
  return (n+0.5)*log(n) - n + Hln2pi;
}

/*
++
double log_gosper(unsigned int n);

        Compute the neperian log of the Gosper approx of n!
--
*/
double log_gosper(unsigned int n)
/*
Approx: n! ~ sqrt( (2n+1/3) Pi ) n^n exp(-n)
        log(n!) ~ 1/2 log( (2n+1/3) Pi ) + n log(n) -n
*/
{
  if(n<1) return 0.;
  return 0.5*log((2.*n+1./3.)*M_PI) + n*log(n) -n;
}

/*
++
double give_log_factoriel(unsigned int n);

        Compute the neperian log of the approx of n!
--
*/
static short give_log_factoriel_OK = 0;
#define give_log_factoriel_LIM 21
static double give_log_factoriel_TAB[give_log_factoriel_LIM];
double give_log_factoriel(unsigned int n)
{
  int i;
  if(give_log_factoriel_OK==0) {
    give_log_factoriel_OK=1;
    for(i=0;i<give_log_factoriel_LIM;i++)
      give_log_factoriel_TAB[i]=log_factoriel(i);
  }

  if(n<give_log_factoriel_LIM) return give_log_factoriel_TAB[n];
  return log_gosper(n);
}
#undef give_log_factoriel_LIM