#include <stdlib.h>
#include <stdio.h>
#include <ctype.h>
#include <math.h>
#include <iostream>
#include <typeinfo>

#include "machdefs.h"

#include "fmath.h"

#include "fsvcache.h"
#include "timing.h"

#include "fsvst.h"

#include "pexceptions.h"  // SOPHYA Exceptions

#include "histats.h"  // SOPHYA HiStats module include files
#include "sopnamsp.h"  // using SOPHYA namespace

#include "machdefs.h"

#include "tvector.h"
#include "srandgen.h"
#include "fioarr.h"
#include "sopemtx.h"
#include "pexceptions.h"
#include "matharr.h"

#include "sambainit.h"

#include "lcproc.h"  // declaration of LightCurveProc class


DataTable  dtclum(int numet, int nmes, TIMEINFO *tminf, MESUREU *mesu);
DataTable dttimeinfo(int nmes, TIMEINFO *tminf);
DataTable createdtmes(int nmes); 
void dtmes(DataTable& dtmes, int nmes, STARINFO *sti, MESUREU *mesu);


struct resolution
  {
  int nbin;
  float xmin, xmax, dxbin;
  int *nst, *nstr;
  float *fmoy, *fmoyr;
  float *dflx, *sflx, *dflxe, *sflxe;
  };

typedef struct resolution  SRESOL;


SRESOL *resolinit(int nmes, int nbin, float xmin, float xmax);
int resolst(SRESOL *sres, int nmes, STARINFO *sti, MESUREU *mesu);
int resolend(SRESOL *sres, int nmes);




/* ------------------------------------------------------------------- */



const char *Pgnm = "hanafsv";

