/*                                                                              *
 *   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: 2013-08-08 13:36:47 $
 * $Revision: 1.103 $
 * $Name: not supported by cvs2svn $
 *
 */

/*-----------------------------------------------------------------------------
                                Includes
 -----------------------------------------------------------------------------*/

#ifdef HAVE_CONFIG_H
#  include <config.h>
#endif
/*----------------------------------------------------------------------------*/
/**
 * @addtogroup uves_wavecal
 */
/*----------------------------------------------------------------------------*/

#include <uves_wavecal_body.h>

/* Definitions */
#include <uves.h>

/* Macro steps */
#include <uves_extract.h>
#include <uves_flatfield.h>
#include <uves_wavecal_search.h>
#include <uves_wavecal_firstsolution.h>
#include <uves_wavecal_identify.h>
#include <uves_rebin.h>
#include <uves_merge.h>

/* Utility functions */
#include <uves_wavecal_utils.h>
#include <uves_utils.h>
#include <uves_utils_wrappers.h>
#include <uves_plot.h>
#include <uves_parameters.h>
#include <uves_dfs.h>
#include <uves_pfits.h>
#include <uves_qclog.h>
#include <uves_recipe.h>
#include <uves_error.h>
#include <flames_reduce_vcorrel.h>

#include <cpl.h>

#include <stdbool.h>
#include <float.h>
#include <string.h>
/*-----------------------------------------------------------------------------
                                Defines
 -----------------------------------------------------------------------------*/
/* threshold values for maximum pixel saturation */
#define DRS_PTHRES_MAX 55000
#define DRS_PTHRES_MIN -20
#define DRS_CVEL_MIN -6.
#define DRS_CVEL_MAX +6.
#define N_FIBRES_MAX 9
/*-----------------------------------------------------------------------------
                            Functions prototypes
 -----------------------------------------------------------------------------*/
static void uves_wavecal_qclog(const cpl_table* table,
                               int firstabs,
                               int lastabs,
                   const cpl_image *arclamp,
                   const uves_propertylist* raw_header,
                   bool flames,
                   int trace_number,
                   int fibre_mask,
                   double offset,
                   enum uves_chip chip,
                   cpl_table* qclog);

static void uves_wavecal_qclog_intmon(cpl_table* table,
                      const cpl_table *line_intmon,
                      const uves_propertylist* raw_header,
                      bool flames,
                      int fibre,
                      enum uves_chip chip,
                      cpl_table* qclog);

/*-----------------------------------------------------------------------------
                            Recipe standard code
 -----------------------------------------------------------------------------*/

const char * const uves_wavecal_desc_short = "Performs the wavelength calibration";

const char * const uves_wavecal_desc =
"The recipe performs a wavelength calibration for each extraction window.\n"
"Conceptually, each chip contains a number of order lines, each of which\n"
"contains a number of fibre traces, each of which contains a number of\n"
"extraction windows. For UVES data, there is only one trace per order and\n"
"three extraction windows (sky, object, sky). For FLAMES/UVES data there\n"
"are multiple traces per order but only one extraction window per trace.\n"
"The number of traces is defined in the order table while the geometry of\n"
"the extraction windows is specified by recipe parameters (see below).\n"
"\n"
"Expected input for this recipe is an arc lamp frame, ARC_LAMP_xxx or\n"
"ECH_ARC_LAMP_xxx (where xxx=BLUE, RED), order table(s) for each chip,\n"
"ORDER_TABLE_xxxx (where xxxx=BLUE, REDL, REDU), 'guess' line table(s)\n"
"for each chip, LINE_TABLE_xxxx, a wavelength catalogue table, \n"
"LINE_REFER_TABLE, and optionally a wavelength table of bright lines,\n"
"LINE_INTMON_TABLE, used only for computing Quality Control parameters.\n"
"\n"
"The output line table(s), LINE_TABLE_xxxx, contains the columns\n"
"X            : Horizontal position (from Gaussian fit) of detected line\n"
"dX           : Uncertainty (one sigma) of X\n"
"Ynew         : Vertical position of detected line\n"
"XWidth       : Width (in pixels) of detected line from Gaussian fit\n"
"Peak         : Intensity of detected line\n"
"Background   : Fitted background (ADU) of detected line\n"
"Slope        : Linear background slope (ADU/pixel) of detected line\n"
"               from Gaussian fit\n"
"Intensity    : Intensity of detected line scaled to unit exposure\n"
"               time. (This column only present if a LINE_INTMON_TABLE\n"
"               is provided.)\n"
"Order        : Absolute order number of detected line\n"
"Y            : Relative order number of detected line\n"
"               (it's not a very descriptive column name)\n"
"WaveC        : Wavelength of this line (computed using the resulting\n"
"               dispersion relation)\n"
"dLambdaC     : Uncertainty (one sigma) of 'WaveC'.\n"
"Pixel        : The width in w.l.u. of a pixel (computed locally).\n"
"Residual     : Residual (in w.l.u.) of this line\n"
"Residual_pix : Residual (in pixels) of this line\n"
"Lambda_candidate : Nearest line in catalogue\n"
"dLambda_cat_sq   : Squared distance to nearest catalogue line\n"
"dLambda_nn_sq    : Squared distance to nearest neighbour multiplied by ALPHA\n"
"Ident        : The wavelength associated with this emission line,\n"
"               or invalid if this line was not identified\n"
"dIdent       : Uncertainty of catalogue wavelength\n"
"Select       : 1 if the line was identified, 0 otherwise\n"
"NLinSol      : 1 if the line was identified and accepted for the\n"
"               polynomial fit, 0 otherwise\n"
"Intensity    : Intensity of detected line scaled to unit exposure\n"
"               time. (This column is present only if a LINE_INTMON_TABLE\n"
"               is provided.)\n"
"\n"
"The 2nd table extension contains the dispersion relation (a 2d polynomial).\n"
"The 3rd table extension contains the map from (pixel, pixel)-space to\n"
" physical order numbers (used internally by the calibration recipe; \n"
"another 2d polynomial).\n"
"\n"
"If there is more than one extraction window, the results of each calibration\n"
"is stored in subsequent table extensions of the same FITS file. For \n"
"example, extensions 4, 5 and 6 would contain the resulting line table \n"
"(and its two associated polynomials) for the second extraction window. \n"
"The results for the calibration of the n'th extraction window is stored \n"
"in extensions (3*n - 2) to 3*n.\n";
/**@{*/

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

/*----------------------------------------------------------------------------*/
/**
  @brief    Setup the recipe options    
  @param    parameters        the parameterlist to fill
  @param    recipe_id         name of calling recipe
  @param    slit              slit length in pixels
  @return   0 if everything is ok
 */
/*----------------------------------------------------------------------------*/
int
uves_wavecal_define_parameters_body(cpl_parameterlist *parameters, 
                    const char *recipe_id, double slit)
{
    const char *subcontext;

    /*****************
     *    General    *
     *****************/
    
    if (uves_define_global_parameters(parameters) != CPL_ERROR_NONE)
    {
        return -1;
    }
    
    /*****************
     *  Extraction   *
     *****************/

    subcontext = NULL;

    /* nwindows */
    uves_par_new_range("nwindows",
               CPL_TYPE_INT,
               "Number of extraction windows per trace. "
               "The windows will be aligned (i.e. no overlap "
               "and no spacing between adjacent windows). "
               "Unless an offset is specified, the middle "
               "window(s) is centered on the trace",
               strcmp(recipe_id, "flames_cal_wavecal") == 0 ?
            1 : 3,      /* FLAMES: 1; UVES: 3 */
            1, INT_MAX);

    /* length */
    uves_par_new_range("length",
               CPL_TYPE_DOUBLE,
               "Length (in pixels) of each extraction window. "
               "This parameter is also equal to the seperation of "
               "adjacent window centers, causing the extraction "
               "windows to always be aligned. The parameter is "
               "automatically adjusted according to the binning "
               "of the input raw frame. If negative, the extraction "
                       "window length is determined automatically "
                       "to cover the full slit",
               slit, -2.0, DBL_MAX);

    /* offset */
    uves_par_new_range("offset",
               CPL_TYPE_DOUBLE,
               "A global offset (in pixels) of all extraction windows",
               0.0, -25., 25.);
    
    /* method */
    if (uves_propagate_parameters_step(UVES_EXTRACT_ID, parameters, 
                       recipe_id, NULL) != 0)
    {
        return -1;
    }
    
    /* Override default optimal extraction profile. Assume constant profile */
#if 0
    /* this should perhaps be enabled but doesn't work properly in the moment.

        ChangeLog:
    uves_cal_wavecal: The arc lamp spectrum is now extracted using
    average extraction and weighting each pixel with its inverse
    variance. This is equivalent to doing an optimal extraction under
    the assumption of a constant spatial profile, and is implemented
    as such. This was a necessary change in order to be robust against
    interorder noisy pixels caused by dividing by the flat-field.

       */
    {
    const char *profile = "constant";
        double kappa = -1;
    
    if (uves_set_parameter_default(parameters, 
                       recipe_id, UVES_EXTRACT_ID ".profile",
                       CPL_TYPE_STRING, &profile)
        != CPL_ERROR_NONE)
        {
        return -1;
        }

        /* Disable cosmic ray rejection, it does not work well for
           this particular profile and very high S/N */
    if (uves_set_parameter_default(parameters,
                       recipe_id, UVES_EXTRACT_ID ".kappa",
                       CPL_TYPE_DOUBLE, &kappa)
        != CPL_ERROR_NONE)
        {
        return -1;
        }
    }
#else
    {
    const char *method = "average";
    
    if (uves_set_parameter_default(parameters, 
                       recipe_id, UVES_EXTRACT_ID ".method",
                       CPL_TYPE_STRING, &method)
        != CPL_ERROR_NONE)
        {
        return -1;
        }
        
    }
#endif
    
    /*****************
     *    Search     *
     *****************/

    subcontext = "search";

    /* range */
    uves_par_new_range("range",
               CPL_TYPE_INT,
               "Width (pix) of search window is 2*range + 1. "
               "This parameter is automatically adjusted "
               "according to binning.",
               8, 1, INT_MAX);

    /* minlines */
    uves_par_new_range("minlines",
               CPL_TYPE_INT,
               "Minimum number of lines to detect. If zero, "
               "the default value (1100 for BLUE/REDL chips; "
               "1000 for REDU chip) is used.",
               0, 0, INT_MAX);

    /* maxlines */
    uves_par_new_range("maxlines",
               CPL_TYPE_INT,
               "Maximum number of lines to detect. If zero, "
               "the default value (1600 for BLUE/REDL chip; "
               "1400 for REDU chip) is used.",
               0, 0, INT_MAX);
    
    /* centeringmethod */
/*  Temporally removed as 'gravity' do not work and it does not make 
    sense a parameter with only one option

    uves_par_new_enum("centeringmethod",
              CPL_TYPE_STRING,
              "Line centering method",
              "gaussian",             // Default 
              1,                      // Number of options 
                      "gaussian");    // List of options 
     
*/
         /* old setting allowed gravity "gaussian", "gravity"); */

    /*******************
     *  First solution *
     *******************/

    subcontext = "first";

    /* shiftmax */
    uves_par_new_range("shiftmax",
               CPL_TYPE_DOUBLE,
               "The maximum shift (pix) in either direction compared to "
               "guess solution. This parameter is automatically "
               "corrected for binning",
               10.0, 0.0, DBL_MAX);

    /* shiftstep */
    uves_par_new_range("shiftstep",
               CPL_TYPE_DOUBLE,
               "The step size (pix) used when searching "
               "for the optimum shift. This parameter is "
               "automatically corrected for binning",
               0.1, 0.0, DBL_MAX);

    /* shifttoler */
    uves_par_new_range("shifttoler",
               CPL_TYPE_DOUBLE,
               "Tolerance (pix) when matching shifted lines. "
               "This parameter is not adjusted according to binning",
               0.05, 0.0, DBL_MAX);


    /*****************
     *  Identify     *
     *****************/

    subcontext = "identify";
    
    /* alpha */
    uves_par_new_range("alpha",
               CPL_TYPE_DOUBLE,
               "The parameter that controls the distance to the "
               "nearest neighbours",
               0.1, 0.0, 1.0);

    /* maxerror */
    uves_par_new_range("maxerror",
               CPL_TYPE_DOUBLE,
               "This parameter controls the graceful exit of "
               "the identification loop. If the RMS of the "
               "global fit exceeds this value (pix) the "
               "iteration stops",
               20.0, 0.0, DBL_MAX);
    
    /* degree */
    uves_par_new_range("degree",
               CPL_TYPE_INT,
               "Degrees of the global 2d dispersion polynomial. If "
                       "a negative number is specified, the polynomial "
                       "degrees are automatically selected by starting from "
                       "(1, 1) and inreasing the degrees as long as the RMS "
                       "residual decreases significantly",
               4, -2, INT_MAX);


    /*****************
     * Calibration   *
     *****************/

    subcontext = "calibrate";

    /* tolerance */
    uves_par_new_value("tolerance",
                       CPL_TYPE_DOUBLE,
                       "Tolerance of fit. If positive, the tolerance "
                       "is in pixel units. If negative, abs(tolerance) "
                       "is in wavelength units. Lines with residuals "
                       "worse than the tolerance are excluded from the "
                       "final fit. Unlike in previous versions, this "
                       "parameter is not corrected for CCD binning. "
                       "This rejection based on the absolute residual in "
                       "pixel can be effectively disabled by setting the "
                       "tolerance to a very large number (e.g. 9999). In "
                       "that case outliers will be rejected using only "
                       "kappa sigma clipping.",
                       0.6); 
/*               0.07); */
/*                       9999.0);*/

    
    /* kappa */
    uves_par_new_range("kappa",
               CPL_TYPE_DOUBLE,
               "Lines with residuals more then kappa stdev "
                       "are rejected from the final fit",
                       4.0,0.,100.);

    /***************
     * Rebinning   *
     ***************/
    if (uves_propagate_parameters_step(UVES_REBIN_ID, parameters, 
                       recipe_id, NULL) != 0)
    {
        return -1;
    }
    
    return (cpl_error_get_code() != CPL_ERROR_NONE);
}

