source: Sophya/trunk/FrEROS/AnaLC/nbsread.cc@ 3891

/* lecture des fichiers de suivi:     cmv 21/12/93 */
#include "machdefs.h"
#include "sopnamsp.h"  // using SOPHYA namespace
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <signal.h>

#include "sophyainit.h"
#include "fsvcache.h"
#include "nbtri.h"
#include "nbmath.h"
#include "strutil.h"
#include "fsvst.h"

// Ne pas changer l'ordre des 3 lignes suivantes!!!
#define  FLAG_MAIN_NBSREAD
#include "nbsread.h"
#include "nbgene.h"

void INIT_NBSR(void);
void OPENSUIV(int ic);
void READST(int ic);
void READET(int ic);
void DATCLEAN(int ic);
void DATSORT(int ic);

void SELSTAR(void);
void GET_STAR_ID(void);

int  DECODESUIVINAME(char* str);
void DECODE_GLOBINFO(GLOBINFO *p);
void DECODE_STARINFO(STARINFO *p, int ic);
void DECODE_TIMEINFO(TIMEINFO *p, int ic, int imes);
void DECODE_MESURE  (MESURE   *p, int ic, int imes);

void DO_MALLOC(void);
void DO_FREE(void);

GLOBINFO     glinf;
STARINFO     star;
TIMEINFO     timinf;
MESURE       mes;

/* Pour catcher un CTRL-C ou des time limites */
static struct sigaction Nbsread_Action;
void Nbsread_SigCatchInit(void);
void Nbsread_SigCatch(int s);

void usage(int);
void usage(int lp) {
 printf("nbsread\n");
 printf("  nom du fichier des fichiers de suivi\n");
 printf("  [-h] help detaille\n");
 printf("  [-et et1,et2,etinc]\n");
 printf("  [-prt niv]\n");
 printf("  [-dbg et1,et2,ccd,niv]\n");
 printf("  [-kill mes1 mes2 mes3 ...]\n");
 printf("  [-mc ngenstar]\n");
 printf("  [-u0 U0MinI,U0MaxI,U0SimI]\n");
 printf("  [-t0 T0MinI,T0MaxI,T0SimI]\n");
 printf("  [-tau TauMinI,TauMaxI,TauSimI]\n");
 printf("  [-tf type(1ou2),UMinI,UMaxI,USimI,Usu0_Cut]\n");
 printf("  [-bl Blending (1ou2)]\n");
 printf("  [-seed seed1,seed2,seed3]\n");
 if(lp<=0) exit(-1);
 printf("\n");
 printf("et1,et2,etinc: on utilise les etoiles de et1 a et2 par pas de etinc\n");
 printf("niv: niveau de print\n");
 printf("et1,et2,ccd,niv: debug etoiles et1 a et2 ccd ccd (-1=tous) niveau niv\n");
 printf("mes1 mes2 mes3...: mesure a ne pas considerer dans les cdl\n"); 
 printf("ngenstar: nombre de generations montecarlo par etoile\n");
 printf("U0MinI,U0MaxI,U0SimI: generation du parametre d impact (def=0,0,0)\n"); 
 printf("   si U0MinI<U0MaxI on tire u0 plat entre U0MinI,U0MaxI\n"); 
 printf("   si U0MinI>=U0MaxI u0 est fixe a U0SimI\n"); 
 printf("TauMinI,TauMaxI,TauSimI: generation du parametre tau (def=0.01,365.)\n"); 
 printf("   ... meme logique que u0 ...\n"); 
 printf("T0MinI,T0MaxI,T0SimI: generation du temps du maximum (def=1,-1,-1)\n"); 
 printf("   si T0MinI<T0MaxI on tire t0 plat entre T0MinI,T0MaxI\n"); 
 printf("   si T0MinI>=T0MaxI ET T0SimI>0 t0 est fixe a T0SimI\n"); 
 printf("                     ET T0SimI<0 t0 tout l intervalle du suivi (def)\n"); 
 printf("type,UMinI,UMaxI,USimI,Usu0_Cut: pour tirage taille finie (def=0,0,0,0.001)\n");
 printf("   type = 1 taille finie tirage plat en U = Rs(proj)/Re\n"); 
 printf("          2 taille finie tirage plat en UC (U = UC*Rs)\n");
 printf("   si UMinI>=UMaxI U est fixe a U0SimI\n"); 
 printf("   si UMinI<UMaxI on tire U plat entre UMinI,UMaxI si type=1\n"); 
 printf("                                 entre UMinI*Rs,UMaxI*Rs si type=2\n"); 
 printf("Blending = 1 sur l'etoile la + brillante du couple\n");
 printf("           2 sur l'etoile la - brillante du couple\n");
 printf("Il y a possibilite de definir ses propres fct de tirage (cf nbgene.c)\n");
 printf("seed1,seed2,seed3: initialisation des aleatoires (def=36117,51106,21478)\n");
exit(-1);
}

