/*                                                                              *
 *   This file is part of the ESO UVES Pipeline                                 *
 *   Copyright (C) 2004,2005 European Southern Observatory                      *
 *                                                                              *
 *   This library 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       *
 *                                                                              */

/*
 * $Author: amodigli $
 * $Date: 2011-12-08 14:06:33 $
 * $Revision: 1.18 $
 * $Name: not supported by cvs2svn $
 */

/*----------------------------------------------------------------------------*/
/**
 * @addtogroup uves_physmod
 */
/*----------------------------------------------------------------------------*/

#ifdef HAVE_CONFIG_H
#  include <config.h>
#endif

/*-----------------------------------------------------------------------------
                                Includes
 -----------------------------------------------------------------------------*/
#include <uves_physmod_regress_echelle.h>

#include <uves_physmod_necregr.h>
#include <uves_utils_polynomial.h>
#include <uves_utils_wrappers.h>
#include <uves_pfits.h>
#include <uves_error.h>
#include <uves_msg.h>

#include <cpl.h>

/*-----------------------------------------------------------------------------
                                Defines
 -----------------------------------------------------------------------------*/
#define SELECT_FUNC_DOUBLE cpl_table_and_selected_double
#define SELECT_FUNC_INT    cpl_table_and_selected_int
/* Selected rows are ignored by many CPL functions (unlike MIDAS)
Therefore SELECT/TABLE t SELECT.and. ...... 
has different effect from cpl_table_and_selected_double(t, ...)

To remove selected rows (which is also not exactly what MIDAS does):
#define SELECT_FUNC_DOUBLE uves_extract_table_rows_local
#define SELECT_FUNC_INT    uves_extract_table_rows_local
(but this did not really make a difference)
*/
/**@{*/
/*-----------------------------------------------------------------------------
                            Functions prototypes
 ----------------------------------------------------------------------------*/
/*-----------------------------------------------------------------------------
                            Static variables
 -----------------------------------------------------------------------------*/

/*-----------------------------------------------------------------------------
                            Functions code
 -----------------------------------------------------------------------------*/

/*----------------------------------------------------------------------------*/
/**
  @brief    Fit mesured position of orders by a bivariate polynomial
  @param    raw_header  An array containing the input image headers. 
                         The ordering must be the same as the ordering 
                         of images in the input image list
  @param    chip         CCD chip
  @param    recipe_id    name of calling recipe
  @param    parameters   The recipe parameter list
  @param    o_tbl        The input order table
  @param    num_outliers The input no of outliers 
  @param    tol          The input tolerance
  @param    kappa        The input value of kappa in kappa-sigma clipping
  @param    s_tbl        The output table. 
  @param    w_tbl        The output table. 

  @return   0 if everything is ok, -1 otherwise
  Fit mesured position of orders by a bivariate polynomial
 */
/*----------------------------------------------------------------------------*/

