/* quelques fonctions mathematiques et constantes utiles */
/*                                         cmv  23/06/94 */
#ifndef NBMATH_H_SEEN
#define NBMATH_H_SEEN

#include "machdefs.h"

#ifdef __cplusplus
extern "C" {
#endif

/* EA conflit entre les diverses copies de nbmath.h, nbconst.h... A nettoyer... */
#undef  Pi   
#undef  Pis2
#undef  DeuxPi 
#undef  SPi   
#undef  S2Pi   
#undef  Rac2   
#undef  Log2   
#undef  LnPi   
#undef  LgPi   
#undef  Ln10   
#undef  DftoDm 
#undef  Hln2pi 
#undef  JourSec 
#undef  AnSec  

#undef  GRAND  
#undef  GRAND2 
#undef  IGRAND 

#define  Pi     (double) ( 3.14159265358979310e+00)  /* c'est Pi */
#define  Pis2   (double) (1.57079632679489655)       /* c'est Pi/2 */
#define  DeuxPi (double) ( 6.28318530717958620e+00)  /* c'est 2*Pi */
#define  SPi    (double) ( 1.77245385090551590e+00)  /* c'est sqrt(Pi) */
#define  S2Pi   (double) ( 2.50662827463100020e+00)  /* c'est sqrt(2*Pi) */
#define  Rac2   (double) ( 1.41421356237309510e+00)  /* c'est sqrt(2) */
#define  Log2   (double) ( 3.01029995663981200e-01)  /* c'est log10(2) */
#define  LnPi   (double) ( 1.14472988584940020e+00)  /* c'est ln(Pi) */
#define  LgPi   (double) ( 4.97149872694133850e-01)  /* c'est log10(Pi) */
#define  Ln10   (double) ( 2.30258509299404590e+00)  /* c'est ln(10) */
#define  DftoDm (double) ( 1.08573620475812940e+00)  /* c'est 2.5/ln(10) */
#define  Hln2pi (double) ( 9.18938533204672670e-01)  /* c'est Ln(2*pi)/2 */

#define  JourSec (float) (86400.000000) /* nombre de secondes dans 24H */
#define  AnSec  (int) (31557600) /* nombre de secondes dans 365.25 jours */

#define  GRAND  (float)  (1.e+35)
#define  GRAND2 (double) (1.e+35)
#define  IGRAND (int_4) (2147483647)

/* pour le fit 1-dim */
#define NPAR_FITFUN 10
#define FITFUN_FLOAT 1
#define FITFUN_DOUBLE 2
union FloatDouble { float  *fx; double *dx; };
typedef union FloatDouble FLOATDOUBLE;

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);

double GausPiv(double *A,int_4 nca,int_4 n,double *B,int_4 ncb,int_4 m,int_4 Inv);

int_4 paramga(float sxin, float syin, float rhoin
             , float *smax, float *axisrat, float *tiltdeg);
int_4 gaparam(float smax, float axisrat, float tiltdeg
             , float *sxin, float *syin, float *rhoin);

double nberfc(double x);
float probnb(float sci2,int_4 inddl,int_4 *ipass);
double nbgammln(double xx);
double nbgammq(double a,double x);
double nbgammp(double a,double x);

void Set_Ihoq(int degre,int *mIhoqN,double *IhoqNX,double *IhoqNW);
int Get_Ihoq(int degre,double *IhoqNX,double *IhoqNW);
void nbgauleg(double x1,double x2,double *x,double *w,int n);
double Integ_Fun(double xmin,double xmax,double (*fonc)(double),int npas);
double Integ_Fun_2D(double (*fonc)(double x,double y)
                ,double xmin,double xmax,double ymin,double ymax
                ,int npasx,int npasy);

void Set_FitFunDPol(int DPol);
double Gauss1DPolF(double x,double *Par,double *DgDpar);
double GaussI1DPol(double x,double *Par,double *DgDpar);
double Polyn1D(double x,double *Par,double *DgDpar); 

double FitProp(double *x,double *y,double *ey,int *n,double *a1);
double FitLin(double *x,double *y,double *ey,int *n
             ,double *a0,double *a1);
double FitPar(double *x,double *y,double *ey,int *n
             ,double *a0,double *a1,double *a2);
double FitParLin(double *xx,double *y,double *ey,int *n,double x0,int Deg_d
                ,double *a0,double *a1,double *b1,double *b2);
double FitPropClean(double *x,double *y,double *ey,int *n
                  ,double per_clean,double *a1);
double FitLinClean(double *x,double *y,double *ey,int *n
                  ,double per_clean,double *a0,double *a1);
double FitParClean(double *x,double *y,double *ey,int *n
                  ,double per_clean,double *a0,double *a1,double *a2);

float interpole(float x,int n,float *t1,float *t2);

int_4 MeanSig(float *mag,float *err,int_4 *ndata,float nsigma,float *mean,float *sigma);
int Most_Probable( float *val, int nval, float binin, int nmoy
                 , float low, float high, float *most, int deb);
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);
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 Ajust_MnSg(int Ns,float *fcfr
           ,float *mean,float *emean,float *sigma,float *esigma
           ,float perc_net,float n_sigma,int deb);

float Int3DCube(float(*fonction)(float, float, float), float step, 
	       float precision);

int IntFaisDr(int n,float *cs,float *sn,float *p,float *x0,float *y0
             ,int npass,float perclean,int lp);

#ifdef __cplusplus
}
#endif

#endif