int main(int narg, char *arg[] )
{
char strg[128],*svf;
SUIVIFIP *sfip;
int numes1,numes2;
int numet1, numet2, incn;
float xmin, xmax;
float ymin, ymax;

int i,k,jt,rc,dnp,fgq;
int  nbstar,nbmes;
int nstar, nmes;

int prtl,prfg,clfg,ntfg;

int nbin;
float minlog,maxlog;


GLOBINFO glinf;
STARINFO  star;
TIMEINFO *tminf, *tmi;

TIMEINFOU *tminu;
MESURE *mes;
MESUREU *mesu;

double *TStart;

SRESOL *sres;

/* ...................................................  */

InitTim();

if (narg > 1)
 if (strcmp(arg[1],"-h") == 0) 
    {
    printf("Usage: %s [-Flag] NomSuivi\n",Pgnm);
    puts  ("  -inq : Marque une pause apres lecture TimeInfo ");
    puts  ("  -act prfg,clfg,ntfg : Flags de controle ");
    puts  ("  -prt prtl : Niveau d'impression (0->3)");
    puts  ("  -his nbin,MinLog,MaxLog : Binning histo resolution");
    puts  ("  -mes m1-m2 : Traitement des mesures m1 a m2 ");
    puts  ("  -numet n1-n2,incn : Traitement des etoiles n1 a n2 (C1)");
    puts  ("  -xpos min-max : Coupure min <= XPosEtoile <= max (C2)");
    puts  ("  -ypos min-max : Coupure min <= YPosEtoile <= max (C3)");
    puts  ("  On fait C1 & C2 & C3, Valeur par defaut = Tout");
    puts  ("  prfg (PrintFlag) : 0 , 1= Glob+TimeInfo"); 
    puts  ("     2= +StarInfo, 3= +MesuresU (=default), 4= +Mesures");
    puts  ("  clfg (Flag Courbe de lumiere) : 0(=defaut),1"); 
    puts  ("  ntfg (Bits Flag Ntuple): B0(=defaut) TimeInfo, B1 Resol, B2 Mesures\n"); 
    return(0);
    }

/*    Decodage des arguments   */
numes1 = 1;  numes2 = 9999999;
numet1 = 1;  numet2 = 9999999;  incn = 1;
xmin = -9.e19;  xmax = 9.e19;
ymin = -9.e19;  ymax = 9.e19;
prtl = 0;  prfg = 3;
clfg = 0;  ntfg = 1;
nbin = 30;
minlog = 2.0;  maxlog = 5.0;
fgq = 0;   /* On continue jusqu'au bout par defaut */

if (narg < 2) goto PbArg;
k = 1;
while (*(arg[k]) == '-')      
  {
  if ( strcmp(arg[k],"-inq") == 0)  fgq = 1;
  else if ( strcmp(arg[k],"-mes") == 0)
    { if (k == narg-1)  goto PbArg;
    sscanf(arg[k+1],"%d-%d",&numes1,&numes2);   k++; }
  else if ( strcmp(arg[k],"-numet") == 0)
    { if (k == narg-1)  goto PbArg;
    sscanf(arg[k+1],"%d-%d,%d",&numet1,&numet2,&incn);   k++; }
  else if ( strcmp(arg[k],"-xpos") == 0)
    { if (k == narg-1)  goto PbArg;
    sscanf(arg[k+1],"%g-%g",&xmin,&xmax);   k++; }
  else if ( strcmp(arg[k],"-ypos") == 0)
    { if (k == narg-1)  goto PbArg;
    sscanf(arg[k+1],"%g-%g",&ymin,&ymax);   k++; }
  else if ( strcmp(arg[k],"-prt") == 0)
    { if (k == narg-1)  goto PbArg;
    sscanf(arg[k+1],"%d",&prtl);   k++; }
  else if ( strcmp(arg[k],"-act") == 0)
    { if (k == narg-1)  goto PbArg;
    sscanf(arg[k+1],"%d,%d,%d",&prfg,&clfg,&ntfg);   k++; }
  else if ( strcmp(arg[k],"-his") == 0)
    { if (k == narg-1)  goto PbArg;
    sscanf(arg[k+1],"%d,%g,%g",&nbin,&minlog,&maxlog);   k++; }

  if (++k >= narg)  break;
  }

if (k < narg)  goto Suite;
PbArg:
printf("%s_Erreur : Pb. arguments ( %s -h pour aide) \n", Pgnm, Pgnm);
return(1);

Suite:
svf = arg[k];

cout << " hanafsv: Args decoded, starting suivi processing ..." << endl;
SophyaInit();
sfip = NULL;  rc = 0;
try {

  POutPersist po("LCInfo.ppf");
  POutPersist po2("TimInfo.ppf");
  POutPersist po3("Mesu.ppf");

  if ( (sfip = SuiviOpen(svf,SUOF_RO_MEM2)) == NULL )   return(5);
  
  nbstar = SuiviGetNbStars(sfip);
  nbmes = SuiviGetNbMesures(sfip);
  
  if (numes1 < 1)  numes1 = 1;
  if (numes2 < 1)  numes2 = 1;
  if (numes2 > nbmes)  numes2 = nbmes;
  if (numes2 < numes1) numes1 = numes2;
  nmes = numes2-numes1+1;
  if (numet1 < 1)  numet1 = 1;
  if (numet2 < 1)  numet2 = 1;
  if (numet2 > nbstar)  numet2 = nbstar;
  if (numet2 < numet1) numet1 = numet2;
  nstar = numet2 - numet1 +1;
  
  //--> new  tminf = (TIMEINFO *)malloc(nmes*sizeof(TIMEINFO));
  tminf = new TIMEINFO[nmes]; 
  //--> new  tminu = (TIMEINFOU *)malloc(nmes*sizeof(TIMEINFOU));
  tminu = new TIMEINFOU[nmes]; 
  //  mes = (MESURE *)malloc(nmes*sizeof(MESURE));
  mes = new MESURE[nmes]; 
  //--> new  mesu = (MESUREU *)malloc(nmes*sizeof(MESUREU));
  mesu = new MESUREU[nmes]; 
  //-> new TStart = (double *) malloc(nmes*sizeof(double));
  TStart = new double[nmes]; 
  if ( (tminf == NULL) || (tminu == NULL) || (mes == NULL) ||
       (mesu == NULL) || (TStart == NULL))
    { puts("hanafsv_Erreur : Pb allocation");  return(7); }

  DataTable dtm = createdtmes(nmes);

  LightCurveProc mylcp("lcproc.ppf");  // mylcp("lcproc.ppf", 512); 
  
  printf("\n ====> hanafsv : Lecture du fichier de suivi %s \n", svf);
  printf("  %d <= NumMes <= %d   %d <= NumEt <= %d (mod %d)\n", 
	 numes1, numes2, numet1, numet2, incn);
  printf("  %g <= XPosRef <= %g  ET  %g <= YPosRef <= %g \n\n",
	 xmin,xmax,ymin,ymax);
  printf("PrtFg= %d (Level=%d)  ClFg= %d  NtFg= %d   Nbin,Min,MaxLog= %d %g %g \n",
	 prfg, prtl, clfg, ntfg, nbin, minlog, maxlog);
  
  
  
  sres = NULL;
  sres =  resolinit(nmes, nbin, minlog, maxlog);
  if (sres == NULL)  { puts("hanafsv_Erreur : Pb resolinit !"); goto Fin; }
  
  
  /*  Lecture GlobInfo   */
  SuiviReadGlobInfo(sfip, (char *)(&glinf));
  if (prfg > 0)  PrtGlobInfo(&glinf, prtl);
  
  /*  Boucle lecture mesure */
  for(jt=0; jt<nmes; jt++)
    {
      tmi = tminf+jt;
      rc = SuiviReadTimeInfo(sfip,jt+numes1,jt+numes1,(char *)tmi);
      if (rc != 0) 
	{ printf("hanafsv_Erreur: Pb SuiviReadTimeInfo( , %d, ...) \n", jt+numes1);
	return(10); }
      DecodeTim(tminf+jt, tminu+jt);
      TStart[jt] = (double)tmi->TStart / (double)86400;
      if (prfg > 0)  PrtTimeInfo(tmi, jt+numes1, prtl);
      if ((tminu+jt)->FgCalib < 1)
	{
	  (mesu+jt)->Xi2 = -1.; 
	  (mesu+jt)->Flux = (mesu+jt)->FluxB = 0.; 
	  (mesu+jt)->ErrFlux = (mesu+jt)->ErrFluxB = 99999.0;
	}
    } // fin de boucle for(jt) - boucle de lecture des TimeInfo
  
  mylcp.TimeInfo(nmes, tminf);

  if (ntfg&1)  { 
    DataTable TimInf =  dttimeinfo(nmes, tminf);
    cout << " Writing TimInf to po2 (TimInfo.ppf) ... " << endl; 
    po2 << PPFNameTag("dttimeinfo") << TimInf ;  
  }
  
  
  if (fgq) {
    puts("\n\n  =====> <CR> pour continuer, <q><CR> pour stop ...");
    gets(strg); 
    if (toupper(strg[0]) == 'Q')  return(0);
  }
  
  dnp = numet2-numet1; 
  if (dnp > 20)  dnp /= 10;
  else dnp = 1;
  
  for (i=numet1; i<=numet2; i+=incn )  /*  Loop over stars in suivi file */
    { 
      
      if (((i-numet1)%dnp) == 0)
	{ sprintf(strg,"LectureStar[%d] ",i);  PrtTim(strg); }
      
      rc = SuiviReadStarInfo(sfip,i,(char *)(&star));

      if (rc != 0)  return(10);
      
      if ( (star.XPos < xmin) || (star.XPos > xmax) ||
	   (star.YPos < ymin) || (star.YPos > ymax) )   continue;
      
      if (prfg > 1)
	{
	  printf(" --------- StarInfo/Mesures Etoile [%d] ----------\n",i);
	  PrtStarInfo(&star, i, prtl);
	}
      if (star.NumEt <  0)  continue;
      if (star.FgRef == 0)  continue;
      
      
      for (jt=0; jt<nmes; jt++)
	{
	  rc = SuiviReadMesures(sfip,i,jt+numes1,jt+numes1,
				(char *)(mes+jt));
	  if (rc != 0)  
	    {
	      printf("Erreur: Star[%d] Pb lecture Mesure[%d]  (Rc=%d) \n",
		     i+1,jt+1,rc);
	      return(15);    
	    }
	  if ((tminu+jt)->FgCalib > 0)
	    {
	      DecodeMes(mes+jt, mesu+jt);
	      Calibre_F_E(mesu+jt, tminu+jt);
	      if (prfg > 3)  PrtMesure(mes+jt, jt+numes1, prtl);
	      if (prfg > 2)  PrtMesureU(mesu+jt, jt+numes1, prtl);
	    }
	  else if (prfg > 2)
	    PrtMesure(mes+jt, jt+numes1, prtl);
	  
	}     /*  Fin de boucle sur les mesures   */

      
      // a star with all its measured fluxes is ready here to be processed 
      // We process our light-curve
      mylcp.ProcessLC(i, nmes, &star, tminf, mesu); 

      if (clfg > 0){
        DataTable LiCu = dtclum(i, nmes, tminf, mesu); 
	char tagbuf[32];
	sprintf(tagbuf,"cl%d",i);
        po << PPFNameTag(tagbuf) << LiCu ;
      }

      if (ntfg&2)  
	resolst(sres, nmes, &star, mesu);
  
      if (ntfg&4)  { 
	dtmes(dtm, nmes, &star, mesu);
      }
     

    }   /*  Fin de boucle sur les etoiles  */

Fin:

  cout << " Writing dtmes to file po3 Mesu.ppf  , dtm: \n" << dtm << endl;
  po3 << PPFNameTag("dtmes") << dtm ;  
  SuiviClose(sfip);

  cout << "---- End Of Loop over stars, LightCurveProc.getDTS() : \n" << mylcp.getDTS() << endl;

  if ((ntfg&2) && (sres != NULL))
    resolend(sres,nmes);
  rc = 0;
}
catch (PThrowable & pex) {
  cerr << " hanafsv/Error - Exception catched " << (string)typeid(pex).name()
       << " - Msg= " << pex.Msg() << endl;
  rc = 98;
}
catch ( ... )
 {
   cerr << " hanafsv/Error - Exception ... catched " << endl;
   rc = 99;
}



return(rc);
}