int uves_physmod_regress_echelle(const uves_propertylist *raw_header,
                 enum uves_chip chip,
                 const char *recipe_id,
                 const cpl_parameterlist* parameters,
                 cpl_table** o_tbl, 
                 int num_outliers, 
                 double tol, 
                 double kappa,
                 cpl_table** s_tbl,
                 cpl_table** w_tbl)
{


  int  CNT =0;
  double THRESH=0.;
  int ORDER=0;
  float RMSMAX=0.;
  float RMS=0.;
  int REJECT[300];
  int NBREJ=0;
  double start1=0;
  double start2=0;
  double step1=0;
  double step2=0;
  int scan1=0;
  int scan2=0;
  int naxis1=0;
  int naxis2=0;

  int def_pol1=0;
  int def_pol2=0;

  double inpr1=0;
  double inpr2=0;
  double inpr3=0;
  char OUTMODE='V';
  int ECHORD_1=0;
  int INPUTI_1=0;
  int INPUTI_2=0;
  int imsize1=0;
  int nraw=0;
  int ord_min=0;
  int ord_max=0;
  int status=0;
  double mean_sq_err_y=0;
  double outputd[10][10];
  polynomial* poly2d_y=NULL;

  uves_propertylist* plist=NULL;




  check (start1 = uves_pfits_get_crval1(raw_header),
     "Could not read start factor from input header");

  check (start2 = uves_pfits_get_crval2(raw_header),
     "Could not read start factor from input header");
  
  check (step1 = uves_pfits_get_cdelt1(raw_header),
     "Could not read step factor from input header");
  check (step2 = uves_pfits_get_cdelt2(raw_header),
     "Could not read step factor from input header");
  
  /* AMO */
  /* here should be arm dependent */
  check (scan2 = uves_pfits_get_nx(raw_header,chip),
     "Could not read scan1 factor from input header");
  check (scan1 = uves_pfits_get_ny(raw_header,chip),
     "Could not read scan2 factor from input header");
  /*
  check (scan2 = uves_pfits_get_out1nx(raw_header), 
       "Could not read scan1 factor from input header");

  check (scan1 = uves_pfits_get_out1ny(raw_header), 
       "Could not read scan1 factor from input header");
  */
  /* scan1 is allways 1 */
  scan1=1;
  check (naxis1 = uves_pfits_get_naxis1(raw_header),
     "Could not read naxis1 from input header");
  check (naxis2 = uves_pfits_get_naxis2(raw_header),
     "Could not read naxis2 from input header");

  uves_msg_debug("start=%f %f step=%f %f naxis=%d %d",
          start1,start2,step1,step2,naxis1,naxis2);

  imsize1=naxis2;


  /* for some reason in MIDAS prg STEP=START */
  step1=start1;
  step2=start2;


  check( uves_get_parameter(parameters, NULL, recipe_id, 
           "def_pol1", CPL_TYPE_INT, &def_pol1 ),
       "Could not read parameter");

  check( uves_get_parameter(parameters, NULL, recipe_id, 
           "def_pol2", CPL_TYPE_INT, &def_pol2 ),
       "Could not read parameter");

  uves_msg_debug("Polynomial %d %d",def_pol1,def_pol2);


  uves_msg_debug("Display orders positions...");

  inpr1 = start1;
  inpr2 = start1+(imsize1-1)*step1;
 

if (inpr2 < inpr1) { 
  inpr3 = inpr1;
  inpr1 = inpr2;
  inpr2 = inpr3;
}


/* select in the table values which are in the range covered by the
   actual formatcheck frame */
  uves_msg_debug("1st select: inputr=%f %f",inpr1,inpr2);
  uves_msg_debug("ord Xmin=%f Xmax=%f",
          cpl_table_get_column_min(*o_tbl,"X"),
      cpl_table_get_column_max(*o_tbl,"X"));

  check(nraw=SELECT_FUNC_DOUBLE(*o_tbl,"X",CPL_NOT_LESS_THAN,inpr1),"Selection on X failed");

  uves_msg_debug("nraw=%d",nraw);

  check(nraw=SELECT_FUNC_DOUBLE(*o_tbl,"X",CPL_NOT_GREATER_THAN,inpr2),"Selection on X failed");
  uves_msg_debug("nraw=%d",nraw);

  /* scan1=1, scan2=naxis2*/
  inpr1 = start2+(scan1-1)*step2;
  inpr2 = start2+(scan2-1)*step2;

  
  if (inpr2 < inpr1) {
    inpr3 = inpr1;
    inpr1 = inpr2;
    inpr2 = inpr3;
  }

  uves_msg_debug("2nd select: inputr=%f %f",inpr1,inpr2);

  check(nraw=SELECT_FUNC_DOUBLE(*o_tbl,"Y",CPL_NOT_LESS_THAN,inpr1),"Selection on Y failed");
  uves_msg_debug("nraw=%d",nraw);

  check(nraw=SELECT_FUNC_DOUBLE(*o_tbl,"Y",CPL_NOT_GREATER_THAN,inpr2),"Selection on Y failed");
  uves_msg_debug("nraw=%d",nraw);

  /* this table has slightly different values with respect to MIDAS */
  
  if(-1 == uves_physmod_necregr(o_tbl,w_tbl)) {
    uves_msg_error("Error using uves_necregr");
    return -1;
  }

*s_tbl=cpl_table_duplicate(*w_tbl);


/* we sort the table */
  plist=uves_propertylist_new();
  uves_propertylist_append_bool(plist,"RMS",0);  /* 0 for ascending order */
  uves_table_sort(*s_tbl,plist);
  uves_free_propertylist(&plist);

  ord_min=cpl_table_get_column_min(*o_tbl,"ORDER");
  ord_max=cpl_table_get_column_max(*o_tbl,"ORDER");
  ECHORD_1=(ord_max-ord_min+1);
  ORDER=(ord_max-ord_min+1)/2;
  ORDER=uves_max_int(1,ORDER);



  RMSMAX=cpl_table_get_double(*s_tbl,"RMS",ORDER,&status);

  RMSMAX=3.5*RMSMAX;
  uves_msg_debug("RMSMAX=%f",RMSMAX);
 if (RMSMAX < 0.05)  RMSMAX = 0.05;
  uves_msg_debug("RMSMAX=%f",RMSMAX);
 




if( OUTMODE != 'S') {
  uves_msg("Maximum admissible rms: %f pixels",RMSMAX);
}

 
 for (ORDER = 0; ORDER<ECHORD_1; ORDER++){
 RMS = cpl_table_get_double(*w_tbl,"RMS",ORDER,&status);
 /*
 uves_msg_debug("RMS=%f",RMS);
 */
 if (RMS > RMSMAX){ 
   if (OUTMODE != 'S') {
     uves_msg_warning("Rejected order number %d RMS = %f pixels",ORDER,RMS);
   }
   REJECT[NBREJ] = ORDER;
   NBREJ = NBREJ + 1;
 }

}


   for(ORDER = 1; ORDER<NBREJ; ORDER++){
       check(SELECT_FUNC_INT(*o_tbl,"ORDER",
         CPL_NOT_EQUAL_TO,REJECT[NBREJ]),"Error1 selecting ORDER");

   }
  INPUTI_1 = def_pol2 + 1;      /* Degree Y + 1 */
  INPUTI_2 = ECHORD_1 - NBREJ;  /* Number of orders minus rejected orders */
if (INPUTI_1 > INPUTI_2) {
  uves_msg("*****************************************************");
  uves_msg("**** Warning : Number of selected orders {INPUTI(2)}");
  uves_msg("**** is too small for the current value of echelle");
  uves_msg("**** parameter DEFPOL(2)=%d", def_pol2);
  uves_msg("*****************************************************");


}

  check(poly2d_y=uves_polynomial_regression_2d(*o_tbl,"X","ORDER","Y",
                         NULL,def_pol1,def_pol2,"YFIT",
                                             NULL,NULL,
                         &mean_sq_err_y,NULL,NULL,-1,-1),
                                         "Fitting YFIT failed");
  cpl_table_set_column_unit(*o_tbl,"YFIT","pix");
  cpl_table_set_column_unit(*o_tbl,"X","pix");
  cpl_table_set_column_unit(*o_tbl,"Y","pix");

  outputd[0][0]=uves_polynomial_get_coeff_2d(poly2d_y,0,0);
  outputd[1][0]=uves_polynomial_get_coeff_2d(poly2d_y,1,0);
  outputd[0][1]=uves_polynomial_get_coeff_2d(poly2d_y,0,1);
  outputd[1][1]=uves_polynomial_get_coeff_2d(poly2d_y,1,1);
  outputd[2][0]=uves_polynomial_get_coeff_2d(poly2d_y,2,0);
  outputd[0][2]=uves_polynomial_get_coeff_2d(poly2d_y,0,2);
  outputd[1][2]=uves_polynomial_get_coeff_2d(poly2d_y,1,2);
  outputd[2][1]=uves_polynomial_get_coeff_2d(poly2d_y,2,1);
  outputd[2][2]=uves_polynomial_get_coeff_2d(poly2d_y,2,2);

  cpl_table_duplicate_column(*o_tbl,"RESIDUAL",*o_tbl,"Y");
  cpl_table_subtract_columns(*o_tbl,"RESIDUAL","YFIT");
  cpl_table_set_column_unit(*o_tbl,"RESIDUAL","pix");

   uves_msg_debug("Error %f",sqrt(mean_sq_err_y));
  THRESH = kappa * sqrt(mean_sq_err_y);
  if (THRESH > tol)   THRESH = tol;

  /* do a selection on the absolute value of RESIDUAL */
  check(SELECT_FUNC_DOUBLE(*o_tbl,"RESIDUAL",
         CPL_LESS_THAN,THRESH),"Error1 selecting RESIDUAL");

  check(SELECT_FUNC_DOUBLE(*o_tbl,"RESIDUAL",
         CPL_GREATER_THAN,-THRESH),"Error2 selecting RESIDUAL");

  if (num_outliers >= 1)
      {
  
      uves_msg_debug("Kappa-Sigma clipping...");
      for(CNT = 1; CNT<num_outliers; CNT++) {
      
          /*
        REGRE/POLY  {ORDTAB} :Y :X,:ORDER {DEFPOL(1)},{DEFPOL(2)}  KEYLONG 
        SAVE/REGR   {ORDTAB} COEFF   KEYLONG
        COMPUTE/REGR   {ORDTAB} :YFIT = COEFF 
        COMPUTE/TABLE    {ORDTAB} :RESIDUAL = :Y - :YFIT;
          */
          
          uves_polynomial_delete(&poly2d_y);
          check(poly2d_y=uves_polynomial_regression_2d(*o_tbl,"X","ORDER","Y",
                               NULL,def_pol1,def_pol2,"YFIT",
                               NULL,NULL,
                               &mean_sq_err_y,NULL,NULL,-1,-1),
            "Fitting YFIT failed");
          
          if (OUTMODE != 'S') {
          uves_msg_debug("Ndata = %" CPL_SIZE_FORMAT " - Rms = %f pixels",
               cpl_table_get_nrow(*o_tbl),sqrt(mean_sq_err_y));
          }
          
          cpl_table_duplicate_column(*o_tbl,"RESIDUAL",*o_tbl,"Y");
          cpl_table_subtract_columns(*o_tbl,"RESIDUAL","YFIT");
          uves_msg_debug("error %f",sqrt(mean_sq_err_y));
          THRESH = kappa * sqrt(mean_sq_err_y);
          if (THRESH > tol)   THRESH = tol;
          /* do a selection on the absolute value of RESIDUAL */
          check(SELECT_FUNC_DOUBLE(*o_tbl,"RESIDUAL",
          CPL_LESS_THAN,THRESH),"Error3 selecting RESIDUAL");
          check(SELECT_FUNC_DOUBLE(*o_tbl,"RESIDUAL",
          CPL_GREATER_THAN,-THRESH),"Error4 selecting RESIDUAL");
          
      }
      }


  cleanup:
  uves_polynomial_delete(&poly2d_y);
  uves_free_propertylist(&plist);
  return 0;

}
/**@}*/