/*----------------------------------------------------------------------------*/
/**
   @brief    Calibrate one extraction window
   @param    arclamp           arc lamp image, possibly debiased and flat-fielded
   @param    arclamp_noise     error bars
   @param    rotated_header    Header describing the geometry of the arc lamp
                               image
   @param    ordertable        Order table describing the order locations on 
                               the arc lamp image
   @param    order_locations   The polynomial describing the order locations
   @param    flat_fielded      was the arc lamp image flat-fielded?
   @param    weights           if non-NULL, weighted extraction is used
   @param    guess             A line table containing line locations and
                               their associated wavelengths. This 'guess' solution 
                               is used as a starting point for the calibration.
   @param    line_refer        Wavelength catalogue of wavelenghts and their
                               uncertainties
   @param    flames            FLAMES reduction?
   @param    tab_in_out_oshift order number shift
   @param    tab_in_out_yshift y-position shift
   @param    chip              CCD chip
   @param    bin_disp          binning in dispersion direction
   @param    trace             Fibre trace ID to calibrate (FLAMES/UVES)
   @param    window            Window number to calibrate
   @param    debug_mode             If set to true, intermediate results are saved
                               to the current directory
   @param    offset            Extraction offset wrt spatial center of order
   @param    slitlength        EXtraction slit length
   @param    parameters        The recipe parameter list
   @param    RANGE             See recipe man page or @c uves_wavecal_create()
   @param    MINLINES          See recipe man page or @c uves_wavecal_create()
   @param    MAXLINES          See recipe man page or @c uves_wavecal_create()
   @param    CENTERING_METHOD  See recipe man page or @c uves_wavecal_create()
   @param    SHIFT_MAX         See recipe man page or @c uves_wavecal_create()
   @param    SHIFT_STEP        See recipe man page or @c uves_wavecal_create()
   @param    SHIFT_TOLERANCE   See recipe man page or @c uves_wavecal_create()
   @param    ALPHA             See recipe man page or @c uves_wavecal_create()
   @param    MAXERROR          See recipe man page or @c uves_wavecal_create()
   @param    DEGREE            See recipe man page or @c uves_wavecal_create()
   @param    TOLERANCE         See recipe man page or @c uves_wavecal_create()
   @param    corvel            Table for velocity correlation
   @param    dispersion_relation  (output) Result of the calibration ('returned')
   @param    absolute_order       (input/output) The map from (pixel, pixel)-space
                                   to absolute order numbering.
                                   This parameter is read (and therefore must
                   contain a 'guess' map) and also updated
                   after the calibration is done.
   @param    first_absolute_order   (output) The absolute order number of the
                                    first relative order
   @param    last_absolute_order    (output) The absolute order number of the
                                    last relative order
   @return   The line table for this trace/window containing the detected
             lines and their associated wavelengths.

   This function executes the macro-steps involved in processing a single 
   extraction window. These are
   - Extract the spectrum using the extraction window defined by the
   arguments @em offset, @em slitlength using the extraction method specified
   by @em EXTRACT_METHOD. (See @c uves_extract() .)
   - Detect emission lines in the extracted spectrum.
   (See @c uves_wavecal_search() .)
   - Get the absolute (physical) numbering of the orders by using the
   provided map @em absolute_order from (pixel, pixel)-space to physical 
   order number. (See @c uves_wavecal_firstsolution() .)
   - Obtain a first solution by estimating the @em x-shift of the detected
   lines with respect to the locations of the lines in the provided
   @em 'guess' line table, then apply this shift to the 'guess' 
   dispersion relation. (See @c uves_wavecal_firstsolution() .)
   - Iteratively identify emission lines and update the dispersion relation, 
   until no new identifications can be made. (See @c uves_wavecal_identify() .)
*/
/*----------------------------------------------------------------------------*/
static cpl_table *
uves_wavecal_process_window(const cpl_image *arclamp, 
                            const cpl_image *arclamp_noise,
                            const uves_propertylist *rotated_header,
                            const cpl_table *ordertable, 
                            const polynomial *order_locations,
                            bool flat_fielded,
                            cpl_image *weights,
                            /* const cpl_table *drs_table,  Not used */
                            const cpl_table *guess,
                            const cpl_table *line_refer,
                            bool flames,
                            int tab_in_out_oshift,
                            double tab_in_out_yshift,
                            enum uves_chip chip, int bin_disp,
                            int trace, int window, int NWINDOWS,
                            /* General */
                            bool debug_mode,
                            /* Extraction */
                            double offset,
                            double slitlength,
                            const cpl_parameterlist *parameters,
                            const char *recipe_id,
                            /* Search */
                            int RANGE,
                            int MINLINES,
                            int MAXLINES,
                            centering_method CENTERING_METHOD,
                            /* First solution */
                            double SHIFT_MAX,
                            double SHIFT_STEP,
                            double SHIFT_TOLERANCE,
                            /* Identify */
                            double ALPHA,
                            double MAXERROR,
                            int DEGREE,
                            /* Calibrate */
                            double TOLERANCE,
                            double kappa,
                            cpl_frame* corvel_frm,
                            cpl_table** flames_qclog,
                            /* Result */
                            polynomial **dispersion_relation,
                            polynomial **absolute_order,
                            int *first_absolute_order,
                            int *last_absolute_order)
{
    cpl_table        *linetable          = NULL;  /* Result */
    cpl_table        *temp               = NULL;
    cpl_image        *spectrum           = NULL;
    cpl_image        *spectrum_noise     = NULL;
    cpl_image        *debug_image        = NULL;
    polynomial       *initial_dispersion = NULL;
    int              *relative_order     = NULL;  /* Map from physical 
                             order to relative order */
    uves_propertylist *spectrum_header    = NULL;

    cpl_image        *rebinned           = NULL;  /* Used for calculating 
                             the instrument resolution */
    cpl_image        *rebinned_noise     = NULL;  /* Used for calculating 
                             the instrument resolution */
    uves_propertylist *rebinned_header    = NULL;
    cpl_image        *merged             = NULL;
    cpl_image        *merged_noise       = NULL;
    uves_propertylist *merged_header      = NULL;
    cpl_table        *info_tbl           = NULL;

    /* Needed for optimal extraction */
    cpl_image *weights_opt = NULL;
    cpl_table *cr_table = NULL;
    cpl_image *cr_image = NULL;
    cpl_table *order_trace = NULL;

    merge_method m_method = flat_fielded ? MERGE_OPTIMAL : MERGE_SUM;

    /* Extract the spectrum */
    uves_free_table(&info_tbl);


    check( spectrum = uves_extract((cpl_image *)arclamp,/* Const-casts are okay,
                                                           the image (bpm) + error bars
                                                           is changed
                                                           only in optimal extraction */
                                   (cpl_image *)arclamp_noise,
                                   NULL,            /* Header (optimal only)    */
                                   ordertable, 
                                   order_locations,
                                   slitlength,      /* Slit length (pixels)     */
                                   offset,          /* Slit center offset       */
                                   parameters,      /* Extraction method        */
                                   recipe_id, 
                                   "",     
                                   true,            /* Extraction partial bins? */
                                   debug_mode,
                                   chip,
                                   &spectrum_header,
                                   &spectrum_noise,
                                   NULL, 
                                   NULL,            /* Optimal:  sky+noise      */
                                   &cr_table, 
                                   &cr_image, 
                                   NULL, 
                                   (weights != NULL) ? &weights : &weights_opt,
                                   &info_tbl,
                                   &order_trace),
           "Error extracting spectrum");
    uves_free_table(&info_tbl);


    check(uves_propertylist_copy_property_regexp(spectrum_header,
                                                 rotated_header,
						 "^ESO ", 0),"error copying hierarch keys");

    /* Set bad pixels to 0, so that the search algorithm doesn't
       fail because of bad pixels (but simply does not detect anything)
    */
    cpl_image_fill_rejected(spectrum, 0);
    cpl_image_accept_all(spectrum);
    cpl_image_fill_rejected(spectrum_noise, 1);
    cpl_image_accept_all(spectrum_noise);

    /* Save spectrum + noise */
    if (debug_mode)
    {


      check(uves_propertylist_copy_property_regexp(spectrum_header, rotated_header,
		  "^ESO ", 0),
       "Error copying hieararch keys");

        check( uves_save_image_local("Extracted spectrum", "spectrum", 
                     spectrum, chip, trace, window, spectrum_header, true),
           "Error saving spectrum");

        check( uves_save_image_local("Extracted spectrum noise", "spectrum_noise",
                     spectrum_noise, chip, trace, window, spectrum_header, true),
           "Error saving spectrum");
    }

    /* Locate lines */
    debug_image = cpl_image_duplicate(arclamp);
    check( linetable = uves_wavecal_search(spectrum, 
                                           spectrum_noise,
                                           spectrum_header,
                                           flat_fielded,
                                           order_locations,
                                           debug_image,
                                           RANGE,
                                           MINLINES, 
                                           MAXLINES, 
                                           CENTERING_METHOD,
                                           bin_disp,trace,window,
                                           flames_qclog[0]),
       "Line search failed");

     /* Read first solution from guess line table */
    {
        int degree = 5; /* for the initial solution only. For the
                           final solution the degree given as recipe
                           parameter is used */
        uves_polynomial_delete(&initial_dispersion);
        cpl_free(relative_order);
        check( initial_dispersion = uves_wavecal_firstsolution(linetable,
                                                               guess, 
                                                               absolute_order,
                                                               ordertable,
                                                               order_locations,
                                                               flames,
                                                               offset,
                                                               &relative_order,
                                                               degree,
                                                               SHIFT_MAX, 
                                                               SHIFT_STEP, 
                                                               SHIFT_TOLERANCE, 
                                                               MAXERROR,
                                                               first_absolute_order,
                                                               last_absolute_order),
               "Could not get first solution");
    }


    if (flames)
    {
        /* !AM: To correct eventual residual shifts between Guess and Final order
           ! (and line) table
           compute/table {LINTAB} :YNEW =
           :YNEW - {{ORDTAB},TAB_IN_OUT_YSHIFT} - {{OLDORDTAB},FIBREPOS({i1})}
           compute/table {LINTAB} :Y = :Y +{{ORDTAB},TAB_IN_OUT_OSHIFT}
        */

        cpl_table_add_scalar(linetable, "Y", tab_in_out_oshift);
        cpl_table_add_scalar(linetable, "Ynew", - tab_in_out_yshift - offset);
    }

    /* Calibrate */
    check( *dispersion_relation = uves_wavecal_identify(linetable, /* Guess solution */
                                                        line_refer,
                                                        initial_dispersion,
                                                        DEGREE,
                                                        TOLERANCE, ALPHA, 
                                                        MAXERROR,
                                                        kappa,trace,window,flames_qclog[0]),
       "Could not calibrate orders");

    if (flames)
    {
        /* AM: To have correct split of fibres in line table:
           > compute/table {LINTAB} :YNEW = :YNEW + {{ORDTAB},TAB_IN_OUT_YSHIFT}
                     + {{OLDORDTAB},FIBREPOS({i1})}
        */

        cpl_table_add_scalar(linetable, "Ynew", + tab_in_out_yshift + offset);
    }

    /* UVES: make plots (resolution + etc.) for the central window,
     * FLAMES: all fibres
     */
    if (flames || (trace == 0 && window == 2)|| (window == 1 && NWINDOWS == 1))
    {
        /* Create x-FWHM column:   FWHM = 2.3548 sigma */
        check(( cpl_table_duplicate_column(linetable, "deltaX", linetable, "Xwidth"),
            cpl_table_multiply_scalar (linetable, "deltaX", TWOSQRT2LN2)),
          "Error creating FWHM column");
        

        check_nomsg( temp = uves_extract_table_rows(
                             linetable, "NLinSol", CPL_NOT_EQUAL_TO, 0) );

        check( uves_plot_table(temp, "Order", LINETAB_RESIDUAL, "Residual of fit"), 
           "Plotting failed");
        
        check( uves_plot_table(temp, "X", "deltaX", "line FWHM (mean = %.2f pixels)",
                   cpl_table_get_column_mean(linetable, "deltaX")),
           "Plotting failed");
        /*
          check( uves_plot_table(linetable, "Y", "deltaX", 
          "line FWHM (mean FWHM = %.2f pixels)",
          cpl_table_get_column_mean(linetable, "deltaX")), "Plotting failed");
        */
        
        
        /* Compute resolution as as lambda / deltalambda   where deltalambda
               is the peak FWHM in wavelength space (after resampling using 
           WAVESTEP = average pixelsize) 
        */
        
        {
        /* Rebin using steps of median pixelsize */
        double wavestep;
        double lambda_start = 0;
        int n_traces = 1;      /* We didn't do a 2d extraction; 
                      there's only 1 trace
                      per order */
        int i, nbins;
                bool threshold_to_positive = true;

        cpl_table_new_column(linetable, "deltaLambda", CPL_TYPE_DOUBLE);
        
        check( rebinned_noise = uves_rebin(spectrum_noise,
                                           parameters, 
                                           recipe_id,
                                           linetable, 
                                           *dispersion_relation, 
                                           *first_absolute_order, 
                                           *last_absolute_order,
                                           n_traces,
                                           threshold_to_positive,
                                           true,
                                           &rebinned_header),
               "Could not rebin noise of arc lamp spectrum");

        threshold_to_positive = false;
        uves_free_propertylist(&rebinned_header);
        check( rebinned = uves_rebin(spectrum,
                                     parameters, 
                                     recipe_id,
                                     linetable, 
                                     *dispersion_relation, 
                                     *first_absolute_order, 
                                     *last_absolute_order,
                                     n_traces,
                                     threshold_to_positive,
                                     false,
                                     &rebinned_header),
               "Could not rebin arc lamp spectrum");

        /* Save order-by-order rebinned spectrum+noise */
        if (debug_mode)
            {
            check( uves_save_image_local("Rebinned spectrum",
                             "wxb", rebinned, chip, 
                             trace, window, rebinned_header, true),
                   "Error saving rebinned spectrum");

            check( uves_save_image_local("Noise of rebinned spectrum",
                             "errwxb", rebinned_noise, chip, 
                             trace, window, rebinned_header, true),
                   "Error saving noise of rebinned spectrum");
            }
        
        check( merged = uves_merge_orders(rebinned, 
					  rebinned_noise,
					  rebinned_header,
					  m_method,
					  n_traces,
					  &merged_header,
					  0,0,chip,
					  &merged_noise),
               "Could not merge arc lamp spectrum");

        check( uves_plot_image_rows(merged, 1, 1, 1, 
                        "Wavelength (arbitrary units)", 
                        "Flux", "Resampled arc lamp spectrum"),
               "Plotting failed");

        /* Save merged arc lamp spectrum */
        if (debug_mode)
            {
            check( uves_save_image_local("Rebinned, merged spectrum",
                             "merged", merged, chip, 
                             trace, window, merged_header, true),
                   "Error saving merged spectrum");
            }
        
        nbins = cpl_image_get_size_x(merged);

        check( wavestep = uves_pfits_get_cdelt1(merged_header),
               "Error reading resampling step size");
        
        check( lambda_start = uves_pfits_get_crval1(merged_header), 
               "Could not read start wavelength of merged spectrum");


	//Begin commented region
	
   if (flames && trace == 0 && corvel_frm != NULL)
    {
       //The following (flames_reduce.VCORREL) calculates 
       //a cross correlation and does some QC


//> !To support simcal Mode
//> if i1 .eq. 1 then
//>    w/o "DRS_CVEL_SWITCH={DRS_CVEL_SWITCH}"
//>    if "{DRS_CVEL_SWITCH}" .eq. "Y" then
//>       w/o "To support simcal Mode"
//>       define/local ord_min/i/1/1 0
//>       define/local ord_max/i/1/1 0
//>       define/local rsample/d/1/1 0
//>       statistic/table {ORDTAB} :ORDER >Null
//>       ord_min = outputr(1)
//>       ord_max = outputr(2)
//>       rsample = {{LINTAB},PIXEL(1)}
//>       rsample = 2./3. * rsample
//>       rebin/echelle  {ofrm} w{ofrm} {rsample} NONL {LINTAB} {SESSOUTV}
//>       mercut/echelle w{ofrm} mw{ofrm} {ord_min},{ord_max} NOAPPEND
//>       !corvel stuff
//>       define/local OLD_CVEL_MAX/d/1/1 {DRS_CVEL_MAX}
//>       define/local OLD_CVEL_MIN/d/1/1 {DRS_CVEL_MIN}
//>       @p flames_reduce,VCORREL x1_rbf_ cvel1 0 {ord_max} {parCorVelTab} _0 _{chip({PATHID})} 0
//>       DRS_CVEL_MAX = DRS_CVEL_MAX + {q1}
//>       DRS_CVEL_MIN = DRS_CVEL_MIN + {q1}
//>
//>       @p flames_reduce,VCORREL x1_rbf_ cvel2 0 {ord_max} {parCorVelTab} _0 _{chip({PATHID})} 0
//>       cvel_0 = {q1}
//>       DRS_CVEL_MAX = OLD_CVEL_MAX
//>       DRS_CVEL_MIN = OLD_CVEL_MIN
//>    endif
//>    write/keyword DRS_CVEL_SWITCH Y
//> endif




       


      const char* drs_base_name=NULL;
      const char* prefid=NULL;
      double ccf_posmax_zero_point=0;
      double ccf_posmax_zero_point_iter0=0;
      double cvel_max=0;
      double cvel_sig=0;
      if(chip == UVES_CHIP_REDL) {
	prefid="l";
      } else {
	prefid="u";
      }
      const char* name=NULL;
      const char* file=NULL;
      cpl_propertylist* plist=NULL;
      double drs_cvel_min=0;
      double drs_cvel_max=0;

 
      check( uves_save_image_local("Rebinned spectrum",
				   "wxb", rebinned, chip, 
				   trace, window, rebinned_header, true),
	     "Error saving rebinned spectrum");

      check( uves_save_image_local("Rebinned, merged spectrum",
				   "mwxb", merged, chip, 
				   trace, window, merged_header, true),
	     "Error saving merged spectrum");
     

      check( file = uves_local_filename("wxb", chip, trace, window),
	     "Error getting filename");
      check_nomsg(plist=cpl_propertylist_load(file,0));

      name=uves_sprintf("wfxb_%s%s%4.4d%s",prefid,"_raw",1,".fits");
      drs_base_name=uves_sprintf("fxb_%s",prefid);

      cpl_image_save(rebinned,name, CPL_BPP_IEEE_FLOAT,plist, 
		     CPL_IO_DEFAULT);


      name=uves_sprintf("mwfxb_%s%s%4.4d%s",prefid,"_raw",1,".fits");
      cpl_image_save(merged,name, CPL_BPP_IEEE_FLOAT,plist, 
		     CPL_IO_DEFAULT);
     
      cpl_propertylist_delete(plist);


      int ord_max=(*first_absolute_order-*last_absolute_order)+1;
     
      uves_msg("cvel max:%g %g",DRS_CVEL_MAX,DRS_CVEL_MIN);
      drs_cvel_max =DRS_CVEL_MAX;
      drs_cvel_min =DRS_CVEL_MIN;

      check_nomsg(flames_reduce_vcorrel(drs_base_name,
					"cvel2", 
					prefid,
					ord_max,
					corvel_frm, 
					"_raw0001", 
					"_raw0001", 
                                        DRS_CVEL_MIN,
                                        DRS_CVEL_MAX,
					&ccf_posmax_zero_point,
					&cvel_max,
					&cvel_sig,
					flames_qclog[0]));
     
     drs_cvel_max +=cvel_max;
     drs_cvel_min +=cvel_max;
     ccf_posmax_zero_point_iter0 =cvel_max;
     uves_msg("cvel max:%g %g",drs_cvel_max,drs_cvel_min);

     check_nomsg(flames_reduce_vcorrel(drs_base_name,
					"cvel2", 
					prefid,
					ord_max,
					corvel_frm, 
					"_raw0001", 
					"_raw0001", 
                                        drs_cvel_min,
                                        drs_cvel_max,
					&ccf_posmax_zero_point,
					&cvel_max,
					&cvel_sig,
					flames_qclog[1]));
      
      drs_cvel_max +=cvel_max;
      drs_cvel_min +=cvel_max;
      ccf_posmax_zero_point =ccf_posmax_zero_point_iter0;
      ck0_nomsg(uves_qclog_add_double(flames_qclog[1],
				   "QC CCF POSOFF",
				   ccf_posmax_zero_point,
				   "CCF pos avg from ThAr calibration",
				   "%f"));


      uves_msg("cvel max:%g min: %g zp: %g",
	       drs_cvel_max,drs_cvel_min,ccf_posmax_zero_point);

   
    }
   //End commented region


        /* For all identified lines fit the line width in the
           merged spectrum to get the resolution  */
         for (i = 0; i < cpl_table_get_nrow(linetable); i++)
            {
            double lambda = cpl_table_get_double(
                linetable, LINETAB_LAMBDAC, i, NULL);
            double width  = 
                cpl_table_get_double(linetable, "Xwidth"   , i, NULL) *
                fabs(cpl_table_get_double(linetable, LINETAB_PIXELSIZE, i, NULL)); 
            /* in wlu */
            
            /* Convert line wavelength and width to 'bin' units */
            int bin       = 1 + 
                uves_round_double((lambda - lambda_start) / wavestep);
            double width_bin = width / wavestep;

            /* Set fitting window to +-5 sigma */
            int first_bin = uves_max_int(    1, uves_round_double(bin - 5*width_bin));
            int last_bin  = uves_min_int(nbins, uves_round_double(bin + 5*width_bin));

            double my, sigma, norm, background;    /* Results of fit */
            double lambda_fwhm;
            
            if (cpl_table_is_valid(linetable, "Ident", i) && first_bin < last_bin)
                {
                /* Fit a gaussian to the merged arc spectrum */
                uves_fit_1d_image(merged, 
#if 1
                                                  merged_noise, 
#else /* Unweighted fit like MIDAS which gives larger sigma */
                                                  NULL,
#endif
                                                  NULL,
                          true,     /* Horizontal?    */
                          false,    /* Fix background?*/
                          false,    /* Fit background?*/
                          first_bin, 
                                                  last_bin, 
                                                  1,  /* xlo, xhi, y    */
                          &my, 
                                                  &sigma, 
                                                  &norm, 
                                                  &background, NULL, /* slope */
                          NULL, 
                                                  NULL, 
                                                  NULL,    /* mse, red_chisq,
                                  covariance     */
                          uves_gauss, 
                                                  uves_gauss_derivative, 
                                                  4);
                                        
                if (cpl_error_get_code() == CPL_ERROR_CONTINUE)
                    {
                    uves_error_reset();
                    uves_msg_debug("Gaussian fitting failed "
                               "at lambda = %f wlu, bins = "
                               "%d - %d, ignoring line",
                               lambda, first_bin, last_bin);
                    
                    cpl_table_set_invalid(linetable, "deltaLambda", i);
                    
                    }
                else
                    {
                        assure(cpl_error_get_code() == CPL_ERROR_NONE,
                               cpl_error_get_code(), "Gaussian fitting failed");
                    
                    /* Convert from bins to wavelength */
                    lambda_fwhm = TWOSQRT2LN2 * sigma * wavestep;
                    
                    cpl_table_set_double(linetable, "deltaLambda",
                                 i, lambda_fwhm);

                    }
                }
            else
                {
                cpl_table_set_invalid(linetable, "deltaLambda", i);
                }
            }
        
        
        /* Create column 'Resol'(ution) = lambda / deltalambda */
        check(( cpl_table_duplicate_column(linetable, "Resol", 
                           linetable, LINETAB_LAMBDAC),
            cpl_table_divide_columns  (linetable, "Resol",
                           "deltaLambda")),
              "Error creating 'Resol' column");
        
        /* Filter out extreme outliers (due to bad gaussian fit) */
        {
         
          int ninvalid=0;
          int nrows=0;
           double resol_avg   = 0;
           double resol_stdev = 0;
           double kappar = 10.0;
           nrows=cpl_table_get_nrow(linetable);
           ninvalid=cpl_table_count_invalid(linetable,"Resol");
           assure(ninvalid < nrows,CPL_ERROR_ILLEGAL_INPUT, 
                 "No valid elements in Resol column. "
                 "You must decrease parameter rebin.wavestep");
           check_nomsg(resol_avg=cpl_table_get_column_median(linetable, "Resol"));
           check_nomsg(resol_stdev=cpl_table_get_column_stdev (linetable, "Resol"));

            for (i = 0; i < cpl_table_get_nrow(linetable); i++)
            {
                double r = cpl_table_get_double(linetable, "Resol", i, NULL);
                if (r < resol_avg - kappar*resol_stdev ||
                r > resol_avg + kappar*resol_stdev)
                {
                    cpl_table_set_invalid(linetable, "Resol", i);
                    cpl_table_set_invalid(linetable, "deltaLambda", i);
                }
            }
        }

        /* check( uves_plot_table(linetable, "X", "Resol", 
           "(x, l / dl)"), "Plotting failed"); 
           check( uves_plot_table(linetable, "Y", "Resol", 
           "(y, l / dl)"), "Plotting failed"); 
        */
        check( uves_plot_table(linetable, LINETAB_LAMBDAC, "Resol",
                       "(l, l / dl)"), "Plotting failed");
        }
        
        /* Plot identifications */
            uves_free_table(&temp);
        check( temp = cpl_table_duplicate(linetable), 
           "Error copying line table");
        check( uves_erase_invalid_table_rows(temp, "Ident"), 
           "Error removing un-identified lines");
        check( uves_plot_table(temp, "X", "Ynew", 
                   "Line identifications"), 
           "Plotting failed");
        uves_free_table(&temp);

    } /* Plots for middle (object) window */

    if (debug_mode) 
    {
        /* Results of uves_wavecal_search are already 
           drawn on debug_image */
        
        /* Draw guess table lines using the initial solution */
        if (0) check( uves_draw_lines(debug_image, initial_dispersion,
                      order_locations, guess,
                      "Ident", "Order", relative_order,
                      -1, -1,
                      false,     /* true = vertical */
                      12), "Error drawing guess solution");
        
        /* Draw catalogue lines using the initial solution */
        check( uves_draw_lines(debug_image, initial_dispersion, order_locations,
                   line_refer, "Wave", NULL, relative_order,
                   uves_min_int(*first_absolute_order, *last_absolute_order),
                   uves_max_int(*first_absolute_order, *last_absolute_order),
                   true,                   /* true = vertical */
                   8), "Error drawing catalogue lines");
        
        /* Draw catalogue lines using the final solution */
        check( uves_draw_lines(debug_image, *dispersion_relation, order_locations, 
                   line_refer, "Wave", NULL, relative_order,
                   uves_min_int(*first_absolute_order, *last_absolute_order),
                   uves_max_int(*first_absolute_order, *last_absolute_order),
                   true,                     /* true = vertical */
                   16), "Error drawing catalogue lines");
        
        /* Draw detected lines using initial solution */
        if (0) check( uves_draw_lines(debug_image, initial_dispersion, 
                      order_locations, linetable, 
                      LINETAB_LAMBDAC, "Order", relative_order,
                      -1, -1,
                      false,            /* true = vertical */
                      -16), "Error drawing detected lines");
        
        /* Draw IDed lines */
            uves_free_table(&temp);
        check(( temp = cpl_table_duplicate(linetable),
            /* Delete rows with invalid 'Ident' */
            uves_erase_invalid_table_rows(temp, "Ident")),
          "Error duplicating table");
        
        check( uves_draw_lines(debug_image, *dispersion_relation, order_locations,
                   temp, LINETAB_LAMBDAC, "Order", relative_order,
                   -1, -1,
                   true,                  /* true = vertical */
                   0), "Error drawing detected lines");
        

        /* Explanation of drawing produced by code above:
           The output frame emission lines will look like this


                     #### |1
                         ####   |2  
                         #|3#
                         -+--
                         #|##
                         ####
                         ####
        

                Legend:

                  ##: The emmission line
                  --: (horizontal line) The line was detected
                  |1: (vertical line)   Predicted position (final solution)
                  |2: (vertical line)   Predicted position (initial solution)
                  |3: (vertical line)   Is drawn iff the line was identified


          
        */

        /* Save the raw arc frame with detected emission lines marked */
        check( uves_save_image_local("Debug image", "rawdebug", 
                     debug_image, chip, trace, window,
                     rotated_header, true),
           "Error saving spectrum");
    }

    if (flames)
    {
        int start = 0;
        int count = cpl_table_get_nrow(linetable);

        check_nomsg( cpl_table_new_column(linetable, "Fibre", CPL_TYPE_INT) );

        cpl_table_fill_column_window(linetable, "Fibre",
                     start, count,
                     trace + 1);   /* Write value in range 1-9 */
    }

  cleanup:

    uves_free_table(&info_tbl);
    uves_free_table(&temp);
    uves_free_image(&weights_opt);
    uves_free_table(&cr_table);
    uves_free_image(&cr_image);
    uves_free_image(&spectrum);
    uves_free_image(&spectrum_noise);
    uves_free_image(&debug_image);
    uves_free_image(&rebinned);
    uves_free_image(&rebinned_noise);
    uves_free_image(&merged);
    uves_free_image(&merged_noise);
    uves_free_propertylist(&spectrum_header);
    uves_free_propertylist(&rebinned_header);
    uves_free_propertylist(&merged_header);
    cpl_free(relative_order);
    uves_polynomial_delete(&initial_dispersion);
    uves_free_table(&order_trace);

    return linetable;
}