/*==========================================================================*/
int main (int narg, char *arg[])
{
int_4 i,ic,iet1,iet2,iet0,ietmax,ietinc=1;
char str[NBSREAD_LSTR];

printf("=====> ");
for(i=0;i<narg;i++) printf("%s ",arg[i]);
printf("\n");

/* printf("main: narg=%d %s\n",narg,arg[0]); */
if (narg < 2) usage(0);
if( !strcmp(arg[1],"-h") != 0 ) usage(1);

/*------------------------------------------*/
/*-------- Initialisation Aleatoire --------*/
/*------------------------------------------*/
 SophyaInit();
 RandGen = new DR48RandGen;
 RandomGeneratorInterface::SetGlobalRandGenP(RandGen);


/*------------------------------------------*/
/*-------- Initialisation variables --------*/
/*------------------------------------------*/

INIT_NBSR();

/*---------------------------------------*/
/*-------- Lecture des arguments --------*/
/*---------------------------------------*/

i=2;
while ( i < narg ) {
  /* printf("arg[%d] %s\n",i,arg[i]); */
  if( !strcmp(arg[i],"-h") != 0 ) {
    usage(1);
  } else if( !strcmp(arg[i],"-et") != 0 ) {
    i++;
    sscanf(arg[i],"%d,%d,%d",&ietlim[0],&ietlim[1],&ietinc);
    if( ietinc<=0 ) ietinc = 1;
    printf("Etoiles traitees de %d a %d par pas de %d\n"
          ,ietlim[0],ietlim[1],ietinc);
  } else if( !strcmp(arg[i],"-prt") != 0 ) {
    i++;
    sscanf(arg[i],"%d",&ffprt);
    printf("Niveau de print %d\n",ffprt);
  } else if( !strcmp(arg[i],"-dbg") != 0 ) {
    i++;
    sscanf(arg[i],"%d,%d,%d,%d",&idebug[0],&idebug[1],&idebug[2],&idebug[3]);
    printf("Debug de l'et %d a l'et %d ccd %d niveau %d\n"
          ,idebug[0],idebug[1],idebug[2],idebug[3]);
  } else if( !strcmp(arg[i],"-mc") != 0 ) {
    i++;
    sscanf(arg[i],"%d",&mc.NGenStar);
    if(mc.NGenStar>0) mc.montecar=1;
      else {mc.montecar=mc.NGenStar; mc.NGenStar=0;}
    printf("MonteCarlo actif, %d generations par etoiles\n",mc.NGenStar);
  } else if( !strcmp(arg[i],"-u0") != 0 ) {
    i++;
    sscanf(arg[i],"%f,%f,%f",&mc.U0MinI,&mc.U0MaxI,&mc.U0SimI);
    printf("U0 generation de %g a %g ou %g\n"
          ,mc.U0MinI,mc.U0MaxI,mc.U0SimI);
  } else if( !strcmp(arg[i],"-tau") != 0 ) {
    i++;
    sscanf(arg[i],"%f,%f,%f",&mc.TauMinI,&mc.TauMaxI,&mc.TauSimI);
    printf("Tau generation de %g a %g ou %g\n"
          ,mc.TauMinI,mc.TauMaxI,mc.TauSimI);
  } else if( !strcmp(arg[i],"-t0") != 0 ) {
    i++;
    sscanf(arg[i],"%lf,%lf,%lf",&mc.T0MinI,&mc.T0MaxI,&mc.T0SimI);
    printf("T0 generation de %g a %g ou %g\n"
          ,mc.T0MinI,mc.T0MaxI,mc.T0SimI);
  } else if( !strcmp(arg[i],"-tf") != 0 ) {
    i++;
    sscanf(arg[i],"%d,%f,%f,%f,%f"
          ,&mc.Taille_Finie,&mc.UMinI,&mc.UMaxI,&mc.USimI,&mc.Usu0_Cut);
    printf("Taille finie %d generation de U %g a %g ou %g (prec=%g)\n"
          ,mc.Taille_Finie,mc.UMinI,mc.UMaxI,mc.USimI,mc.Usu0_Cut);
  } else if( !strcmp(arg[i],"-bl") != 0 ) {
    i++;
    sscanf(arg[i],"%d",&mc.Blending);
    printf("Blending %d\n",mc.Blending);
  } else if( !strcmp(arg[i],"-seed") != 0 ) {
    i++;
    sscanf(arg[i],"%d,%d,%d",&mc.iseed1,&mc.iseed2,&mc.iseed3);
    printf("Seed %d %d %d\n",mc.iseed1,mc.iseed2,mc.iseed3);
  } else if( !strcmp(arg[i],"-kill") != 0 ) {
    i++;
    while ( i<narg ) {
      if ( arg[i][0] != '-' ) {
        if ( nkillms<NKILLMSMX ) {
          sscanf(arg[i],"%d",&killms[nkillms]);
          nkillms++;
          i++;
        } else {
          printf("Trop de mesures a tuer! nkillms=%d / %d\n",nkillms,NKILLMSMX);
          exit(-1);
        }
      } else { i--; break;}
    }
    printf("Kill: %d mesures a tuer:",nkillms);
    if(nkillms>0) for(ic=0;ic<nkillms;ic++) printf(" %d",killms[ic]);
    printf("\n");
  } else {
    printf("arguments: cas non prevu: %s\n",arg[i]); usage(0);
  }
  i++;
}
ietlim[0] = ( ietlim[0] <= 0 ) ? 1 : ietlim[0] ;
ietlim[1] = ( ietlim[1] < ietlim[0] ) ? IGRAND : ietlim[1] ;
printf("Etoiles traitees de %d a %d\n",ietlim[0],ietlim[1]);

/*----------------------------------------------------------------*/
/*----- Lecture nom de fichier data et selection des etoiles -----*/
/*----------------------------------------------------------------*/

printf("Fichier de donnees suivi: %s\n",arg[1]);
if( (dataccd = fopen(arg[1],"r") ) == 0 ) {
  printf("impossible d'ouvrir .dataccd\n"); exit(-1);
}

GET_STAR_ID();

/*-----------------------------------------*/
/*-------- Initialisation generale --------*/
/*-----------------------------------------*/

/* initialisation des compteurs et des variables */

numccd = nstdeb = netoiles = 0;

rewind(dataccd);
if( fgets(str,NBSREAD_LSTR,dataccd) == 0 )
    {printf("fichier dataccd illisible\n"); exit(-1);}
nbccd = -1;
sscanf(str,"%d",&nbccd);
printf("on veut traiter %d ccd\n",nbccd);
if( nbccd <= 0 ) exit(-1);

/* initialisation USER */

Calibration_Absolue = 0;
Calibration_Absolue_UParms = 0;
UINIT();

/* initialisation montecarlo */

AMPLML( 1 );

/* gestion des interrupts */

Nbsread_SigCatchInit();

/*-------------------------------------------------------*/
/*-------- Boucle sur les CCD (ou paves Plaques) --------*/
/*-------------------------------------------------------*/

while ( fgets(str,NBSREAD_LSTR,dataccd) != 0 ) {   /* While... */

strip(str,'B',' ');
printf("%s\n",str);

if ( strstr(str,"@END") == str ) break;

if( strstr(str,"@CCD") == str) {

  numccd++;
  sscanf(&str[4],"%d",&ccdnum);
  printf("Nouveau CCD: c'est le %d ieme, de numero %d\n",numccd,ccdnum);
  szMXBYTE=0;
  nbmsMX=0;
  for (ic=0;ic<NCOULMX;ic++) {
    nbfil[ic] = iet[ic] = net[ic] = nmes[ic]= 0;
    for (i=0;i<NBFILMX;i++) nbms[ic][i] = nms1fil[ic][i] = 0;
  }

} else if( strstr(str,"@FINCCD") == str) {

  for(ic=0;ic<NCOULMX;ic++) {
    if( nbfil[ic] <= 0 ) continue;
    /* printf("Main_Debug: lunfil[%d][0]=%ld\n",ic,lunfil[ic][0]); */
    sztmi[ic] = lunfil[ic][0]->Head.RecSize[3] / lunfil[ic][0]->Head.NItem;
    szmes[ic] = lunfil[ic][0]->Head.RecSize[4] / lunfil[ic][0]->Head.NItem;
    if( szMXBYTE < sztmi[ic] ) szMXBYTE = sztmi[ic];
    if( szMXBYTE < szmes[ic] ) szMXBYTE = szmes[ic];
    if(ffprt>0) printf("coul=%3d sztmi=%d szmes=%d\n",ic+1,sztmi[ic],szmes[ic]);
  }
  if(ffprt>0) printf("szMXBYTE=%d\n",szMXBYTE);

  DO_MALLOC();
  for (ic=0;ic<NCOULMX;ic++) OPENSUIV(ic);
  for (ic=0;ic<NCOULMX;ic++) DATCLEAN(ic);
  for (ic=0;ic<NCOULMX;ic++) DATSORT (ic); /* fill indexu */
  TFIRST = Tfirst[0]; TLAST  = Tlast[0];
  for (ic=1;ic<NCOULMX;ic++) {
    if( Tfirst[ic] != -1. && Tfirst[ic]<TFIRST ) TFIRST = Tfirst[ic];
    if( Tlast[ic]  != -1. && Tlast[ic] >TLAST  ) TLAST  = Tlast[ic];
  }
  if(ffprt>0) printf("Tfirst=%.2f (1=.2%f,2=%.2f)  Tlast=%.2f (1=%.2f,2=%.2f)\n"
                    ,TFIRST,Tfirst[0],Tfirst[1],TLAST,Tlast[0],Tlast[1]);
  UINITCCD();
  AMPLML( 2 );

  iet1 = ( ietlim[0] > 0 ) ? ietlim[0] : 1;
  iet2 = ( ietlim[1] > net[0] ) ? net[0] : ietlim[1];

  for(iet0=iet1;iet0<=iet2; iet0+=ietinc) {   /* boucle sur les etoiles */

    netoiles++;

    for (ic=0;ic<NCOULMX;ic++) iet[ic]=iet0;

    debug=0;
    if( iet0>=idebug[0] && iet0<=idebug[1] && ( ccdnum==idebug[2] || idebug[2]<0) ) {
      nstdeb++;
      debug=idebug[3];
    }
    if(debug!=0) printf("\n\n************** DEBUG CCD %d ET %d LUES %d \n"
                       ,ccdnum,iet0,netoiles);

    /* lecture des starinfos */
    for (ic=0;ic<NCOULMX;ic++) {
      if( ic == 0 ) {   /*fichier non-composite ou couleur 1*/
        READST(ic);
      } else {
        iet[ic] = staru[0].XRef;
        READST(ic);
      }
    }

    /* calibration absolue */
    starcal.Mag_B= starcal.Mag_V= starcal.Mag_R
                 = starcal.Mag_I= starcal.Rstar= -GRAND;
    starcal.Rtype = -1;
    if( iet[0]>0 && iet[1]>0 && Calibration_Absolue!=0 )
      Calibration_Absolue(staru[0].FluxRef,staru[1].FluxRef,
                          &starcal,Calibration_Absolue_UParms);

    /* on saute l'etoile et la lecture des mesures ? */
    if(starcut.Mag_B_cut>0 &&
      (starcal.Mag_B<starcut.Mag_B_min || starcal.Mag_B>starcut.Mag_B_max) ) continue;
    if(starcut.Mag_V_cut>0 &&
      (starcal.Mag_V<starcut.Mag_V_min || starcal.Mag_V>starcut.Mag_V_max) ) continue;
    if(starcut.Mag_R_cut>0 &&
      (starcal.Mag_R<starcut.Mag_R_min || starcal.Mag_R>starcut.Mag_R_max) ) continue;
    if(starcut.Mag_I_cut>0 &&
      (starcal.Mag_I<starcut.Mag_I_min || starcal.Mag_I>starcut.Mag_I_max) ) continue;
    if(starcut.Rstar_cut>0 &&
      (starcal.Rstar<starcut.Rstar_min || starcal.Rstar>starcut.Rstar_max) ) continue;

    if(starcut.Rtype_cut>0 &&
      (starcal.Rtype<starcut.Rtype_min || starcal.Rtype>starcut.Rtype_max) ) continue;

    if(starcut.Mag_BV_cut>0 &&
      (starcal.Mag_B-starcal.Mag_V<starcut.Mag_BV_min ||
       starcal.Mag_B-starcal.Mag_V>starcut.Mag_BV_max) ) continue;
    if(starcut.Mag_BR_cut>0 &&
      (starcal.Mag_B-starcal.Mag_R<starcut.Mag_BR_min ||
       starcal.Mag_B-starcal.Mag_R>starcut.Mag_BR_max) ) continue;
    if(starcut.Mag_VR_cut>0 &&
      (starcal.Mag_V-starcal.Mag_R<starcut.Mag_VR_min ||
       starcal.Mag_V-starcal.Mag_R>starcut.Mag_VR_max) ) continue;
    if(starcut.Mag_RI_cut>0 &&
      (starcal.Mag_R-starcal.Mag_I<starcut.Mag_RI_min ||
       starcal.Mag_R-starcal.Mag_I>starcut.Mag_RI_max) ) continue;

    /* l'etoile a t-elle ete selectionnee ? */
    SELSTAR();

    /* lecture des mesures */
    for (ic=0;ic<NCOULMX;ic++) READET(ic);

    ietmax = (mc.montecar>0 && mc.NGenStar>0) ? mc.NGenStar : 1;
    ietgen=0;
    while( ietgen < ietmax ) {
      if(mc.montecar>0) AMPLML( 3 );
      UEVT();
      ietgen++;
    }

  }      /* fin de boucle sur les etoiles */

  UENDCCD();

  printf("on ferme les fichiers de suivi du ccd %d\n",ccdnum);
  DO_FREE();
  for(ic=0;ic<NCOULMX;ic++) {
    printf("Couleur: %d\n",ic+1);
    for(i=0;i<nbfil[ic];i++) {
      printf("fichier %d\n",i);
      SuiviClose(lunfil[ic][i]);
  } }

} else if( strstr(str,"@COUL") == str) {

  sscanf(&str[5],"%d",&ic);
  printf("Couleur: %d\n",ic);
  ic--;

} else {
  str[strlen(str)-1] = ' ';   /*  lu par fgets() */
  strip(str,'B',' ');
  DECODESUIVINAME(str);
  printf("%s\n",str);
  strcpy(lunnam[ic][nbfil[ic]], str);
  printf("lunnam[%d][nbfil[%d]] == lunnam[%d][%d] = %s\n"
        ,ic,ic,ic,nbfil[ic],lunnam[ic][nbfil[ic]]);
  if( (lunfil[ic][nbfil[ic]] = SuiviOpen(str,OpenType)) == 0 ) {
    printf("Impossible d'ouvrir le fichier %s\n",str);
    exit(-1);
  }
  /* lecture du nombre de mesures */
  nbms[ic][nbfil[ic]] = SuiviGetNbMesures(lunfil[ic][nbfil[ic]]);
  if(ffprt>1) printf("Ouverture lunfil[%3d][%5d] = %d mes\n"
                    ,ic,nbfil[ic],nbms[ic][nbfil[ic]]);
  if( nbms[ic][nbfil[ic]] > nbmsMX ) nbmsMX = nbms[ic][nbfil[ic]];
  /* numero du 1er fichier de l'unite nbfil[ic]: va de 0->nmes */
  nms1fil[ic][nbfil[ic]] = nmes[ic];
  nmes[ic] += nbms[ic][nbfil[ic]];
  nbfil[ic]++;

}
}   /* Fin du While... */

UEND();
AMPLML( 5 );
if( StId1 != 0 ) free(StId1);
if( StId2 != 0 ) free(StId2);
fclose(dataccd);

printf("Nombre total d'etoiles lues %d\n",netoiles);
delete RandGen;

return(0);
}

