/* 
 * This file is part of the ESO UVES Pipeline
 * Copyright (C) 2004,2005 European Southern Observatory
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, 51 Franklin St, Fifth Floor, Boston, MA  02111-1307  USA
 */
/*
  ============================================================================
  flames_corvel:
  Purpose: 
  to cross correlate in velocity space a wavelength calibrated spectra with a 
  reference mask to get eventual velocity shift of one with respect to the 
  other.

  This code implements the Geneva alghorithm as for HARPS. Information and 
  reference alghorithms where provided from Claudio Melo, ESO-Paranal.
  ============================================================================
*/

/* 
  ----------------------------------------------------------------------------
  INCLUDES
  ----------------------------------------------------------------------------
*/ 

#ifdef HAVE_CONFIG_H
#  include <config.h>
#endif
//#include <flames_lfit.h>
#include <flames_midas_def.h>   /* MIDAS environment interface functions */
#include <flames_corvel.h>        /* FLAMES-UVES functions */
#include <flames_newmatrix.h>   /* FLAMES-UVES functions for array manipolation */
#include <uves_utils.h>         /* M_PI */
#include <uves_msg.h>        
#include <stdio.h> 
#include <math.h> 
#include <stdlib.h> 
#include <irplib_utils.h>
#include <string.h>

/* 
  ----------------------------------------------------------------------------
  LOCAL DEFINITIONS
  ----------------------------------------------------------------------------
*/ 
#define MAX_LEN 512
#define MAX_DIM 2
#define MAX_ORD 4
#define MAX_PIX 10000
#define MAX_DEG 4
#define FLAMES_SPEED_OF_LIGHT 299792.458 //is defined also in uves_utils.h

/*
static void 
fpoly(double x,double p[],int np);
*/

/*
static void 
get_mask(char*  tpl_name,double in_msk_wgt_min, double in_msk_hole_wid, 
         char* log_opt, double** msk_hole_width, double** msk_hole_center,
         double** msk_hole_wgt);
*/

static void 
do_cor_vel(double* wcal_sol,float** sp_flux,
           double* rv_ccf,double* msk_hole_siz, double* msk_hole_cen,
           double* msk_hole_wgt, double bar_v,double bar_v_max,
           int fit_type,int in_ima_nrow,
           int in_msk_nrow,int rv_ccf_size, double* ccf,double* ccf_max,
           double* pix_passed_ord,int* tot_line,double* ll_range_ord,
           int in_ima_id);

static void 
fit_ccf(double* rv_ccf,double* ccf_nor,int type,double* ccf_res,
        double* ccf_fit);





/* void fgauss(double x,double g[],int ng); */
/*
static double 
fgauss(double x,double a[],double y,double dyda[],int na);
*/

static void 
gaussian_fit(const double *   xfit, const double * yfit,int size,
                    double * norm, double * xcen, double * sig_x,
                    double * fwhm_x); 
static void 
correl_bin(int sp_flux_sz, float** sp_flux,double* sp_ll,double* sp_dll,
                  int *in_msk_nrow,double* msk_blu,double* msk_red,
                  double* msk_w, int* i_blue_masques,int* i_red_masques,
                  double* intensity_s,double* pix,double* ll_range);
    
static int 
hunt(double* xx, int n, double x, int jlo);


static void 
do_ccf_f(double* mask_ll,double* mask_d,double* mask_w,double* sp_ll,
         float** sp_flux,double* sp_dll,double* rv_ccf,double* ccf_o,
         double* pix_passed_ord,double* wcal_range_ord,int in_msk_nrow, 
         int in_ima_ncol, int rv_ccf_size, int in_ima_id);



/**
  @brief    find offset to be applied to wavecal solution to be appropriate
            to night calibrations
  @param    IN_A input image (merged spectrum)
  @param    IN_B input correlation velocity reference table
  @param    IN_N order number
  @param    OU_A output corvel table
  @param    OU_B output total image
  @param    OU_C output normalized image
  @param    rv_ccf_min radial velocity cross correlation function min
  @param    rv_ccf_max radial velocity cross correlation function max
  @param    rv_ccf_step radial velocity cross correlation function step
  @return   0 if successfull
 */