///////////////////////// DataTabale dtclum ///////////////////////////
DataTable  dtclum(int numet, int nmes, TIMEINFO *tminf, MESUREU *mesu)
{
int idc,jt;

 DataTable dt(100);
 dt.AddFloatColumn("TStart");
 dt.AddFloatColumn("Flux");
 dt.AddFloatColumn("Xi2");
 dt.AddFloatColumn("Fond");
 dt.AddFloatColumn("ErrFlux");
 dt.AddFloatColumn("SigX");
 dt.AddFloatColumn("SugY");
 dt.AddFloatColumn("Rho");
 dt.AddFloatColumn("DelX");
 dt.AddFloatColumn("DelY");
 dt.AddFloatColumn("X");
 dt.AddFloatColumn("Y");

 double xnt[15];

idc = numet;

for(jt=0; jt<nmes; jt++)
  { 
  xnt[0] = tminf[jt].TStart/86400.;
  xnt[1] = (double) mesu[jt].Flux;
  xnt[2] = (double) mesu[jt].Xi2;
  xnt[3] = (double) mesu[jt].Fond;
  xnt[4] = (double) mesu[jt].ErrFlux;
  xnt[5] = (double) tminf[jt].SigX;
  xnt[6] = (double) tminf[jt].SigY;
  xnt[7] = (double) tminf[jt].Rho;
  xnt[8] = (double) tminf[jt].DelX;
  xnt[9] = (double) tminf[jt].DelY;
  xnt[10] = (double) mesu[jt].X ;
  xnt[11] = (double) mesu[jt].Y ;
  dt.AddRow(xnt);
  }

return(dt);
}