/*==========================================================================*/
void INIT_NBSR(void)
{
int i;

OpenType= SUOF_RO_MEM2;
GRAND2_KILL = 5. * GRAND2;

ietlim[0]= -IGRAND;
ietlim[1]= IGRAND;
idebug[0]=idebug[1]=idebug[2]=idebug[3]=0;
NStId = NSelStId = FgSelSt = 0;
StId1 = StId2 = 0;
nkillms = 0; for(i=0;i<NKILLMSMX;i++) killms[i]=0;
ffprt=0;

starcut.Mag_B_cut  = 0;   starcut.Mag_B_min  = starcut.Mag_B_max = 0.;
starcut.Mag_V_cut  = 0;   starcut.Mag_V_min  = starcut.Mag_V_max = 0.;
starcut.Mag_R_cut  = 0;   starcut.Mag_R_min  = starcut.Mag_R_max = 0.;
starcut.Mag_I_cut  = 0;   starcut.Mag_I_min  = starcut.Mag_I_max = 0.;
starcut.Mag_BV_cut = 0;   starcut.Mag_BV_min = starcut.Mag_BV_max = 0.;
starcut.Mag_BR_cut = 0;   starcut.Mag_BR_min = starcut.Mag_BR_max = 0.;
starcut.Mag_VR_cut = 0;   starcut.Mag_VR_min = starcut.Mag_VR_max = 0.;
starcut.Mag_RI_cut = 0;   starcut.Mag_RI_min = starcut.Mag_RI_max = 0.;
starcut.Rstar_cut  = 0;   starcut.Rstar_min  = starcut.Rstar_max = 0.;
starcut.Rtype_cut  = 0;   starcut.Rtype_min  = starcut.Rtype_max = 0;

mc.montecar=0;
mc.Taille_Finie=0;
mc.Blending=0;
mc.NGenStar=1;
mc.Usu0_Cut = 0.001;
mc.U0SimI=0.; mc.U0MinI=0.; mc.U0MaxI=0.;
mc.TireU0 = 0;
mc.T0SimI= -1.; mc.T0MinI=1.; mc.T0MaxI= -1.;
mc.TauSimI=0.; mc.TauMinI=0.01; mc.TauMaxI=365.;
mc.TireTau = 0;
mc.USimI=0.; mc.UMinI=0.; mc.UMaxI=0.;
mc.TireU = 0;
mc.Get_Par_Blending = 0;
mc.coeff_Arec1[0] = mc.coeff_Arec1[1] = 0.;
mc.coeff_Arec2[0] = mc.coeff_Arec2[1] = 0.;
mc.starcal1.Rstar = mc.starcal2.Rstar = 0.;
mc.iseed1 = 36117;
mc.iseed2 = 51106;
mc.iseed3 = 21478;

}

/*==========================================================================*/
void OPENSUIV(int ic)
/* ic = couleur (0 a 1) */
{
int_4 i,j,rc;
int_4 l, i1, i2;
int lp;

if(ffprt>0) printf("===> OPENSUIV couleur %3d\n",ic+1);

if( nbfil[ic] <= 0 ) return;

if(nmes[ic]>0 && ffprt>0) {
  printf("Coul %2d nmes %5d on lit %4d fichiers, nbmsMAX %5d\n"
        ,ic+1,nmes[ic],nbfil[ic],nbmsMX);
  if(ffprt>1) {
    printf("nbms:\n");
    for(i=0;i<nbfil[ic];i++) {printf(" %6d",nbms[ic][i]); if(i%10==9) printf("\n");}
    if((nbfil[ic]-1)%10!=9) printf("\n");

    printf("nms1fil:\n");
    for(i=0;i<nbfil[ic];i++) {printf(" %6d",nms1fil[ic][i]); if(i%10==9) printf("\n");}
    if((nbfil[ic]-1)%10!=9) printf("\n");
} }

for(i=0;i<nbfil[ic];i++) {

  l = SuiviGetType(lunfil[ic][i]);

  /* type de fichier de suivi */
  if( i==0 && ic==0 ) {
    typsuivi=l;
    printf("Type du fichier de Suivi: %d\n",typsuivi);
  } else if ( typsuivi != l ) {
    printf("couleur %3d fichier de suivi %5d mauvais type %d %d\n"
          ,ic+1,i+1,l,typsuivi);
    exit(-1);
  }

  /* lecture de global info */
  if( ( rc=SuiviReadGlobInfo( lunfil[ic][i],(char *) (&glinf) ) ) != 0 ) {
    printf("Problemes lecture global info coul %3d fichier %d, rc=%d\n",ic+1,i+1,rc);
    exit(-1);
  }
  DECODE_GLOBINFO(&glinf);

  //if( ffprt>3 || (ffprt>0 && i==0) ) PrtGlobInfo(&glinf,ffprt-3);
  if( ffprt>3 || (ffprt>0 && i==0) ) PrtGlobInfo(&glinf,99);
    
  /* lecture du nombre d'etoiles */
  l= SuiviGetNbStars(lunfil[ic][i]);
  if( i==0 ) {
    net[ic] = SuiviGetNbStars(lunfil[ic][i]);
  } else if( l !=net[ic] ) {
    printf("Probleme dans nombre d'etoiles coul %3d fichier %d: %d %d\n"
          ,ic+1,i+1,l,net[ic]);
    exit(-1);
  }

}

/* REMARQUE
   lecture de timeinfo:
   attention: on passe par un buffer intermediaire car quand on definit
   une structure, la longueur est un nombre entier de fois l'element
   le plus long de la structure (ici (double)) et donc on a, pour TIMEINFO
   sizeof(TIMEINFO)=416 au lieu de 412 effectivement ecrit 
   Depuis la version avec gestion du Byte swap entre machine il n'est
   plus possible de lire directement une serie de timeinfo (meme pb pour mesure),
   il faut boucler sur toutes les photos! (debug 13/03/98 Reza)
   Voici le commentaire exact de Reza:  
     "On ne peut plus lire un groupe de mesure ou de timeinfo s'il faut
      ByteSwaper et si Taille RecSuivi <> sizeof(structure)
      il faut lire les mesures/timeinfo un a un"
   Versions CVS ou il y avait encore la lecture par buffer:
     nbsread.h 1.14 , nbsread.c 1.30
*/

lp = ffprt-2;

if(nmes[ic]>0) {
  for(i=0;i<nbfil[ic];i++) {
    i2=nbms[ic][i];
    i1=nms1fil[ic][i];
    for(j=0;j<i2;j++) {
      if( (rc=SuiviReadTimeInfo(lunfil[ic][i], j+1, j+1, (char *)(&timinf))) != 0 ) {
        printf("Problemes lecture timinf info coul %3d fichier %5d, rc=%d\n"
              ,ic+1,i+1,rc);
        exit(-1);
      }
      l = j+i1;
      DECODE_TIMEINFO(&timinf,ic,l);
      indexu[ic][l] = 0;
      date  [ic][l] = (double) timeu[ic][l].TStart / JourSec;
      ampli [ic][l] = 1.;
      meslun[ic][l] = i;
      if( ffprt>1 ) {
        PrtTimeInfo (&timinf,l,ffprt-2);
        if( ffprt>5 ) PrtTimeInfoU (&timeu[ic][l],l,1);
      }
    }
  }
}

/* some prints ? */
if(nmes[ic]>0 && ffprt>0) {
  printf("date:\n");
  for(i=0;i<nmes[ic];i++) {printf(" %9.2f",date[ic][i]); if(i%10==9) printf("\n");}
  if((nmes[ic]-1)%10!=9) printf("\n");
  if(lp>=3) {
    printf("meslun:\n");
    for(i=0;i<nmes[ic];i++) {printf(" %9d",meslun[ic][i]); if(i%10==9) printf("\n");}
    if((nmes[ic]-1)%10!=9) printf("\n");
} }

}