/*----------------------------------------------------------------------------*/
/**
   @brief    Get the command line options and execute the data reduction
   @param    frames      The input frames list
   @param    flames      FLAMES mode?
   @param    recipe_id   recipe name
   @param    parameters  The recipe parameters list
   @return   CPL_ERROR_NONE if everything is ok

   This function is the executor function for the wavecal recipe. It 
   - reads the parameter list,
   - adjusts recipe parameters for binning, 
   - loads the first raw frame and performs the calibration by
   calling @c uves_wavecal_process_window() on each extraction window of each
   fibre trace.
   - saves the linetables to disk.

   The fibre trace information is read from the provided order table, while
   the information on the number and width of extraction windows are defined by
   recipe parameters.

   The line tables for all traces and windows are stored in a single output FITS
   file, however line tables for different chips (REDL, REDU) are stored in separate
   output files.
*/
/*----------------------------------------------------------------------------*/
void
uves_wavecal_exe_body(cpl_frameset *frames,
              bool flames,
              const char *recipe_id,
              const cpl_parameterlist *parameters,
              const char *starttime)
{
    /*
     * Variables containg the values of recipe parameters 
     */

    /* General */
    bool debug_mode;

    /* Extraction */
    int NWINDOWS;
    double OFFSET;
    double SLITLENGTH_par;  /* slitlength given by user */

    /* Search */
    int    RANGE;
    int    MAXLINES;
    int    MINLINES;
    centering_method CENTERING_METHOD;

    /* First solution */
    double SHIFT_MAX;
    double SHIFT_STEP;
    double SHIFT_TOLERANCE;

    /* Identify */
    double  ALPHA;
    double  MAXERROR;
    int     DEGREE;
    /* Calibrate */
    double  TOLERANCE;
    double  KAPPA;

    /* Input */
    cpl_image        *arclamp[2]          = {NULL, NULL};
    cpl_image        *arclamp_noise       = NULL;
    uves_propertylist *arclamp_header[2]   = {NULL, NULL};
    uves_propertylist *rotated_header[2]   = {NULL, NULL};
    
    /* Order table */
    cpl_table        *ordertable            = NULL;
    uves_propertylist *ordertable_header     = NULL;
    polynomial       *order_locations       = NULL;
    cpl_table        *traces                = NULL;
    
    /* Bias */
    cpl_image        *master_bias        = NULL;
    uves_propertylist *master_bias_header = NULL;

    /* Flat field */
    cpl_image        *master_flat        = NULL;
    cpl_image        *mflat_noise        = NULL;
    uves_propertylist *master_flat_header = NULL;

    /* Weight map */
    cpl_image        *weights            = NULL;

    /* DRS guess table is not used */
    /* 
       cpl_table        *drs_table             = NULL;
       uves_propertylist *drs_header            = NULL;
    */

    //FLAMES-DRS specific descriptors
    int* fibres_mask=NULL;
    double* fibres_pos=NULL;

    /* Guess line table */
    cpl_table        *guess              = NULL;
    polynomial       *absolute_order     = NULL;
    
    /* Velocity correction table */
    cpl_table        *corvel             = NULL;
    cpl_frame        *corvel_frm         = NULL;
    uves_propertylist *corvel_header      = NULL;

    /* Reference catalogue */
    cpl_table        *line_refer            = NULL;
    cpl_table        *line_intmon           = NULL;
   
    /* Output */
    lt_type          *linetable             = NULL;

    uves_propertylist *primary_header        = NULL;
    uves_propertylist *table_header          = NULL;
    /* QC for resolution + intmon + NULL */
    cpl_table        *qclog[3]              = {NULL, NULL, NULL};

    /* Local variables */
    cpl_image        *absorder_image     = NULL;
    const char *arclamp_filename = "";
    const char *line_refer_filename = "";
    const char *line_intmon_filename = "";
    char *product_filename = NULL;
    char *temp = NULL;
    bool blue  = false;
    bool sim_cal = false;
    enum uves_chip chip;
    int binx = 0;
    int biny = 0;
    bool drs_cvel_sw=false;
    const char* PROCESS_CHIP=NULL;
    extract_method em;

    /* Read recipe parameters */
    {
    const char *centering_m = "gaussian";
        const char *profile = "";

 

    /* General */
    check( uves_get_parameter(parameters, NULL, "uves", "debug", 
                  CPL_TYPE_BOOL, &debug_mode), "Could not read parameter");
    
    check( uves_get_parameter(parameters, NULL, "uves", "process_chip", CPL_TYPE_STRING, &PROCESS_CHIP),
               "Could not read parameter");
    uves_string_toupper((char*)PROCESS_CHIP);

        /* Extraction */
    check( uves_get_parameter(parameters, NULL, recipe_id, "nwindows",
                  CPL_TYPE_INT   , &NWINDOWS    ), "Could not read parameter");
    check( uves_get_parameter(parameters, NULL, recipe_id, "length",
                  CPL_TYPE_DOUBLE, &SLITLENGTH_par), "Could not read parameter");
    check( uves_get_parameter(parameters, NULL, recipe_id, "offset",
                  CPL_TYPE_DOUBLE, &OFFSET      ), "Could not read parameter");
    
    /* Don't allow optimal extraction. This requires that
       additional arguments (weight image, ..) are passed to uves_extract */
    temp = uves_sprintf("%s.%s", recipe_id, UVES_EXTRACT_ID);
    check( em = uves_get_extract_method(parameters, NULL, temp),
           "Could not read extraction method");
    
    check( uves_get_parameter(parameters, NULL, recipe_id, UVES_EXTRACT_ID ".profile",
                  CPL_TYPE_STRING, &profile), "Could not read parameter");

    assure( em == EXTRACT_LINEAR || em == EXTRACT_AVERAGE || em == EXTRACT_WEIGHTED ||
            (em == EXTRACT_OPTIMAL && strcmp(profile, "constant") == 0),
            CPL_ERROR_UNSUPPORTED_MODE, 
            "Only linear/average/weighted/optimal(constant profile) extraction "
            "methods are supported by this recipe");
    
    /* Search */
    check( uves_get_parameter(parameters, NULL, recipe_id, "search.range",
                  CPL_TYPE_INT   , &RANGE       ), "Could not read parameter");
    check( uves_get_parameter(parameters, NULL, recipe_id, "search.minlines",
                  CPL_TYPE_INT   , &MINLINES    ), "Could not read parameter");
    check( uves_get_parameter(parameters, NULL, recipe_id, "search.maxlines",
                  CPL_TYPE_INT   , &MAXLINES    ), "Could not read parameter");
/*
    check( uves_get_parameter(parameters, NULL, recipe_id, "search.centeringmethod",
                  CPL_TYPE_STRING, &centering_m ), "Could not read parameter");
*/
    if      (strcmp(centering_m, "gravity" ) == 0) CENTERING_METHOD = CENTERING_GRAVITY;
    else if (strcmp(centering_m, "gaussian") == 0) CENTERING_METHOD = CENTERING_GAUSSIAN;
    else
        { 
        /* Impossible */ assure(false, CPL_ERROR_ILLEGAL_INPUT,
                    "Unrecognized parameter value '%s'", centering_m);
        }
    
    /* First solution */
    check( uves_get_parameter(parameters, NULL, recipe_id, "first.shiftmax",
                  CPL_TYPE_DOUBLE , &SHIFT_MAX      ),
           "Could not read parameter");
    check( uves_get_parameter(parameters, NULL, recipe_id, "first.shiftstep", 
                  CPL_TYPE_DOUBLE , &SHIFT_STEP     ), 
           "Could not read parameter");
    check( uves_get_parameter(parameters, NULL, recipe_id, "first.shifttoler",
                  CPL_TYPE_DOUBLE , &SHIFT_TOLERANCE), 
           "Could not read parameter");
    
    /* Identify */
    check( uves_get_parameter(parameters, NULL, recipe_id, "identify.alpha",
                  CPL_TYPE_DOUBLE , &ALPHA      ), "Could not read parameter");
    check( uves_get_parameter(parameters, NULL, recipe_id, "identify.maxerror",
                  CPL_TYPE_DOUBLE , &MAXERROR   ), "Could not read parameter");
    check( uves_get_parameter(parameters, NULL, recipe_id, "identify.degree",
                  CPL_TYPE_INT    , &DEGREE     ), "Could not read parameter");
    
    /* Calibrate */
    check( uves_get_parameter(parameters, NULL, recipe_id, "calibrate.tolerance",
                  CPL_TYPE_DOUBLE, &TOLERANCE   ), "Could not read parameter");
    check( uves_get_parameter(parameters, NULL, recipe_id, "calibrate.kappa",
                  CPL_TYPE_DOUBLE, &KAPPA   ), "Could not read parameter");
    
    /* Additional checks */
    if (CENTERING_METHOD == CENTERING_GRAVITY)
        {
        uves_msg_warning("Centering method 'gravity' might lead to inaccurate "
                 "results. Recommended is 'gaussian'");
        }
    }

    /* Load raw image and header, and identify input frame as red or blue */
    check( uves_load_arclamp(frames, flames, &arclamp_filename, arclamp, arclamp_header, 
                 rotated_header, &blue, &sim_cal), "Error loading raw frame");

    /* Load reference line table */
    check( uves_load_linerefertable(frames, &line_refer_filename, &line_refer, NULL),
       "Could not load line reference table");
    uves_msg("Using line reference table '%s'", line_refer_filename);
    
    /* Load INTensity MONitoring table if present */
    if (cpl_frameset_find(frames, UVES_LINE_INTMON_TABLE) != NULL)
    {
        uves_free_table(&line_intmon);
        check( uves_load_lineintmon(frames, &line_intmon_filename,
                    &line_intmon),
           "Error loading line reference table");
        
        uves_msg("Using bright line table '%s'", line_intmon_filename);
    }
  
    /*
     * Adjust parameters according to binning 
     */
    check (binx = uves_pfits_get_binx(arclamp_header[0]),
       "Could not read x binning factor from input header");
    check (biny = uves_pfits_get_biny(arclamp_header[0]),
       "Could not read y binning factor from input header");
    SLITLENGTH_par /= (1.0*binx);  /* Extraction slit length */
    RANGE /= (1.0*biny);           /* Search window */
    SHIFT_MAX /= (1.0*binx);       /* Max shift compared to guess solution */
    SHIFT_STEP /= (1.0*binx);

    /* After the default tolerance was lowered to 0.07, do not adjust it according to binning,
       which would cause too many lines to be rejected
       
       TOLERANCE /= (1.0*biny); 
    */


    /* Loop over one or two chips, over traces and
       over extraction windows */
    for (chip = uves_chip_get_first(blue); 
         chip != UVES_CHIP_INVALID; 
         chip = uves_chip_get_next(chip)) {


        if(strcmp(PROCESS_CHIP,"REDU") == 0) {
            chip = uves_chip_get_next(chip);
        }
      
        const char *ordertable_filename = "";
        const char *corvel_filename = "";
        const char *master_flat_filename = "";
        const char *master_bias_filename = "";
        const char *weights_filename = "";
        /* const char *drs_filename        = "";    not used */
        const char *guess_filename      = "";
        const char *chip_name = "";
        int ntraces;
        int tracerow;                       /* Index of table row */
        int raw_index = uves_chip_get_index(chip);
        int current_linetable_extension;
        int tab_in_out_oshift = -1;
        double tab_in_out_yshift = -1;
        double slitlength;
      
      uves_msg("Processing %s chip in '%s'",
               uves_chip_tostring_upper(chip), arclamp_filename);
      
      check_nomsg( chip_name = uves_pfits_get_chipid(arclamp_header[raw_index], chip));
      
      
      uves_msg_debug("binning = %dx%d", binx, biny);
      
      
      /* Load the order table for this chip */
      uves_free_table       (&ordertable);
      uves_free_propertylist(&ordertable_header);
      uves_polynomial_delete(&order_locations);
      uves_free_table       (&traces);
      
      
      check( uves_load_ordertable(frames, 
                                  flames,
                                  chip_name, 
                                  &ordertable_filename,
                                  &ordertable, 
                                  &ordertable_header,
                                  NULL,
                                  &order_locations, 
                                  &traces, 
                                  (flames) ? &tab_in_out_oshift : NULL,
                                  (flames) ? &tab_in_out_yshift : NULL,
                                  &fibres_mask, 
                                  &fibres_pos, /* fibre_pos,fibre_mask */
                                  chip, 
                                  false  /* load guess table? */),
             "Could not load order table");
      uves_msg("Using order table in '%s'", ordertable_filename);
      ntraces = cpl_table_get_nrow(traces);
      uves_free_double(&fibres_pos);
      uves_free_int(&fibres_mask);
      
      /* Load master bias if present */
      uves_free_image(&master_bias);
      uves_free_propertylist(&master_bias_header);
      if (cpl_frameset_find(frames, UVES_MASTER_BIAS(chip)) != NULL)
          {
              check( uves_load_mbias(frames, chip_name, &master_bias_filename, &master_bias, 
                                     &master_bias_header, chip),
                     "Error loading master bias");
              
              uves_msg_low("Using master bias in '%s'", master_bias_filename);
          }
      else
          {
              uves_msg_warning("Master bias not provided. Bias subtraction not done");
          }
      
      
      /* Load master flat if present */
      uves_free_image(&master_flat);
      uves_free_propertylist(&master_flat_header);
      if ((cpl_frameset_find(frames, UVES_MASTER_FLAT(chip)) != NULL ||
           cpl_frameset_find(frames, UVES_MASTER_DFLAT(chip)) != NULL ||
           cpl_frameset_find(frames, UVES_MASTER_IFLAT(chip)) != NULL ||
           cpl_frameset_find(frames, UVES_MASTER_TFLAT(chip)) != NULL))
          {
              check( uves_load_mflat(frames, chip_name, &master_flat_filename, &master_flat, 
                                     &master_flat_header, chip, NULL),
                     "Error loading master flat");
              
              uves_msg_low("Using master flat in '%s'", master_flat_filename);
          }
      else
          {
              uves_msg_warning("Master flat not provided. Flat-fielding not done");
          }

      /* Load weight map if present */
      if (em == EXTRACT_WEIGHTED) {
          uves_free_image(&weights);
          check( weights = uves_load_weights(frames, &weights_filename, chip),
                 "Error loading weight map");

          uves_msg_low("Using weight map %s", weights_filename);
      }

      if (flames)
          /* Load CORVEL table */
          {
             if ((corvel_frm=cpl_frameset_find(frames, FLAMES_CORVEL_MASK)))
                  {
                      check( uves_load_corvel(frames,
                                              &corvel, &corvel_header,
                                              &corvel_filename),
                             "Could not load velocity correction table");
                      
                      uves_msg("Using velocity correction table %s", corvel_filename);
                      drs_cvel_sw=true;
                      
                      
                  }
              else
                  {
                      uves_msg("No corvel table found. Switch off corvel");
                      corvel = NULL;
                  }
          }

      /* Allocate all line tables for this chip */
      uves_lt_delete(&linetable);
      linetable = uves_lt_new(NWINDOWS, ntraces);


      /* Init QC tables for this chip */
      uves_qclog_delete(&qclog[0]); qclog[0] = uves_qclog_init(arclamp_header[raw_index], chip);
      uves_qclog_delete(&qclog[1]); qclog[1] = uves_qclog_init(arclamp_header[raw_index], chip);

      /* Saving the rotated raw arc frame */
      if (debug_mode) check( uves_save_image_local("Arc lamp frame", "raw", 
                                              arclamp[raw_index], 
                                              chip, -1, -1, rotated_header[raw_index], true),
                        "Error saving arc lamp frame");

      if (master_bias != NULL)
          {
              uves_msg("Subtracting master bias");
                
              check( uves_subtract_bias(arclamp[raw_index], master_bias),
                     "Error during bias subtraction");
          }
      else {
          /* In lack of a real master bias frame, one might subtract the bias by estimating it
             as the median value across the chip (which should be okay for arc lamp frames)
          */
            
          /* Disabled. Would need to be tested. probably doesn't make any big difference anyway
             double bias = cpl_image_get_median(ff);
             uves_msg_debug("Estimated bias level is %f ADU", bias);
             cpl_image_subtract_scalar(ff, bias);
          */
      }


      /* Define arc lamp noise */
      uves_free_image(&arclamp_noise);
      check( arclamp_noise = uves_define_noise(arclamp[raw_index], 
                                               arclamp_header[raw_index], 1, chip),
             "Could not set arc lamp noise");

      if (master_flat != NULL)
          {
              uves_msg("Dividing by master flat");
                
              uves_free_image(&mflat_noise);
              check( mflat_noise = 
                     uves_define_noise(master_flat, master_flat_header, 
                                       uves_pfits_get_datancom(master_flat_header),
                                       chip),
                     "Could not set master flat error bars");
                
              check( uves_flatfielding(arclamp[raw_index], arclamp_noise,
                                       master_flat, mflat_noise),
                     "Error while dividing by flat field");
          }

      if (debug_mode) check( uves_save_image_local("Pre-processed arc lamp frame",
                                              "preproc", 
                                              arclamp[raw_index], chip, -1, -1,
                                              rotated_header[raw_index], true),
                        "Error saving arc lamp frame");

      /* Set appropriate slitlength if user did not */
      if (SLITLENGTH_par < 0) {
            
          double header_full_slit;
            
          check( header_full_slit = 
                 uves_pfits_get_slitlength_pixels(arclamp_header[raw_index], chip),
                 "Could not read slit length");
            
          /* Avoid pixels at the edge of the slit
           *  which are likely to be noisy 
           */
          slitlength = uves_max_double(1.0, (header_full_slit - 2)/NWINDOWS);

          uves_msg("Full slit = %.2f pixels", header_full_slit);
      }
      else {
          slitlength = SLITLENGTH_par;
      }
                
   
      /* Loop over traces */
      for(tracerow = 0; tracerow < ntraces; tracerow++) {
        double trace_offset;
        int trace_number;
        int trace_enabled;
        
        trace_offset  = cpl_table_get_double(traces, "Offset"   , tracerow, NULL);
        trace_number  = cpl_table_get_int   (traces, "TraceID"  , tracerow, NULL);
        trace_enabled = cpl_table_get_int   (traces, "Tracemask", tracerow, NULL);

            if (ntraces > 1) {
                uves_msg("Processing trace %d", trace_number);
            }

        if (flames && sim_cal) 
        {
            /* Only calibrate SIMCAL fibre in SIMCAL mode */
            trace_enabled = (trace_number == 1) ? 1 : 0;
        }

        uves_msg_low("Trace offset = %.2f pixels ; enabled = %d",
             trace_offset, trace_enabled);
        
        assure( flames || trace_number == 0, CPL_ERROR_ILLEGAL_INPUT,
            "%s: UVES trace number must be 0, it is %d", 
            ordertable_filename, trace_number );
        
        
        if (trace_enabled != 0) {
        int window;          /* window number */

        /* Load guess line table for this trace, any window */
        uves_free_table       (&guess);
        uves_polynomial_delete(&absolute_order);

        check( uves_load_linetable(
               frames, flames, chip_name, order_locations,
               cpl_table_get_column_min(ordertable, "Order"),
               cpl_table_get_column_max(ordertable, "Order"),
               &guess_filename, &guess, 
               NULL, NULL,
               &absolute_order, chip, trace_number, -1),
               "Could not load guess line table for trace number %d", 
               trace_number);
        uves_msg("Using guess line table '%s'", guess_filename);
 
        if (debug_mode)
            {
            /* Create an image showing the polynomial   m = f(x,y)
               where m is the absolute order number
            */
            int x, y;

            absorder_image = cpl_image_new(cpl_image_get_size_x(arclamp[raw_index]),
                               cpl_image_get_size_y(arclamp[raw_index]),
                               CPL_TYPE_FLOAT);
            assure_mem(absorder_image);
            
            for (y = 1; y <= cpl_image_get_size_y(arclamp[raw_index]); y++)
                {
                for (x = 1; x <= cpl_image_get_size_x(arclamp[raw_index]); x++)
                    {
                    double absorder = 
                        uves_polynomial_evaluate_2d(absolute_order, x, y);
                    
                    cpl_image_set(absorder_image, x, y, absorder);
                    }
                }
            
            check( uves_save_image_local("Absolute order image", "absord",
                             absorder_image, chip, trace_number,
                             1, rotated_header[raw_index], true),
                   "Error saving absolute order image");
            
            uves_free_image(&absorder_image);
            }

        /* Loop over sky windows */
        for (window = 1; window <= NWINDOWS; window ++) {
            /*      
                | -trace_offs- | 
                |              |  -win_offs-  |
                |              |              |
                |              | | -glb_offs- | 
                        order_loc.=0       |
                |
                window_loc.
            */
                
            /* Adjacent windows are separated by slitlength,
               and offset is zero when window_number = (NWINDOWS+1)/2,
               so the formula is */
            double window_offset =
            slitlength * (window - (NWINDOWS+1) / 2.0);

            /* Total offset (see sketch) */
            double offset = trace_offset + window_offset + OFFSET;

            /* Number of lines to detect. Use defaults if 
               parameter values are negative */
#if 0
            int lines_min = (MINLINES >= 1) ? MINLINES : 
            (chip == UVES_CHIP_REDU) ? 1000 : 2000;
            int lines_max = (MAXLINES >= 1) ? MAXLINES : 
            (chip == UVES_CHIP_REDU) ? 1400 : 2400;
#else
/* like MIDAS (fewer lines): */
            int lines_min = (MINLINES >= 1) ? MINLINES : 
            (chip == UVES_CHIP_REDU) ? 1000 : 1100;
            int lines_max = (MAXLINES >= 1) ? MAXLINES : 
            (chip == UVES_CHIP_REDU) ? 1400 : 1600;
#endif

            assure( lines_min <= lines_max , CPL_ERROR_ILLEGAL_INPUT, 
                "Minimum and maximum number of requested line "
                "detections don't make sense (min = %d; max = %d)",
                lines_min, lines_max);

            if (NWINDOWS > 1) {
                uves_msg("Processing window %d of %d", window, NWINDOWS);
            }
                    
            passure( *(uves_lt_get_disprel(linetable, window, trace_number))
                 == NULL, "%d %d", window, trace_number);
            passure( *(uves_lt_get_absord (linetable, window, trace_number))
                 == NULL, "%d %d", window, trace_number);


            if (weights != NULL) {
                /* Object weighted extraction, use offset = 0 three times */
                offset = 0;
            }
                    
            /* Set absord guess solution */
            *uves_lt_get_absord(linetable, window, trace_number) = 
            uves_polynomial_duplicate(absolute_order);

            /* Execute macrosteps */
            check( *uves_lt_get_table(linetable, window, trace_number) =
                   uves_wavecal_process_window(
                       /* Raw */
                       arclamp[raw_index], 
                       arclamp_noise,
                       rotated_header[raw_index],
                       /* Calibration */
                       ordertable, order_locations,
                       master_flat != NULL,
                       weights,
                       /* drs_table,  not used */
                       guess,
                       line_refer,
                       flames,
                       tab_in_out_oshift,
                       tab_in_out_yshift,
                   /* Which window? */
                       chip, biny, trace_number, window,NWINDOWS,
                   /* Parameters    */
                   debug_mode,
                   /* Extract       */
                   offset, slitlength, parameters, recipe_id,
                   /* Search        */
                   RANGE, lines_min, lines_max, CENTERING_METHOD,
                   /* First solution */
                   SHIFT_MAX, SHIFT_STEP, SHIFT_TOLERANCE,
                   /* Identify      */
                   ALPHA, MAXERROR, DEGREE,
                   /* Calibrate     */
		       TOLERANCE, KAPPA,
                       corvel_frm,qclog,
                   /* Returned */
                   uves_lt_get_disprel(linetable, window, trace_number),
                   uves_lt_get_absord (linetable, window, trace_number),
                   uves_lt_get_firstabs(linetable, window, trace_number),
                   uves_lt_get_lastabs(linetable, window, trace_number)),
               "Wavelength calibration failed");


            //If CORVEL map is provided we perform the 
                    //corresponding analysis
            if(drs_cvel_sw) {
              /*
              define/local ord_min/i/1/1 0
              define/local ord_max/i/1/1 0
              define/local rsample/d/1/1 0
              statistic/table {ORDTAB} :ORDER >Null
              ord_min = outputr(1)
              ord_max = outputr(2)
              rsample = {{LINTAB},PIXEL(1)}
              rsample = 2./3. * rsample
              rebin/echelle  {ofrm} w{ofrm} {rsample} NONL {LINTAB} {SESSOUTV}
              mercut/echelle w{ofrm} mw{ofrm} {ord_min},{ord_max} NOAPPEND 
              !corvel stuff
              define/local OLD_CVEL_MAX/d/1/1 {DRS_CVEL_MAX}
              define/local OLD_CVEL_MIN/d/1/1 {DRS_CVEL_MIN}
              @p flames_reduce,VCORREL x1_rbf_ cvel1 0 {ord_max} {parCorVelTab} _0 _{chip({PATHID})} 0
              DRS_CVEL_MAX = DRS_CVEL_MAX + {q1}
              DRS_CVEL_MIN = DRS_CVEL_MIN + {q1}

              @p flames_reduce,VCORREL x1_rbf_ cvel2 0 {ord_max} {parCorVelTab} _0 _{chip({PATHID})} 0
              cvel_0 = {q1}
              DRS_CVEL_MAX = OLD_CVEL_MAX
              DRS_CVEL_MIN = OLD_CVEL_MIN
              */


            }



            if (flames ||                             /* FLAMES: all fibres */
            (window == 2 && trace_number == 0)|| /* UVES: central window */
            (window == 1 && NWINDOWS == 1)) { /* UVES: special user setting */

            check( uves_wavecal_qclog(
                   *uves_lt_get_table(
                       linetable,
                       window,
                       trace_number),
                                   *uves_lt_get_firstabs(linetable, window, trace_number),
                                   *uves_lt_get_lastabs(linetable, window, trace_number),
                   arclamp[raw_index], 
                   arclamp_header[raw_index],
                   flames, 
                   trace_number, trace_enabled, trace_offset,
                   chip,
                                   qclog[0]),
                   "Could not calculate resolution QC parameters");

            if (line_intmon != NULL) {
                check( uves_wavecal_qclog_intmon(
                       *uves_lt_get_table(
                       linetable,
                       window,
                       trace_number), 
                       line_intmon,
                       arclamp_header[raw_index],
                       flames, trace_number,
                       chip,
                       qclog[1]),
                   "Could not calculate int.mon. QC parameters");
            }
            else
                {
                /* Kill initialization and set pointer to NULL */
                uves_qclog_delete(&qclog[1]);
                }
            }
                    
            /* Finished processing. Save later (because
               all QC parameters must be available
               when the product file is first created). */

        }/* for each window... */
            
        }/* if trace enabled? */
        else
        {
            uves_msg("Skipping trace number %d", trace_number);
        }
    }/* for each trace... */

    /* Finished calculating all line tables for current chip. Now save. */

    /* Prepare product filename 
       (which need not be calculated for each trace and window) */
    cpl_free(product_filename);
    check( product_filename = uves_line_table_filename(chip), "Error getting filename");
    current_linetable_extension = 1;
        
    /* Loop over traces */
    char extname[80];
    for(tracerow = 0; tracerow < cpl_table_get_nrow(traces); tracerow++)
        {
        int trace_number;
        double trace_offset;
        int trace_enabled;
            
        trace_offset  = cpl_table_get_double(traces, "Offset"    , tracerow, NULL);
        trace_number  = cpl_table_get_int   (traces, "TraceID"   , tracerow, NULL);
        trace_enabled = cpl_table_get_int   (traces, "Tracemask" , tracerow, NULL);
            
        if (trace_enabled != 0)
                    {
            int window;
                
            /* Loop over sky windows */
            for (window = 1; window <= NWINDOWS; window ++)
                {
                double window_offset =
                    slitlength * (window - (NWINDOWS+1) / 2.0);
            
                double offset = trace_offset + window_offset + OFFSET;
        
                /* Table header */
                uves_free_propertylist(&table_header);
                table_header   = uves_propertylist_new();
                check( uves_pfits_set_traceid     ( table_header, trace_number),
                       "Error writing trace ID to product header");
                check( uves_pfits_set_offset      ( table_header, offset),
                       "Error writing trace offset to product header");
                check( uves_pfits_set_windownumber( table_header, window),
                       "Error writing window number to product header");

                sprintf(extname,"LINE_T%d_W%d_X%d",trace_number,window,
                          current_linetable_extension);
                 uves_pfits_set_extname(table_header,extname);

                check( uves_pfits_set_firstabsorder( table_header, 
                                     *uves_lt_get_firstabs(
                                     linetable, 
                                     window,
                                     trace_number)),
                       "Error writing order number to product header");
                check( uves_pfits_set_lastabsorder( table_header,
                                    *uves_lt_get_lastabs(
                                    linetable, 
                                    window,
                                    trace_number)),
                       "Error writing order number to product header");

                /* Save line table + 2 polynomials (in 3 extensions) */

                if (current_linetable_extension == 1) {
                    uves_free_propertylist(&primary_header);
                    primary_header = uves_propertylist_new();
                    
                    if (flames)
		      {
			char values[80];
			/* The MIDAS pipeline writes this QC to the
			   header (always zero), but not as part of
			   the QC logging */
			check_nomsg( uves_flames_pfits_set_ccfposmax(
								     primary_header, 0.0) );
			//Add descriptors needed for science reduction
			/* FIBREMASK */

			uves_propertylist_append_string(primary_header, 
							"HISTORY", 
							"'FIBREMASK','I*4'");

			{
			  int i;
			  for (i = 0; i < N_FIBRES_MAX; i++) {
			    snprintf(values, 80, "%1.1d ",
				     cpl_table_get_int(traces,"Tracemask",
						       i,NULL));
			    uves_propertylist_append_string(primary_header,
							    "HISTORY", values);
			    uves_msg_debug("value=%d",
					     cpl_table_get_int(traces,
							       "Tracemask",
							       i,NULL));
			  }
			}
			uves_propertylist_append_string(primary_header,
							"HISTORY", " ");


			/* PIXEL */        
                        double pixelsize;
                        double wavestep;

			check( pixelsize = 
			       cpl_table_get_column_mean(
							 *uves_lt_get_table(
									    linetable,
									    window,
									    trace_number),
						     LINETAB_PIXELSIZE),
			       "Error reading mean pixelsize");
			uves_msg_warning("Average pixelsize = %f w.l.u.", 
				       pixelsize);

			wavestep = pixelsize*2.0/3;

			uves_propertylist_append_string(primary_header, 
							"HISTORY", 
							"'PIXEL','R*4'");
			snprintf(values,80,"%14.7g %14.7g",pixelsize,pixelsize);
			uves_propertylist_append_string(primary_header,
							"HISTORY", values);
			uves_propertylist_append_string(primary_header,
							"HISTORY", " ");
                      


		      }

                    uves_msg("Creating line table '%s'", product_filename);
                    uves_pfits_set_extname(primary_header,"wavelength calibration solution");
                    sprintf(extname,"LINE_T%d_W%d_X%d",trace_number,window,
                             current_linetable_extension);
                    uves_pfits_set_extname(table_header,extname);

                    check( uves_frameset_insert(
                           frames,
                           *uves_lt_get_table(
                               linetable,
                               window,
                               trace_number),
                           CPL_FRAME_GROUP_PRODUCT,
                           CPL_FRAME_TYPE_TABLE,
                           CPL_FRAME_LEVEL_INTERMEDIATE,
                           product_filename,
                           UVES_LINE_TABLE(flames, chip),
                           arclamp_header[raw_index],     
                           primary_header,
                           table_header,
                           parameters,
                           recipe_id,
                           PACKAGE "/" PACKAGE_VERSION,
                           qclog, starttime, true, 0),
                           "Could not add line table '%s' (%s) to frameset",
                           product_filename, UVES_LINE_TABLE(flames, chip));

                    uves_msg("Line table '%s' (%s) added to frameset",
                         product_filename, UVES_LINE_TABLE(flames, chip));
                    }
                else   /* If this is not the first line table, 
                      append to the existing file */
                    {
                    check( uves_table_save(
                           *uves_lt_get_table(linetable,
                                      window,
                                      trace_number),
                           NULL,            /* Primary header,
                                       ignored when mode
                                       is IO_EXTEND */
                           table_header,    /* Extension header */
                           product_filename,/* This file 
                                       already exists */
                           CPL_IO_EXTEND),   /* Append to
                                    existing file */
                           "Error appending table to file '%s'", 
                           product_filename);
                    }
                current_linetable_extension += 1;
                sprintf(extname,"LINE_T%d_W%d_X%d",trace_number,window,
                         current_linetable_extension);
                uves_pfits_set_extname(table_header,extname);
                /* Save in next extension */
                check( uves_save_polynomial(*uves_lt_get_disprel(
                                linetable, 
                                window,
                                trace_number),
                                product_filename,
                                table_header), 
                       "Could not write polynomial to file '%s'", 
                       product_filename);
                current_linetable_extension += 1;
                sprintf(extname,"LINE_T%d_W%d_X%d",trace_number,window,
                         current_linetable_extension);
                uves_pfits_set_extname(table_header,extname);
                /* Save in next extension */
                check( uves_save_polynomial(*uves_lt_get_absord(
                                linetable, 
                                window,
                                trace_number), 
                                product_filename,
                                table_header), 
                       "Could not write polynomial to file '%s'", 
                       product_filename);
                current_linetable_extension += 1;

                uves_msg("Line table for trace %d, window #%d "
                     "saved to extensions %d-%d of '%s'",
                     trace_number, window, 
                     current_linetable_extension - 3,
                     current_linetable_extension - 1, 
                     product_filename);

                } /* for each window */
            } /* if trace enabled */
        } /* for each trace */

      if(strcmp(PROCESS_CHIP,"REDL") == 0) {
	chip = uves_chip_get_next(chip);
      }


    }/* For each chip */

  cleanup:

    /* Input */
    uves_free_image(&arclamp[0]);
    uves_free_image(&arclamp[1]);
    uves_free_image(&arclamp_noise);
    uves_free_image(&absorder_image);
    uves_free_propertylist(&arclamp_header[0]);
    uves_free_propertylist(&arclamp_header[1]);
    uves_free_propertylist(&rotated_header[0]);
    uves_free_propertylist(&rotated_header[1]);


    uves_free_table(&ordertable);
    uves_free_propertylist(&ordertable_header);
    uves_free_table(&corvel);
    uves_free_propertylist(&corvel_header);
    uves_polynomial_delete(&order_locations);
    uves_polynomial_delete(&absolute_order);
    uves_free_table(&traces);

    
    uves_free_image(&master_bias);
    uves_free_propertylist(&master_bias_header);
    uves_free_image(&master_flat);
    uves_free_image(&mflat_noise);
    uves_free_propertylist(&master_flat_header);
    uves_free_image(&weights);

    
    /* DRS not used
       uves_free_table(&drs_table);
       uves_free_propertylist(&drs_header);
    */
    
    uves_free_table(&guess);

    uves_free_table(&line_refer);
    uves_free_table(&line_intmon);


    /* Output */
    uves_lt_delete(&linetable);
    uves_free_propertylist(&primary_header);
    uves_free_propertylist(&table_header);
    uves_qclog_delete(&qclog[0]);
    uves_qclog_delete(&qclog[1]);

    cpl_free(product_filename);
    cpl_free(temp);
     
    return;
}