//////////////////// DataTable dttimeinfo //////////////////////////
DataTable dttimeinfo(int nmes, TIMEINFO *tminf)
{
int idn,jt,i;

 DataTable dt(10000);
 dt.AddFloatColumn("NumPhoto");
 dt.AddFloatColumn("TStart");
 dt.AddFloatColumn("Expose");
 dt.AddFloatColumn("Nb0kGF");
 dt.AddFloatColumn("Nb0kPS");
 dt.AddFloatColumn("Fond");
 dt.AddFloatColumn("SigFond");
 dt.AddFloatColumn("SigX");
 dt.AddFloatColumn("SigY");
 dt.AddFloatColumn("Rho");
 dt.AddFloatColumn("SgX");
 dt.AddFloatColumn("SgY");
 dt.AddFloatColumn("DelX");
 dt.AddFloatColumn("DelY");
 dt.AddFloatColumn("Absorption");
 dt.AddFloatColumn("Airmass");
 dt.AddFloatColumn("FgCalib");
 
  for(i=0; i<6; i++){
   char calname[32];
   sprintf(calname,"Calib%d",i);
   dt.AddFloatColumn(calname);
  }
 
 for(i=0; i<6; i++)
   for(int j =0; j<3;j++){
     char name[32];
     sprintf(name,"Resol%d%d",j,i);
      dt.AddFloatColumn(name);
   }
 

 double xnt[48];

 idn = 300;
 for(i=0; i<36; i++)  xnt[i] = 0.;

for(jt=0; jt<nmes; jt++)
  {
  xnt[0] = tminf[jt].NumPhoto;
  xnt[1] = tminf[jt].TStart/86400.;
  xnt[2] = tminf[jt].Expose;
  xnt[3] = tminf[jt].NbOkGF;
  xnt[4] = tminf[jt].NbOkPS;
  xnt[5] = tminf[jt].Fond;
  xnt[6] = tminf[jt].SigFond;

  xnt[7] = tminf[jt].SigX;
  xnt[8] = tminf[jt].SigY;
  xnt[9] = tminf[jt].Rho;
  xnt[10] = tminf[jt].SgX;
  xnt[11] = tminf[jt].SgY;  
  xnt[12] = tminf[jt].DelX;
  xnt[13] = tminf[jt].DelY;
  xnt[14] = tminf[jt].Absorption;
  xnt[15] = tminf[jt].AirMass;

  xnt[16] = tminf[jt].FgCalib;
  for(i=0; i<6; i++) 
    {
    xnt[17+i] = tminf[jt].Calib[i];
    xnt[23+i*3] = tminf[jt].Resol[0][i];
    xnt[24+i*3] = tminf[jt].Resol[1][i]; 
    xnt[25+i*3] = tminf[jt].Resol[2][i]; 
    }
  dt.AddRow(xnt);
  }

return(dt);
}