/*==========================================================================*/
void READST(int ic)
/* ic = couleur (0 a 1) */
{
int_4 i,rc;

staru[ic].XRef=0;

if(nmes[ic]>0)
  for(i=0;i<nmes[ic];i++) mesu[ic][i].Flux=mesu[ic][i].Fond=0.;

if(debug>0) printf("READST: coul %3d et %d\n",ic+1,iet[ic]);

if( iet[ic]<=0 || iet[ic]>net[ic] ) return;

/*lecture de star info sur premier fichier */
if((rc=SuiviReadStarInfo(lunfil[ic][0],iet[ic],(char *)(&star))) != 0) {
  printf("Problemes lecture star info coul %3d etoile %d, rc=%d\n"
        ,ic+1,iet[ic],rc);
  exit(-1);
}

DECODE_STARINFO(&star,ic);

if(debug>1) {
  printf("starinfo: coul %3d nbfil=%5d\n",ic+1,nbfil[ic]);
  PrtStarInfo(&staru[ic],ic,5);
}

}

/*==========================================================================*/
void READET(int ic)
/* ic = couleur (0 a 1) */
{
int_4 i,j,rc;
int_4 i1,i2;

if( iet[ic]<=0 || iet[ic]>net[ic] ) return;
if( nmes[ic]<=0 ) return;

/* lecture des mesures: cf lecture expliquee pour timeinfo */
for(i=0;i<nbfil[ic];i++) {
  if(debug>1) printf("    Fichier %5d nbms %5d\n",i,nbms[ic][i]);
  i2=nbms[ic][i];
  i1=nms1fil[ic][i];
  for(j=0;j<i2;j++) {
    if((rc=SuiviReadMesures(lunfil[ic][i],iet[ic],j+1,j+1,(char *)(&mes))) != 0) {
      printf("Problemes lecture mes coul %3d etoile %d lim %5d %5d, rc=%d\n"
            ,ic+1,iet[ic],i1,i2,rc);
      exit(-1);
    }
    DECODE_MESURE(&mes,ic,j+i1);
    if(debug>5) PrtMesure(&mes,j+i1,5);
  }
}


if(debug>3) {
  if(debug>4) for(i=0;i<nmes[ic];i++) PrtMesureU(&mesu[ic][i],i,5);
  printf("\nmesu.Flux couleur=%d:\n",ic);
  for(i=0;i<nmesure[ic];i++) {
    j=indexu[ic][i];
    printf(" %12.5g",mesu[ic][j].Flux);
    if(i%10==9 || i==nmesure[ic]-1) printf("\n");
  }
  if(debug>4) {
    printf("mesu.ErrFlux couleur=%d:\n",ic);
    for(i=0;i<nmesure[ic];i++) {
       j=indexu[ic][i]; printf(" %12.5g",mesu[ic][j].ErrFlux);
       if(i%10==9 || i==nmesure[ic]-1) printf("\n");
     }
    printf("mesu.FluxB couleur=%d:\n",ic);
    for(i=0;i<nmesure[ic];i++) {
       j=indexu[ic][i]; printf(" %12.5g",mesu[ic][j].FluxB);
       if(i%10==9 || i==nmesure[ic]-1) printf("\n");
    }
  }
}

}

/*==========================================================================*/
void DATCLEAN(int ic)
/* ic = couleur (0 a 1) */
{
int_4 i,j,k;
double d1,d2;

if(ffprt>0) printf("=====> DATCLEAN coul=%3d ccd %5d\n",ic+1,ccdnum);

if(nmes[ic]<=0) return;

/* on tue selon les criteres utilisateur */
UDATCLEAN(ic+1);
/* a ce niveau les dates=GRAND2_KILL ne doivent absolument plus servir */

/* on tue les mesures 'n' pour chaque fichier de suivi */
if( nbfil[ic] > 0 ) {
  for(i=0;i<nbfil[ic];i++) {
    if( nbms[ic][i]<= 0 ) continue;
    for(j=0;j<nbms[ic][i];j++) {
      if( date[ic][nms1fil[ic][i]+j] >= GRAND2_KILL) continue;
      for(k=0;k<nkillms;k++) if ( j+1 == killms[k] ) {
        if(ffprt>1)
          printf(" Kill: mesure %6d (seq=%6d) du fichier suivi %4d\n"
                         ,j,nms1fil[ic][i]+j,i);
        date[ic][nms1fil[ic][i]+j] = GRAND2;
      }
    }
  }
}


/* on tue les temps en double et on prend le dernier */
if( nmes[ic]>1 ) {
  for(i=0;i<nmes[ic]-1;i++) {
    d1=date[ic][i];
    if( d1 >= GRAND2 ) continue;
    for(j=i+1;j<nmes[ic];j++) {
      d2=date[ic][j];
      if( d2 >= GRAND2 ) continue;
      if( fabs(d1-d2) < (double) 0.00000001 ) {
        if(ffprt>1)
          printf(" Kill: %5d t=%f car dupliquee avec %5d t=%f soit %f\n"
                ,i,d1,j,d2,fabs(d1-d2));
        date[ic][i] = GRAND2;
        break;
      }
    }
  }
}

}

/*==========================================================================*/
void DATSORT(int ic)
/* ic = couleur (0 a 1) */
{
int_4 i,j;

if(ffprt>0) printf("DATSORT: coul %3d ccd %5d\n",ic+1,ccdnum);

nmesure[ic]=0;
Tfirst[ic] = Tlast[ic] = -1.;

if(nmes[ic]<=0) return;

tri_double (date[ic],indexu[ic],nmes[ic]);

for(i=0;i<nmes[ic];i++) if(date[ic][i]<GRAND2) nmesure[ic]++;

if(nmesure[ic]>0) {
  Tfirst[ic]=date[ic][indexu[ic][0]];
  Tlast[ic]=date[ic][indexu[ic][nmesure[ic]-1]];
}

if(ffprt>0) {
  printf("coul %3d nmesure %6d nmes %6d Tint %.3f %.3f\n"
                  ,ic+1,nmesure[ic],nmes[ic],Tfirst[ic],Tlast[ic]);
  if(ffprt>1) {
    printf("date:\n");
    for(i=0;i<nmes[ic];i++)
      {
      j=indexu[ic][i];
      if(date[ic][j]<GRAND2) printf(" %9.2f",date[ic][j]);
        else printf(" %9.2f",-1.);
      if(i%10==9) printf("\n");
      }
    if((nmes[ic]-1)%10!=9) printf("\n");
  
    printf("timeu.NumPhoto:\n");
    for(i=0;i<nmes[ic];i++) {
       j=indexu[ic][i]; printf(" %9d",timeu[ic][j].NumPhoto); if(i%10==9) printf("\n");}
    if((nmes[ic]-1)%10!=9) printf("\n");

    if(ffprt>2) {
      printf("indexu:\n");
      for(i=0;i<nmes[ic];i++) {
         printf(" %9d",indexu[ic][i]); if(i%10==9) printf("\n");}
      if((nmes[ic]-1)%10!=9) printf("\n");
} } }

}