int flames_corvel(const char *IN_A,
                  const char *IN_B,
                  const int  IN_N,
                  const char *OU_A,
                  const char *OU_B,
                  const char *OU_C,
                  const double rv_ccf_min,
                  const double rv_ccf_max,
                  const double rv_ccf_step)
{



  char in_ima[MAX_LEN];   /* char array for input ima */
  char ou_ima[MAX_LEN];   /* char array for output ima */
  char ou_tab[MAX_LEN];   /* char array for output ima */
  char in_msk[MAX_LEN];   /* char array for input mask */

  /* MIDAS stuff */
  int midas_unit = 0;
  int midas_null = 0;
  int midas_nval = 0;
  int midas_status = 0;

  /* tmp variable used in MIDAS env calls */
  int in_ima_id =0;
  int ou_ima_id =0;
  int in_msk_id =0;
  int in_ima_naxis =0;
  int ou_ima_naxis =1;
  int ou_ima_npix[2] ={0,0};
  double ou_ima_start[2] ={0.,0.};
  double ou_ima_step[2] ={0.,0.};

  float cuts[4]={0.,0.,0.,0.};
  int tid=0;
  int ccf_pos_col=0;
  int ccf_nrm_col=0;
  int ccf_out_col=0;
  
  int in_ima_npix[MAX_DIM];
  int in_msk_ncol=0; 
  int in_msk_nrow=0;  
  char ident[73];
  char cunit[3][16];


  /* Other useful variables */
  int in_ima_nx = 0;      /* No of columns */
  int in_ima_ny = 0;      /* No of rows */
  int in_ima_ord=0;       /* order number of input image */
  
  float ** m_in_ima=NULL;     /* input image array */

  double* in_ima_wcal_sol=NULL;

  double* msk_hole_sta=NULL;
  double* msk_hole_end=NULL;
  double* msk_hole_cen=NULL;
  double* msk_hole_siz=NULL;
  double* msk_hole_wgt=NULL;

  double* msk_hole_cen_selw=NULL;
  double* msk_hole_siz_selw=NULL;
  double* msk_hole_wgt_selw=NULL;
  

  double in_msk_wgt_min=0.9;  /*1 */
  double in_msk_hole_wid=0.; /*0 */
  double tmp_double=0;
  double in_ima_wstart =0.;
  double in_ima_wstep  =0.;


  double* rv_ccf=NULL;
  int rv_ccf_size=0;
  //double rv_ccf_par[3] ={0.,0.,0.};

  int    wstart_id=0;
  int    wend_id=0;
  int    weight_id=0;

  int i=0;


  double  tmp_dbl=0;

  double ccf_max=0;
  double ccf_avg=0;

  double* ccf_nrm=NULL;
  double pix_passed_ord=0;

  int tot_line=0;
  double ll_range_ord=0;
  double* ccf_res=NULL;
  double* ccf_fit=NULL;
  double* ccf_o=NULL;
  char wstart_key[80];

  /* Program's Id */
  SCSPRO("flames_corvel");

 
  memset(ident, '\0', 73);
  memset(cunit[0], '\0', 48);
  strncpy(cunit[1], "PIXEL           ", 16);
  strncpy(cunit[2], "PIXEL           ", 16);
  /* ================================================================ */
  /* GET INPUT DATA                                                   */
  /* ================================================================ */
  /* get input ima name */

  in_ima_ord=IN_N;
  if((midas_status = SCKGETC(IN_A,1,MAX_LEN,&midas_nval,in_ima)) !=0) {
     uves_msg_warning("Error reading char keyword %s",IN_A);
     return flames_midas_error(MAREMMA);
  }

  //sprintf(in_ima,IN_A);
  /* get input ima order number */
  //midas_status = SCKRDI(IN_N,1,1,&midas_nval,&in_ima_ord, 
  //          &midas_unit, &midas_null);

  /* ================================================================ */
  /* Read 2D extracted input spectra */
  /* ================================================================ */
  /* get input ima frame */

  
  if( (midas_status = SCFOPN(in_ima,D_R4_FORMAT,0,F_IMA_TYPE,&in_ima_id))!=0) {
     uves_msg_warning("Error opening input image %s",IN_A);
     return flames_midas_error(MAREMMA);
  }

  /* get input ima dimension */
  if((midas_status = SCDRDI(in_ima_id,"NAXIS",1,1,&midas_nval,
                            &in_ima_naxis,&midas_unit,&midas_null)) !=0)
  {
     uves_msg_warning("Error reading NAXIS from image %s",IN_A);
     return flames_midas_error(MAREMMA);
  } 
  /* get input ima no of columns and rows */

  if((midas_status = SCDRDI(in_ima_id,"NPIX",1,in_ima_naxis,&midas_nval,
                            in_ima_npix,&midas_unit,&midas_null))!=0) {
     uves_msg_warning("Error reading NPIX from image %s",IN_A);
     return flames_midas_error(MAREMMA);
  }



  if (in_ima_naxis > 1) {
     in_ima_nx = in_ima_npix[0];
     in_ima_ny = in_ima_npix[1];
  }
  else {
     in_ima_nx = in_ima_npix[0];
     in_ima_ny = 1;            /* input image is one extracted order */
  }

  /* Prepare memory area to hold input image */
  m_in_ima = matrix(0, in_ima_ny-1, 0, in_ima_nx-1);
 
  memset(&m_in_ima[0][0], '\0', in_ima_nx*in_ima_ny*sizeof(float));

  /* get input ima in prepared area */
  if((midas_status = SCFGET(in_ima_id,1,in_ima_nx*in_ima_ny,&midas_nval, 
                            (char *)&m_in_ima[0][0])) != 0) {
    uves_msg_warning("Error mapping image %s",IN_A);
     return flames_midas_error(MAREMMA);

  }

  /* ================================================================ */
  /* PREPARE WCAL SOLUTION                                            */
  /* ================================================================ */
  /* get WSTART and WSTEP values to calculate array of wcal pix values
     in_ima_wcal_sol stores the wavelength calibration solution */


  sprintf(wstart_key,"%s%d","WSTART",in_ima_ord);
  if((midas_status = SCDRDD(in_ima_id,wstart_key,1,1,
                        &midas_nval,&tmp_double,
                            &midas_unit,&midas_null)) != 0) {
     uves_msg_warning("Error reading %s from input image %s",wstart_key,IN_A);
     return flames_midas_error(MAREMMA);
  }

  
 
  in_ima_wstart=(float)tmp_double;

 
  if((midas_status = SCDRDD(in_ima_id,"CDELT1",1,1,&midas_nval,&tmp_double,
                            &midas_unit,&midas_null))!=0) {
     uves_msg_warning("Error reading CDELT1 from input image %s",IN_A);
     return flames_midas_error(MAREMMA);
  }
 
  
  in_ima_wstep=(float)tmp_double;
  in_ima_wcal_sol=dvector(0,in_ima_nx);


  for (i=0; i< in_ima_nx; i++){
    in_ima_wcal_sol[i]=(double)(in_ima_wstart+in_ima_wstep*i);
  }
  /* get input mask table name */
  if((midas_status = SCKGETC(IN_B,1,MAX_LEN,&midas_nval,in_msk))!=0) {
     uves_msg_warning("Error reading input table %s",IN_B);
     return flames_midas_error(MAREMMA);
  }

  /* ================================================================ */
  /* GET INPUT MASK                                                   */
  /* ================================================================ */
  /* ================================================================ */
  /*  
   The input mask is as follows. 
   First column tells you where the hole begins, 
   The second one where the hole ends and the third is the weight of each hole
   (this last value is important for the stellar case where one may want to 
    give more importance to stellar lines of a given type)
   We get from the input mask the following parameters:
   1) the minimum weight of the holes of the mask used in the CCF 
   2) the width of the holes 
   3) the weight of the holes 

   After this operation the parameters which counts are:
       msk_hole_siz_selw[i] 
       msk_hole_cen_selw[i]  
       msk_hole_wgt_selw[i]  

   in_msk_wgt_min=1;
   in_msk_hole_wid=1.;
   strcpy(log_opt," ");

   get_mask(in_msk,in_msk_wgt_min,in_msk_hole_wid,log_opt,
           &msk_hole_width,&msk_hole_center,&msk_hole_wgt);
  */
  /* ================================================================ */
  /* get input mask table frame */


  if((midas_status = TCTOPN(in_msk,F_I_MODE,&in_msk_id))!=0) {
     uves_msg_warning("Error reading input mask %s",in_msk);
     return flames_midas_error(MAREMMA);
  }


  TCIGET (in_msk_id, &in_msk_ncol, &in_msk_nrow);
  /* get input mask table column id */
  if((midas_status = TCCSER(in_msk_id,"WSTART",&wstart_id))!=0) {
    uves_msg_warning("Error reading WSTART from input mask %s",in_msk);
    return flames_midas_error(MAREMMA);
  }

  /* get input mask table column id */
  if((midas_status = TCCSER(in_msk_id,"WEND",&wend_id))!=0) {
    uves_msg_warning("Error reading WEND from input mask %s",in_msk);
    return flames_midas_error(MAREMMA);
  }

  /* get input mask table column id */
  if((midas_status = TCCSER(in_msk_id,"WEIGHT",&weight_id))!=0) {
    uves_msg_warning("Error reading WEIGHT from input mask %s",in_msk);
    return flames_midas_error(MAREMMA);
  }

  /* Defines and initializes all necessary vectors */
  msk_hole_sta=dvector(0,in_msk_nrow);  
  msk_hole_end=dvector(0,in_msk_nrow);    
  msk_hole_siz=dvector(0,in_msk_nrow);
  msk_hole_wgt=dvector(0,in_msk_nrow);
  msk_hole_cen=dvector(0,in_msk_nrow);

  /* selected values...*/
  msk_hole_siz_selw=dvector(0,in_msk_nrow);
  msk_hole_wgt_selw=dvector(0,in_msk_nrow);
  msk_hole_cen_selw=dvector(0,in_msk_nrow);


  for(i=1;i<in_msk_nrow;i++) {
     TCERDD(in_msk_id,i,wstart_id,&tmp_dbl,&midas_null); 
     msk_hole_sta[i-1]=tmp_dbl;
     TCERDD(in_msk_id,i,wend_id,&tmp_dbl,&midas_null);
     msk_hole_end[i-1]=tmp_dbl;
     TCERDD(in_msk_id,i,weight_id,&tmp_dbl,&midas_null);
     msk_hole_wgt[i-1]=tmp_dbl;
     msk_hole_siz[i-1]=msk_hole_end[i-1]-msk_hole_sta[i-1];
     msk_hole_cen[i-1]=msk_hole_sta[i-1]+msk_hole_siz[i-1]*0.5;
     /*
     uves_msg_debug("sta=%f end=%f wgt=%f siz=%f cen=%f",
              msk_hole_sta[i-1],
              msk_hole_end[i-1],
              msk_hole_wgt[i-1],
              msk_hole_siz[i-1],
              msk_hole_cen[i-1]);
     */
  }
  TCTCLO(in_msk_id);
  /*ADAPTED*****/
  /* 
     If a fixed width is given as input parameter in_msk_hole_wid then
     is calculated msk_hole_siz
     in_msk_hole_wid is the fixed width given in km/s
     In our case in_msk_hole_wid=0 and the following if is not entered
  */
  
  if (in_msk_hole_wid > 0) {

     for(i=1;i<in_msk_nrow;i++) {
        msk_hole_siz[i-1]=
           in_msk_hole_wid*msk_hole_siz[i-1]/FLAMES_SPEED_OF_LIGHT;
     }

  }


  /* selects mask on in_msk_wgt_min of force weight=1
     in our case in_msk_wgt_min =1 and the following if is not entered
     is executed instead the else part
  */
  if (in_msk_wgt_min < 1) {
    /* If a lower limit of the weight of the holes is specified as input 
       parameter in_msk_wgt_min, then selects values of 
       wsize,wcenter,weight
       If no condition is given keep the vectors intact
    */
     int counter=0;
     for(i=1;i<in_msk_nrow;i++) {
       if (msk_hole_wgt[counter] > in_msk_wgt_min) {
          msk_hole_siz_selw[counter] = msk_hole_siz[i];  
          msk_hole_cen_selw[counter] = msk_hole_cen[i];  
          msk_hole_wgt_selw[counter] = msk_hole_wgt[i]; 
          counter++; 
       }
     }
  }
  else {
     for(i=1;i<in_msk_nrow;i++) {
       if (msk_hole_wgt[i] > in_msk_wgt_min) {
          msk_hole_siz_selw[i] = msk_hole_siz[i];  
          msk_hole_cen_selw[i] = msk_hole_cen[i];  
          msk_hole_wgt_selw[i] = msk_hole_wgt[i]; 
       }
     }
  }

 
  /* ================================================================ */
  /* END GET INPUT MASK                                               */
  /* ================================================================ */
  /* ================================================================ */
  /* COMPUTE CCF                                                      */
  /* ================================================================ */
  /* we allocate memory and define the vector to be used to evaluate CCF */
  /* this vector defines the points at which the CCF is computed */
  //midas_status = SCKRDD(IN_C,1,3,&midas_nval,rv_ccf_par, 
  //        &midas_unit, &midas_null);

  rv_ccf_size=(int)((rv_ccf_max-rv_ccf_min)/rv_ccf_step+1);
  rv_ccf=dvector(0,rv_ccf_size);
  ccf_o=dvector(0,rv_ccf_size);

 
  rv_ccf[0]=rv_ccf_min;
  for(i=1;i<rv_ccf_size;i++){
    rv_ccf[i]=rv_ccf[i-1]+rv_ccf_step;
  }
  /* 
     ======================================================================= 
     Do correlation. Values calculated by this subriutine are:
       ccf:            ccf matrix containing the ccf for each order (ccf_i)

       ccf_max:        vector containing the highest value of each ccf_i, 

       pix_passed_all: number of pixels of the input spectrum used for the
                       computation of each ccf_i,
       pix_passed_ord is the currespondent order value

       tot_line:       number of holes used in the computation of each ccf_i, 

       ll_range_all:   wavelength interval of each order the input spectrum 
                       used in the computation of each ccf_i
       ll_range_ord is the correspondent order value
     ======================================================================= 
   */


  do_cor_vel(in_ima_wcal_sol,   /* wave calibration solution */
         m_in_ima,          /* extracted spectrum */
         rv_ccf,            /* points at which the CCF is computed */
         msk_hole_siz_selw, /* hole size   selected on weight criteria */
         msk_hole_cen_selw, /* hole center selected on weight criteria */
         msk_hole_wgt_selw, /* hole weight selected on weight criteria */
             0,                 /* barv     :Baricentric Velocity Corr */
             0,                 /* barv_max :Its maximum               */
             0,                 /* fit_type (Gaussian): 0/1 emis/absorb */
             in_ima_nx,         /* X sise of input spectra */
             in_msk_nrow,       /* size of input mask */
         rv_ccf_size,       /* size of CCF */
         ccf_o,               /* out: ccf for each order (ccf_i) */     
             &ccf_max,          /* out: max(ccf) for each order (ccf_i) */
             &pix_passed_ord,   /* out: each order in sp's no of pix to 
                                             get ccf_i */
         &tot_line,         /* out: no of holes used to get ccf_i */
         &ll_range_ord,     /* out: each order's wav interval to get ccf_i */
             in_ima_id);        /* input ima id (to write descriptors) */

 
  /* Sum the individual ccf_i for each bin and normalize the final ccf */

  SCFCLO(in_ima_id); //not needed anymore
  ccf_nrm=dvector(0,rv_ccf_size);
  for(i=0;i<rv_ccf_size;i++){
     ccf_avg +=ccf_o[i];
     if(!irplib_isinf(ccf_o[i])) {
        if(ccf_o[i] > ccf_max) {
           ccf_max=ccf_o[i];
        }
     }
  }
 

  /* Creating a new table for offline plotting of peaks */
  SCKGETC(OU_A,1,MAX_LEN,&midas_nval,ou_tab);
  /* jmlarsen: use F_O_MODE for new table
     old code: TCTINI(ou_tab,F_IO_MODE,rv_ccf_size,&tid);*/
  TCTINI(ou_tab,F_O_MODE,rv_ccf_size,&tid);
  
  /* Creating a new column */
  TCCINI(tid, D_R8_FORMAT, 1, "F8.4", " ", "ccf_pos", &ccf_pos_col);
  TCCINI(tid, D_R8_FORMAT, 1, "F8.4", " ", "ccf_nrm", &ccf_nrm_col);
  TCCINI(tid, D_R8_FORMAT, 1, "F8.4", " ", "ccf_out", &ccf_out_col);

  /* Writing table values */
  /*
  if (abs(ccf_max) >= FEPSILON) { 
     for(i=0;i<rv_ccf_size;i++){
        ccf_nrm[i]=ccf_o[i]/ccf_max;
        TCEWRD(tid, i+1, ccf_pos_col, &rv_ccf[i]);
        TCEWRD(tid, i+1, ccf_nrm_col, &ccf_nrm[i]);
        TCEWRD(tid, i+1, ccf_out_col, &ccf_o[i]);
     }
  } else {
     for(i=0;i<rv_ccf_size;i++){
        ccf_nrm[i]=0.;
        TCEWRD(tid, i+1, ccf_pos_col, &rv_ccf[i]);
        TCEWRD(tid, i+1, ccf_nrm_col, &ccf_nrm[i]);
        TCEWRD(tid, i+1, ccf_out_col, &ccf_o[i]);
     }
  }
  */
 
     for(i=0;i<rv_ccf_size;i++){
        ccf_nrm[i]=ccf_o[i]/ccf_max;
        TCEWRD(tid, i+1, ccf_pos_col, &rv_ccf[i]);
        TCEWRD(tid, i+1, ccf_nrm_col, &ccf_nrm[i]);
        TCEWRD(tid, i+1, ccf_out_col, &ccf_o[i]);
     }

  SCDWRD(tid,"CCF_MAX",&ccf_max,1,1,&midas_unit); 
  SCDWRD(tid,"WAV_RNG",&ll_range_ord,1,1,&midas_unit); 
  SCDWRD(tid,"PIX_TOT",&pix_passed_ord,1,1,&midas_unit); 
  SCDWRI(tid,"LIN_TOT",&tot_line,1,1,&midas_unit); 

  TCTCLO(tid);
 
 
  /* TO BE IMPLEMENTED */
  /* Gaussian Fit of the normalized CCF */
  /* 
     one fit normalized_ccf as a function of rv_ccf using as fit type an
     emission Gaussian. Output of the fit are the Gaussian fit coefficients 
     ccf_res and ccf_fit is the fitted Gaussian computed on the rv_ccf 
     velocity bins
  */

  /* ccf_res[0]=ccf_res[0]/(1.-ccf_res[3]); */
  fit_ccf(rv_ccf,ccf_nrm,1,ccf_res,ccf_fit);


 
  /* dump results in ouput image*/

  ou_ima_npix[0]=rv_ccf_size;
  ou_ima_npix[1]=1;
  ou_ima_start[0]=rv_ccf[0];
  ou_ima_start[1]=ccf_nrm[0];
  ou_ima_step[0]=ccf_max;
  ou_ima_step[1]=1;
  cuts[0] = 0;
  cuts[1] = 0;
  cuts[2] = 0;
  cuts[3] = 1;
 
 
  SCKGETC(OU_B,1,MAX_LEN,&midas_nval,ou_ima);

  SCFCRE(ou_ima,D_R8_FORMAT,F_O_MODE,F_IMA_TYPE,rv_ccf_size,&ou_ima_id);
  SCDWRC(ou_ima_id,"IDENT", 1, ident, 1, 72, &midas_unit);
  SCDWRI(ou_ima_id,"NAXIS",&ou_ima_naxis,1,1,&midas_unit); 
  SCDWRI(ou_ima_id,"NPIX",ou_ima_npix,1,2,&midas_unit); 
  SCDWRD(ou_ima_id,"START",ou_ima_start, 1, 2, &midas_unit);
  SCDWRD(ou_ima_id,"STEP", ou_ima_step, 1, 2, &midas_unit);
  SCDWRC(ou_ima_id,"CUNIT", 1, cunit[0], 1, 48, &midas_unit);
  SCDWRR(ou_ima_id,"LHCUTS", cuts, 1, 4, &midas_unit);
  SCFPUT(ou_ima_id,1,rv_ccf_size,(char *)ccf_o);
  SCDWRD(ou_ima_id,"CCF_MAX",&ccf_max,1,1,&midas_unit); 
  SCDWRD(ou_ima_id,"WAV_RNG",&ll_range_ord,1,1,&midas_unit); 
  SCDWRD(ou_ima_id,"PIX_TOT",&pix_passed_ord,1,1,&midas_unit); 
  SCDWRI(ou_ima_id,"LIN_TOT",&tot_line,1,1,&midas_unit); 
  SCFCLO(ou_ima_id);

  cuts[3] = ccf_max;
  SCKGETC(OU_C,1,MAX_LEN,&midas_nval,ou_ima);
  SCFCRE(ou_ima,D_R8_FORMAT,F_O_MODE,F_IMA_TYPE,rv_ccf_size,&ou_ima_id);
  
  SCDWRC(ou_ima_id,"IDENT", 1, ident, 1, 72, &midas_unit);
  SCDWRI(ou_ima_id,"NAXIS",&ou_ima_naxis,1,1,&midas_unit); 
  SCDWRI(ou_ima_id,"NPIX",ou_ima_npix,1,2,&midas_unit); 
  SCDWRD(ou_ima_id,"START",ou_ima_start, 1, 2, &midas_unit);
  SCDWRD(ou_ima_id,"STEP", ou_ima_step, 1, 2, &midas_unit);
  SCDWRC(ou_ima_id,"CUNIT", 1, cunit[0], 1, 48, &midas_unit);
  SCDWRR(ou_ima_id,"LHCUTS", cuts, 1, 4, &midas_unit);
  SCFPUT(ou_ima_id,1,rv_ccf_size,(char *)ccf_nrm);
  SCDWRD(ou_ima_id,"CCF_MAX",&ccf_max,1,1,&midas_unit); 
  SCDWRD(ou_ima_id,"WAV_RNG",&ll_range_ord,1,1,&midas_unit); 
  SCDWRD(ou_ima_id,"PIX_TOT",&pix_passed_ord,1,1,&midas_unit); 
  SCDWRI(ou_ima_id,"LIN_TOT",&tot_line,1,1,&midas_unit); 
  SCFCLO(ou_ima_id);
  /* free allocated memory */
  /* free_matrix(m_in_ima,0,in_ima_ny-1,0,in_ima_nx-1); */
  free_dvector(msk_hole_sta,0,in_msk_nrow);
  free_dvector(msk_hole_end,0,in_msk_nrow);
  free_dvector(msk_hole_siz,0,in_msk_nrow);
  free_dvector(msk_hole_wgt,0,in_msk_nrow);
  free_dvector(msk_hole_cen,0,in_msk_nrow);
  free_dvector(msk_hole_siz_selw,0,in_msk_nrow);
  free_dvector(msk_hole_wgt_selw,0,in_msk_nrow);
  free_dvector(msk_hole_cen_selw,0,in_msk_nrow);
  free_dvector(rv_ccf,0,rv_ccf_size);
  free_dvector(ccf_nrm,0,rv_ccf_size);
  free_dvector(in_ima_wcal_sol,0,in_ima_nx);
  free_dvector(ccf_o,0,rv_ccf_size);


  SCSEPI();
  return 0;

}