/**
@brief computes QC log
@param linetable        line table on which QC log is computed
@param firstabs     first absolute order number (not necessarily
                    in line table but corresponding to orders
                    present in order table)
@param lastabs
@param raw_header   FITS header
@param chip         UVES chip ID
@param tolerance    tolerance of calibration (recipe parameter)
@param qclog        ouput QC log table
*/
static void uves_wavecal_qclog(const cpl_table* linetable,
                               int firstabs,
                               int lastabs,
                   const cpl_image *arclamp,
                   const uves_propertylist* raw_header,
                   bool flames,
                   int trace_number,
                   int fibre_mask,
                   double offset,
                   enum uves_chip chip,
                   cpl_table* qclog)
{

    const char *qc_fib_drsno_name= uves_qclog_get_qc_name("DRSNO", flames, trace_number);
    const char *qc_fib_seq_name  = uves_qclog_get_qc_name("SEQ", flames, trace_number);
    const char *qc_fib_pos_name  = uves_qclog_get_qc_name("POS", flames, trace_number);
    const char *qc_fib_msk_name  = uves_qclog_get_qc_name("MSK", flames, trace_number);
    const char *qc_fwhmavg_name  = uves_qclog_get_qc_name("FWHMAVG", flames, trace_number);
    const char *qc_fwhmrms_name  = uves_qclog_get_qc_name("FWHMRMS", flames, trace_number);
    const char *qc_fwhmmed_name  = uves_qclog_get_qc_name("FWHMMED", flames, trace_number);
    const char *qc_resolavg_name = uves_qclog_get_qc_name("RESOLAVG", flames, trace_number);
    const char *qc_resolrms_name = uves_qclog_get_qc_name("RESOLRMS", flames, trace_number);
    const char *qc_resolmed_name = uves_qclog_get_qc_name("RESOLMED", flames, trace_number);
    const char *qc_wlenmin_name  = uves_qclog_get_qc_name("WLENMIN", flames, trace_number);
    const char *qc_wlenmax_name  = uves_qclog_get_qc_name("WLENMAX", flames, trace_number);
    const char *qc_ordmin_name   = uves_qclog_get_qc_name("ORDMIN", flames, trace_number);
    const char *qc_ordmax_name   = uves_qclog_get_qc_name("ORDMAX", flames, trace_number);
    const char *qc_detected_ordmin_name   = uves_qclog_get_qc_name("ORDMIN DETECTED", flames, trace_number);
    const char *qc_detected_ordmax_name   = uves_qclog_get_qc_name("ORDMAX DETECTED", flames, trace_number);
    const char *qc_nlintot_name  = uves_qclog_get_qc_name("NLINTOT", flames, trace_number);
    const char *qc_nlinsel_name  = uves_qclog_get_qc_name("NLINSEL", flames, trace_number);
    const char *qc_nlinsol_name  = uves_qclog_get_qc_name("NLINSOL", flames, trace_number);
    const char *qc_line_werr_name  = uves_qclog_get_qc_name("LINE WAVEERR", flames, trace_number);
    const char *qc_line_wsys_name  = uves_qclog_get_qc_name("LINE SYSERR", flames, trace_number);
    const char *qc_nlinres1_name = uves_qclog_get_qc_name("NLINRES1", flames, trace_number);
    const char *qc_lineresidavg_name = 
    uves_qclog_get_qc_name("LINE RESIDAVG", flames, trace_number);
    const char *qc_lineresidrms_name = 
    uves_qclog_get_qc_name("LINE RESIDRMS", flames, trace_number);
    char comment[80];
    cpl_table *selected = NULL;
    double wmin=0;
    double wmax=0;
    double wcen=0;
    int nfinal=0;

    char test_id[80];
    sprintf(test_id,"%sResolution-Test-Results",flames ? "Fibre-" : "");

    check_nomsg(uves_qclog_add_string(qclog,
                      "QC TEST1 ID",
                      test_id,
                      "Name of QC test",
                      "%s"));

    check_nomsg( uves_qclog_add_common_wave(raw_header, chip, qclog) );

    if (flames)
    {
        /* Fibre ID */
        check_nomsg(uves_qclog_add_int(qclog,
                       qc_fib_drsno_name,
                       trace_number + 1,
                       "DRS det. fibre seq. pos.",
                       "%d"));
        
        /* Index */
        check_nomsg(uves_qclog_add_int(qclog,
                       qc_fib_seq_name,
                       trace_number + 1,
                       "det. fibre seq. no.",
                       "%d"));

        check_nomsg(uves_qclog_add_double(qclog,
                          qc_fib_pos_name,
                          offset,
                          "det. fibre seq. rel. pos.",
                          "%.4f"));

        check_nomsg(uves_qclog_add_int(qclog,
                       qc_fib_msk_name,
                       fibre_mask,
                       "DRS det. fibre mask value",
                       "%d"));

        {

        double exptime;
      
        check( exptime = uves_flames_pfits_get_dit(raw_header),
               "Error reading exposure time");
      
        check_nomsg(uves_qclog_add_double(qclog,
                          "QC FIB ABSTRANS",
                          cpl_image_get_flux(arclamp) / exptime,
                          "abs. trans. countrate",
                          "%.4f"));
        }
        {
        int n_hpix;
        int x, y;
      
        n_hpix = 0;
        for (y = 1; y <= cpl_image_get_size_y(arclamp); y++)
            for (x = 1; x <= cpl_image_get_size_x(arclamp); x++)
            {
                int pis_rejected;
                int value = cpl_image_get(arclamp, x, y, &pis_rejected);
          
                if (!pis_rejected &&
                (value < DRS_PTHRES_MIN || value > DRS_PTHRES_MAX))
                {
                    n_hpix += 1;
                }
            }

        check_nomsg(uves_qclog_add_int(qclog,
                           "QC NHOTPIX",
                           n_hpix,
                           "no. of hot pixels",
                           "%d"));
      
        }
  
        {
        int plate_id;
        check( plate_id = uves_flames_pfits_get_plateid(raw_header),
               "Error reading plate ID");
        check_nomsg(uves_qclog_add_int(qclog,
                           "QC PLATENO",
                           plate_id,
                           "Plate Id.",
                           "%d"));
        }

    } /* if flames */

    /* FLAMES + UVES common QC params */
    selected = uves_extract_table_rows(linetable, "NLinSol", CPL_NOT_EQUAL_TO, 0);
    /* FWHM in pixels */
    sprintf(comment,"average FWHM in X of sel lines on TRACE%d WIN2 [pix]",trace_number);
    check_nomsg(uves_qclog_add_double(qclog,
                      qc_fwhmavg_name,
                      cpl_table_get_column_mean(selected,"Xwidth")*TWOSQRT2LN2,
                      comment,
                      "%.2f"));

    sprintf(comment,"stdev FWHM in X of sel lines on TRACE%d WIN2 [pix]",trace_number);
    check_nomsg(uves_qclog_add_double(qclog,
                      qc_fwhmrms_name,
                      cpl_table_get_column_stdev(selected,"Xwidth")*TWOSQRT2LN2,
                      comment,
                      "%.4f"));

    sprintf(comment,"median FWHM in X of sel lines on TRACE%d WIN2 [pix]",trace_number);
    check_nomsg(uves_qclog_add_double(qclog,
                      qc_fwhmmed_name,
                      cpl_table_get_column_median(selected,"Xwidth")*TWOSQRT2LN2,
                      comment,
                      "%.4f"));
   
    sprintf(comment,"mean R of sel lines on TRACE%d WIN2",trace_number);
    check_nomsg(uves_qclog_add_double(qclog,
                      qc_resolavg_name,
                      cpl_table_get_column_mean(selected,"Resol"),
                      comment,
                      "%.4f"));

    sprintf(comment,"stdev R of sel lines on TRACE%d WIN2",trace_number);
    check_nomsg(uves_qclog_add_double(qclog,
                      qc_resolrms_name,
                      cpl_table_get_column_stdev(selected,"Resol"),
                      comment,
                      "%.4f"));

    sprintf(comment,"median R of sel lines on TRACE%d WIN2",trace_number);
    check_nomsg(uves_qclog_add_double(qclog,
                      qc_resolmed_name,
                      cpl_table_get_column_median(selected,"Resol"),
                      comment,
                      "%.4f"));

    /* Convert A -> picometers */
    sprintf(comment,"mean line pos resid on TRACE%d WIN2 [pm]",trace_number);
    check_nomsg(uves_qclog_add_double(qclog,
                      qc_lineresidavg_name,
                      cpl_table_get_column_mean(selected, LINETAB_RESIDUAL)*100,
                      comment,
                      "%.4f"));

    sprintf(comment,"sigma line pos resid on TRACE%d WIN2 [pm]",trace_number);
    check_nomsg(uves_qclog_add_double(qclog,
                      qc_lineresidrms_name,
                      cpl_table_get_column_stdev(selected, LINETAB_RESIDUAL)*100,
                      comment,
                      "%.4f"));

    /* Convert A -> nm */
    wmin=cpl_table_get_column_min(linetable,LINETAB_LAMBDAC)/10.0;
    sprintf(comment,"minimum wavelength on TRACE%d WIN2 [nm]",trace_number);
    check_nomsg(uves_qclog_add_double(qclog,
                      qc_wlenmin_name,
                      wmin,
                      comment,
                      "%.4f"));

    wmax=cpl_table_get_column_max(linetable,LINETAB_LAMBDAC)/10.0;
    sprintf(comment,"maximum wavelength on TRACE%d WIN2 [nm]",trace_number);
    check_nomsg(uves_qclog_add_double(qclog,
                      qc_wlenmax_name,
                      wmax,
                      comment,
                      "%.4f"));

    sprintf(comment,"minimum order number expected on TRACE%d WIN2",trace_number);
    check_nomsg(uves_qclog_add_int(qclog,
                   qc_ordmin_name,
                   uves_min_int(firstabs, lastabs),
                   comment,
                   "%d"));

    sprintf(comment,"maximum order number expected on TRACE%d WIN2",trace_number);
    check_nomsg(uves_qclog_add_int(qclog,
                   qc_ordmax_name,
                   uves_max_int(firstabs, lastabs),
                   comment,
                   "%d"));


    sprintf(comment,"minimum order number detected on TRACE%d WIN2",trace_number);
    check_nomsg(uves_qclog_add_int(qclog,
                   qc_detected_ordmin_name,
                   cpl_table_get_column_min(linetable,"Order"),
                   comment,
                   "%d"));

    sprintf(comment,"maximum order number detected on TRACE%d WIN2",trace_number);
    check_nomsg(uves_qclog_add_int(qclog,
                   qc_detected_ordmax_name,
                   cpl_table_get_column_max(linetable,"Order"),
                   comment,
                   "%d"));

    sprintf(comment,"No. of lines found on TRACE%d WIN2",trace_number);
    check_nomsg(uves_qclog_add_int(qclog,
                   qc_nlintot_name,
                   cpl_table_get_nrow(linetable),
                   comment,
                   "%d"));

    sprintf(comment,"No. of lines selected on TRACE%d WIN2",trace_number);
    check_nomsg(uves_qclog_add_int(qclog,
                   qc_nlinsel_name,
                   cpl_table_get_nrow(linetable) -
                                   cpl_table_count_invalid(linetable, "Ident"),
                   comment,
                   "%d"));

    nfinal=cpl_table_get_nrow(selected);
    sprintf(comment,"Final No. of lines used on TRACE%d WIN2",trace_number);
    check_nomsg(uves_qclog_add_int(qclog,
                   qc_nlinsol_name,
                   nfinal,
                   comment,
                   "%d"));

    cpl_table* extracted=NULL;
    double rms_wlu=0;
    const char* rms_wlu_alpha=NULL;
    double lines_sqrt=0;
    double lines_werr=0;
    int next=0;
    wcen=0.5*(wmin+wmax);
    check_nomsg(next=cpl_table_and_selected_string(qclog,"key_name",
                                                   CPL_EQUAL_TO,
                                                   "QC LINE RESIDRMS WLU"));
    check_nomsg(extracted=cpl_table_extract_selected(qclog));
  
    check_nomsg(rms_wlu_alpha=cpl_table_get_string(extracted,"key_value",0));
    rms_wlu=atof(rms_wlu_alpha);
    lines_sqrt=sqrt(nfinal);
    lines_werr=rms_wlu/lines_sqrt;
    sprintf(comment,"Wavelength error on TRACE%d [Angstrom]",trace_number);
    check_nomsg(uves_qclog_add_double(qclog,
                                   qc_line_werr_name,
                                   lines_werr,
                                   comment,
                                   "%.3g"));

    uves_free_table(&extracted);
    uves_free_table(&selected);
    selected = cpl_table_duplicate(linetable);
    assure_mem( selected );


    sprintf(comment,"Wavelength systematic error on TRACE%d [Angstrom]",trace_number);
    check_nomsg(uves_qclog_add_double(qclog,
                                   qc_line_wsys_name,
                                   wcen*100./SPEED_OF_LIGHT,
                                   comment,
                                   "%f"));


    /* Remove unidentified lines and
       lines with residual > 1 A      */
    check_nomsg( uves_erase_invalid_table_rows(selected, "Ident") );
    check_nomsg( uves_erase_table_rows(selected, LINETAB_RESIDUAL, 
                       CPL_NOT_LESS_THAN,
                       1.0) );
  
    sprintf(comment,"No. of lines with residuals < 0.1 nm on TRACE%d",trace_number);
    check_nomsg(uves_qclog_add_int(qclog,
                   qc_nlinres1_name,
                   cpl_table_get_nrow(selected),
                   comment,
                   "%d"));


  cleanup:

    uves_free_string_const(&qc_fib_drsno_name);
    uves_free_string_const(&qc_fib_seq_name);
    uves_free_string_const(&qc_fib_pos_name);
    uves_free_string_const(&qc_fib_msk_name);
    uves_free_string_const(&qc_fwhmavg_name);
    uves_free_string_const(&qc_fwhmrms_name);
    uves_free_string_const(&qc_fwhmmed_name);

    uves_free_string_const(&qc_resolavg_name);
    uves_free_string_const(&qc_resolrms_name);
    uves_free_string_const(&qc_resolmed_name);
    uves_free_string_const(&qc_wlenmin_name);
    uves_free_string_const(&qc_wlenmax_name);
    uves_free_string_const(&qc_ordmin_name);
    uves_free_string_const(&qc_ordmax_name);
    uves_free_string_const(&qc_nlintot_name);
    uves_free_string_const(&qc_nlinsel_name);
    uves_free_string_const(&qc_line_werr_name);
    uves_free_string_const(&qc_line_wsys_name);
    uves_free_string_const(&qc_nlinsol_name);
    uves_free_string_const(&qc_nlinres1_name);
    uves_free_string_const(&qc_lineresidavg_name);
    uves_free_string_const(&qc_lineresidrms_name);
    uves_free_string_const(&qc_detected_ordmin_name);
    uves_free_string_const(&qc_detected_ordmax_name);

    uves_free_table(&selected);

    return;
 
}