/*==========================================================================*/
int DECODESUIVINAME(char* str)
/*
Pour decoder les noms des variables d'environnement dans une chaine de characteres
- La chaine "str" est mofifiee (et doit donc etre modifiable)
- Variable d'environnement = /home/.../toto/ ou /home/.../toto (avec ou sans / final)
- La chaine "str" peut etre:
  $AAA ${AAA} $(AAA) $AAA/blabla ${AAA}/blabla $(AAA)/blabla /.../$AAA/... etc...
  Le cas ${AAA}blabla $(AAA)blabla n'est pas traite (/ impose en fin de nom de variable)
- Retourne: -1 si echec, 0 si pas de variable a decoder, sinon nombre de variables decodees
*/
{
char *str1,*vnam,*cenv;
int ls=0,lce,rc=0,i,ip,jp;

ls = strlen(str);
/* Pas de $ ou chaine trop courte */
if(posc(str,'$')<0 || ls<1) return rc;
str1 = (char *) malloc((NBSREAD_LSTR+5)*sizeof(char));
vnam = (char *) malloc(NBSREAD_LSTR*sizeof(char));

i= ip = rc = 0;
while( i<ls ) {
  if(str[i] == ' ') {               /* on vire les blancs */
    i++;
  } else if(str[i] == '/') {        /* on vire les doubles "/" */
    if(str1[ip]=='/') i++; else {str1[ip] = str[i]; ip++; i++;}
  } else if(str[i] != '$') {        /* pas de variable d'environnement */
    str1[ip] = str[i]; ip++; i++;
  } else {                          /* variable d'environnement */
    jp = 0;
    while( i<ls ) {
      if(str[i]=='/') {i++; break;}
      else if(str[i]=='$'||str[i]=='{'||str[i]=='}'||str[i]=='('||str[i]==')')  i++;
      else {vnam[jp]=str[i]; jp++; i++;}
    }
    vnam[jp] = '\0';
    if( (cenv=getenv(vnam)) == NULL )
      {printf("variable inconnue=%s.\n",vnam); free(str1); free(vnam); return -1;}
    lce = strlen(cenv);
    rc++;
    if(lce>0) {
      jp = 0;
      while( jp<lce ) {
        /* on evite TOTO=/home/toto/ et blabla/$TOTO -> blabla//home/toto */
        if(ip>0) if(jp==0 && cenv[jp]=='/' && str1[ip-1]=='/') {jp++; continue;}
        str1[ip] = cenv[jp]; ip++; jp++;
      }
      /* on evite TOTO=/home/toto/ et $TOTO/blabla -> /home/toto//blabla */
      if(cenv[lce-1]!='/') {str1[ip] = '/'; ip++;}
    }
  }
}
str1[ip] = '\0';
strcpy(str,str1);

free(str1); free(vnam);
return rc;
}

/*==========================================================================*/
void DECODE_GLOBINFO(GLOBINFO *p)
{
globu = *p;
}

/*==========================================================================*/
void DECODE_TIMEINFO(TIMEINFO *p, int ic, int imes)
/* ic = couleur (0 a 1), imes = mesure (0 a nmes-1) */
{
DecodeTim( p , &timeu[ic][imes]);
}

/*==========================================================================*/
void DECODE_STARINFO(STARINFO *p, int ic)
/* ic = couleur (0 a 1) */
{
staru[ic] = *p;
}

/*==========================================================================*/
void DECODE_MESURE(MESURE *p, int ic, int imes)
/* ic = couleur (0 a 1), imes = mesure (0 a nmes-1) */
{
if(date[ic][imes]<GRAND2_KILL) {
  DecodeMes( p , &mesu[ic][imes] );
  Calibre_F_E( &mesu[ic][imes], &timeu[ic][imes]);
} else {
  Mes_a_zero(&mesu[ic][imes]);
}
}

/*==========================================================================*/
void DO_MALLOC(void)
{
int n,ic;
size_t ll,sof, sofM, sofT;

if(ffprt>1)
  printf("DO_MALLOC: nbmsMX=%d  szMXBYTE=%d\n",nbmsMX,szMXBYTE);

for(ic=0;ic<NCOULMX;ic++) {          /* BOUCLE SUR LES COULEURS */

  n = (nmes[ic] > 0) ? nmes[ic] : 1;

  sofT = sizeof(TIMEINFOU);
  ll= n * sofT;
  if(ffprt>0)
    printf("timeu: allocation de %5d unites de %d bytes soit %d\n"
          ,n,(int)sofT,(int)ll);
  if( ( timeu[ic] = (TIMEINFOU *) malloc(ll) ) == 0 ) {
    printf("impossible d'allouer *timeu par malloc coul %3d\n",ic+1);
    exit(-1);
  }

  sofM = sizeof(MESUREU);
  ll= n * sofM;
  if(ffprt>0)
    printf("mesu: allocation de %5d unites de %d bytes soit %d\n"
          ,n,(int)sofM,(int)ll);
  if( (mesu[ic] = (MESUREU *) malloc(ll) ) == 0 ) {
    printf("impossible d'allouer *mesu par malloc coul %3d\n",ic+1);
    exit(-1);
  }

  sof = sizeof(int_4);
  ll= n * sof;
  if(ffprt>0)
    printf("indexu: allocation de %5d unites de %d bytes soit %d\n"
          ,n,(int)sof,(int)ll);
  if( (indexu[ic] = (int_4 *) malloc(ll) ) == 0 ) {
    printf("impossible d'allouer *indexu par malloc coul %3d\n",ic+1);
    exit(-1);
  }

  sof = sizeof(unsigned short);
  ll= n * sof;
  if(ffprt>0)
    printf("meslun: allocation de %5d unites de %d bytes soit %d\n"
          ,n,(int)sof,(int)ll);
  if( (meslun[ic] = (unsigned short *) malloc(ll) ) == 0 ) {
    printf("impossible d'allouer *meslun par malloc coul %3d\n",ic+1);
    exit(-1);
  }

  sof = sizeof(double);
  ll= n * sof;
  if(ffprt>0)
    printf("date: allocation de %5d unites de %d bytes soit %d\n"
          ,n,(int)sof,(int)ll);
  if( (date[ic] = (double *) malloc(ll) ) == 0 ) {
    printf("impossible d'allouer *date par malloc coul %3d\n",ic+1);
    exit(-1);
  }

  sof = sizeof(float);
  ll= n * sof;
  if(ffprt>0)
    printf("ampli: allocation de %5d unites de %d bytes soit %d\n"
          ,n,(int)sof,(int)ll);
  if( (ampli[ic] = (float *) malloc(ll) ) == 0 ) {
    printf("impossible d'allouer *ampli par malloc coul %3d\n",ic+1);
    exit(-1);
  }

}             /* BOUCLE SUR LES COULEURS */

}

/*==========================================================================*/
void DO_FREE(void)
{
int ic;
for(ic=0;ic<NCOULMX;ic++) {
  free(timeu[ic]);
  free(mesu[ic]);
  free(indexu[ic]);
  free(date[ic]);
  free(ampli[ic]);
}
}

/*==========================================================================*/
int read_timeinfo(int ic,int imes,TIMEINFO *tim)
/*
Lecture de timeinfo pour la couleur ic (0->1), mesure imes (0->nmes[ic]-1)
  return 0 si OK, <0 si problemes
*/
{
int rc,ifil,ime;

if ( imes<0 || imes>=nmes[ic] ) return(-1);
ifil = meslun[ic][imes];
if ( ifil<0 || ifil>=nbfil[ic] ) return(-2);
ime = imes-nms1fil[ic][ifil];
if ( ime<0 || ime>=nbms[ic][ifil] ) return(-3);

rc=SuiviReadTimeInfo(lunfil[ic][ifil], ime+1, ime+1, (char *)(tim));
if( rc != 0 ) {
  printf("read_timeinfo: problemes lecture timinf info coul %3d fichier %5d, rc=%d\n"
        ,ic+1,imes+1,rc);
  return(-4);
}

/*
printf("ic=%2d imes=%5d ifil=%4d ime=%5d NumPhoto=%d TStart=%d Expose=%d\n"
      ,ic,imes,ifil,ime,tim->NumPhoto,tim->TStart,tim->Expose);
printf("          Calib=%f %f sig=%f %f %f\n"
      ,tim->Calib[0][0],tim->Calib[0][1],tim->SigX,tim->SigY,tim->Rho);
*/

return(0);
}

/*==========================================================================*/
void SELSTAR(void)
/* selection de l etoile courante ? */
{
int i;
FgSelSt = 0;
if( NSelStId <= 0 || NStId <= 0 || StId1 == 0 || StId2 == 0 ) return;
for(i=0;i<NSelStId;i++)
  if( StId1[i]<=iet[0] && iet[0]<=StId2[i] ) {FgSelSt = iet[0]; break;}
}

