/*
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
*/
#include "machdefs.h"
#include "sopnamsp.h"  // using SOPHYA namespace
#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#include <string.h>

#include "fsvcache.h"
#include "fsvst.h"
#include "fmath.h"
#include "nbmath.h"
#include "nbrandom.h"
#include "nbsread.h"
#include "nbgene.h"

#ifdef __cplusplus
extern "C" {
#endif
/* Prototype des routines locales en C ou C++ */
void mycalib(float f1,float f2, STARCALU* stcal, void *user_parm);
void myblending(void);
float mytiretau(float);
#ifdef __cplusplus
}
#endif

/*================================================================*/
void UINIT(void)
{
static char CcdParm[20];

/* fonction utilisateur de calibration absolue */
Calibration_Absolue = mycalib;
strcpy(CcdParm,"CCD");
Calibration_Absolue_UParms = CcdParm;
/* fonction utilisateur de generation des parametres du blending */
mc.Get_Par_Blending = myblending;
/* fonction utilisateur de tirage du parametre Tau */
mc.TireTau = mytiretau;
}

/*=================== ============================================*/
void UDATCLEAN(int coul)
/*
  Attention ici coul = 1,2
    donc pour tuer une image pour la couleur coul:
    date[coul-1][...] = GRAND2;
*/
{
int ic = coul - 1;
printf("UDATCLEAN: couleur %d, indice tableau %d\n",coul,ic);
}

/*================================================================*/
void UINITCCD(void)
{
printf("UINITCCD: Debut de l'etude du CCD %d %d\n",numccd,ccdnum);
}

/*================================================================*/
void UEVT(void)
{
int ic,i,j;
float xnt[30];
double fball[2], fbgood[2], w;

if( ffprt > 0 || iet[0]%100==0 ) {
  printf("****> Et= %d %d   Sel=%d  fref=%g,%g\n"
        ,iet[0],iet[1],FgSelSt,staru[0].FluxRef,staru[1].FluxRef);
  printf("      B=%g V=%g R=%g I=%g Rs=%g t=%d\n"
        ,starcal.Mag_B,starcal.Mag_V,starcal.Mag_R,starcal.Mag_I
        ,starcal.Rstar,starcal.Rtype);
  printf("      u0=%g t0=%g tau=%g U=%g\n",mc.U0Sim,mc.T0Sim,mc.TauSim,mc.USim);
  printf("      Arec1=%g,%g Arec2=%g,%g Rstar1=%g Rstar2=%g\n"
        ,mc.coeff_Arec1[0],mc.coeff_Arec1[1],mc.coeff_Arec2[0],mc.coeff_Arec2[1]
        ,mc.starcal1.Rstar,mc.starcal2.Rstar);
  printf("      A0_pct=%g A0_tf=%g A0_blend=%g,%g\n"
        ,mc.A0MaxPt,mc.A0Max,mc.A0Max_blend[0],mc.A0Max_blend[1]);
  fflush(stdout);
}

/* moyenne des flux bruts de toutes les images et des images valides */
for(ic=0;ic<NCOULMX;ic++) {

  fball[ic] = w = 0.;
  if( nmes[ic] > 0 ) {
    for(i=0;i<nmes[ic];i++) {
      if( mesu[ic][i].Xi2 < 0. ) continue;
      if( mesu[ic][i].ErrFluxB <= 0. ) continue;
      fball[ic] += mesu[ic][i].FluxB*mesu[ic][i].ErrFluxB;
      w += mesu[ic][i].ErrFluxB;
    }
    if( w > 0. ) fball[ic] /= w;
  }

  fbgood[ic] = w = 0.;
  if( nmesure[ic] > 0 ) {
    for(j=0;j<nmesure[ic];j++) {
      i = indexu[ic][j];
      if( mesu[ic][i].Xi2 < 0. ) continue;
      if( mesu[ic][i].ErrFluxB <= 0. ) continue;
      fbgood[ic] += mesu[ic][i].FluxB*mesu[ic][i].ErrFluxB;
      w += mesu[ic][i].ErrFluxB;
    }
    if( w > 0. ) fbgood[ic] /= w;
  }

}

xnt[ 0] = iet[0];
xnt[ 1] = iet[1];
xnt[ 2] = staru[0].FluxRef;
xnt[ 3] = staru[1].FluxRef;
xnt[ 4] = starcal.Mag_B;
xnt[ 5] = starcal.Mag_V;
xnt[ 6] = starcal.Mag_R;
xnt[ 7] = starcal.Mag_I;
xnt[ 8] = starcal.Rstar;
xnt[ 9] = starcal.Rtype;
xnt[10] = mc.U0Sim;
xnt[11] = mc.T0Sim;
xnt[12] = mc.TauSim;
xnt[13] = mc.USim;
xnt[14] = mc.A0MaxPt;
xnt[15] = mc.A0Max;
xnt[16] = mc.A0Max_blend[0];
xnt[17] = mc.A0Max_blend[1];
xnt[18] = mc.starcal1.Rstar;
xnt[19] = mc.starcal2.Rstar;
xnt[20] = fball[0];
xnt[21] = fball[1];
xnt[22] = fbgood[0];
xnt[23] = fbgood[1];
// a stoquer dans votre fichier prefere

}

/*================================================================*/
void UENDCCD(void)
{
printf("UENDCCD: Fin de l'etude du CCD %d %d\n",numccd,ccdnum);
}

/*================================================================*/
void UEND(void)
{
}

/*================================================================*/
void mycalib(float f1,float f2, STARCALU* stcal, void *user_parm)
/* fonction de calibration absolue completement idiote */
{
if(f1<=0. || f2<0.) return;
if( strcmp((char *) user_parm,"CCD") ) return;
stcal->Mag_B = stcal->Mag_V = -2.5*log10((double) f2);
stcal->Mag_R = stcal->Mag_I = -2.5*log10((double) f1);
stcal->Rtype = (int) (drand01()+1.5);
stcal->Rstar = 3.*drand01();
}

/*================================================================*/
void myblending(void)
/* fonction de generation de blending completement idiote */
{
mc.coeff_Arec1[0] = 0.5+drand01()/2.;
mc.coeff_Arec1[1] = 0.5+drand01()/2.;
mc.starcal1.Rstar = 3.*drand01();

mc.coeff_Arec2[0] = 0.5-drand01()/2.;
mc.coeff_Arec2[1] = 0.5-drand01()/2.;
mc.starcal2.Rstar = 3.*drand01()/2.;
}

/*================================================================*/
float mytiretau(float alea)
/* fonction de tirage aleatoire de Tau personnelle */
/* (tirage logarithmique en Tau)                   */
{
if( mc.TauMinI<=0. || mc.TauMaxI<=0. || mc.TauMaxI<=mc.TauMinI ) {
  mc.TauMinI = 1.; mc.TauMaxI = 100.;
}
alea *= log(mc.TauMaxI/mc.TauMinI);
return( (float) mc.TauMinI*exp((double) alea) );
}