////////////////// DataTable creatdtmes ////////////////////////
DataTable createdtmes(int nmes) 

{

 DataTable dt(100);
 dt.AddFloatColumn("NumEt");
 dt.AddFloatColumn("FluxRef");
 dt.AddFloatColumn("XPos");
 dt.AddFloatColumn("YPos");
 dt.AddFloatColumn("FgRef");
 dt.AddFloatColumn("NbVois");
 dt.AddFloatColumn("DisMin");
 dt.AddFloatColumn("DisM2");

 if (nmes > 10) nmes = 10;

 char nflux[32]; char nerrflux[32]; char nfluxb[32]; char nerrfluxb[32];
 char nfond[32]; char nxi2[32]; char nx[32]; char ny[32]; char nsigx[32];
 char nsigy[32]; char nflx[32]; char nfnd[32]; char nxi[32];

 for (int k=0; k<nmes; k++) {
   sprintf(nflux,"Flux%d",k); sprintf(nerrflux,"ErrFlux%d",k);
   sprintf(nfluxb,"FluxB%d",k); sprintf(nerrfluxb,"ErrFluxB%d",k);
   sprintf(nfond,"Fond%d",k); sprintf(nxi2,"Xi2%d",k); sprintf(nx,"X%d",k);
   sprintf(ny,"Y%d",k); sprintf(nsigx,"SigX%d",k); sprintf(nsigy,"SigY%d",k);
   sprintf(nflx,"Flx%d",k); sprintf(nfnd,"Fnd%d",k); sprintf(nxi,"Xi%d",k); 

   dt.AddFloatColumn(nflux);
   dt.AddFloatColumn(nerrflux);
   dt.AddFloatColumn(nfluxb);
   dt.AddFloatColumn(nerrfluxb);
   dt.AddFloatColumn(nfond);
   dt.AddFloatColumn(nxi2);
   dt.AddFloatColumn(nx);
   dt.AddFloatColumn(ny);
   dt.AddFloatColumn(nsigx);
   dt.AddFloatColumn(nsigy);
   dt.AddFloatColumn(nflx);
   dt.AddFloatColumn(nfnd);
   dt.AddFloatColumn(nxi);
 }
 return dt;
}