/**
@brief computes intmon QC log
@param table          line table on which QC log is computed
@param line_intmon    table of bright lines
@param raw_header     FITS header
@param chip           UVES chip ID
@param qclog          ouput QC log table
*/
static void uves_wavecal_qclog_intmon(cpl_table* table,
                      const cpl_table *line_intmon,
                      const uves_propertylist* raw_header,
                      bool flames,
                      int fibre,
                      enum uves_chip chip,
                      cpl_table* qclog)
{
    const char *qc_intavg_name = NULL;
    const char *qc_nlinint_name = NULL;

    cpl_table *temp = NULL;
        
    check_nomsg(uves_qclog_add_string(qclog,
                                      "QC TEST2 ID",
                                      flames ? "Fibre-Line-Intensity-Test-Results"
                                      : "Line-Intensity-Test-Results",
                                      "Name of QC test",
                                      "%s"));

    check_nomsg( uves_qclog_add_common_wave(raw_header,
                                            chip, qclog) );

    {
    double tolerance = 0.001; /* (A) 
                     The lines in the line table
                     and intmon table are considered
                     identical if the difference
                     is less than this number.
                  */
                      
    double exptime;

    int N_bright = cpl_table_get_nrow(line_intmon);
    int i;

    check( exptime = uves_pfits_get_exptime(raw_header),
           "Could not get exposure time");
    
    cpl_table_new_column(table, "Intensity", CPL_TYPE_DOUBLE); 
    for (i = 0; i < cpl_table_get_nrow(table); i++)
    {
        int is_null;
        double ident = cpl_table_get_double(table, "Ident", i, &is_null);

        if (!is_null)
        {
            int bright_index = uves_wavecal_find_nearest(
            line_intmon, ident, 0, N_bright-1);

            double bright = cpl_table_get_double(
            line_intmon, "Wave", bright_index, NULL);

            if (fabs(bright - ident) < tolerance)
            {
                double peak = cpl_table_get_double(table, "Peak", i, NULL);
                double pixelsize = 
                fabs(cpl_table_get_double(table, LINETAB_PIXELSIZE, i, NULL));

                double lambda_fwhm = cpl_table_get_double(table, "Xwidth", i, NULL)
                * TWOSQRT2LN2 * pixelsize;
                /* Line FWHM in wlu */

                double intensity = peak * lambda_fwhm / exptime;
                /* Same formula as in MIDAS */

                cpl_table_set_double(table, "Intensity", i, intensity);
            }
            else
            {
                cpl_table_set_invalid(table, "Intensity", i);
            }
        }
        else
        {
            cpl_table_set_invalid(table, "Intensity", i);
        }
    }
    }

    uves_free_table(&temp);
    temp = cpl_table_duplicate(table);
    uves_erase_invalid_table_rows(temp, "Intensity");

    {
    double mean;
    if (cpl_table_get_nrow(temp) == 0)
        {
        uves_msg_warning("No bright lines found!");
        mean = 0;
        }
    else
        {
        mean = cpl_table_get_column_mean(temp, "Intensity");
        }

    if (flames)
        {
        qc_intavg_name  = uves_sprintf("QC FIB%d INTAVG", fibre+1); /* Count 1-9 */
        qc_nlinint_name = uves_sprintf("QC FIB%d NLININT", fibre+1);
        }
    else
        {
        qc_intavg_name  = uves_sprintf("QC INTAVG");
        qc_nlinint_name = uves_sprintf("QC NLININT");
        }
        
    check_nomsg(uves_qclog_add_double(qclog,
                      qc_intavg_name,
                      mean,
                      "average intensity of line list on TRACE0 WIN2",
                      "%.4f"));
    
    check_nomsg(uves_qclog_add_int(qclog,
                       qc_nlinint_name,
                       cpl_table_get_nrow(temp),
                       "No. of lines to measure INTAVG on TRACE0 WIN2",
                       "%d"));
    }

 cleanup:
    uves_free_string_const(&qc_intavg_name);
    uves_free_string_const(&qc_nlinint_name);
  uves_free_table(&temp);
  return;
}


/**@}*/