/*==========================================================================*/
void GET_STAR_ID(void)
/*
Lecture du fichier dataccd pour identifier les cartes de selection des etoiles.
*/
{
char str[NBSREAD_LSTR];
int nallmin = 50, l;

NStId = NSelStId = FgSelSt = 0;
if( StId1 != 0 ) { free(StId1); StId1 = 0;}
if( StId2 != 0 ) { free(StId2); StId2 = 0;}

StId1 = (int *)malloc(nallmin*sizeof(int));
StId2 = (int *)malloc(nallmin*sizeof(int));
if( StId1==0 || StId2==0 ) {
  printf("GET_STAR_ID: Impossible d allouer StId1/2 pour la premiere fois %d\n"
        ,nallmin);
  return;
} else NStId = nallmin;

rewind(dataccd);
while ( fgets(str,NBSREAD_LSTR,dataccd) != 0 ) {
  strip(str,'B',' ');
  if ( strstr(str,"@SELET") != str ) continue;
  if( NSelStId >= NStId ) {
    l = NStId + nallmin;
    StId1 = (int *)realloc(StId1,l*sizeof(int));
    StId2 = (int *)realloc(StId2,l*sizeof(int));
    if( StId1==0 || StId2==0 ) {
      printf("GET_STAR_ID: Impossible de re-allouer %d pour StId1/2\n",l);
      return;
    } else {
      NStId = l;
      /* printf("GET_STAR_ID: reallocation de %d pour StId1/2 reussie\n",NStId); */
    }
  }
  StId1[NSelStId]=StId2[NSelStId]=-1;
  sscanf(&str[6],"%d %d",&StId1[NSelStId],&StId2[NSelStId]);
  if( StId2[NSelStId]<StId1[NSelStId] ) StId2[NSelStId]=StId1[NSelStId];
  NSelStId++;
}

if(ffprt>0) {
  printf("GET_STAR_ID: %d selections d etoiles (lim=%d)\n",NSelStId,NStId);
  if(ffprt>5)
    for(l=0;l<NSelStId;l++) {
      printf("    sel[%d]   et %d",l+1,StId1[l]);
      if( StId1[l]!=StId2[l] ) printf("  a %d\n",StId2[l]);
        else printf("\n");
    }
}

}

/*==========================================================================*/
void Nbsread_SigCatchInit(void)
{
int rc;
#ifdef OSF1
Nbsread_Action.sa_mask = 0;
Nbsread_Action.sa_flags = SA_RESTART;
#endif
#ifdef HPUX
memset( &(Nbsread_Action.sa_mask), 0, sizeof(sigset_t) ); 
Nbsread_Action.sa_flags = 0;
#endif
#ifdef Linux
memset( &(Nbsread_Action.sa_mask), 0, sizeof(sigset_t) ); 
Nbsread_Action.sa_flags = 0;
#endif
Nbsread_Action.sa_handler = Nbsread_SigCatch;

if( (rc=sigaction(SIGINT, &Nbsread_Action, NULL)) != 0 )
  printf("Nbsread_SigCatchInit_Error: signal SIGINT=%d rc=%d\n",SIGINT,rc);
if( (rc=sigaction(SIGUSR1, &Nbsread_Action, NULL)) != 0 )
  printf("Nbsread_SigCatchInit_Error: signal SIGUSR1=%d rc=%d\n",SIGUSR1,rc);
if( (rc=sigaction(SIGUSR2, &Nbsread_Action, NULL)) != 0 )
  printf("Nbsread_SigCatchInit_Error: signal SIGUSR2=%d rc=%d\n",SIGUSR2,rc);
}

/*==========================================================================*/
void Nbsread_SigCatch(int s)
{
switch(s) {
  case SIGINT :
    printf("\n\nNbsread_SigCatch: signal SIGINT=%d catched\n\n",s);
    break;
  case SIGUSR1 :
    printf("n\nNbsread_SigCatch: signal SIGUSR1=%d catched\n\n",s);
    break;
  case SIGUSR2 :
    printf("n\nNbsread_SigCatch: signal SIGUSR2=%d catched\n\n",s);
    break;
  default :
    printf("n\nNbsread_SigCatch: signal 1=%d catched et non traite\n\n",s);
    break;
}
printf("UENDCCD is called\n");
UENDCCD();
printf("UEND is called\n");
UEND();
printf("exit is called\n");
exit(-1);
}