void dtmes(DataTable& dt, int nmes, STARINFO *sti, MESUREU *mesu)
{
int idn,jt,off,sz;


double xnt[200];


//idn = 400;
if (nmes > 10)  nmes=10;

xnt[0] = sti->NumEt;
xnt[1] = sti->FluxRef;
xnt[2] = sti->XPos;
xnt[3] = sti->YPos;
xnt[4] = sti->FgRef;
xnt[5] = sti->NbVois;
xnt[6] = sti->DisMin; 
xnt[7] = sti->DisM2;

off = 7;  sz = 13;
if (nmes > 10) nmes = 10;
for(jt=0; jt<nmes; jt++)
  {
  xnt[off+sz*jt+1] = mesu[jt].Flux;
  xnt[off+sz*jt+2] = mesu[jt].ErrFlux;
  xnt[off+sz*jt+3] = mesu[jt].FluxB;
  xnt[off+sz*jt+4] = mesu[jt].ErrFluxB;
  xnt[off+sz*jt+5] = mesu[jt].Fond;
  xnt[off+sz*jt+6] = mesu[jt].Xi2;
  xnt[off+sz*jt+7] = mesu[jt].X;
  xnt[off+sz*jt+8] = mesu[jt].Y;
  xnt[off+sz*jt+9] = mesu[jt].SigX;
  xnt[off+sz*jt+10] = mesu[jt].SigY;
  xnt[off+sz*jt+11] = mesu[jt].Flx0;
  xnt[off+sz*jt+12] = mesu[jt].Fnd0;
  xnt[off+sz*jt+13] = mesu[jt].Xi20;
  }
 dt.AddRow(xnt);

return;
}



/* Nouvelle-Fonction */

SRESOL *resolinit(int nmes, int nbin, float xmin, float xmax)