void
do_cor_vel(double* wcal_sol,float** sp_flux,double* rv_ccf,
           double* msk_hole_siz,double* msk_hole_cen,
           double* msk_hole_wgt,double bar_v,double bar_v_max,
           int fit_type,int in_ima_ncol,int in_msk_nrow,
           int rv_ccf_size,
       double* ccf_o,            /* matrix with ccf_i */
       double* ccf_max,        /* vector with max(ccf_i) */
       double* pix_passed_ord, /* no of in spct pixels used to get ccf_i */
       int*    tot_line,       /* no of holes used to get ccf_i */
       double* wcal_range_ord, /* wave range of each order in spct used to get ccf_i */
           int in_ima_id)
{

  /* Local variables */
  double* dw_map=NULL;
  double* ccf_all=NULL;
  double* ccf_all_fit=NULL;
  double* msk_hole_cen_selr=NULL;
  double* msk_hole_siz_selr=NULL;
  double* msk_hole_wgt_selr=NULL;

  double* ccf_o_results=NULL;
  /* double* ccf_o_fit=NULL; */
  /* ccf_o_fit is commented out as not really used */
  double* rv_ccf_cor=NULL;

  double wcal_min=0;
  double wcal_max=0;
  double d_secular_red=0;
  double d_secular_blu=0;

  int i=0;
  int sel_no=0;

  /* Local Functions */

/* 
   ==========================================================================
   Subroutine body 
   ==========================================================================
*/

/* The following 2 lines has de facto no effect as bar_v and bar__max are 0 */
  d_secular_red=bar_v_max-bar_v;
  d_secular_blu=bar_v_max-bar_v;


  dw_map=dvector(0,in_ima_ncol);
  ccf_all=dvector(0,rv_ccf_size);
  ccf_all_fit=dvector(0,rv_ccf_size);
  rv_ccf_cor=dvector(0,rv_ccf_size);
  /* ccf_o_fit=dvector(0,in_ima_ncol); */
  /* ccf_o_fit is commented out as not really used*/
  ccf_o_results=dvector(0,4);


  msk_hole_cen_selr=dvector(0,in_msk_nrow);
  msk_hole_siz_selr=dvector(0,in_msk_nrow);
  msk_hole_wgt_selr=dvector(0,in_msk_nrow);
  /* defines delta_lambda vector as delta_lambda=lambda(i+1)-lambda(i) */
  for(i=0;i<in_ima_ncol-1;i++){
    dw_map[i]=wcal_sol[i+1]-wcal_sol[i];
  }
  /* Not relevant for the ThAr correlation.
     This computes the minimum and the maximum wavelengths given the velocity
     point extremes in which the CCF is going to be computed
     (rv_ccf[0] is the first velocity bin and rv_ccf[-1] is the last) and the
     max BAR_V velocity (baricentric velocity) possible
  */

  
  /* Here should start a loop over orders: we do not do it as we assume
     to have in input the spectra relative to each order */

  /* The following two lines are not relevant in case of ThAr spectra */
  /* They are to compute the min and max wavelength being given the velocity
     point extremes in which the CCF is going to be computed and the max
     baricentric velocity possible */

  wcal_min=wcal_sol[0]-(rv_ccf[0]-bar_v-d_secular_blu)*
           wcal_sol[0]/FLAMES_SPEED_OF_LIGHT; 

  wcal_max=wcal_sol[in_ima_ncol-1]-(rv_ccf[rv_ccf_size-1]-bar_v+d_secular_red)*           
           wcal_sol[in_ima_ncol-1]/FLAMES_SPEED_OF_LIGHT;


/*
>From the python version:

ll_max=ll_map[order,-1]-(RV_CCF[-1]-berv+D_secular_red)*ll_map[order,-1]/speed_of_light
*/ 
  /* Filter wcenter,wsize,weight to include holes whose center is within the
     limits wcal_min and wcal_max
  */

  for(i=0;i<in_msk_nrow;i++){
    if((msk_hole_cen[i]>wcal_min) && (msk_hole_cen[i]<wcal_max)) {
      msk_hole_cen_selr[sel_no]=msk_hole_cen[i];
      msk_hole_siz_selr[sel_no]=msk_hole_siz[i];
      msk_hole_wgt_selr[sel_no]=msk_hole_wgt[i];
      sel_no++;
    }
  }
  *tot_line=sel_no;

  if(sel_no) {
    /* If at least one is left after filtering the mask */
    *wcal_range_ord=0.;
  /* we get the velocity bins were the CCF is going to be computed 
     corrected for bar_v */
    for(i=0;i<rv_ccf_size;i++){
        rv_ccf_cor[i]=rv_ccf[i]-bar_v; 
    }

    /* computes the ccf on the order order. 

       The input arguments are:
       msk_hole_cen_selr, centers of each hole selected on wave range criteria
       msk_hole_siz_selr, widths  of each hole selected on wave range criteria
       msk_hole_wgt_selr, weights of each hole selected on wave range criteria
       wcal_sol, the vector containing the correspondence pixel to 
                 lambda for the order order
       sp_flux[order] is the vector containing the intensity of each pixel
                  for the order order
       dw_map is the delta lambda between consecutive pixels
       rv_ccf-bar_v is the velocity bin where the CCF is going to be 
                    computed corrected for the BAR_V.

       OUTPUT arguments are:
       ccf_o, the ccf of the order order,
       pix_passed tells you how many pixels have participated in the ccf,
       wcal_range is the length (in Angstroms) of the region covered by the
                  holes which participated in the CCF 
                  (i.e., the sum of the vector wcal_msk_size_selr);

    */
 
  do_ccf_f(msk_hole_cen_selr, msk_hole_siz_selr, msk_hole_wgt_selr,  
       wcal_sol, sp_flux, dw_map, rv_ccf_cor, ccf_o, pix_passed_ord,
           wcal_range_ord, sel_no, in_ima_ncol, rv_ccf_size, in_ima_id);

  }
  else {
    /* there is no mas holes in the wavelength interval wcal_min, wcal_max
       then everything is set to zero */
    printf("No hole between wcal_min=%f and wcal_max=%f all set to 0. \n",
            wcal_min,wcal_max);
     for(i=0;i<rv_ccf_size;i++){
       /* rv_ccf[i]=0.; */
         ccf_o[i]=rv_ccf[i]*0.;
         /* ccf_o_fit[i]=ccf_o[i]; */
         /* ccf_o_fit is commented out as not really used */
     }
     *pix_passed_ord=0.;
     *wcal_range_ord=0.;

    ccf_o_results[0]=0.;
    ccf_o_results[1]=0.;
    ccf_o_results[2]=0.;
    ccf_o_results[3]=0.;

  }
 

  /* write results on output table */

  /* Free memory */
  free_dvector(rv_ccf_cor,0,rv_ccf_size);
  free_dvector(dw_map,0,in_ima_ncol);
  /* free_dvector(ccf_o_fit,0,in_ima_ncol); */
  /* ccf_o_fit is commented out as not really used */
  free_dvector(ccf_o_results,0,4);


  free_dvector(ccf_all,0,rv_ccf_size);
  free_dvector(ccf_all_fit,0,rv_ccf_size);

  free_dvector(msk_hole_cen_selr,0,in_msk_nrow);
  free_dvector(msk_hole_siz_selr,0,in_msk_nrow);
  free_dvector(msk_hole_wgt_selr,0,in_msk_nrow);

  return;
 
} /* end function do_corvel */