/*************************************************************************/
/*************************************************************************/
/*******************  Commentaire nbsread et nbgene **********************/
/*************************************************************************/
/*************************************************************************/
/* 
++ 
  Module        nbsread (C et C++)
  Lib   libnbsread++.a et libnbsread.a
  include       nbsread.h nbgene.h
|
        *- nbsread*
        
        C'est un programme qui permet de lire des fichiers de suivi
        pour 2 couleurs associees. Il fournit des points d'entrees permettant
        a l'utilisateur de piloter une analyse sans devoir s'occuper de
        la machinerie relative a la lecture de fichiers de suivis.
|   nbsread  nbsread.dataccd.Exemple
        Les noms des fichiers de suivi a lire doivent etre mis dans un fichier
        dont un descriptif est donne dans le fichier
        `AnaSuiv/nbsread.dataccd.Exemple'. La lecture de plusieurs fichiers
        de suivi pour plusieurs CCD dans 2 couleurs est permise.
        Ce fichier permet aussi de selectionner des etoiles par leur
        numero ou un intervalle de numeros, le programme renvoie `FgSelSt'
        non nulle si l'etoile a ete selectionnee par les datacards.
|
        *- Structure des donnees*
        
        Le programme lit les etoiles les unes apres les autres.
        Pour chaque etoile il remplit diverses structure ou tableau
        de structures permettant l'analyse de sa courbe de lumiere
        (la description des structures se trouve dans `fsvst.h'):
|  EXTERN GLOBINFO globu : decodage de GLOBINFO
|  EXTERN STARINFO staru[NCOULMX] : decodage de STARINFO
|  EXTERN int_4 nmes[NCOULMX] : nombre de mesures pour le CCD courant
|                               et les couleurs 0 et 1
|  EXTERN TIMEINFOU *timeu[NCOULMX] : tableau de structures TIMEINFOU
|  EXTERN MESUREU *mesu[NCOULMX] : tableau de structures MESUREU
|  EXTERN double *date[NCOULMX] : tableau des dates en jours
|  EXTERN float *ampli[NCOULMX] : tableau des amplifications montecarlo
|  EXTERN int_4 nmesure[NCOULMX] : nombre de mesure valides
|                                  pour le CCD courant
|  EXTERN int_4 *indexu[NCOULMX] : tableau des images valides
|                                  rangees chronologiquement
|  EXTERN STARCALU starcal : strucure de magnitude absolue
|  EXTERN STARCUTU starcut : structure de coupure sur les magnitudes
|                            et les couleurs (permet de ne pas lire
|                            certaines etoiles pour gagner du temps)

        Les tableaux de structures ont `nmes[0] / nmes[1]' elements correspondant
        au nombre total d'images pour les couleurs 0 et 1 dans les fichiers
        de suivi d'un CCD.
|
        *- Les points d'entree utilisateurs et structure de nbsread*
        
        Ce sont les routine commencant par la lettre `U...()'
        ayant un nom en majuscule.

|  Initialisation du job {
|
|    void UINIT(void);  <*********** USER
|    Initialisation montecarlo
|
|    Boucle sur les CCD a lire {
|
|      ---> Lecture de GLOBINFO et de TIMEINFOU pour un CCD
|      void UDATCLEAN(int coul);  <*********** USER
|      Classement des images valides par dates croissantes
|      void UINITCCD(void);  <*********** USER 
|      Initialisation montecarlo pour le CCD courant
|
|      Boucle sur les etoiles a lire {
|
|        ---> lecture de MESUREU
|        Calcul des magnitudes absolues
|        Coupure sur les magnitudes absolues et la couleur des etoiles
|        Generation montecarlo pour chaque etoile
|          (eventuellement plusieurs generation sur la meme etoile)
|        void UEVT(void);   <*********** USER
|
|      } Fin de la boucle sur les etoiles
|
|      void UENDCCD(void);  <*********** USER
|
|    } Fin de la boucle sur les CCD
|
|    void UEND(void);  <*********** USER
|    Fin du montecarlo, impression des statistiques des generations
|
|  } Fin du job
|
        *- void UINIT(void)*
        
        Tout ce qui ne doit etre fait qu'une fois par job
        definitions de parametres generaux, ouverture hbook etc...
|
        *- UDATCLEAN(int icoul)*
        
        Permet de tuer des images
        Si date[icoul][imes] est mise a `GRAND2' l'image correspondante
        ne sera pas prise en compte dans la sequence en temps.
        Si date[icoul][imes] est mise a `GRAND2_KILL', l'image correspondante
        ne sera ni decodee ni prise en compte (structure `MESUREU' mis a zero).
        Ceci permet d'eliminer des images ou il y a eu des erreurs de codage
        qui conduisent a un core dump.
|
        *- Classement et elimination des images*
        
        Les images sont ensuite rangees chronologiquement en ne tenant pas
        compte des images tuees. Pour les images ayant la meme date, seules
        celles de plus grand indice sont prises en compte.
        Le tableau `indexu[icoul][]' est fournit: il contient les indices des images
        valides rangees chronologiquement et sans doublons. L'indice de indexu
        varie de `0 a nmesure[icoul]' valeur qui represente le nombre d'images
        utilisables pour l'analyse.
        Par exemple, l'indice de la 5ieme image valide de la couleur bleue (reference
        secondaire) est:
|   j = indexu[1][5];
        et l'acces au flux brut est donne par:
|   fb = mesu[1][j].FluxB;
        Toutes les images sont encore accessibles directement, par exemple
        le flux brut de la 7ieme image (la 7ieme des fichiers de suivi
        sans aucun choix ni classement chronologique) est, pour la reference
        primaire (couleur rouge):
|   fb = mesu[0][7].FluxB;
|
        *- UINITCCD(void)*
        
        C'est un point d'entree a la fin de l'initialisation
        du CCD courant. A ce niveau toutes les informations de `TIMEINFOU',
        les notions d'images valides et le classement chronologique sont
        connus. Les variables `TFIRST' et `TLAST' representent respectivement
        la valeur de `date[][]' la plus petit et la plus grande pour les images
        valides. A ce niveau peuvent etre faites des initialisation montecarlo
        propres a un CCD.
|
        *- Calcul des magnitudes absolues*
                
        Par defaut, il n'y a pas
        de calcul de magnitudes absolues. L'utilisateur peut fournir une
        routine de prototype:
|   void mycalib(float,float,STARCALU *);
        puis la connecter (dans `UINIT' ou `UINITCCD') au prototype:
|   void (*Calibration_Absolue)(float,float,STARCALU *, void*);
        en ecrivant: `Calibration_Absolue = mycalib;'.
        Cette routine sera appelee automatiquement pour chaque etoile
        avec les arguments:
|   Calibration_Absolue(staru[0].FluxRef,staru[1].FluxRef,&starcal,uparms);
        et remplira la structure `starcal'. Dand sa propre routine,
        l'utilisateur peut remplir des magnitudes `B,V,R,I' ainsi
        que le type de l'etoile (`SP=1 ou GR=2 ou indefini=0') et
        son rayon en unites de rayon solaire. Le remplissage de
        toutes ces valeurs n'est pas obligatoire, il depend de ce qui
        doit etre fait apres, soit par l'utilisateur dans son analyse,
        soit par la simulation montecarlo (l'effet de taille finie par
        exemple a besoin du rayon de l'etoile).
        Pour passer des parametres supplementaires a la fonction,
        il faut initialiser le pointeur par l'assignation:
|   Calibration_Absolue_UParms = &user_data_block
        et l'adresse sera passee comme 4eme argument de Calibration_Absolue.
        `user_data_block' doit etre declaree comme variable `globale'
        dans le fichier ou comme variable `static' dans la routine
        ou est faite l'assignation.
|
        *- Coupure sur les magnitudes absolues et la couleur des etoiles*
        
        Pour des raisons de rapidite, l'utilisateur peut couper sur la
        magnitude absolue ou la couleur de l'etoile en remplissant la
        structure `starcut'. Le temps de lecture de l'etoile est
        economise, et le point d'entree `UEVT()' n'est pas appele.
|
        *- Generation montecarlo pour chaque etoile*
        
        L'effet de microlentille est simule sur chaque etoile
        (ou plusieurs fois par etoile). Une description detaillee des
        possibilites de simulation est donnee plus loin.
|
        *- UEVT(void)*
        
        C'est l'endroit ou l'utilisateur met le code
        d'analyse et de traitement des courbes de lumieres. A ce niveau
        les etoiles selectionnees dans le fichier datacard par une carte
        `@SELET' sont identifier par le flag `FgSelSt>0'. Les variables
        de niveau d'impression, debug, etc... sont accessibles:
|  iet[NCOULMX] : numero de l'etoile (ref ou ref sec)
|  FgSelSt : flag de selection
|  idebug[4] : niveau de debug [iet1,iet2,ccd,niveau de debug]
|  ffprt : niveau d'impression
|  etc...
|
        *- void UENDCCD(void)*
                
        Point d'entree a la fin du traitement
        d'un CCD. Il permet de fermer ou d'imprimer toute structure
        relative a un CCD.
|       
        *- void UEND(void)*
                
        Point d'entree a la fin du job. C'est l'endroit
        pour desallouer les tableaux et pour imprimer les conclusions du job
|
        *- Arguments de l'executable*
        
|     [-h] :  help tres tres detaille
|     [-et et1,et2,etinc] :  numeros des etoiles a traiter
|     [-prt niv] :  niveau de print
|     [-dbg et1,et2,ccd,niv] : debug de certaine etoiles/ccd
|     [-kill mes1 mes2 mes3 ...] : pour tuer les images mes1...
        -kill permet de tuer les images mes1,mes2,mes3,... de `chaque'
        fichier de suivi pour un CCD et une couleur. C'est particulierement
        utile si, par exemple, la premiere image de chaque fichier de suivi
        est l'images de reference: on s'en debarasse en indiquant:
|  -kill 1
|
        Pour les arguments de nbsread, conformenent aux conventions prises
        pour les routines de lecture de fichiers de suivi,
        les numeros d'etoiles et d'images vont de `1 a nmes[]'
        (et non pas de `0 a nmes[]-1').
|
        *- Comment creer ses propres routines utilisateur*
        
        Il faut ecrire un fichier dans lequel sont definies
        les point d'entree utilisateurs:
        `UINIT(), UDATCLEAN(int), UINITCCD(), UEVT(), UENDCCD() et UEND()'.
        Le nom du fichier doit obligatoirement s'appeler
|  nbsreadu_myroutine.cc  (meme si on compile en C)
|  nbsreaduF_myroutine.f  pour la partie Fortran (si il y a des hbook)
        myroutine est un nom choisi librement par l'utilisateur.
        La compilation et le link s'effectue en executant:
|  make -f Makeus++ nbsread_myroutine   (si on travaille en C++)
|  make -f Makeus   nbsread_myroutine   (si on travaille en C  )
        Un exemple detaille a ete sauve dans la base CVS:
|  nbsreadu_squel.cc , nbsreaduF_squel.f
        Pour l'executer et regarder ses impressions,
        copiez `nbsreadu_squel.cc' et `nbsreaduF_squel.f' dans un de
        vos repertoires, et recuperez dans la base CVS `AnaSuiv/Makeus++'
        (ou `AnaSuiv/Makeus' si vous voulez compiler et linker en C).
        Compilez et linkez:
|  make -f Makeus++ nbsread_squel
|  (ou make -f Makeus nbsread_squel)
        Creez un fichier `fichier.dataccd' contenant des fichiers de suivi.
        Executez `nbsreadu_squel' en tapant la ligne de commande suivante:
|  nbsreadu_squel fichier.dataccd -prt 1 \
|          -mc 2 -u0 0.01,0.5 -t0 1,-1,-1 -tau 2,10,-1 -tf 2,0,2 -bl 1
        (cette procedure est egalement decrite aussi dans `nbsreadu_squel.cc')
|
        *- Initialisation du montecarlo et generation*
        
        L'initialisatuon peut etre effectuee, soit par arguments de
        l'executable:
|     [-mc ngenstar]
|     [-u0 U0MinI,U0MaxI,U0SimI]
|     [-t0 T0MinI,T0MaxI,T0SimI]
|     [-tau TauMinI,TauMaxI,TauSimI]
|     [-tf type(1ou2),UMinI,UMaxI,USimI,Usu0_Cut]
|     [-bl Blending (1ou2)]
|     [-seed seed1,seed2,seed3]
        soit en remplissant directement la structure `MONTECARLO mc' dans
        UINIT() et eventuellement UINITCCD() pour l'initialisation du `t0'.
        `nbsread -h' donne un help online tres detaille sur la facon de piloter
        la generation. La description des variables de la structure `mc'
        se trouve dans le fichier `nbgene.h'.
        La description des possibilites du montecarlo est donnee
        ci-apres en fonction des variables de la structure `mc':

| ----- CHOIX DU TYPE DE SIMULATION -----
| mc.Taille_Finie = 0 : generation taille ponctuelle
|                   1 :            taille finie tirage plat en U = Rs(proj)/Re
|                   2 :            taille finie tirage plat en UC (U = UC*Rs)
| mc.Blending = 1 : simulation du "blending" sur l'etoile principale (la + brillante)
|               2 : simulation du "blending" sur l'etoile cachee (la - brillante)
| 
| ----- GENERATION PLATE -----
| si max<=min
|   - TauSim=TauSimI Tau du phenomene
|   - T0Sim=T0SimI   instant du maximum si T0SimI>0
|   - U0Sim=U0SimI   parametre d'impact (en Re) si U0SimI>0
|   - USim=USimI     projection rayon etoile dans plan naine (en Re) si USimI>0
| si max>min
|   - TauSim est engendre avec une loi plate entre TauMinI ET TauMaxI
|   - T0Sim  est engendre avec une loi plate entre T0MinI  ET T0MaxI
|            si T0MaxI<T0MinI et T0SimI<0 on prend plat sur toute la periode
|   - U0Sim  est engendre avec une loi plate entre U0MinI  ET U0MaxI
|   - USim   est engendre avec une loi plate entre UMinI ET UMaxI si mc.Taille_Finie=1
|                                      entre UMinI*Rs ET UMaxI*Rs si mc.Taille_Finie=2
| 
| ----- GENERATION SELON UNE LOI DETERMINEE PAR L'UTILISATEUR -----
| l'utilisateur peut fournir la loi de tirage aleatoire pour u0,tau et u:
|   - mc.TireU0(float), mc.TireT0(float), mc.TireTau(float) et 
|        mc.TireU(float) sont 3 pointeurs sur des fonctions 
|        "float (*TireU0)(float alea)"
|        Il sont initialises a 0 dans nbsread.c
|   - mise en oeuvre (ex: pour le tirage de tau plat en log):
|     1-/ dans nbsreadu.c, l'utilisateur ecrit les routines de tirage
|         ( alea est uniformement distribue entre [0.,1.[ )
|       float my_tiretau(float alea)
|         {
|         float tau;
|         tau = log(mc.TauMinI)
|             + alea * ( log(mc.TauMaxI) - log(mc.TauMinI) );
|         return( (float) exp( (double) tau ) );
|         }
|      2-/ dans UINIT (par exemple)
|         mc.TireTau = my_tiretau;

        Evidemment dans l'hypothese d'un tirage avec effet de taille finie
        le rayon de l'etoile (starcal.Rstar) doit etre rempli
        dans la calibration absolue donnee par l'utilisateur.

        Pour une generation avec effet de blending, l'utilisateur doit
        une routine de simulation de l'effet de blending et la connecter
        au prototype mc.Get_Par_Blending dans UINIT() ou UINITCCD().
|  mc.Get_Par_Blending = my_blending_gene

        Pour une simulation avec effet de blending sur l'etoile la plus
        brillante du couple, il faut remplir
|  mc.coeff_Arec1[ic] pour les 2 couleurs (ic=0,1)
        Arec est le coefficient de blending:
|  Ampli(blend) = (Ampli-1.)*Arec +1

        Pour une simulation avec effet de blending sur l'etoile la moins
        brillante du couple, il faut remplir
|  mc.coeff_Arec2[ic] pour les 2 couleurs (ic=0,1)

        Pour une simulation avec effet de blending sur l'etoile la plus
        brillante du couple et avec effet de taille finie, il faut remplir
|  mc.coeff_Arec1[ic] pour les 2 couleurs (ic=0,1)
|  mc.starcal1.Rstar :   rayon de l'etoile la plus brillante

        Pour une simulation avec effet de blending sur l'etoile la moins
        brillante du couple et avec effet de taille finie, il faut remplir
|  mc.coeff_Arec2[ic] pour les 2 couleurs (ic=0,1)
|  mc.starcal2.Rstar :   rayon de l'etoile la moins brillante
|
        *- Initialisation des aleatoires*
        
        Le generateur utilise est "drand48" standard unix.
        Il est prevu 3 variables permettant d'initialiser le generateur:
|  mc.iseed1, mc.iseed2, mc.iseed3
        A la fin du job l'etat des graines est automatiquement imprime
        et permet de redemarrer un autre job pour la suite de la
        sequence d'aleatoires.
|
        *- Exemple de fichier datacard*
        
         Cf AnaSuiv/nbsread.dataccd.Exemple
|4                  nombre de CCD a traiter
|@CCD 03            numero du 1er CCD a traiter
|@COUL 1            numero de la couleur de reference
|/home/suivi/file1_coul1_ccd03.suivi
|/home/suivi/file2_coul1_ccd03.suivi
|@COUL 2            numero de la couleur secondaire (facultatif)
|/home/suivi/file1_coul2_ccd03.suivi
|/home/suivi/file2_coul2_ccd03.suivi
|/home/suivi/file3_coul2_ccd03.suivi
|@FINCCD            carte de fin de CCD
|@CCD 01            numero du 2sd CCD a traiter
|@COUL 1            numero de la couleur de reference
|/home/suivi/file1_coul1_ccd01.suivi
|/home/suivi/file2_coul1_ccd01.suivi
|@COUL 2            numero de la couleur secondaire (facultatif)
|/home/suivi/file1_coul2_ccd01.suivi
|@FINCCD            carte de fin de CCD
|@CCD 05            numero du 3ieme CCD a traiter
|@COUL 1            numero de la couleur de reference
|/home/suivi/file1_coul1_ccd05.suivi
|/home/suivi/file2_coul1_ccd05.suivi
|/home/suivi/file3_coul1_ccd05.suivi
|/home/suivi/file4_coul1_ccd05.suivi
|@COUL 2            numero de la couleur secondaire (facultatif)
|/home/suivi/file1_coul2_ccd05.suivi
|/home/suivi/file2_coul2_ccd05.suivi
|/home/suivi/file3_coul2_ccd05.suivi
|@FINCCD            carte de fin de CCD
|@CCD 07            numero du 4ieme CCD a traiter
|@COUL 1            numero de la couleur de reference
|/home/suivi/file1_coul1_ccd07.suivi
|/home/suivi/file2_coul1_ccd07.suivi
|/home/suivi/file3_coul1_ccd07.suivi
|@FINCCD            carte de fin de CCD
|@END               carte de fin de fichier datacard
|
|.... liste d etoiles a selectionner: FgSelSt != 0 dans nbsread.h
|@SELET 1001
|@SELET 1234
|@SELET 2032 2050
|@SELET 5044
|@SELET 10001 10100
|@SELET 6022

        A l'exception de la 1ere carte et des noms de fichiers de suivi
        toute carte valide commence par un caractere @.
        Il est fournit un moyen de selectionner des etoiles dans le
        fichier datacard et de les identifier dans `UEVT()'.
        grace au flag `FgSelSt'.
|  @SELET 1001        selection de l'etoile 1001 (attention et=[1,...])
|  @SELET 1001 2001   selection des etoiles de 1001 a 2001 inclusif
|
        *- Commentaire sur les structures des fichers de suivi et de nbsread*
        
        Dans le fichier de suivi sont codees les structures:
|     GLOBINFO STARINFO MESURE TIMEINFO
        La description des diverses structures et variables qu'elles
        contiennent se trouve dans `fsvst.h'.
        Pour gagner de la place, certaines d'entre elles (`TIMEINFO,MESURE')
        ont des variables codees de facon non triviale dans les mots
        de 16 bits. A la relecture du fichier de suivi, `nbsread'
        decode automatiquement les variables pour les ranger en clair
        dans des structures utilisateurs (dont le nom se termine par un `U').
        Notamment la structure `MESUREU' decode la structure `MESURE' et
        renvoit le flux et les erreurs calibrees et tenant compte des
        informations codees dans la structure `TIMEINFO'.
|     TIMEINFO  -> TIMEINFOU (qui ne contient qu'une partie des infos)
|     MESURE    -> MESUREU
        Les routines de decodage sont dans fsvst.h mais ne doivent
        en pratique pas etre utilisees directement par
        l'utilisateur de `nbsread':
|     TIMEINFOU * DecodeTim( TIMEINFO *time, TIMEINFOU *timeu);
|     MESUREU * DecodeMes( MESURE *mesc, MESUREU *mesu);
|
        *- Routines pratiques utilisables dans nbsread*
        
|  int read_timeinfo(int ic,int imes,TIMEINFO *tim)
        Comme `TIMEINFOU' ne contient qu'une partie des informations
        de la structure `TIMEINFO' du fichier de suivi, il peut etre pratique
        de pouvoir recuprer `TIMEINFO'. `read_timeinfo' permet de
        recuperer `TIMEINFO' pour la couleur `ic (0,1)' et la mesure
        `imes (0,nmes[ic]-1)'.
|    void PrtGlobInfo (GLOBINFO *glinf, int lp);
|    void PrtStarInfo (STARINFO *sti, int n, int lp);
|    void PrtMesure (MESURE *mes, int n, int lp);
|    void PrtMesureU (MESUREU *mesu, int n, int lp);
|    void PrtTimeInfo (TIMEINFO *tim, int n, int lp);
|    void PrtTimeInfoU (TIMEINFOU *tim, int n, int lp);
        Ces routines permettent d'imprimer les diverses structures
        d'un fichiers de suivi. Elles sont dans `fsvst.h'
|
        *- Routines de travail avec les transformations*
        
        Elles permettent, par exemple, de calculer les coordonnees
        sur l'image courante a partir des coordonnees de references de STARINFO.
        Les rputines se trouvent dans `fsvst.h' et `transfost.h'.
|    void GlobInfoToTransf(GLOBINFO *gli, TRANSFO *t1, TRANSFO *t2);
|    void TimeInfoToTransf(&timloc,&t1,&t2);
|    void CordTransf(double xS, double yS, double *xD, double *yD, 
|                    TRANSFO *transf)
        Voici un exemple d'utilisation a mettre dans `UEVT()' pour l'etoile
        courante.
|   TIMEINFO timloc;
|   TRANSFO t1,t2;
|   double Xcur[2],Ycur[2],x,y;
|   for(ic=0;ic<NCOULMX;ic++) {
|     for(k=0;k<nmes[ic];k++) {
|       read_timeinfo(ic,k,&timloc);
|       TimeInfoToTransf(&timloc,&t1,&t2);
|       -- des coordonnees de reference vers les courantes
|       x = staru[ic].XPos; y = staru[ic].YPos;
|       CordTransf(x,y,&Xcur[ic],&Ycur[ic],&t1);
|       PrintTransfo(&t1);
|       PrintTransfo(&t2);
|       -- des coordonnees courantes vers celles de reference
|       CordTransf(Xcur[ic],Ycur[ic],&x,&y,&t2);
|     }
|   }

--
*/