{
int i,j;
SRESOL *sres;

if (nbin < 1)  return(NULL);
if (xmax <= xmin)   return(NULL);

//-> new sres = malloc(nmes*sizeof(SRESOL));
sres = new SRESOL[nmes]; 
if (sres == NULL)  return(NULL);
for(i=0; i<nmes; i++)
  {
  sres[i].nbin = nbin;
  sres[i].xmin = xmin;
  sres[i].xmax = xmax;
  sres[i].dxbin = (xmax-xmin)/(float)nbin;
  //-> new  sres[i].nst = malloc(nbin*sizeof(int));
  sres[i].nst = new int[nbin]; 
  if (sres[i].nst == NULL)   break;
  //-> new  sres[i].nstr = malloc((nbin+2)*sizeof(int));
  sres[i].nstr = new int[nbin+2]; 
  if (sres[i].nstr == NULL)   
    { delete[] sres[i].nst;  break; }
  //->new  sres[i].fmoy =  malloc(nbin*sizeof(float));
  sres[i].fmoy =  new float[nbin]; 
  if (sres[i].fmoy == NULL)   
    { delete[] sres[i].nst; delete[] sres[i].nstr;  break; }
  //->new  sres[i].fmoyr =  malloc(nbin*sizeof(float));
  sres[i].fmoyr =  new float[nbin] ; 
  if (sres[i].fmoyr == NULL)   
    { delete[] sres[i].nst; delete[] sres[i].nstr; 
    delete[] sres[i].fmoy; break; }
  //-> new  sres[i].dflx =  malloc(nbin*sizeof(float));
  sres[i].dflx =  new float[nbin]; 
  if (sres[i].dflx == NULL)   
    { delete[] sres[i].nst;  delete[] sres[i].nstr;
    delete[] sres[i].fmoy;  delete[] sres[i].fmoyr;  break; }
  //->new  sres[i].sflx =  malloc(nbin*sizeof(float));
  sres[i].sflx =  new float[nbin]; 
  if (sres[i].sflx == NULL)   
    { delete[] sres[i].nst; delete[] sres[i].nstr; 
      delete[] sres[i].fmoy;  delete[] sres[i].fmoyr;
      delete[] sres[i].dflx;  break; }
  //->new  sres[i].dflxe =  malloc(nbin*sizeof(float));
  sres[i].dflxe =  new float[nbin]; 
  if (sres[i].dflxe == NULL)   
    { delete[] sres[i].nst;  delete[] sres[i].nstr; 
    delete[] sres[i].fmoy; delete[] sres[i].fmoyr; 
    delete[] sres[i].dflx;  delete[] sres[i].sflx;  break; }
  //->new  sres[i].sflxe =  malloc(nbin*sizeof(float));
  sres[i].sflxe =  new float[nbin]; 
  if (sres[i].sflxe == NULL)   
    { delete[] sres[i].nst;  delete[] sres[i].nstr;
    delete[] sres[i].fmoy;  delete[] sres[i].fmoyr; 
    delete[] sres[i].dflx;  delete[] sres[i].sflx;  
    delete[] sres[i].dflxe;  break; }
  }

if (i != nmes)
  {
/* BUG REZA ???  puts("resolinit_Erreur : Pb malloc/new  %d %d \n", i, nmes); */
  printf("resolinit_Erreur : Pb malloc  %d %d \n", i, nmes);
  for(j=0; j<i; j++)
    {
    delete[] sres[j].nst; 
    delete[] sres[j].nstr; 
    delete[] sres[j].fmoy;
    delete[] sres[j].fmoyr;
    delete[] sres[j].dflx;
    delete[] sres[j].sflx;
    delete[] sres[j].dflxe;
    delete[] sres[j].sflxe;    
    }
  delete[] sres; 
  return(NULL);
  }
for(i=0; i<nmes; i++)
  {
  for(j=0; j<nbin; j++)
    {
    sres[i].nst[j] = sres[i].nstr[j] = 0 ; 
    sres[i].fmoy[j] = sres[i].fmoyr[j] = 0. ;
    sres[i].dflx[j] = 0. ;
    sres[i].sflx[j] = 0. ;
    sres[i].dflxe[j] = 0. ;
    sres[i].sflxe[j] = 0. ;    
    }
  }

return(sres);
}


/* Nouvelle-Fonction */