void 
do_ccf_f(double* mask_ll,double* mask_d,double* mask_w,double* sp_ll,
         float** sp_flux,double* sp_dll,double* rv_ccf,double* ccf_o,
         double* pix_tot,double* ll_range_tot,int in_msk_nrow, 
         int in_ima_ncol, int rv_ccf_size, int in_ima_id)
{

  /* This routine should evaluate and return:
     ccf_o[rv_ccf_size]-the resulting CCF for a given order (not normalized)
     pix_passed-a double scalar
     ll_range-a double scalar
  */ 
     /* iter for v */ 

  /* at rest the mask holes are centered on the vector mask_ll.
     at a velocity rv, they will be centered on 
     mask_ll+rv*mask_ll/FLAMES_SPEED_OF_LIGHT
     The blue edge of the holes (Mask_blue) is then this new center minus
     half of the size of the hole. The same is valid for the red edge of 
     the hole.

  */

  /* local variable definition-initializzation */

  double** covar;
  double** alpha;

  double* msk_blu=NULL;
  double* msk_red=NULL;
  double* sp_ll_prime=NULL;
  double* sfit=NULL;
  double* xfit=NULL;
  double* yfit=NULL;
  double* aa=NULL;
  double* erraa=NULL;


  int* i_blu_masques=NULL;
  int* i_red_masques=NULL;
  int* ia=NULL;

  double intensity_s=0;
  double pix=0;
  double ll_range=0;
  double norm=0;
  double cen=0;
  double sig=0;
  double fwhm=0;
  double rv=0;

  int i=0;
  int j=0;

  int midas_unit = 0;
  int sp_ll_sz = in_ima_ncol;

  int ndeg=6;

  /* Function prototype */


  xfit=dvector(1,rv_ccf_size);
  yfit=dvector(1,rv_ccf_size);
  sfit=dvector(1,rv_ccf_size);

  covar = dmatrix(1,ndeg,1,ndeg);
  alpha = dmatrix(1,ndeg,1,ndeg);

  aa=dvector(1,ndeg);
  erraa=dvector(1,ndeg);
  ia=ivector(1,ndeg);
  sp_ll_prime=dvector(0,in_ima_ncol);
  msk_blu = dvector(0,in_msk_nrow);
  msk_red = dvector(0,in_msk_nrow);
  i_blu_masques = ivector(0,in_msk_nrow);
  i_red_masques = ivector(0,in_msk_nrow);

  for(i=0;i<rv_ccf_size;i++) {
    rv=rv_ccf[i];
    sfit[i]=1.0;
  }

  for(i=0;i<rv_ccf_size;i++) {
    rv=rv_ccf[i];
    /* 
       we define the 1st derivative: sp_ll_prime[j]=sp_ll[j]+sp_dll[j]*0.5; 
       j is a counter variable of values up to sp_ll_sz equal to the No of
       extracted spectra definition points
    */
    for(j=0; j<sp_ll_sz; j++) {
      sp_ll_prime[j]=sp_ll[j]+sp_dll[j]*0.5;
    }
    for(j=0;j<in_msk_nrow;j++) {
       /* shift the mask holes for a velocity RV[i] */
      msk_blu[j]=mask_ll[j]+rv*mask_ll[j]/FLAMES_SPEED_OF_LIGHT-0.5*mask_d[j];
      msk_red[j]=mask_ll[j]+rv*mask_ll[j]/FLAMES_SPEED_OF_LIGHT+0.5*mask_d[j];

      /*  
      The idea is to know where (i.e. in which pixel) a given hole will start 
      because we won't want to scan through the vector wave to find the pixel 
      i where lambda(i-1) < mask_start <lambda(i). The command search_sorted 
      does it (see below).
      It returns the position where the element mask_blue will fit in the 
      vector lamda+delta_lambda/2.
      This is done for the blue edge of the mask and for the red edge. 
      The +1 in the end is because phython vectors starts at 0 and F77 at 1.
      */
      
    }
    /*
    Look for the first and the last holes available for the crooss-correlation
    assuming the spectrum has a dimension nspec and sp_ll(nspec) and 
    flux(nspec) are the wavelength and spectral flux vectors

    Then finds the first hole such as
    wave[0]<=mask_blu[first_hole] && mask_red[first_hole-1]<wave[0]
    find last_hole such as 
    wave[nspec]>=mask_red[first_hole] && mask_red[first_hole+1]>wave[nspec]
    
    This search is done using
    find_pos_d(vector,len(vector),x,i,j,guess)
    which returns the index of the element in the vector such as
    vector[i]<=x<vector[i+1]
    The search is carried out between the elements:
    vector[i] and vector[j]
    and using "guess" and first "guess" for the position of "x" within "vector"

    (see NR F77 chapter 3.4)
    */
    
    int first_hole=hunt(msk_blu-1, in_msk_nrow, sp_ll[0],0);
    int guess=first_hole;

    for(j=0;j<in_msk_nrow;j++) {
      //for(j=0;j<3;j++) {

       i_blu_masques[j]=hunt(sp_ll_prime-1,sp_ll_sz,msk_blu[j],0)+1;

       guess=i_blu_masques[j];
       i_red_masques[j]=hunt(sp_ll_prime-1,sp_ll_sz,msk_red[j],guess)+1;
       guess=i_red_masques[j];
       //uves_msg_debug("masques: %d %d",i_blu_masques[j],i_red_masques[j]);

    }
  correl_bin(sp_ll_sz,sp_flux,sp_ll,sp_dll,
            &in_msk_nrow,msk_blu,msk_red,mask_w,i_blu_masques,
            i_red_masques,&intensity_s,&pix,&ll_range);

  ccf_o[i]=intensity_s;

  }
  *pix_tot+=pix;
  *ll_range_tot+=ll_range;
  for(i=0;i<rv_ccf_size;i++) {
    j=i+1;
    xfit[j]=rv_ccf[i];
    yfit[j]=ccf_o[i];
    sfit[j]=1;
  }

  aa[1]=300;
  aa[2]=0;
  aa[3]=1;
  aa[4]=1.;


  ia[1]=1;
  ia[2]=1;
  ia[3]=1;
  ia[4]=0;

  gaussian_fit(rv_ccf,ccf_o,rv_ccf_size,&norm,&cen,&sig,&fwhm);

  /* write output in descriptor */
  uves_msg_debug("Position max corvel=%f",cen);
  SCDWRD(in_ima_id,"CORVEL_MAX",&cen,1,1,&midas_unit);



  /* Free allocated memory */
  free_dmatrix(covar,1,ndeg,1,ndeg);
  free_dmatrix(alpha,1,ndeg,1,ndeg);

  free_dvector(aa,1,ndeg);
  free_dvector(erraa,1,ndeg);
  free_ivector(ia,1,ndeg);
 
  free_dvector(xfit,1,rv_ccf_size);
  free_dvector(yfit,1,rv_ccf_size);
  free_dvector(sfit,1,rv_ccf_size);

  free_dvector(msk_blu,0,in_msk_nrow);
  free_dvector(msk_red,0,in_msk_nrow);
  free_ivector(i_blu_masques,0,in_msk_nrow);
  free_ivector(i_red_masques,0,in_msk_nrow);

  free_dvector(sp_ll_prime,0,in_ima_ncol);
  
} /* end function do_ccf_f */

void 
correl_bin(int nx,        /* in: dimension of flux (is it necessary?) */
        float** flux,   /* in: Spectral flux (dim nx) */
        double *ll,        /* in: wavelength (dim nx) */
        double *dll,    /* in: Delta(lambda) D_ll (dim nx) */
        int *nbr_trou,  /* in: Number of holes read from the mask file */
        double *ll_s,   /* in: Mask hole start wavelength (dim nbr_trou) */
        double *ll_e,   /* in: Mask hole end wavelength (dim nbr_trou) */
        double *ll_wei, /* in: Mask hole weight wavelength (dim nbr_trou) */
        int *i_start,   /* in: see python code (page 9 line 76-77)  */
        int *i_end,        /* in: see python code (page 9 line 76-77)  */
        double *out_ccf, /* out: Value of the CCF for a given velocity 
                                point */
        double *pix,     /* out: number of pixelx used in the 
                                computation of the CCF */
        double *llrange) /* out: wavelenght interval covered by the 
                                pixels used in computation of the CCF */
{

    /* pointers */
    
    float *pflux=NULL;
    double *pll=NULL;
    double *pdll=NULL;
    double *pll_s=NULL;
    double *pll_e=NULL;
    double *pll_wei=NULL;
    int *pi_start=NULL;
    int *pi_end=NULL;
    int trou=0;
    int i=0;

    pflux = *flux;
    pll   = ll;
    pdll  = dll;
    pll_s = ll_s;
    pll_e = ll_e;
    pll_wei = ll_wei;
    pi_start = i_start;
    pi_end   = i_end;
    

    /*local param */
    
    
    *out_ccf=0.0;
    *pix=0.0;
    *llrange=0.0;

    
    for (trou=0;trou < *nbr_trou;trou++) {
        
      if (pi_start[trou] == pi_end[trou]) {
         *out_ccf=*out_ccf+(pll_e[trou]-pll_s[trou])/pdll[pi_start[trou]]*
               pflux[pi_start[trou]]*(pll_wei[trou]);
                 
         *pix=*pix+(pll_e[trou]-pll_s[trou])*pll_wei[trou]/
                   pdll[pi_start[trou]];
            
         *llrange=*llrange+(pll_e[trou]-pll_s[trou])*pll_wei[trou];
            
      } else if (pi_start[trou]+1 == pi_end[trou]) {
            
         *out_ccf=*out_ccf+
                   ((pll[pi_start[trou]]+pdll[pi_start[trou]]*.5-pll_s[trou])*
                 pflux[pi_start[trou]]/pdll[pi_start[trou]]+
           (pll_e[trou]-(pll[pi_start[trou]]+pdll[pi_start[trou]]*.5))*
              pflux[pi_end[trou]]/pdll[pi_start[trou]])*pll_wei[trou];
                
             *pix=*pix+((pll[pi_start[trou]]+pdll[pi_start[trou]]*.5-
                  pll_s[trou])/pdll[pi_start[trou]]+
        (pll_e[trou]-(pll[pi_start[trou]]+pdll[pi_start[trou]]*.5))/
        pdll[pi_end[trou]])*pll_wei[trou];
                
 
         *llrange=*llrange+((pll[pi_start[trou]]+pdll[pi_start[trou]]*.5-
                      pll_s[trou])+
        (pll_e[trou]-(pll[pi_start[trou]]+pdll[pi_start[trou]]*.5)))*
                      pll_wei[trou];


      } else {
        
         *out_ccf=*out_ccf+((pll[pi_start[trou]]+pdll[pi_start[trou]]*0.5-
                      pll_s[trou])*pflux[pi_start[trou]]/pdll[pi_start[trou]]+
               (pll_e[trou]-(pll[pi_end[trou]]-pdll[pi_end[trou]]*.5))*
            pflux[pi_end[trou]]/pdll[pi_end[trou]])*pll_wei[trou];

              *pix=*pix+
           ((pll[pi_start[trou]]+pdll[pi_start[trou]]*0.5-pll_s[trou])/
            pdll[pi_start[trou]]+
            (pll_e[trou]-(pll[pi_end[trou]]-pdll[pi_end[trou]]*.5))/
            pdll[pi_end[trou]])*pll_wei[trou];

               *llrange=*llrange+
          ((pll[pi_start[trou]]+pdll[pi_start[trou]]*0.5-pll_s[trou])+
           (pll_e[trou]-(pll[pi_end[trou]]-pdll[pi_end[trou]]*.5)))
            *pll_wei[trou];

          for (i=pi_start[trou]+1;i<=pi_end[trou]-1;i++) {

                  *out_ccf=*out_ccf+pflux[i]*pll_wei[trou];
                  *pix=*pix+pll_wei[trou];
                  *llrange=*llrange+pdll[i]*pll_wei[trou];
         }
      }
        }
} /* end function correl_bin */
            
            
void
fit_ccf(double* rv_ccf,double* ccf_nor,int type,double* ccf_res,
        double* ccf_fit)
{
  /* Gaussian Fit either in emission or in absorbtion depending on the flag, 
     emission for the ThAr 
     It first computes a single fit in order to find the first guess 
     parameters. Then it does the fit again now putting more weight on the 
     core of the Gaussian. It returns the fit coefficients and the fitted 
     function.

*/

} /* end function fit_ccf */




  
int 
hunt(double* xx, int n, double x, int jlo) {

    int jhi;

    
    int ascnd=(xx[n] >= xx[1]);
    if (jlo <= 0 || jlo >n) {
        jlo=0;
        jhi=n+1;
    } else {
        int inc=1;
        if ((x>=xx[jlo]) == ascnd) {
            if (jlo == n) return jlo-1;
            jhi=jlo+1;
            while ((x>=xx[jhi]) == ascnd) {
                jlo=jhi;
                inc +=inc;
                jhi=jlo+inc;
                if (jhi>n) {
                    jhi=n+1;
                    break;
                }
            }
        } else {
            if (jlo==1) {
                jlo=0;
                return jlo-1;
            }
            jhi=jlo--;
            while ((x<xx[jlo])==ascnd) {
                jhi=jlo;
                inc *=2;
                if (inc >= jhi) {
                    jlo=0;
                    break;
                }
                else jlo=jhi-inc;
            }
        }
    }
    while ((jhi-jlo) != 1) {
        int jm=(jhi+jlo) >> 1;
        if ( (x >= xx[jm]) ==ascnd)
            jlo=jm;
        else
            jhi=jm;
    }
    if (x == xx[n]) jlo=n-1;
    if (x == xx[1]) jlo=1;

    return jlo-1;
} /* end function hunt */