int resolst(SRESOL *sres, int nmes, STARINFO *sti, MESUREU *mesu)
{
int i,j,ib;
float err, flx, flxr;
float xlg, dx;

flxr = sti->FluxRef;
if (flxr < 1.)  return(-1);
xlg = (float)log10((double)flxr);

for(j=0; j<nmes; j++)
  {
  ib = (int)((xlg-sres[j].xmin)/sres[j].dxbin);
  if (ib < 0)  sres[j].nstr[sres[j].nbin]++;
  if (ib >= sres[j].nbin)  sres[j].nstr[sres[j].nbin+1]++;
  if ( (ib < 0) || (ib >= sres[j].nbin) )   continue;
  sres[j].nstr[ib]++;
  sres[j].fmoyr[ib] += flxr;
  if (mesu[j].Xi2 < 0.)  continue;
  err = mesu[j].ErrFlux;
  flx = mesu[j].Flux;
  if (err/flxr < 0.01)  err = 0.01*flxr;
  /* if (fabs((double)(err/flx)) > (1./xlg))  continue; */
  if (fabs((double)((flx-flxr)/err)) > 8.)  continue;
  sres[j].nst[ib]++;
  sres[j].fmoy[ib] += flx;
  dx = flxr-flx;
  sres[j].dflx[ib] += dx;
  sres[j].sflx[ib] += (dx*dx);
  dx /= err;
  sres[j].dflxe[ib] += dx;
  sres[j].sflxe[ib] += (dx*dx);  
  }

return(ib);
}

/* Nouvelle-Fonction */
int resolend(SRESOL *sres, int nmes)

{
int i,j,ib,idn;
float ns;
float xnt[10];

printf("Resolend_Info: Under/OverFlow  %d %d \n", 
       sres[0].nstr[sres[0].nbin], sres[0].nstr[sres[0].nbin+1]);
for(j=0; j<nmes; j++)
  {
  idn = 501+j;
  for(ib=0; ib<sres[j].nbin; ib++)
    {
    xnt[0] = ib;   
    xnt[1] = sres[j].xmin+ib*sres[j].dxbin; 
    for(i=2; i<10; i++)  xnt[i] = 0.;
    if (sres[j].nstr[ib] < 1)   goto  Fill;
    sres[j].fmoyr[ib] /= (float)sres[j].nstr[ib];
    xnt[8] = sres[j].nstr[ib];
    xnt[9] = sres[j].fmoyr[ib];
    if (sres[j].nst[ib] < 1)   goto  Fill;
    xnt[2] = ns = sres[j].nst[ib];
    sres[j].fmoy[ib] /= ns;
    xnt[3] = sres[j].fmoy[ib];
    sres[j].dflx[ib] /= ns;
    sres[j].sflx[ib] /= ns;
    sres[j].sflx[ib] -= sres[j].dflx[ib]*sres[j].dflx[ib];
    if (sres[j].sflx[ib] >= 0.) 
      sres[j].sflx[ib]  = (float)sqrt((double)sres[j].sflx[ib]);
    xnt[4] = sres[j].dflx[ib];
    xnt[5] = sres[j].sflx[ib];
    sres[j].dflxe[ib] /= ns;
    sres[j].sflxe[ib] /= ns;
    sres[j].sflxe[ib] -= sres[j].dflxe[ib]*sres[j].dflxe[ib];
    if (sres[j].sflxe[ib] >= 0.) 
      sres[j].sflxe[ib]  = (float)sqrt((double)sres[j].sflxe[ib]);
    xnt[6] = sres[j].dflxe[ib];
    xnt[7] = sres[j].sflxe[ib];
    Fill:
    // --- SOPHYA NTuple fill    hfn_(&idn, xnt);
    printf("Resol[%d,%d]  ns,nsr %d %d fm,df,sf= %g (%g) %g %g \n",j,ib, 
           sres[j].nst[ib], sres[j].nstr[ib],
           sres[j].fmoy[ib], sres[j].fmoyr[ib],
           sres[j].dflx[ib], sres[j].sflx[ib]);  
    }
  }

for(j=0; j<nmes; j++)
  {
  delete[] sres[j].nst; 
  delete[] sres[j].nstr; 
  delete[] sres[j].fmoy;
  delete[] sres[j].fmoyr;
  delete[] sres[j].dflx;
  delete[] sres[j].sflx;
  delete[] sres[j].dflxe;
  delete[] sres[j].sflxe;    
  }
delete[] sres; 

return (0);

}