/*
void fgauss(double x,double g[],int ng)
{
    int i=0;
    double arg=0.0;
    double ex=0.0;
    double fac=0.0;
    arg=(x-g[2])/g[3];
    ex=exp(-arg*arg);
    fac=g[4]+g[1]*ex*2.0*arg;
    return fac;
}
*/


 /*
static double 
fgauss(double x,double a[],double y,double dyda[],int na)
{
    double arg=0.0;
    double ex=0.0;
    double fac=0.0;
   
    arg=(x-a[2])/a[3];
    ex=exp(-arg*arg);
    fac=a[4]+a[1]*ex*2.0*arg;
    y = a[4]+fac;

    dyda[1]=ex;
    dyda[2]=fac/a[2];
    dyda[3]=fac*arg/a[2];
    dyda[4]=0;
    return fac;
}
 */

/*
static void 
fpoly(double x,double p[],int np)
{
    int j=0;
    p[1]=1;
    for (j=2; j<=np;j++) p[j]=p[j-1]*x;
}
*/


static void
gaussian_fit(const double *   x,
             const double *   y,
                   int        size,
                   double *   norm,
                   double *   xcen,
                   double *   sig_x,
                   double *   fwhm_x) 
{
    double          u0, ux, uxx;
    double          max_val ;
    int             i;

    /* Check entries */
    /* Extraction zone */
    
    /* Extract the image zone to fit */
    /* Check if there are enough good pixels */
    /* Convert the image to double */
    /* Compute xcen  */
    u0 = ux = 0.0 ;
    for (i=0 ; i<size ; i++) {
        u0 += y[i] ;
        ux += x[i] * y[i] ;
    }
    /* Compute sig_x  */
    uxx = 0.0 ;
    for (i=0 ; i<size ; i++) {
        uxx += (x[i]-(ux/u0)) * (x[i]-(ux/u0)) * y[i] ;
    }
    if (sig_x) *sig_x = sqrt(fabs(uxx/u0)) ;
    if (fwhm_x) *fwhm_x = 2 * sqrt(2 * log(2.0)) * sqrt(fabs(uxx/u0)) ;

    max_val=y[0];
    for (i=1 ; i<size ; i++) {
      if(y[i] > max_val) max_val=y[i];
    }
    /* Compute norm */
    if (norm) *norm = max_val*2*M_PI*sqrt(fabs(uxx/u0)) ; 
    
    /* Shift xcen and ycen to coordinates in the input big image */
    if (xcen) *xcen = ux/u0;
    
}