/*                                                                              *
 *   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:46 $
 * $Revision: 1.125 $
 * $Name: not supported by cvs2svn $
 *
 */

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

/*----------------------------------------------------------------------------*/
/**
 * @defgroup uves_reduce        Substep: Reduce
 *
 * This module is for reducing a science frame.
 */
/*----------------------------------------------------------------------------*/
/**@{*/


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

#include "uves_reduce.h"

#include <uves.h>
#include <uves_extract.h>
#include <uves_backsub.h>
#include <uves_parameters.h>
#include <uves_flatfield.h>
#include <uves_rebin.h>
#include <uves_merge.h>
#include <uves_utils_cpl.h>
#include <uves_utils_wrappers.h>
#include <uves_pfits.h>
#include <uves_dfs.h>
#include <uves_dump.h>
#include <uves_qclog.h>
#include <uves_plot.h>
#include <uves_error.h>

#include <cpl.h>
#include <float.h>
#include <string.h>
/*-----------------------------------------------------------------------------
                            Functions prototypes
  -----------------------------------------------------------------------------*/
#define UVES_MIN_LINE_ROWS_TO_MAKE_FIT 5
#define UVES_BLAZE_DUMMY_VAL 999.
static cpl_error_code
extract_ff_rebin_merge(cpl_image *back_subbed, 
                       cpl_image *backsubbed_noise,
               const uves_propertylist *backsubbed_header,
               const cpl_image *master_flat, 
                       cpl_image *mflat_noise,
               const cpl_table *ordertable, 
                       const polynomial *order_locations,
               const cpl_table *linetable, 
                       const uves_propertylist *linetable_header[3],
                       const polynomial *dispersion_relation[3],
               double slit_length, 
                       double slit_offset, 
                       int window,
               enum uves_chip chip,
               bool blaze_correct,
                       bool tilt_corr,
                       bool debug_mode,
               const cpl_parameterlist *parameters, 
                       const char *context,
                       const char *mode,
               flatfielding_method ff_method,
               extract_method ee_method,
               merge_method m_method,
               /* Output */
               cpl_image **x, 
                       uves_propertylist **x_header,
               cpl_image **fx,
               cpl_table **cosmic_mask,
               cpl_image **wave_map,
               cpl_image **flatfielded_variance,
               uves_propertylist **flatfielded_variance_header,
               cpl_image **resampled_spectrum,
               cpl_image **resampled_mf,
               cpl_image **merged_sky,
               cpl_image **rebinned_spectrum, 
                       cpl_image **rebinned_noise, 
                       uves_propertylist **rebinned_header,
               cpl_image **merged_spectrum, 
                       cpl_image **merged_noise, 
                       uves_propertylist **merged_header,
                       cpl_table** info_tbl,
                       cpl_table **order_trace);

static cpl_image *
subtract_sky(cpl_image *rebinned_obj, 
             cpl_image *rebinned_obj_noise, 
             uves_propertylist *rebinned_obj_header,
         const cpl_image *rebinned_sky1, 
             const cpl_image *rebinned_sky1_noise, 
             const uves_propertylist *rebinned_sky1_header,
         const cpl_image *rebinned_sky2, 
             const cpl_image *rebinned_sky2_noise, 
             const uves_propertylist *rebinned_sky2_header,
         cpl_image **merged_obj, 
             cpl_image **merged_obj_noise, 
             uves_propertylist *merged_obj_header,
         const cpl_image *merged_sky1, 
             const cpl_image *merged_sky1_noise, 
             const uves_propertylist *merged_sky1_header,
         const cpl_image *merged_sky2, 
             const cpl_image *merged_sky2_noise, 
             const uves_propertylist *merged_sky2_header,
         double obj_slit, 
             double sky1_slit, 
             double sky2_slit);


static cpl_image *
subtract_sky_row(cpl_image *obj, 
                 cpl_image *obj_noise, 
                 double obj_start, 
                 double obj_end, 
                 double obj_slit,
         const cpl_image *sky1, 
                 const cpl_image *sky1_noise, 
                 double sky1_start, 
                 double sky1_end, 
                 double sky1_slit,
         const cpl_image *sky2, 
                 const cpl_image *sky2_noise, 
                 double sky2_start, 
                 double sky2_end, 
                 double sky2_slit,
         int row, 
                 double wavestep, 
                 double *common_start, 
                 double *common_end);

static double get_offset(const cpl_image *back_subbed, 
             const cpl_table *ordertable, 
                         const polynomial *order_locations,
             double search_range, 
                         int nsamples, 
                         double *doffset);

static cpl_image *
uves_get_blaze_ratio(const cpl_image *spectrum,
             const cpl_image *spectrum_noise);

/*-----------------------------------------------------------------------------
                              Implementation
  -----------------------------------------------------------------------------*/



/*----------------------------------------------------------------------------*/
/**
   @brief    Define recipe parameters used for reducing a frame
   @return   Newly allocated parameter list for this sub-step

   See source code or 'esorex --man-page' for a description of each parameter.
*/
/*----------------------------------------------------------------------------*/

cpl_parameterlist *
uves_reduce_define_parameters(void)
{
    const char *name = "";
    char *full_name = NULL;
    cpl_parameterlist *parameters = NULL;
    cpl_parameter *p = NULL;

    parameters = cpl_parameterlist_new();
    
    /**************
     *  Backsub   *
     **************/
    if (cpl_error_get_code() == CPL_ERROR_NONE)
    {
        uves_propagate_parameters_step(UVES_BACKSUB_ID, parameters, 
                       UVES_REDUCE_ID, NULL);
    }
    
    
    /*****************
     *  Extraction   *
     *****************/
    if (cpl_error_get_code() == CPL_ERROR_NONE)
    {
        uves_propagate_parameters_step(UVES_EXTRACT_ID, parameters, 
                       UVES_REDUCE_ID, NULL);
    }


    /******************
     *  Slit geometry *
     ******************/
    if (cpl_error_get_code() == CPL_ERROR_NONE)
    {
        name = "slitlength";
        full_name = uves_sprintf("%s.%s", UVES_REDUCE_ID, name);
        
        uves_parameter_new_range(p, full_name,
                    CPL_TYPE_DOUBLE,
                    "Extraction slit length (in pixels). "
                    "If negative, the value "
                    "inferred from the raw frame header is used",
                    UVES_REDUCE_ID,
                    -1.0,
                    -2.0, DBL_MAX);
        
        cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
        cpl_parameterlist_append(parameters, p);
        cpl_free(full_name);
    }

    if (cpl_error_get_code() == CPL_ERROR_NONE)
    {
        name = "skysub";
        full_name = uves_sprintf("%s.%s", UVES_REDUCE_ID, name);
        
        uves_parameter_new_value(p, full_name,
                     CPL_TYPE_BOOL,
                     "Do sky-subtraction (only applicable to linear "
                     "and average extractions)?",
                     UVES_REDUCE_ID,
                     true);
        
        cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
        cpl_parameterlist_append(parameters, p);
        cpl_free(full_name);
    }

    if (cpl_error_get_code() == CPL_ERROR_NONE)
    {
        name = "objoffset";
        full_name = uves_sprintf("%s.%s", UVES_REDUCE_ID, name);
        
        uves_parameter_new_value(p, full_name,
                     CPL_TYPE_DOUBLE,
                     "Offset (in pixels) of extraction slit "
                                     "with respect to center of order. "
                                     "This parameter applies to linear/average/"
                                     "optimal extraction. "
                                     "For linear/average extraction, if the related "
                                     "parameter objslit is negative, the offset is "
                     "automatically determined by measuring the "
                     "actual object position. ",
                     UVES_REDUCE_ID,
                     0.0);
        
        cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
        cpl_parameterlist_append(parameters, p);
        cpl_free(full_name);
    }

    if (cpl_error_get_code() == CPL_ERROR_NONE)
        {
            name = "objslit";
            full_name = uves_sprintf("%s.%s", UVES_REDUCE_ID, name);
            
            uves_parameter_new_range(p, full_name,
                                     CPL_TYPE_DOUBLE,
                                     "Object window size (in pixels). This must "
                                     "be less than the total slit length. If "
                                     "negative, the default value (half of full "
                                     "slit length) is used. The upper and lower "
                                     "sky windows are defined as the part of the "
                                     "full slit (if any) outside the object "
                                     "window. The center of the object window "
                                     "is determined by the offset parameter. "
                                     "This parameter does not apply to optimal "
                                     "extraction.",
                                     UVES_REDUCE_ID,
                                     -1.0,
                                     -2.0, DBL_MAX);
        
        cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
        cpl_parameterlist_append(parameters, p);
        cpl_free(full_name);
    }

    if (cpl_error_get_code() == CPL_ERROR_NONE)
    {
        name = "tiltcorr";
        full_name = uves_sprintf("%s.%s", UVES_REDUCE_ID, name);
        
        uves_parameter_new_value(p, full_name,
                     CPL_TYPE_BOOL,
                     "If enabled (recommended), the provided "
                                     "dispersion solutions "
                                     "obtained at different slit positions are "
                                     "interpolated linearly at the actually "
                                     "measured position of the object/sky. "
                                     "Line tilt correction is currently not supported "
                                     "for 2d extraction, in which case the "
                                     "dispersion solution obtained at the middle of "
                                     "the slit is always used.",
                     UVES_REDUCE_ID,
                     true);
        
        cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
        cpl_parameterlist_append(parameters, p);
        cpl_free(full_name);
    }



    /*****************
     *  Flatfielding *
     *****************/

    if (cpl_error_get_code() == CPL_ERROR_NONE)
    {
        name = "ffmethod";
        full_name = uves_sprintf("%s.%s", UVES_REDUCE_ID, name);
        
        uves_parameter_new_enum(p, full_name,
                    CPL_TYPE_STRING,
                    "Flat-fielding method. If set to 'pixel', "
                    "flat-fielding is done in pixel-pixel space "
                    "(before extraction); if set to 'extract', "
                    "flat-fielding is performed in pixel-order "
                    "space (i.e. after extraction). If set to "
                    "'no', no flat-field correction is done",
                    UVES_REDUCE_ID,
                    "extract",    /* 'Pixel' method is usually preferred,
                             but do like UVES/MIDAS */
                    3, 
                    "pixel", "extract", "no");
        
        cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
        cpl_parameterlist_append(parameters, p);
        cpl_free(full_name);
    }

    /*****************
     *  Blaze corr.  *
     *****************/

    if (cpl_error_get_code() == CPL_ERROR_NONE)
    {
/*
        name = "blazecorr";
        full_name = uves_sprintf("%s.%s", UVES_REDUCE_ID, name);
        
        uves_parameter_new_value(p, full_name,
                    CPL_TYPE_BOOL,
                    "(highly experimental, recommended=false) "
                    "Apply a correction for the different shapes "
                    "of flat-field and science blaze functions? "
                    "For this to be possible, flat-fielding method "
                    "must be different from 'no'.",
                    UVES_REDUCE_ID,
                    false);
        
        cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
        cpl_parameterlist_append(parameters, p);
        cpl_free(full_name);
*/
    }

    /*****************
     *  Rebinning    *
     *****************/
    if (cpl_error_get_code() == CPL_ERROR_NONE)
    {
        uves_propagate_parameters_step(UVES_REBIN_ID, parameters, 
                       UVES_REDUCE_ID, NULL);
    }
    

    /*****************
     *   Merging     *
     *****************/
    if (cpl_error_get_code() == CPL_ERROR_NONE)
    {
      name = "merge";
      full_name = uves_sprintf("%s.%s", UVES_REDUCE_ID, name);
        
      uves_parameter_new_enum(p, full_name,
			      CPL_TYPE_STRING,
                              "Order merging method. If 'optimal', the "
                              "flux in the overlapping region is set "
                              "to the (optimally computed, using the "
                              "uncertainties) average of single order "
                              "spectra. If 'sum', the flux in the "
                              "overlapping region is computed as the "
                              "sum of the single order spectra. If 'noappend' "
                              "the spectrum is simply rebinned but not merged."
                              "If flat-fielding is done, method 'optimal' "
                              "is recommended, otherwise 'sum'.",
                              UVES_REDUCE_ID,
                              "optimal",
                              3, 
			      "optimal", "sum", "noappend");
        
      cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
      cpl_parameterlist_append(parameters, p);
      cpl_free(full_name);


      name = "merge_delt1";
      full_name = uves_sprintf("%s.%s", UVES_REDUCE_ID, name);
        
      uves_parameter_new_range(p, full_name,
			       CPL_TYPE_DOUBLE,
                               "Order merging left hand (short wavelength) "
                               "cut. To reduce the amount of order "
                               "overlapping regions we allow to cut short and "
                               "long wavelength ranges. "
                               "This may reduce the ripple possibly "
                               "introduced by the order merging. "
                               "Suggested values are: "
                               "10 (W<=390), 12 (390<W<=437, 520<W<=564), "
			       "14 (437<W<=520, 564<W) ",
                               UVES_REDUCE_ID,
                               0.,0.,20.);
        
      cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
      cpl_parameterlist_append(parameters, p);
      cpl_free(full_name);


      name = "merge_delt2";
      full_name = uves_sprintf("%s.%s", UVES_REDUCE_ID, name);
        
      uves_parameter_new_range(p, full_name,
			       CPL_TYPE_DOUBLE,
                               "Order merging right hand (long wavelength) "
                               "cut. To reduce the amount of order "
                               "overlapping regions we allow to cut short and "
                               "long wavelength ranges. "
                               "This may reduce the ripple possibly "
                               "introduced by the order merging. "
                               "Suggested values is 4",
                               UVES_REDUCE_ID,
                               0.,0.,20.);
        
      cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
      cpl_parameterlist_append(parameters, p);
      cpl_free(full_name);



    }

    if (cpl_error_get_code() != CPL_ERROR_NONE)
    {
        cpl_msg_error(__func__, "Creation of reduction parameters failed: '%s'", 
              cpl_error_get_where());
        cpl_parameterlist_delete(parameters);
        return NULL;
    }

    
    return parameters;
}

/*----------------------------------------------------------------------------*/
/**
  @brief    Reduce a science frame
  @param    raw_image          the image to reduce (pre and overscan 
                               are already removed)
  @param    raw_header         FITS header of raw image
  @param    rotated_header     Header describing the geometry of the 
                               raw image after
                               rotation and removal of pre- and overscan areas
  @param    master_bias        The master bias image. If NULL, no bias subtraction is done.
  @param    master_dark        The master dark image. If NULL, no dark subtraction is done.
  @param    mdark_header       The master dark header containing the master
                               dark exposure time.
  @param    master_flat        Master flat image. May be NULL, but only if no
                               flatfielding is requested.
  @param    mflat_header       The master flat header. Contains the number of 
                               input flat frames which
                               is needed to generate the master flat noise image.
  @param    ordertable         Order table describing the order locations on the raw image
  @param    order_locations    The polynomial describing the order positions
  @param    linetable          Length 3 array of linetable for sky, object, sky.
  @param    linetable_header   Length 3 array of linetable headers for sky, object, sky.
  @param    dispersion_relation Length 3 array of dispersion relations for sky, object, sky.
  @param    chip               The CCD chip id
  @param    debug_mode              If true, intermediate results are saved to the
                               current directory
  @param    parameters         Parameters used for background subtraction, flat-fielding,
                               extraction and rebinning.
                               Add parameters to the caller recipe by calling
                               @c uves_propagate_parameters_step() from
                   the @c recipe_create() function
  @param    rec_id             Use @em parameters belonging to this recipe
  @param    mode               mode to set different default params (efficiencY)
  @param    x                  (Output) In 2d extraction mode, the extracted spectrum
  @param    x_header           (Output) In 2d extraction mode, header of extracted spectrum
  @param    fx                 (Output) In 2d extraction mode, the flat-fielded, 
                               extracted spectrum
  @param    cosmic_mask        (Output) In optimal extraction mode, the list of hot pixels
  @param    background         (Output) The background that was subtracted from the raw image
  @param    flatfielded_variance (Output) If non-NULL variance of extracted+flatfielded
                                 (not resampled) object
  @param    flatfielded_variance_header (Output) If non-NULL header of @em flatfielded_variance
  @param    resampled_spectrum (Output) Extracted + resampled but not FF object. 
                               If input is NULL, this
                               is not produced
  @param    resampled_mf       (Output) Extracted + resampled master flat. If input is NULL, 
                               this is not produced
  @param    merged_sky         (Output) Input pointer must be non-NULL. Returned is the sky
                               spectrum, or NULL if no sky was extracted (image
                               slicer or 2d modes)
  @param    rebinned_spectrum  (Output) The resampled 2d spectrum (non sky-subtracted)
  @param    rebinned_noise     (Output) Noise (1 sigma) of @em rebinned_spectrum
  @param    rebinned_header    (Output) Header describing geometry of @em rebinned_spectrum
                               (and of course @em rebinned_noise)
  @param    merged_spectrum    (Output) The merged spectrum (non sky-subtracted)
  @param    merged_noise       (Output) Noise (1 sigma) of @em merged_spectrum
  @param    merged_header      (Output) Header belonging to @em merged_spectrum
                                and @em merged_noise
  @param    reduced_rebinned_spectrum (Output) The sky subtracted spectrum rebinned to
                                      (wavelength, order) space (i.e. non-merged).
  @param    reduced_rebinned_noise (Output) Noise (1 sigma) of @em reduced_rebinned_spectrum
  @param    reduced_spectrum   (Output) Merged, sky-subtracted spectrum
  @param    reduced_noise      (Output) Noise (1 sigma) of @em reduced_spectrum
  @param    info_tbl           (output) table with information on object location, FWHM, and S/N
  @param    extraction_slit    (output) extraction slit actually used
  @param    order_trace        (Output) In optimal extraction mode, table describing the
                               measured spatial profile
  @return   CPL_ERROR_NONE iff okay.

  This function pre-processes the input image:
  - subtracts master bias if provided,
  - subtracts master dark if provided,
  - subtracts background (see @c uves_backsub_spline()),
  - gets extraction slit length
  - divides by master flat if @em FF_METHOD is FF_PIXEL.
  
  Then, if optimal extraction:
  - extracts (including sky subtraction), optionally flat-field corrects and merges
    the object window (see @c extract_ff_rebin_merge()).

  If linear/average extraction,
  - extracts, optionally flat-field corrects and merges the object and both 
    sky windows, then
  - subtracts the sky spectra from the object spectrum (after normalizing 
    to the object slit length)

  If 2d extraction,
  - extracts, optionally flat-field corrects and merges each of the spatial traces
  (see @c extract_ff_rebin_merge()). In 2d extraction mode, no sky subtraction can
  be performed as the object position is not defined.
**/
/*----------------------------------------------------------------------------*/

cpl_error_code uves_reduce(const cpl_image *raw_image, 
                           const uves_propertylist *raw_header, 
                           const uves_propertylist *rotated_header,
                           const cpl_image *master_bias,
                           const uves_propertylist *mbias_header,
                           const cpl_image *master_dark,
                           const uves_propertylist *mdark_header, 
                           const cpl_image *master_flat,
                           const uves_propertylist *mflat_header,
                           const cpl_table *ordertable,
                           const polynomial *order_locations,
                           const cpl_table *linetable[3],
                           const uves_propertylist *linetable_header[3], 
                           const polynomial *dispersion_relation[3],
                           enum uves_chip chip,
               /* General */
               bool   debug_mode,
               /* Backsub */
               /* Flat fielding */
               /* Extraction */
               /* Rebinning  */
               const cpl_parameterlist *parameters, 
                           const char *rec_id,
                  const char *mode,
               /* Output */
               cpl_image **x, uves_propertylist **x_header,
               cpl_image **fx,
               cpl_table **cosmic_mask,
               cpl_image **wave_map,
               cpl_image **background,
               cpl_image **flatfielded_variance,
               uves_propertylist **flatfielded_variance_header,
               cpl_image **resampled_spectrum,
               cpl_image **resampled_mf,
               cpl_image **merged_sky,
               /* Before sky-subtraction */
               cpl_image **rebinned_spectrum, 
                           cpl_image **rebinned_noise, 
                           uves_propertylist **rebinned_header,
               cpl_image **merged_spectrum  , 
                           cpl_image **merged_noise, 
                           uves_propertylist **merged_header,
               /* After sky-subtraction */
               cpl_image **reduced_rebinned_spectrum, 
                           cpl_image **reduced_rebinned_noise,
               cpl_image **reduced_spectrum         , 
                           cpl_image **reduced_noise, 
                           cpl_table **info_tbl,
               double *extraction_slit,
                           cpl_table **order_trace)
{
    /* Recipe parameters */
  char context[80];
    flatfielding_method ff_method;
    merge_method m_method;
    extract_method ex_method;
    bool blaze_corr=false;
    bool sky_sub;
    bool tilt_corr;
    double full_slit;
    double obj_slit;
    double obj_offset;
    double exptime=0;
    cpl_image *back_subbed         = NULL;         /* Image before extraction */
    cpl_image *backsubbed_noise    = NULL;

    cpl_image *mflat_noise         = NULL;         /* Master flat noise */

    cpl_image *simple_extracted    = NULL;         /* Needed only for blaze-correction */
    cpl_image *simple_extracted_mf = NULL; 

    cpl_image *sky_lo                = NULL;       /* Merged sky spectra */
    cpl_image *sky_lo_noise          = NULL;
    cpl_image *sky_hi                = NULL;
    cpl_image *sky_hi_noise          = NULL;
    uves_propertylist *sky_lo_header  = NULL;
    uves_propertylist *sky_hi_header  = NULL;

    cpl_image *sky_lo_rebinned       = NULL;       /* Rebinned sky spectra */
    cpl_image *sky_lo_rebinned_noise = NULL;
    cpl_image *sky_hi_rebinned       = NULL;
    cpl_image *sky_hi_rebinned_noise = NULL;
    uves_propertylist *sky_lo_rebinned_header = NULL;
    uves_propertylist *sky_hi_rebinned_header = NULL;

    char *subcontext                 = NULL;
    double header_full_slit;                       /* Slit length in pixels
                              from FITS header */
    char *ex_method_string           = NULL;
    double dnoise=-999;
    double bnoise=-999;
    bool has_fnoise=false;
    double fnoise=0;
    uves_propertylist* local_raw_header=NULL;
    /* Check, initialize input */
    passure( background                 != NULL, " "); *background                = NULL;
    /* resampled_spectrum, resampled_mf may be NULL */
    passure( rebinned_spectrum          != NULL, " "); *rebinned_spectrum         = NULL;
    passure( rebinned_noise             != NULL, " "); *rebinned_noise            = NULL;
    passure( rebinned_header            != NULL, " "); *rebinned_header           = NULL;
    passure( merged_spectrum            != NULL, " "); *merged_spectrum           = NULL;
    passure( merged_sky                 != NULL, " "); *merged_sky                = NULL;
    passure( merged_noise               != NULL, " "); *merged_noise              = NULL;
    passure( merged_header              != NULL, " "); *merged_header             = NULL;
    passure( reduced_rebinned_spectrum  != NULL, " "); *reduced_rebinned_spectrum = NULL;
    passure( reduced_rebinned_noise     != NULL, " "); *reduced_rebinned_noise    = NULL;
    passure( reduced_spectrum           != NULL, " "); *reduced_spectrum          = NULL;
    passure( reduced_noise              != NULL, " "); *reduced_noise             = NULL;

    passure( (flatfielded_variance == NULL) == (flatfielded_variance_header == NULL),
         "%d %d", flatfielded_variance == NULL, flatfielded_variance_header == NULL);

    assure_nomsg( extraction_slit != NULL, CPL_ERROR_NULL_INPUT );

    /* Get flat-fielding/extract method (recipe parameters) 
       These parameters determine the overall reduction strategy. */
    {
     if(strcmp(mode,".efficiency")==0) {
       sprintf(context,"%s%s",rec_id,mode);
     } else {
       sprintf(context,"%s",rec_id);
     }
    check( ff_method = uves_get_flatfield_method(parameters, context, UVES_REDUCE_ID),
           "Could not read flat-fielding method");
     
    assure( ff_method == FF_NO || master_flat != NULL, CPL_ERROR_NULL_INPUT,
        "Flat-fielding requested, but no flat field provided");
    
    /* Read extract method from  <context>.<uves_reduce>.<extract>.method  */
    check( ex_method = uves_get_extract_method(parameters, context,
                           UVES_REDUCE_ID "." UVES_EXTRACT_ID),
           "Could not get extraction method");
           
    assure( ex_method != EXTRACT_WEIGHTED, CPL_ERROR_ILLEGAL_INPUT, 
        "Weighted extraction is used only internally, "
        "as a part of optimal extraction");

    check( m_method = uves_get_merge_method(parameters, context, UVES_REDUCE_ID),
           "Could not get merging method");

/*
    check( uves_get_parameter(parameters, rec_id, UVES_REDUCE_ID,
                  "blazecorr", CPL_TYPE_BOOL, &blaze_corr),
           "Could not read parameter");
*/
    /*Forcing blazecorection to be off */
    blaze_corr=false;
    check( uves_get_parameter(parameters, rec_id, UVES_REDUCE_ID,
                  "skysub", CPL_TYPE_BOOL, &sky_sub),
           "Could not read parameter");
    
    check( uves_get_parameter(parameters, rec_id, UVES_REDUCE_ID,
                  "tiltcorr", CPL_TYPE_BOOL, &tilt_corr),
           "Could not read parameter");
    
    assure( !blaze_corr || ff_method != FF_NO, CPL_ERROR_INCOMPATIBLE_INPUT,
        "Sorry, cannot apply blaze function "
        "correction when no flatfielding is done");

    if (blaze_corr && ex_method == EXTRACT_2D)
        {
        uves_msg_warning("There will be no blaze function correction "
                 "for 2d extraction");
        }

    if (ff_method == FF_NO && m_method == MERGE_OPTIMAL)
        {
        uves_msg_warning("No flat-fielding done, but merge method = optimal. "
                 "Is this really what you want?");
        }
    if (ff_method != FF_NO && m_method == MERGE_SUM)
        {
        uves_msg_warning("Flat-fielding will be done, but merge method = sum. "
                 "Is this really what you want?");
        }

    check( uves_get_parameter(parameters, rec_id, UVES_REDUCE_ID, "slitlength", 
                  CPL_TYPE_DOUBLE, &full_slit), "Could not read parameter");

        check( uves_get_parameter(parameters, rec_id, UVES_REDUCE_ID, "objoffset",
                                  CPL_TYPE_DOUBLE, &obj_offset), 
               "Could not read parameter");
        check( uves_get_parameter(parameters, rec_id, UVES_REDUCE_ID, "objslit",
                                  CPL_TYPE_DOUBLE, &obj_slit), 
               "Could not read parameter");
    }

    /* Append '.uves_reduce' to context */
    subcontext = uves_sprintf("%s.%s", rec_id, UVES_REDUCE_ID);

    /* Subtract bias */
    check( back_subbed = cpl_image_duplicate(raw_image),
       "Error copying raw image");

    local_raw_header=(uves_propertylist* )raw_header;
    if (master_bias != NULL)
    {
        uves_msg("Subtracting master bias");
        check( uves_subtract_bias(back_subbed, master_bias),
           "Error subtracting master bias");


        int pn = PORT_ID(chip);
        char key_name[80];
        sprintf(key_name, "ESO QC OUT%d RON MASTER", pn);

        if(uves_propertylist_has(mbias_header,key_name)) {
           bnoise=uves_propertylist_get_double(mbias_header,key_name);
        } else {
           uves_msg("key %s not found, compute it on the fly",key_name);




           int kappa=5;
           int niter=3;


           for(int i=0;i<niter;i++){
                double bmin=cpl_image_get_min(master_bias);
                double bmax=cpl_image_get_max(master_bias);
                double bavg=cpl_image_get_mean(master_bias);
                double brms=cpl_image_get_stdev(master_bias);
                cpl_mask* msk=cpl_mask_threshold_image_create(master_bias,bavg-kappa*brms,bavg+kappa*brms);
                cpl_mask_not(msk);
                cpl_image_reject_from_mask(master_bias, msk);
                cpl_mask_delete(msk);
           }

           bnoise=cpl_image_get_stdev(master_bias);

        }

        uves_propertylist_append_c_double(local_raw_header,UVES_BNOISE,bnoise,
                                          "Master bias RMS on frame");




    }
    else
    {
        uves_msg("Skipping bias subtraction");
    }
    
    /* Subtract dark if available */
    if (master_dark != NULL)
    {
        uves_msg("Subtracting master dark");
        check( uves_subtract_dark(back_subbed, raw_header,
                      master_dark, mdark_header),
           "Error subtracting master dark");





        int kappa=5;
        int niter=3;


        for(int i=0;i<niter;i++){
             double dmin=cpl_image_get_min(master_dark);
             double dmax=cpl_image_get_max(master_dark);
             double davg=cpl_image_get_mean(master_dark);
             double drms=cpl_image_get_stdev(master_dark);
             cpl_mask* msk=cpl_mask_threshold_image_create(master_dark,davg-kappa*drms,davg+kappa*drms);
             cpl_mask_not(msk);
             cpl_image_reject_from_mask(master_dark, msk);
             cpl_mask_delete(msk);
        }


        dnoise=cpl_image_get_stdev(master_dark);
        uves_msg_warning("noise master dark (not rescaled) %g",dnoise);
        exptime=uves_pfits_get_exptime(mdark_header);
        uves_msg_warning("Exptime master dark %g",exptime);

        uves_propertylist_append_c_double(local_raw_header,UVES_DNOISE,dnoise,
                                          "Master dark RMS on frame");
        uves_propertylist_append_c_double(local_raw_header,UVES_DTIME,exptime,
                                          "Master dark RMS on frame");
    }
    else
    {
        uves_msg("Skipping dark subtraction");
    }

    if (master_flat != NULL)
    {
       has_fnoise=uves_propertylist_contains(mflat_header,UVES_FNOISE);
       if(has_fnoise) {
         fnoise=uves_propertylist_get_double(mflat_header,UVES_FNOISE);
       }
    }
    if (debug_mode) check( uves_save_image_local("Bias/dark subtracted raw image", "pre",
                        back_subbed, chip, -1, -1, rotated_header, true), 
              "Error saving image");

    uves_msg("Creating noise image");

    /* Define/initialize input image noise (r.o.n. and photonic) */
    check( backsubbed_noise = uves_define_noise(back_subbed, raw_header, 
                                                1, chip),
       "Could not calculate noise image");    

    /* Save noise image */
    if (debug_mode) check( uves_save_image_local("Background subtracted raw image noise",
                        "errb", backsubbed_noise,
                        chip, -1, -1, rotated_header, true),
              "Error saving image"); 

    /* Subtract background  */
    uves_msg("Subtracting inter-order background");
    
    check( uves_backsub_spline(back_subbed, raw_header,
                   ordertable, order_locations,
                   parameters, subcontext,
                   chip,
                   false,           /* Use flat-field parameters? */
                   background),
       "Error subtracting background");

    /* Save bias, dark, background subtracted frame */
    if (debug_mode) check( uves_save_image_local("Background subtracted raw image", "b",
                        back_subbed, chip, -1, -1, rotated_header, true),
              "Error saving image");


    /* 
     * Initialize flat-field noise (if necessary)
     */
    if (ff_method == FF_NO)
    {
        uves_msg("Skipping flat-field correction");
    }

    if (ff_method != FF_NO || resampled_mf != NULL)
    {
        int mflat_datancom;

        /* Save master flat image */
        if (debug_mode)
        {
            check( uves_save_image_local("Master flat image", "mf", 
                         master_flat,
                         chip, -1, -1, rotated_header, true), 
               "Error saving master flat image");
        }
        
        
        /* Define master flat noise */
        check( mflat_datancom  = uves_pfits_get_datancom(mflat_header),
           "Error reading number of raw flat field frames "
           "used for master flat image");

        uves_msg("Creating master flat noise image");
        
        check( mflat_noise = uves_define_noise(master_flat, mflat_header, 
                                               mflat_datancom,chip),
           "Could not define master flat noise");


        /* Save master flat noise image */
        if (debug_mode)
        {
            check( uves_save_image_local("Master flat noise", "errmf", mflat_noise,
                         chip, -1, -1, rotated_header, true), 
               "Error saving master flat image");
        }
    }


    /*
     * Get full slit length 
     */
    check( header_full_slit = uves_pfits_get_slitlength_pixels(raw_header, chip),
       "Could not read slit length");
    
    /* If user didn't specify slit length, use header value */
    if (full_slit < 0)
    {
        /* Avoid pixels at the edge of the slit
         *  which are likely to be noisy
         */
        full_slit = uves_max_double(1.0, header_full_slit - 2);
    }
    else
    {
        /* Warn if user specified value is larger than header value */
        if (full_slit > header_full_slit)
        {
            uves_msg_warning("Specified full slit length (%e pixels) "
                     "is larger than input header slit "
                     "length (%e pixels)",
                     full_slit, header_full_slit);
        }
    }
    
    uves_msg("Slit length = %.2f pixels", full_slit);
    *extraction_slit = full_slit;

    if (ff_method == FF_PIXEL)
    {
        uves_msg("Dividing by normalized master flat-field (method = pixel)");
        
        check( uves_flatfielding(back_subbed, backsubbed_noise,
                     master_flat, mflat_noise),
           "Could not perform flat-fielding");
        
            /* Save flat-fielded image + noise */
        if (debug_mode) 
        {
            check( uves_save_image_local("Flat-fielded image", "fb", 
                         back_subbed, chip, -1, -1, 
                         rotated_header, true),
               "Error saving flat-fielded image");
            
            check( uves_save_image_local("Flat-fielded image noise", "errfb", 
                         backsubbed_noise, chip, -1, -1,
                         rotated_header, true),
               "Error saving noise of flat-fielded image");
        }
    }
    
    /* Extract the object window (+ sky windows depending on method) */
    switch(ex_method)
    {
    case EXTRACT_OPTIMAL:
    {
        int window_number = 2;

        check( extract_ff_rebin_merge(back_subbed, 
                                      backsubbed_noise,
                                      raw_header,
                                      master_flat, 
                                      mflat_noise,
                                      ordertable, 
                                      order_locations,
                                      linetable[window_number-1],
                                      linetable_header,
                                      dispersion_relation,
                                      full_slit, 
                                      obj_offset,
                                      window_number,
                                      chip,
                                      blaze_corr,
                                      tilt_corr,
                                      debug_mode,
                                      parameters, 
                                      subcontext,
				      mode,
                                      ff_method,
                                      ex_method,
                                      m_method,
                                      NULL, 
                                      NULL, 
                                      NULL,
                                      cosmic_mask,
				      wave_map, 
                                      flatfielded_variance,
                                      flatfielded_variance_header,
                                      resampled_spectrum,
                                      resampled_mf,
                                      merged_sky,
                                      /* merged_sky will be computed 
                                         during optimal extraction */
                                      rebinned_spectrum, 
                                      rebinned_noise, 
                                      rebinned_header,
                                      merged_spectrum, 
                                      merged_noise, 
                                      merged_header,
                                      info_tbl,
                                      order_trace),
               "Error during reduction");
        
        /* The sky-subtracted spectra are just the optimally extracted spectra
         * (since sky-subtraction is done during extraction)
         */
        check(( *reduced_spectrum = cpl_image_duplicate(*merged_spectrum),
            *reduced_noise    = cpl_image_duplicate(*merged_noise),
            *reduced_rebinned_spectrum = cpl_image_duplicate(*rebinned_spectrum),
            *reduced_rebinned_noise    = cpl_image_duplicate(*rebinned_noise)),
          "Error creating sky-subtracted spectra");
    }
    break;
    case EXTRACT_LINEAR: /* Same as average (pass ex-method to uves_extract) */
    case EXTRACT_AVERAGE:
    {
        /* Average/linear extraction.
         * Define sky+object+sky windows, 
         * extract, rebin, merge, subtract
         */
        
        const char *slicer_name;
        double doffset = 0;
        double obj_hi, obj_lo;
        double sky_lo_slit, sky_hi_slit;
        int window_number;
        
        /*
         *  See if there's an image slicer
         *  Extract sky only if not
         */
        
        check( slicer_name = uves_pfits_get_slit1_name(raw_header),
           "Could not read slicer id");
          
        uves_msg("Slicer name = '%s'%s", slicer_name, 
             (strcmp(slicer_name, "FREE") == 0) ? " (no slicer)" : "");
          
        if ( strncmp(slicer_name, "SLIC", 4) == 0)
        {
            /*
             *    Use full slit for object, no sky
             */
              
                    obj_hi = uves_min_double(+full_slit/2, obj_offset + full_slit/2);
                    obj_lo = uves_max_double(-full_slit/2, obj_offset - full_slit/2);
                    
                    obj_slit = obj_hi - obj_lo;
              
            sky_lo_slit = -1;        /* Don't extract sky */
            sky_hi_slit = -1;
              
            uves_msg("Extraction slits (full slit = %.2f pixels)", full_slit);
            uves_msg("|* Sky 1 *|******** Obj ********|* Sky 2 *|");
            uves_msg("|* %-5.1f *|******* %-5.1f *******|* %-5.1f *|",
                 0.0, obj_slit, 0.0);
        }
        else
        {
            /* There's no slicer */
            assure( strncmp(slicer_name, "FREE", 4) == 0, CPL_ERROR_UNSUPPORTED_MODE,
                "Unrecognized slicer name: '%s'. "
                "Recognized names include 'FREE', 'SLIC#1', 'SLIC#2', 'SLIC#3'.",
                slicer_name);
            
            /* Measure offset if user didn't specify */
            if (obj_slit < 0)
            {
                check( obj_offset = 
                   get_offset(back_subbed, 
                          ordertable, order_locations,
                          full_slit/2, /* Offset search range */
                          10,          /* Samples per order */
                          &doffset),
                   "Could not find object offset");
                
                uves_msg("Measured object position = %.2f +- %.2f pixels", 
                     obj_offset, doffset);

                if (sky_sub)
                {
                    /* Define object extraction slit length 
                       as half of full slit. */
                    obj_hi = uves_min_double(+full_slit/2, 
                                 obj_offset + full_slit/4.0);
                    obj_lo = uves_max_double(-full_slit/2, 
                                 obj_offset - full_slit/4.0);
                }
                else
                /* No sky subtraction. Object = full slit */
                {
                    obj_hi = uves_min_double(+full_slit/2, 
                                 obj_offset + full_slit/2.0);
                    obj_lo = uves_max_double(-full_slit/2, 
                                 obj_offset - full_slit/2.0);
                }
                obj_slit = obj_hi - obj_lo;
            }
            else
            /* User specified object slit */
            {
                uves_msg("Object offset = %.2f pixels", obj_offset);

                obj_hi = obj_offset + obj_slit / 2;
                obj_lo = obj_offset - obj_slit / 2;
            }
            
            uves_msg("Object slit = %.2f pixels", obj_slit);
            
            assure( -full_slit / 2 < obj_offset && obj_offset < full_slit / 2, 
                CPL_ERROR_ILLEGAL_INPUT,
                "Object is outside slit! Offset = %f, Slit length = %f",
                obj_offset, full_slit);
            
            /* Sky slits (might be negative if object has large offset) */
            if (sky_sub)
            {
                sky_lo_slit = obj_lo - (-full_slit/2);
                sky_hi_slit = full_slit/2 - obj_hi;
                
                assure( sky_lo_slit > 0 || sky_hi_slit > 0, CPL_ERROR_ILLEGAL_INPUT,
                    "At least one sky slit length must be positive. "
                    "They are %f and %f pixels", sky_lo_slit, sky_hi_slit);
            }
            else
            {
                sky_lo_slit = -1; /* Don't extract sky */
                sky_hi_slit = -1;
            }
            
            uves_msg("Extraction slits (full slit = %.2f pixels)", full_slit);
            uves_msg("|*** Sky 1 **%s|**** Obj ****|%s** Sky 2 ***|", 
                 (obj_lo > -obj_hi) ? "*" : "",
                 (obj_lo > -obj_hi) ? ""  : "*");
            uves_msg("|*** %-5.1f **%s|*** %-5.1f ***|%s** %-5.1f ***|",
                 sky_lo_slit, (obj_lo > -obj_hi) ? "*" : "",
                 obj_slit   , (obj_lo > -obj_hi) ? ""  : "*",
                 sky_hi_slit);
        }
        
        /* The window geometry has now been deermined. Extract spectra.

           It is important to use the same rebinning step size,
           for sky and object (otherwise the sky spectrum cannot (easily)
           be subtracted). If this step size is not specified (i.e. is negative)
           in the parameter list, it is determined from the 
           average pixelsize, which is read from the line table.
           Therefore pass the object's line table also for the sky windows
           (but still use different dispersion relations for sky/object)
        */
        
        /* Extract sky 1 */
        window_number = 1;
        if ( sky_lo_slit > 0 )
        {
            uves_msg("Processing sky 1 window");
            check( extract_ff_rebin_merge(back_subbed, 
                                          backsubbed_noise,
                                          raw_header,
                                          master_flat, 
                                          mflat_noise,
                                          ordertable, 
                                          order_locations,
                                          linetable[2-1],   /* Object linetable */
                                          linetable_header,
                                          dispersion_relation,
                                          sky_lo_slit,     /* Slit length (pixels) */
                                          -full_slit/2 + sky_lo_slit/2,
                                          /* Slit center offset */
                                          window_number,
                                          chip,
                                          blaze_corr,
                                          tilt_corr,
                                          debug_mode,
                                          parameters, 
                                          subcontext,
					  mode,
                                          ff_method,
                                          ex_method,
                                          m_method,
                                          NULL, 
                                          NULL, 
                                          NULL, 
                                          NULL,
                                          wave_map, 
                                          NULL,
					  NULL,
                                          NULL, 
                                          NULL, 
                                          NULL,
                                          &sky_lo_rebinned, 
                                          &sky_lo_rebinned_noise, 
                                          &sky_lo_rebinned_header,
                                          &sky_lo, 
                                          &sky_lo_noise, 
                                          &sky_lo_header,
                                          NULL, 
                                          NULL),
               "Error processing lower sky window");
        }
        else
        {
            uves_msg("Skipping sky 1 window");
            sky_lo_rebinned = NULL;
            sky_lo = NULL;
        }
        
        /* Extract sky 2 */
        window_number = 3;
        if ( sky_hi_slit > 0 )
        {
            uves_msg("Processing sky 2 window");
              
            uves_free_propertylist(rebinned_header);
            check( extract_ff_rebin_merge(back_subbed, 
                                          backsubbed_noise,
                                          raw_header,
                                          master_flat, 
                                          mflat_noise,
                                          ordertable, 
                                          order_locations,
                                          linetable[2-1], /* Object linetable */
                                          linetable_header, 
                                          dispersion_relation,
                                          sky_hi_slit,    /* Slit length (pixels) */
                                          full_slit/2 - sky_hi_slit/2,
                                          /* Slit center offset      */
                                          window_number,
                                          chip,
                                          blaze_corr,
                                          tilt_corr,
                                          debug_mode,
                                          parameters, 
                                          subcontext,
					  mode,
                                          ff_method,
                                          ex_method,
                                          m_method,
                                          NULL, 
                                          NULL, 
                                          NULL,
                                          NULL,
                                          wave_map, 
                                          NULL,
                                          NULL,
					  NULL, 
                                          NULL, 
                                          NULL,
                                          &sky_hi_rebinned, 
                                          &sky_hi_rebinned_noise, 
                                          &sky_hi_rebinned_header,
                                          &sky_hi, 
                                          &sky_hi_noise, 
                                          &sky_hi_header,
                                          NULL,
                                          NULL),
                   "Error processing upper sky window");
        }
        else
        {
            uves_msg("Skipping sky 2 window");
            sky_hi_rebinned = NULL;
            sky_hi = NULL;
        }
          
        /* Extract object */
        window_number = 2;
        uves_msg("Processing object window");
        uves_free_propertylist(rebinned_header);
        check( extract_ff_rebin_merge(back_subbed, backsubbed_noise, raw_header,
                                      master_flat, mflat_noise,
                                      ordertable, order_locations,
                                      linetable[window_number-1], 
                                      linetable_header, 
                                      dispersion_relation,
                                      obj_slit,   /* Slit length (pixels) */
                                      obj_offset,   
                                      /* Slit center offset */
                                      window_number,
                                      chip,
                                      blaze_corr,
                                      tilt_corr,
                                      debug_mode,
                                      parameters, 
                                      subcontext,
				      mode,
                                      ff_method,
                                      ex_method, 
                                      m_method,
                                      NULL, 
                                      NULL, 
				      NULL, 
                                      NULL,
                                      wave_map,
                                      flatfielded_variance,
                                      flatfielded_variance_header,
                                      resampled_spectrum, 
                                      resampled_mf, 
                                      NULL,
                                      rebinned_spectrum, 
                                      rebinned_noise, 
                                      rebinned_header,
                                      merged_spectrum, 
                                      merged_noise, 
                                      merged_header,
                                      info_tbl,
                                      NULL),
               "Error processing object window");

        if (info_tbl != NULL && *info_tbl != NULL)
        {
            /* Compute obj. position from sky_lo_slit
               for consistency with optimal extraction */
            int i;
            for (i = 0; i < cpl_table_get_nrow(*info_tbl); i++)
            {
                cpl_table_set_double(*info_tbl, "ObjPosOnSlit", i,
                         cpl_table_get_double(*info_tbl, "ObjPosOnSlit", i, NULL)
                         + 
                         ((sky_lo_slit >= 0) ? sky_lo_slit : 0));

            }
        }
                
          
        /* Now subtract sky from both rebinned spectrum and merged spectrum */
          
        /* Duplicate, then subtract */
        
        /* 1d spectrum */
        check(( *reduced_spectrum = cpl_image_duplicate(*merged_spectrum),
            *reduced_noise    = cpl_image_duplicate(*merged_noise)),
          "Error allocating sky-subtracted spectra");
        
        /* 2d (wavelength, order) spectrum */
        check(( *reduced_rebinned_spectrum = 
            cpl_image_duplicate(*rebinned_spectrum),
            *reduced_rebinned_noise    = 
            cpl_image_duplicate(*rebinned_noise)),
          "Error allocating sky-subtracted spectra");
        
        if (sky_lo != NULL || sky_hi != NULL)
        {
            uves_msg("Subtracting sky");
            
            check( *merged_sky = 
               subtract_sky(*reduced_rebinned_spectrum, 
                    *reduced_rebinned_noise, *rebinned_header,
                    sky_lo_rebinned, sky_lo_rebinned_noise, 
                    sky_lo_rebinned_header,
                    sky_hi_rebinned, sky_hi_rebinned_noise,
                    sky_hi_rebinned_header,
                    reduced_spectrum, reduced_noise, *merged_header,
                    sky_lo, sky_lo_noise, sky_lo_header,
                    sky_hi, sky_hi_noise, sky_hi_header,
                    (ex_method == EXTRACT_AVERAGE) ? 1.0 : obj_slit,
                    (ex_method == EXTRACT_AVERAGE) ? 1.0 : sky_lo_slit,
                    (ex_method == EXTRACT_AVERAGE) ? 1.0 : sky_hi_slit),
               "Could not subtract sky");

            if (*merged_sky == NULL)
            {
                uves_msg_warning("Could not subtract sky");
            }
        }
        else
        {
            uves_msg_low("Skipping sky subtraction");
              
            /* Return no sky spectrum */
            *merged_sky = NULL;
        }
    } /* Simple extraction */
    break;
    case EXTRACT_2D:
    {
        int window_number = 2;    /* Use middle line table for entire
                     slit length (like MIDAS) */

        int half_slit_length;     /* The slit length is
                     2*half_slit_length = an even number */
        
        /* Round to nearest integer, remove (noisy) edge (~2 pixels) */
        half_slit_length = 
        uves_max_int(1, uves_round_double(header_full_slit/2) - 1);

        check( extract_ff_rebin_merge(back_subbed, 
                                          backsubbed_noise,
                      raw_header,
                      master_flat, 
                                          mflat_noise,
                      ordertable, 
                                          order_locations,
                      linetable[window_number-1],
                      linetable_header,
                      dispersion_relation,
                      2*half_slit_length, 
                                          0, /* offset is not used when method=2d */
                      window_number,
                      chip,
                      blaze_corr,
                                          tilt_corr,
                                          debug_mode,
                      parameters, 
                                          subcontext,
				      mode,
                      ff_method,
                      ex_method,
                      m_method,
                      x, 
                                          x_header, 
                                          fx,      /* 2d-ex. output      */
				      NULL,  NULL,
                                          NULL,    /* Optimal-ex. output */
                      NULL, 
                                          NULL,    /* Don't want resampled_spectrum,
                              resampled_mf */
                      NULL,    /* Don't want sky spectrum */
                                          NULL,    /* Don't want 
                              flatfielded_variance+header */
                      rebinned_spectrum, 
                                          rebinned_noise, 
                                          rebinned_header,
                      merged_spectrum, 
                                          merged_noise, 
                                          merged_header, 
                                          info_tbl,
                                          NULL),
           "Error during reduction");

            if (x_header != NULL)
                {
                    uves_pfits_set_hs(*x_header, 
                                      uves_round_double(2*half_slit_length));
                }
            if (rebinned_header != NULL)
                {
                    uves_pfits_set_hs(*rebinned_header, 
                                      uves_round_double(2*half_slit_length));
                }
            if (merged_header != NULL)
                {
                    uves_pfits_set_hs(*merged_header, 
                                      uves_round_double(2*half_slit_length));
                }


        /* No sky-subtraction done. Just copy the merged spectra
         * to get the 'reduced' (i.e. sky-subtracted) spectra
         */
        check(( *reduced_spectrum = cpl_image_duplicate(*merged_spectrum),
            *reduced_noise    = cpl_image_duplicate(*merged_noise),
            *reduced_rebinned_spectrum = cpl_image_duplicate(*rebinned_spectrum),
            *reduced_rebinned_noise    = cpl_image_duplicate(*rebinned_noise)),
          "Error allocating reduced spectra");
    }
    break;
    default:
        assure( false, CPL_ERROR_ILLEGAL_INPUT, 
            "Unknown extraction method: %d", ex_method);
        break;
    } /* switch extraction method optimal/simple/2d */
    
  cleanup:
    uves_free_image(&back_subbed);
    uves_free_image(&backsubbed_noise);
    uves_free_image(&mflat_noise);
    uves_free_image(&simple_extracted);
    uves_free_image(&simple_extracted_mf);
    uves_free_image(&sky_lo);
    uves_free_image(&sky_lo_noise);
    uves_free_image(&sky_hi);
    uves_free_image(&sky_hi_noise);
    uves_free_image(&sky_lo_rebinned);
    uves_free_image(&sky_lo_rebinned_noise);
    uves_free_image(&sky_hi_rebinned);
    uves_free_image(&sky_hi_rebinned_noise);
    uves_free_propertylist(&sky_lo_header);
    uves_free_propertylist(&sky_hi_header);
    uves_free_propertylist(&sky_lo_rebinned_header);
    uves_free_propertylist(&sky_hi_rebinned_header);

    cpl_free(subcontext);
    cpl_free(ex_method_string);

    if (cpl_error_get_code() != CPL_ERROR_NONE)
    {
        uves_free_image(background);
        uves_free_image(flatfielded_variance);
        uves_free_propertylist(flatfielded_variance_header);
        uves_free_image(resampled_spectrum);
        uves_free_image(resampled_mf);
        uves_free_image(merged_sky);
        uves_free_image(rebinned_spectrum);
        uves_free_image(rebinned_noise);
        uves_free_propertylist(rebinned_header);
        
        uves_free_image(merged_noise);
        uves_free_image(merged_spectrum);
        uves_free_propertylist(merged_header);
    }

    return cpl_error_get_code();
}


/*----------------------------------------------------------------------------*/
/**
  @brief    1st order line tilt correction
  @param    dispersion_solution       calibration solutions
  @param    linetable_header          describing where dispersion 
                                      solutions were obtained
  @param    objoffset                 object slit position
  @return   interpolated dispersion solution at the specified object position

**/
/*----------------------------------------------------------------------------*/
static polynomial *
interpolate_wave(const polynomial *dispersion_relation[3],
                 const uves_propertylist *linetable_header[3],
                 double objoffset)
{
    polynomial *dispersion = NULL;
    polynomial *q1 = NULL;
    polynomial *q2 = NULL;
    cpl_table *offset = cpl_table_new(3);
    int ilow, ihigh;
    double offset1, offset2;
    
    /* We need the sort pattern. Use a table for that */
    cpl_table_new_column(offset, "Index", CPL_TYPE_INT);
    cpl_table_new_column(offset, "Offset", CPL_TYPE_DOUBLE);
    
    {
        int i;
        bool reverse;
        for (i = 0; i < 3; i++) {
            cpl_table_set_int(offset, "Index", i, i);
            cpl_table_set_double(offset, "Offset", i,
                                 uves_pfits_get_offset(linetable_header[i]));
            
            uves_msg_debug("Wavecal %d offset is %f pixels", i, 
                           cpl_table_get_double(offset, "Offset", i, NULL));
        }

        reverse = false;
        uves_sort_table_1(offset, "Offset", reverse);
    }
    
    /* Find indices of the two dispersion solutions neares to the object position */
    if (objoffset <= cpl_table_get_double(offset, "Offset", 1, NULL))
        {
            ilow    = cpl_table_get_int(offset, "Index", 0, NULL);
            ihigh   = cpl_table_get_int(offset, "Index", 1, NULL);
            offset1 = cpl_table_get_double(offset, "Offset", 0, NULL);
            offset2 = cpl_table_get_double(offset, "Offset", 1, NULL);
        }
    else
        {
            ilow  = cpl_table_get_int(offset, "Index", 1, NULL);
            ihigh = cpl_table_get_int(offset, "Index", 2, NULL);
            offset1 = cpl_table_get_double(offset, "Offset", 1, NULL);
            offset2 = cpl_table_get_double(offset, "Offset", 2, NULL);
        }

    uves_msg("Interpolating dispersion relation at offset = %.2f",
             objoffset);

    uves_msg_debug("Using previous solutions at %.2f and %.2f pixels",
                   offset1, offset2);
    
    /* Fail cleanly if 2 dispersion solution were obtained at the same offset
       (rather than silently dividing by zero) */
    assure( offset1 < offset2,
            CPL_ERROR_DIVISION_BY_ZERO,
            "Dispersion solution %d offset = %.2f pixels; "
            "dispersion solution %d offset = %.2f pixels; cannot extrapolate",
            ilow, offset1,
            ihigh, offset2);
    
    /* Do simple linear interpolation =
           p = a p1 + b p2
       where
          a = (offset2 - offset) / (offset2 - offset1)
          b = (offset1 - offset) / (offset1 - offset2)

       which corrects for any line tilt to 1st order.

       A 2nd order line tilt correction (distortions) is probably overkill
       because of UVES' short slit and very straight arclines.
    */
    {
        double a = (offset2 - objoffset) / (offset2 - offset1);
        double b = (offset1 - objoffset) / (offset1 - offset2);
        
        q1 = uves_polynomial_duplicate(dispersion_relation[ilow]); 
        uves_polynomial_rescale(q1, 0, a); 
        /* q1 = a p1 */
        
        q2 = uves_polynomial_duplicate(dispersion_relation[ihigh]); 
        uves_polynomial_rescale(q2, 0, b); 
        /* q2 = b p2 */
        
        dispersion = uves_polynomial_add_2d(q1, q2);
    }
    
  cleanup:
    uves_free_table(&offset);
    uves_polynomial_delete(&q1);
    uves_polynomial_delete(&q2);
    return dispersion;
}


/*----------------------------------------------------------------------------*/
/**
  @brief    Reduce one extraction window
  @param    back_subbed              The de-biased, dark subtracted, 
                                     background subtracted, possibly already 
                     flat-fielded image in (pixel, pixel)-space
  @param    backsubbed_noise         Noise of @em back_subbed
  @param    backsubbed_header        FITS header of @em back_subbed
  @param    master_flat              The master flat frame. May be NULL if
                                     flatfielding is not performed.
  @param    mflat_noise              Noise of @em master_flat
  @param    ordertable               Order table describing the order locations 
                                     on the raw image
  @param    order_locations          The polynomial describing the order positions
  @param    linetable                Line table used for the current extraction window
  @param    linetable_header         Line table header containing offsets for previously
                                     obtained dispersion relations
  @param    dispersion_relation      The dispersion relations, will be interpolated
                                     depending on offset
  @param    slit_length              Length (in spatial direction, in pixels)
                                     of current extraction window.      
  @param    slit_offset              Offset of current extraction window
                                     (applies only to linear/average,
                                     not optimal/2d extraction)
  @param    window                   Extraction window number. Must be 1, 2, 3, 
                                     for lower sky, object, upper sky, respectively. 
                     Not used in optimal/2d extraction modes.
  @param    chip                     The CCD chip id
  @param    blaze_correct            Correct (empirically) for the different 
                                     shapes of object/flat-field blaze function?
  @param    tilt_corr                Correct for line tilt? If not, the appropriate
                                     dispersion solutions are used for object/sky
                                     (without interpolation)
  @param    debug_mode                    If true, intermediate results are saved to 
                                     the current directory
  @param    parameters               Parameters used for extraction and rebinning.
  @param    context                  Use @em parameters belonging to this context
  @param    ff_method                The flat-fielding method.
  @param    ex_method                Extraction method.
  @param    m_method                 The order merging method.
  @param    x                        (Output) If non-NULL, this will contain the
                                     extracted spectrum (recipe product
                     only in 2d mode)
  @param    x_header                 (Output) If non-NULL, the header
                                     belonging to @em x
  @param    fx                       (Output) If non-NULL, the extracted, 
                                     flat-fielded spectrum 
                                     (recipe product only in 2d mode)
  @param    cosmic_mask              (Output) Used only for optimal extraction
  @param    flatfielded_variance     (Output) If non-NULL, variance of ff. spectrum
  @param    flatfielded_variance_header (Output) Header of flatfielded_variance
  @param    resampled_spectrum       (Output) Extracted + rebinned. 
                                     Not ff. May be NULL
  @param    resampled_mf             (Output) Extracted + rebinned master flat-field.
                                     May be NULL
  @param    merged_sky               (Output) For optimal extraction, 
                                     the merged sky spectrum
  @param    rebinned_spectrum        (Output) Extracted, flat-fielded, 
                                     rebinned spectrum
  @param    rebinned_noise           (Output) Error (1 sigma) of @em rebinned_spectrum
  @param    rebinned_header          (Output) Header describing geometry 
                                     of @em rebinned_spectrum and @em rebinned_noise
  @param    merged_spectrum          (Output) The merged spectrum
  @param    merged_noise             (Output) Noise (1 sigma) of
                                     @em merged_spectrum
  @param    merged_header            (Output) Header belonging to @em merged_spectrum
                                     and @em merged_noise
  @param    info_tbl                 table with information on object location, FWHM, and S/N
  @return   CPL_ERROR_NONE iff okay

  As its name sort of implies, this function
  - extracts the spectrum (optimally, 2d or simple). If optimal extraction,
  also the sky is extracted.
  - divides by flat-field iff method is FF_EXTRACT. If method is FF_PIXEL
  the input is assumed to have been already flat-fielded.
  - rebins to (order, wavelength)-space the extracted spectrum/spectra.
  - merges the rebinned spectrum/spectra.

**/
/*----------------------------------------------------------------------------*/

static cpl_error_code
extract_ff_rebin_merge(cpl_image *back_subbed, 
                       cpl_image *backsubbed_noise,
               const uves_propertylist *backsubbed_header,
               const cpl_image *master_flat, 
                       cpl_image *mflat_noise,
               const cpl_table *ordertable, 
                       const polynomial *order_locations,
               const cpl_table *linetable, 
                       const uves_propertylist *linetable_header[3], 
                       const polynomial *dispersion_relation[3],
               double slit_length, 
                       double slit_offset,
               int window, 
                       enum uves_chip chip,
               bool blaze_correct,
                       bool tilt_corr,
                       bool debug_mode,
               const cpl_parameterlist *parameters, 
                       const char *context,
                       const char *mode,
               flatfielding_method ff_method,
               extract_method ex_method,
               merge_method m_method,
               /* Output */
               cpl_image **x, uves_propertylist **x_header,
               cpl_image **fx,
               cpl_table **cosmic_mask,
               cpl_image **wave_map,
               cpl_image **flatfielded_variance,
               uves_propertylist **flatfielded_variance_header,
               cpl_image **resampled_spectrum,
               cpl_image **resampled_mf,
               cpl_image **merged_sky,
               cpl_image **rebinned_spectrum, 
                       cpl_image **rebinned_noise, 
                       uves_propertylist **rebinned_header,
               cpl_image **merged_spectrum, 
                       cpl_image **merged_noise, 
                       uves_propertylist **merged_header,
                       cpl_table **info_tbl,
                       cpl_table **order_trace)
{
    cpl_image *extracted           = NULL;
    cpl_image *extracted_noff      = NULL;
    cpl_image *extracted_noise     = NULL;
    cpl_image *extracted_sky       = NULL;   /* For optimal extraction */
    cpl_image *extracted_sky_noise = NULL;
    cpl_image *blaze_ratio         = NULL;   /* The (per-order normalized)
                                                ratio of blaze functions */

    cpl_image *cosmic_image        = NULL;
    cpl_image *weights             = NULL;
    cpl_table *profile_table       = NULL;
    uves_propertylist *extracted_header = NULL;

    cpl_image *extracted_mf        = NULL;
    cpl_image *extracted_mf_noise  = NULL;

    cpl_image *rebinned_sky        = NULL;
    cpl_image *rebinned_sky_noise  = NULL;
    cpl_image *merged_sky_noise    = NULL;

    polynomial *dispersion_int     = NULL;    /* interpolated at object position */
    polynomial *dispersion_int_sky = NULL;    /* if sky was extracted simultaneously
                                                 with the object (optimal extraction)
                                                 this is the dispersion at the average
                                                 sky position */
    cpl_table *poly_table          = NULL;
    cpl_image *temp_image          = NULL;

    int n_traces;                             /* Number of traces. Equal to 1,
                         unless for 2d reduction */
    int first_abs_order, last_abs_order;
    int filename_window;                      /* The window number appended to 
                         the filenames of local products */
    //const char* recipe_id=make_str(UVES_SCIRED_ID);
   double delt1=0;
   double delt2=0;


    /* If ff_method is FF_EXTRACT, or if resampled_mf is requested,
       then master flat must be provided */
    passure((ff_method != FF_EXTRACT && resampled_mf == NULL) 
        ||
        master_flat != NULL, " ");

    /* Blaze correction only makes sense if flat-fielding,
     */
    passure( !blaze_correct || ff_method != FF_NO, " ");

    passure( ex_method != EXTRACT_OPTIMAL || merged_sky != NULL, " ");
    passure( ex_method != EXTRACT_OPTIMAL || cosmic_mask != NULL, " ");

    passure(1 <= window && window <= 3, "Illegal window: %d", window);

    passure( (x == NULL) == (x_header == NULL) &&
	     (x == NULL) == (fx == NULL), " ");





    check( uves_get_parameter(parameters, NULL, 
			      context, "merge_delt1",
			      CPL_TYPE_DOUBLE, &delt1),
	   "Could not read parameter delt1");
    
    check( uves_get_parameter(parameters, NULL, 
			      context, "merge_delt2",
			      CPL_TYPE_DOUBLE, &delt2),
	   "Could not read parameter delt2");






    if (ex_method == EXTRACT_OPTIMAL || ex_method == EXTRACT_2D)
    {
        /* Don't append window number if optimal/2d extraction.
           There's only one window in these cases, and
           it allows the response/efficiency recipe to save 
           both optimally and linearly extracted spectra
           (without overwriting). */

        filename_window = -1;  /* -1 means don't append window
                      number to filename */
    }
    else
    {
        filename_window = window;
    }

    n_traces = (ex_method == EXTRACT_2D) ? uves_round_double(slit_length) : 1;
    
    check( first_abs_order = uves_pfits_get_firstabsorder(linetable_header[0]), 
       "Could not read order numbers from line table header");
    check( last_abs_order  = uves_pfits_get_lastabsorder (linetable_header[0]), 
       "Could not read order numbers from line table header");

    if (window == 2)
    {
        uves_msg("Extracting object");
    }
    if( ff_method == FF_EXTRACT ) {
      check_nomsg(cosmic_image=uves_image_mflat_detect_blemishes(master_flat,
								 backsubbed_header));
    }

    check( extracted = 
       uves_extract(back_subbed, 
            backsubbed_noise,
            backsubbed_header,
            ordertable, 
            order_locations,
            slit_length,               /* Slit length (pixels)     */
            slit_offset,               /* Slit center offset       */
            parameters, 
            context,                   /* Extraction method, kappa */
	    mode,
            false,                     /* Don't extract partial bins */
                        debug_mode,
            chip,
            &extracted_header,         /* Spectrum header          */
            &extracted_noise,          /* Spectrum noise           */
            /* Optimal extraction parameters: */
            &extracted_sky,            /* Sky                      */
            &extracted_sky_noise,      /* Sky                      */
            cosmic_mask, 
            &cosmic_image,
            (debug_mode) ? 
            &profile_table : NULL,
            &weights,                  /* If optimal, weights
                              are defined              */
            info_tbl,
                        order_trace),
       "Error extracting spectrum");

    if (x != NULL) {
        *x = cpl_image_duplicate(extracted);
        *x_header = uves_propertylist_duplicate(extracted_header);
    }

    if (ex_method == EXTRACT_OPTIMAL) {
        uves_msg_low("%" CPL_SIZE_FORMAT " hot pixels were detected during optimal extraction",
                     cpl_table_get_nrow(*cosmic_mask));
        
        if (cpl_table_get_nrow(*cosmic_mask) > 0) {
            check( uves_plot_table(*cosmic_mask, "X", "Y", 
                                   "%" CPL_SIZE_FORMAT "hot/cold pixels",
                                   cpl_table_get_nrow(*cosmic_mask)),
                   "Plotting failed");
        }
    }
    
    /* Save extracted spectrum + noise + sky + noise, and
       if optimal: weightmap + crmask + order_trace */
    if (debug_mode) {
        /* This (bit ugly) code selects filename
         *  and description string depending on
         *  whether flat-fielding was already done
         */



      check(uves_propertylist_copy_property_regexp(extracted_header,
						   backsubbed_header,
						   "^ESO ",0),"error copying hierarch keys");


        check( uves_save_image_local((ff_method == FF_PIXEL) ? 
                                     "Extracted, flatfielded spectrum" :
                                     "Extracted spectrum",
                                     (ff_method == FF_PIXEL) ? 
                                     "xfb" : "xb",
                                     extracted, chip, -1,
                                     filename_window, extracted_header, true),
               "Error saving extracted%s spectrum", 
               (ff_method == FF_PIXEL) ? ", flatfielded" : "");
        
        check( uves_save_image_local((ff_method == FF_PIXEL) ?
                                     "Extracted, flatfielded spectrum noise" :
                                     "Extracted spectrum noise",
                                     (ff_method == FF_PIXEL) ? 
                                     "errxfb" : "errxb",
                                     extracted_noise, chip, -1, 
                                     filename_window, extracted_header, true),
               "Error saving noise of extracted%s spectrum", 
               (ff_method == FF_PIXEL) ? ", flatfielded" : "");
        
        if (extracted_sky != NULL)
            {
                check( uves_save_image_local((ff_method == FF_PIXEL) ? 
                                             "Extracted, flatfielded sky" : 
                                             "Extracted sky",
                                             (ff_method == FF_PIXEL) ? 
                                             "xfsky" : "xsky",
                                             extracted_sky, chip, -1,
                                             filename_window, extracted_header, true),
                       "Error saving extracted%s sky", 
                       (ff_method == FF_PIXEL) ? ", flatfielded" : "");
                
                check( uves_save_image_local((ff_method == FF_PIXEL) ? 
                                             "Noise of extracted, flatfielded sky" : 
                                             "Noise of extracted sky",
                                             (ff_method == FF_PIXEL) ? 
                                             "errxfsky" : "errxsky",
                                             extracted_sky_noise, chip, -1,
                                             filename_window, extracted_header, true),
                       "Error saving extracted%s sky noise", 
                       (ff_method == FF_PIXEL) ? ", flatfielded" : "");
            }
        
        if (ex_method == EXTRACT_OPTIMAL)
            {
                check( uves_save_image_local("Optimal extraction weights", 
                                             "weights",
                                             weights, chip, -1,
                                             filename_window, 
                                             backsubbed_header, true),
                       "Error saving weights map");
                
                check( uves_save_table_local("Cosmic ray table", "crmask", 
                                             *cosmic_mask, chip, -1,
                                             filename_window, 
                                             backsubbed_header, NULL),
                       "Error saving cosmic ray mask");

                check( uves_save_image_local("Cosmic ray image", "crimage",
                                             cosmic_image, chip, -1,
                                             filename_window, 
                                             backsubbed_header, true),
                       "Error saving cosmic ray mask");

                if (profile_table != NULL)
                    {
                        check( uves_save_table_local("Profile table", "profile", 
                                                     profile_table, chip, -1,
                                                     filename_window, 
                                                     backsubbed_header, NULL),
                               "Error saving profile table");
                    }
            }
    }

    /* Extract + resample master flat, only if necessary */
    if (master_flat != NULL && (ff_method == FF_EXTRACT || resampled_mf != NULL))
    {
        uves_msg("Extracting master flat field");
        
        /* Extract the master flat spectrum.
           If object was extracted with method=optimal,
           then temporarily set method=weighted */
        
        if (ex_method == EXTRACT_OPTIMAL)
        {
            const char *temp_method = "weighted";
            
            /* Cast to non-const is okay. After extraction, the 
               parameter is set to 'optimal' (see below), so there
               is not net change (unless the extraction fails, in
               which case parameter list will change).
            */
            check( uves_set_parameter((cpl_parameterlist *) parameters,
                          context, UVES_EXTRACT_ID ".method", 
                          CPL_TYPE_STRING, &temp_method),
               "Error setting extraction method to '%s'", temp_method);
        }
        
        check( extracted_mf = 
           uves_extract((cpl_image *)master_flat,  /* const-cast ok,
                                  outlier pixels are
                                  flagged bad only in
                                  optimal extraction */
                mflat_noise,
                NULL,           /* FITS header */
                ordertable, 
                order_locations,
                slit_length,    /* Slit length (pixels)     */
                slit_offset,    /* Slit center offset       */
                parameters, context, /* Extraction method, 
                            kappa */
			mode,
                false,          /* Extraction partial bins? */
                                debug_mode,
                chip,
                NULL,           /* Spectrum header          */
                &extracted_mf_noise,  /* Spectrum noise     */
                NULL, 
                NULL,           /* Sky                      */
                NULL, 
                NULL, 
                NULL,    /* Cosmic ray table/image,
                         profile     */
                &weights, /* Weights are used unchanged */
                NULL,
                                NULL),  
           "Error extracting master flat spectrum");
        
        /* Reset parameter to previous value 
           (parameter list is declared const) */
        if (ex_method == EXTRACT_OPTIMAL) {
                const char *method = "optimal";
                
                /* Cast to non-const is okay. On successful termination,
                   there is no net change in the parameter list (see above). */
                check( uves_set_parameter((cpl_parameterlist *) parameters,
                                          context, UVES_EXTRACT_ID ".method", 
                                          CPL_TYPE_STRING, &method),
                       "Error resetting extraction method to '%s'", method);
            }
        
        if (debug_mode) {
                double ff_mean;
                
                /* Save normalized master flat spectrum + noise */
                uves_free_image(&temp_image);
                
                ff_mean    = cpl_image_get_mean(extracted_mf);
                
                check( temp_image = 
                       cpl_image_divide_scalar_create(extracted_mf, ff_mean),
                       "Could not normalize master flat spectrum");
                
                check( uves_save_image_local("Normalized master flat spectrum",
                                             "xmf",
                                             temp_image, chip, -1,
                                             filename_window, extracted_header, true),
                       "Error saving image");
                
                /* Also rescale noise before saving */
                uves_free_image(&temp_image);
                check( temp_image = 
                       cpl_image_divide_scalar_create(extracted_mf_noise,
                                                      ff_mean),
                       "Could not rescale master flat noise spectrum");
                
                check( uves_save_image_local("Noise of normalized "
                                             "master flat spectrum", "errxmf",
                                             temp_image, chip, -1,
                                             filename_window, extracted_header, true),
                       "Error saving image");
                
                uves_free_image(&temp_image);
            }
        
        /* Rebin master flat */
        if (resampled_mf != NULL) {
                uves_msg("Rebinning master flat spectrum");


                /* Use dispersion solution obtained at slit center 
                 * (high accuracy is non-essential here, the resampled
                 * flat-field is not used in further processing
                 */

                check( *resampled_mf = uves_rebin(extracted_mf,
                                                  parameters, context,
                                                  linetable,
                                                  dispersion_relation[1],
                                                  first_abs_order,
                                                  last_abs_order,
                                                  n_traces,
                                                  false,
                                                  false,
                                                  rebinned_header),
                       "Error resampling master flat");
	       
 
                check( *wave_map = uves_get_wave_map(back_subbed,
						    context,
						    parameters,
						    ordertable, 
						    linetable,
						    order_locations,
						    dispersion_relation[1],
						    first_abs_order,
						    last_abs_order,
						    slit_length),
                       "Error generating wave map");
 
                check( uves_save_image_local("Wave map",
                                             "wave_map",
                                             *wave_map, chip, -1,
                                             filename_window, 
					     extracted_header, true),
                       "Error saving wavemap image");
                
 

                if (debug_mode) {

      check(uves_propertylist_copy_property_regexp(*rebinned_header,
						   backsubbed_header,
						   "^ESO ", 0),"error copying hierarch keys");


                    check( uves_save_image_local("Resampled master flat spectrum",
                                                 "wxmf", *resampled_mf, chip, -1,
                                                 filename_window, *rebinned_header, true),
                           "Error saving image");
                }
            }
    } /* Extract, rebin master flat */
 
    /* If we didn't already, divide by the flat field */
    if (ff_method == FF_EXTRACT)
    {

    uves_msg("Dividing by normalized master flat-field (method = extract)");
        
    /* Remember this for later */
    extracted_noff = cpl_image_duplicate(extracted);

    check( uves_flatfielding(extracted   , extracted_noise,
                             extracted_mf, extracted_mf_noise),
           "Could not perform flat-fielding");

    if (extracted_sky != NULL)
    {
       check( uves_flatfielding(extracted_sky, extracted_sky_noise,
                                extracted_mf, extracted_mf_noise),
              "Could not perform flat-fielding");
    }

    /* Save flat-fielded spectrum + noise */
    if (debug_mode) 
    {
       check( uves_save_image_local("Flat-fielded spectrum", "fxb",
                                    extracted, chip, -1, filename_window,
                                    extracted_header, true),
              "Error saving image");
            
       check( uves_save_image_local("Flat-fielded spectrum noise", 
                                    "errfxb", extracted_noise, chip,
                                    -1, filename_window, extracted_header, true),
              "Error saving noise of flat-fielded image");
    }

    if (debug_mode && extracted_sky != NULL) 
    {
       check( uves_save_image_local("Flat-fielded sky", "fxsky", 
                                    extracted_sky, chip, -1,
                                    filename_window, extracted_header, true),
              "Error saving image");
            
       check( uves_save_image_local("Flat-fielded sky noise", "errfxsky",
                                    extracted_sky_noise, chip, -1,
                                    filename_window, extracted_header, true),
              "Error saving noise of flat-fielded image");
    }
    }

    if (fx != NULL)
    {
       *fx = cpl_image_duplicate(extracted);
    }

    /* Variance of flat-fielded, pre-rebinned spectrum
       is a product of science recipe (for whatever reason...) */
    if (flatfielded_variance != NULL)
    {

       check( *flatfielded_variance = 
              cpl_image_multiply_create(extracted_noise, 
                                        extracted_noise),
              "Error creating variance of flatfielded spectrum");

       passure(flatfielded_variance_header != NULL, " ");
       check( *flatfielded_variance_header = 
              uves_propertylist_duplicate(extracted_header),
              "Could not copy extracted spectrum header");
    }
    
    if (blaze_correct)
    {
       if (ex_method == EXTRACT_2D)
       {
          /* It requires an extracted spectrum which we don't have in 2d mode */
          uves_msg_low("Skipping blaze function correction for 2d extraction mode");
       }
       else
       {
          uves_msg("Calculating blaze function correction");
            
          check( blaze_ratio = uves_get_blaze_ratio(extracted, extracted_noise),
                 "Error calculating blaze function correction");
            
          uves_msg("Applying blaze function correction");
            
          check(( cpl_image_divide(extracted      , blaze_ratio),
                  cpl_image_divide(extracted_noise, blaze_ratio)),
                "Error applying blaze function correction");
            
          if (extracted_sky != NULL)   /* If sky was extracted (optimal) */
          {
             check(( cpl_image_multiply(extracted_sky, blaze_ratio),
                     cpl_image_multiply(extracted_sky_noise, blaze_ratio)),
                   "Error applying blaze function correction");
          }
       }
    }

    /* Rebin from (x, order) to (wavelength, order) */
    uves_msg("Rebinning spectrum");
    if (ex_method == EXTRACT_2D) {
       if (tilt_corr) {
          uves_msg_warning("Line tilt correction in rebinning "
                           "of 2d spectra is unsupported");
       }
       dispersion_int = uves_polynomial_duplicate(dispersion_relation[window-1]);
    }
    else if (tilt_corr) { 
       double objoffset;
        
       if (info_tbl != NULL) {
          objoffset = cpl_table_get_column_median(*info_tbl, "ObjPosOnSlit");
          /* This is the object position measured from the bottom of
             the of specified extraction window.
             Need to convert this to the same coordinates as used in the wavecal.
          */
            
          uves_msg_debug("Object position (from bottom of extraction window) = %.2f pixels",
                         objoffset);
            
          objoffset -= slit_length / 2;
          /* Now wrt middle of specified window */
            
          objoffset += slit_offset;
          /* Now wrt order trace center */
            
          uves_msg_debug("Object position (from slit center) = %.2f pixels",
                         objoffset);
       }
       else {
          /* Sky windows */
          uves_msg_debug("Object offset not measured during extraction, "
                         "using %.2f pixels", slit_offset);
          objoffset = slit_offset;
       }
        
       check( dispersion_int = interpolate_wave(dispersion_relation,
                                                linetable_header,
                                                objoffset),
              "Could not interpolate dispersion solutions");
        
       if (debug_mode) {
          check( poly_table = uves_polynomial_convert_to_table(dispersion_int), 
                 "Error converting polynomial to table");
            
          check( uves_save_table_local("Interpolated dispersion relation",
                                       "intdisp", 
                                       poly_table, chip, -1,
                                       filename_window, backsubbed_header, NULL),
                 "Error saving interpolated dispersion solution");
       }
    }
    else {
       dispersion_int = uves_polynomial_duplicate(dispersion_relation[window-1]);
    }
    
    uves_free_propertylist(rebinned_header);
    check( *rebinned_spectrum = uves_rebin(extracted,
                                           parameters, context, 
                                           linetable, dispersion_int,
                                           first_abs_order,
                                           last_abs_order,
                                           n_traces,
                                           false,
                                           false,
                                           rebinned_header),
           "Could not rebin spectrum");
    
    uves_msg("Rebinning spectrum noise");

    /* As in UVES/MIDAS the noise spectrum is rebinned to the same
     * level. It is not propagated using error propagation formula.
     * In other words, after this step, the noise level no longer
     *  describes the empirical noise actually observed in the spectrum
     * (which does change during rebinning depending on the bin width)
     */

    {
       bool threshold_to_positive = true;

       uves_free_propertylist(rebinned_header);
       check( *rebinned_noise = uves_rebin(extracted_noise,
                                           parameters, context, 
                                           linetable, dispersion_int,
                                           first_abs_order,
                                           last_abs_order,
                                           n_traces,
                                           threshold_to_positive,
                                           true,
                                           rebinned_header),
              "Could not rebin spectrum noise");
    }

    if (extracted_sky != NULL) {
       uves_msg("Rebinning sky spectrum");

       if (tilt_corr) {
          /* Optimal extraction extracts an average of the sky
             in the entire extraction window.
               
             Calibrate the sky spectrum using the dispersion solution
             at the extraction window center, i.e. at offset = slit_offset
          */
            
          check( dispersion_int_sky = interpolate_wave(dispersion_relation,
                                                       linetable_header,
                                                       slit_offset),
                 "Could not interpolate dispersion solutions");
       }
       else {
          /* Use middle solution */
          dispersion_int_sky = uves_polynomial_duplicate(dispersion_relation[1]);
       }

       /* Re-use the same rebinned_header */
       uves_free_propertylist(rebinned_header);
       check( rebinned_sky = uves_rebin(extracted_sky,
                                        parameters, context, 
                                        linetable, dispersion_int_sky,
                                        first_abs_order,
                                        last_abs_order,
                                        n_traces,
                                        false,
                                        false,
                                        rebinned_header),
              "Could not rebin sky noise");
        
       uves_msg("Rebinning sky spectrum noise");
        
       uves_free_propertylist(rebinned_header);
       check( rebinned_sky_noise = uves_rebin(extracted_sky_noise,
                                              parameters, context, 
                                              linetable, dispersion_int_sky,
                                              first_abs_order,
                                              last_abs_order,
                                              n_traces,
                                              true,
                                              true,
                                              rebinned_header),
              "Could not rebin sky noise");
    }
    
    /* Save rebinned spectrum + noise */
    if (debug_mode) 
    {
       const char *filename         = "";
       const char *filename_err     = "";
       const char *filename_sky     = "";
       const char *filename_sky_err = "";
       if (ff_method == FF_PIXEL)
       {
          filename         = "wxfb";
          filename_err     = "errwxfb";
          filename_sky     = "wxfsky";
          filename_sky_err = "errwxfsky";
       }
       else if (ff_method == FF_EXTRACT)
       {
          filename         =    "wfxb";
          filename_err     = "errwfxb";
          filename_sky     =    "wfxsky";
          filename_sky_err = "errwfxsky";
       }
       else if (ff_method == FF_NO)
       {
          filename         =    "wxb";
          filename_err     = "errwxb";
          filename_sky     =    "wxsky";
          filename_sky_err = "errwxsky";
       }
       else
       {
          passure( false, "Unknown ff_method: %d", ff_method);
       }

       check(uves_propertylist_copy_property_regexp(*rebinned_header,
                                                    backsubbed_header,
                                                    "^ESO ", 1),"error copying hierarch keys");

        
       check( uves_save_image_local("Rebinned spectrum",
                                    filename, *rebinned_spectrum,
                                    chip, -1, filename_window, *rebinned_header, true),
              "Error saving image");
        
       check( uves_save_image_local("Noise of rebinned spectrum", filename_err,
                                    *rebinned_noise, chip, -1, filename_window, 
                                    *rebinned_header, true),
              "Error saving image");

       if (extracted_sky != NULL)
       {
          check( uves_save_image_local("Rebinned sky", filename_sky,
                                       rebinned_sky, chip, -1,
                                       filename_window, *rebinned_header, true),
                 "Error saving image");
            
          check( uves_save_image_local("Noise of rebinned sky",
                                       filename_sky_err, 
                                       rebinned_sky_noise, chip, -1,
                                       filename_window, *rebinned_header, true),
                 "Error saving image");
       }
    }

    /* We also need to produce the  rebinned-immediately-after-extraction
       (but non flat-fielded) spectrum,
       which is a product of the science recipe.
       This is trivial unless ff_method is FF_EXTRACT 
    */
    if (resampled_spectrum != NULL)   /* Not for sky windows */
    {
       if (ff_method == FF_EXTRACT)
       {
          /* Rebin the extracted spectrum (before flatfielding) */
          uves_msg("Rebinning pre-flatfielded spectrum");    
            
          uves_free_propertylist(rebinned_header);
          check( *resampled_spectrum = 
                 uves_rebin(extracted_noff,
                            parameters, context, 
                            linetable, dispersion_int,
                            first_abs_order,
                            last_abs_order,
                            n_traces,
                            false,
                            false,
                            rebinned_header),
                 "Could not rebin spectrum");
            
          if (debug_mode) {

             check( uves_save_image_local("Rebinned, extracted spectrum", 
                                          "wx", *resampled_spectrum, 
                                          chip, -1, filename_window,
                                          *rebinned_header, true),
                    "Error saving image");
          }

       }
       else
       {
          check( *resampled_spectrum = cpl_image_duplicate(*rebinned_spectrum),
                 "Error copying rebinned spectrum");
       }
    }

    /* Merge orders to 1D spectrum */
    if (extracted_sky != NULL)
    {
       uves_msg("Merging sky");
       check( *merged_sky = uves_merge_orders(rebinned_sky, 
                                              rebinned_sky_noise,
                                              *rebinned_header,
                                              m_method,
                                              n_traces,
                                              merged_header,
                                              delt1,delt2,chip,
                                              &merged_sky_noise),
              "Error merging sky");
    }
    
    uves_msg("Merging spectrum");
    uves_free_propertylist(merged_header);
    check( *merged_spectrum = uves_merge_orders(*rebinned_spectrum, 
                                                *rebinned_noise,
						*rebinned_header,
						m_method,
						n_traces,
						merged_header,
						delt1,delt2,chip,
						merged_noise),
           "Error merging orders");

    check(uves_propertylist_copy_property_regexp(*merged_header,
						 backsubbed_header,
						 "^ESO ", 0),"error copying hierarch keys");

    if (debug_mode)
    {
       check( uves_save_image_local("Merged spectrum", "m", *merged_spectrum, 
                                    chip, -1, filename_window, *merged_header, true),
              "Error saving image");
        
       check( uves_save_image_local("Noise of merged spectrum", "errm",
                                    *merged_noise, chip, -1,
                                    filename_window, *merged_header, true),
              "Error saving image");
    }
    
    if (debug_mode && extracted_sky != NULL)
    {
       check( uves_save_image_local("Merged sky", "msky", *merged_sky, 
                                    chip, -1,
                                    filename_window, *merged_header, true),
              "Error saving image");
        
       check( uves_save_image_local("Noise of merged sky", "errmsky", 
                                    merged_sky_noise, chip, -1,
                                    filename_window, *merged_header, true),
              "Error saving image");
    }

  cleanup:
    uves_free_image(&extracted);
    uves_free_image(&extracted_noff);
    uves_free_image(&extracted_noise);
    uves_free_image(&extracted_sky);
    uves_free_image(&extracted_sky_noise);
    uves_free_image(&cosmic_image);
    uves_free_image(&blaze_ratio);
    uves_free_image(&weights);
    uves_polynomial_delete(&dispersion_int);
    uves_polynomial_delete(&dispersion_int_sky);
    uves_free_table(&poly_table);
    uves_free_propertylist(&extracted_header);
    uves_free_table(&profile_table);
    uves_free_image(&extracted_mf);
    uves_free_image(&extracted_mf_noise);
    uves_free_image(&rebinned_sky);
    uves_free_image(&rebinned_sky_noise);
    uves_free_image(&merged_sky_noise);
    
    uves_free_image(&temp_image);
    
    return cpl_error_get_code();
}

                      
/*----------------------------------------------------------------------------*/
/**
  @brief    Subtract sky from extracted spectrum

  @param    rebinned_obj                 The rebinned object spectrum
  @param    rebinned_obj_noise           Noise of rebinned object spectrum
  @param    rebinned_obj_header          Header of rebinned object spectrum
  @param    rebinned_sky1                The rebinned lower sky spectrum
  @param    rebinned_sky1_noise          Noise of rebinned lower sky spectrum
  @param    rebinned_sky1_header         Header of rebinned lower sky spectrum
  @param    rebinned_sky2                The rebinned upper sky spectrum
  @param    rebinned_sky2_noise          Noise of rebinned upper sky spectrum
  @param    rebinned_sky2_header         Header of rebinned upper sky spectrum

  @param    merged_obj                 The merged object spectrum
  @param    merged_obj_noise           Noise of merged object spectrum
  @param    merged_obj_header          Header of merged object spectrum
  @param    merged_sky1                The merged lower sky spectrum
  @param    merged_sky1_noise          Noise of merged lower sky spectrum
  @param    merged_sky1_header         Header of merged lower sky spectrum
  @param    merged_sky2                The merged upper sky spectrum
  @param    merged_sky2_noise          Noise of merged upper sky spectrum
  @param    merged_sky2_header         Header of merged upper sky spectrum

  @param    obj_slit                   Extraction slit length (object)
  @param    sky1_slit                  Extraction slit length (lower sky)
  @param    sky2_slit                  Extraction slit length (upper sky)

  @return   The merged sky spectrum which was subtracted, or NULL if the provided
  sky windows had no good pixels

  Due to (slightly) different dispersion relations, the ranges of sky and object spectra
  might not be exactly identical. This function takes care of subtracting the sky from
  rebinned and merged spectra while aligning the proper wavelengths bin with each other.

**/
/*----------------------------------------------------------------------------*/

static cpl_image *
subtract_sky(cpl_image *rebinned_obj, cpl_image *rebinned_obj_noise,
         uves_propertylist *rebinned_obj_header,
         const cpl_image *rebinned_sky1, const cpl_image *rebinned_sky1_noise,
         const uves_propertylist *rebinned_sky1_header,
         const cpl_image *rebinned_sky2, const cpl_image *rebinned_sky2_noise,
         const uves_propertylist *rebinned_sky2_header,
         cpl_image **merged_obj, cpl_image **merged_obj_noise, 
         uves_propertylist *merged_obj_header,
         const cpl_image *merged_sky1, const cpl_image *merged_sky1_noise, 
         const uves_propertylist *merged_sky1_header,
         const cpl_image *merged_sky2, const cpl_image *merged_sky2_noise, 
         const uves_propertylist *merged_sky2_header,
         double obj_slit, double sky1_slit, double sky2_slit)
{
    double wavestep;
    int norders;

    cpl_image *merged_sky = NULL;               /* Result */
    
    passure( rebinned_obj != NULL, " ");
    passure( rebinned_obj_noise != NULL, " ");
    passure( rebinned_obj_header != NULL, " ");
    passure( merged_obj != NULL, " ");
    passure( merged_obj_noise != NULL, " ");
    passure( merged_obj_header != NULL, " ");
    passure( *merged_obj != NULL, " ");
    passure( *merged_obj_noise != NULL, " ");
    /* Sky spectra may be NULL (if not extracted) */

    check( wavestep = uves_pfits_get_cdelt1(rebinned_obj_header), 
       "Error reading wavelength step");
    norders = cpl_image_get_size_y(rebinned_obj);

    /* Do some consistency checking
       (that 'wavestep' and 'norders' is same for all spectra) */
    assure((rebinned_sky1 == NULL || norders == cpl_image_get_size_y(rebinned_sky1)) &&
       (rebinned_sky2 == NULL || norders == cpl_image_get_size_y(rebinned_sky2)), 
       CPL_ERROR_ILLEGAL_INPUT,
       "Different number of orders in object/sky spectra: obj = %d, "
       "sky1 = %" CPL_SIZE_FORMAT ", sky3 = %" CPL_SIZE_FORMAT "",
       norders,
       cpl_image_get_size_y(rebinned_sky1),
       cpl_image_get_size_y(rebinned_sky2));

    if (rebinned_sky1 != NULL)
    {
        double wavestep1;
        check( wavestep1 = uves_pfits_get_cdelt1(rebinned_sky1_header), 
           "Error reading wavelength step");
        assure( fabs(wavestep1 - wavestep) / wavestep < 0.01, 
            CPL_ERROR_ILLEGAL_INPUT,
            "Different bin widths: sky1 = %f ; obj = %f", 
            wavestep1, wavestep);
    }
    if (rebinned_sky2 != NULL)
    {
        double wavestep2;
        check( wavestep2 = uves_pfits_get_cdelt1(rebinned_sky2_header),
           "Error reading wavelength step");
        assure( fabs(wavestep2 - wavestep) / wavestep < 0.01,
            CPL_ERROR_ILLEGAL_INPUT,
            "Different bin widths: sky3 = %f ; obj = %f",
            wavestep2, wavestep);
    }

    /* Subtract sky (rebinned spectrum) */
    {
    int order;
    for (order = 1; order <= norders; order++)
        {
        double obj_start   , obj_end;
        double sky1_start  , sky1_end;
        double sky2_start  , sky2_end;
        double common_start, common_end;
        
        check( obj_start = uves_pfits_get_wstart(rebinned_obj_header, order),
               "Error reading start wavelength for order #%d", order);
        check( obj_end   = uves_pfits_get_wend  (rebinned_obj_header, order),
               "Error reading end wavelength for order #%d", order);
        
        if (rebinned_sky1 != NULL)
            {
            check( sky1_start = 
                   uves_pfits_get_wstart(rebinned_sky1_header, order), 
                   "Error reading start wavelength for order #%d", order);
            check( sky1_end   = 
                   uves_pfits_get_wend  (rebinned_sky1_header, order), 
                   "Error reading end wavelength for order #%d", order);
            }
        else
            {
            sky1_start = obj_start;
            sky1_end   = obj_end;
            }
        
        if (rebinned_sky2 != NULL)
            {
            check( sky2_start = 
                   uves_pfits_get_wstart(rebinned_sky2_header, order), 
                   "Error reading start wavelength for order #%d", order);
            check( sky2_end   = 
                   uves_pfits_get_wend  (rebinned_sky2_header, order),
                   "Error reading end wavelength for order #%d", order);
            }
        else
            {
            sky2_start = obj_start;
            sky2_end   = obj_end;
            }
        
        check( merged_sky =
               subtract_sky_row(rebinned_obj , rebinned_obj_noise , 
                    obj_start , obj_end,  obj_slit,
                    rebinned_sky1, rebinned_sky1_noise,
                    sky1_start, sky1_end, sky1_slit,
                    rebinned_sky2, rebinned_sky2_noise, 
                    sky2_start, sky2_end, sky2_slit,
                    order, wavestep, &common_start, 
                    &common_end),
               "Could not subtract sky for rebinned spectrum order #%d", order);
        uves_free_image(&merged_sky);
        
        check( uves_pfits_set_wstart(rebinned_obj_header, order, common_start),
               "Error updating start wavelength for order #%d", order);
        check( uves_pfits_set_wend  (rebinned_obj_header, order, common_end  ), 
               "Error updating start wavelength for order #%d", order);
        }
    }

    /* Subtract sky (merged spectrum) */
    {
    double obj_start   , obj_end;
    double sky1_start  , sky1_end;
    double sky2_start  , sky2_end;
    double common_start, common_end;
    
    obj_start = uves_pfits_get_crval1(merged_obj_header);
    obj_end   = obj_start + wavestep * (cpl_image_get_size_x(*merged_obj) - 1);
    
    if (merged_sky1 != NULL)
        {
        sky1_start = uves_pfits_get_crval1(merged_sky1_header);
        sky1_end   = sky1_start +
            wavestep * (cpl_image_get_size_x(merged_sky1) - 1);
        }
    else
        {
        sky1_start = obj_start;
        sky1_end   = obj_end;
        }
    
    if (merged_sky2 != NULL)
        {
        sky2_start = uves_pfits_get_crval1(merged_sky2_header);
        sky2_end   = sky2_start + 
            wavestep * (cpl_image_get_size_x(merged_sky2) - 1);
        }
    else
        {
        sky2_start = obj_start;
        sky2_end   = obj_end;
        }
    
    /* Subtract sky for image row 1 (the only row in the image) */
    check(     merged_sky = subtract_sky_row(*merged_obj, *merged_obj_noise, 
                          obj_start , obj_end,  obj_slit,
                          merged_sky1, merged_sky1_noise,
                          sky1_start, sky1_end, sky1_slit,
                          merged_sky2, merged_sky2_noise,
                          sky2_start, sky2_end, sky2_slit,
                          1, wavestep, &common_start,
                          &common_end),
        "Error subtracting sky of merged spectrum");
    
    check( uves_pfits_set_crval1(merged_obj_header, common_start),
           "Could not update start wavelength");
    
    /* Make sure that the last bin corresponds to 'common_end' wavelength */
    check( uves_crop_image(merged_obj,
                   1, 1,
                   1 + uves_round_double((common_end - 
                              common_start)/wavestep), 
                   1),
           "Error cropping merged spectrum");
    
    check( uves_crop_image(merged_obj_noise,
                   1, 1,
                   1 + uves_round_double((common_end - 
                              common_start)/wavestep), 
                   1),
           "Error cropping merged spectrum noise");
    
    if (merged_sky != NULL)
        {
        /* The image header also applies for the sky */
        assure( cpl_image_get_size_x(merged_sky) == 
            cpl_image_get_size_x(*merged_obj), CPL_ERROR_ILLEGAL_OUTPUT,
            "Sky and object spectrum sizes differ, "
            "sky = %" CPL_SIZE_FORMAT " bins, obj = %" CPL_SIZE_FORMAT " bins", 
            cpl_image_get_size_x(merged_sky),
            cpl_image_get_size_x(*merged_obj));
        }
    }
    
  cleanup:
    if (cpl_error_get_code() != CPL_ERROR_NONE)
    {
        uves_free_image(&merged_sky);
    }
    return merged_sky;
}    

/*----------------------------------------------------------------------------*/
/**
  @brief    Subtract sky from 1 row of spectrum
  @param    obj                 The object spectrum
  @param    obj_noise           Noise of object spectrum
  @param    obj_start           Wavelength of first bin of row
  @param    obj_end             Wavelength of last bin to consider
  @param    obj_slit            Object slit length

  @param    sky1                The lower sky spectrum (or NULL if not extracted)
  @param    sky1_noise          Noise of lower sky spectrum
  @param    sky1_start          Wavelength of first bin of row
  @param    sky1_end            Wavelength of last bin to consider
  @param    sky1_slit           Lower sky slit length

  @param    sky2                The upper sky spectrum (or NULL if not extracted)
  @param    sky2_noise          Noise of upper sky spectrum
  @param    sky2_start          Wavelength of first bin of row
  @param    sky2_end            Wavelength of last bin to consider
  @param    sky2_slit           Upper sky slit length
  @param    row                 Row (order number) of spectrum images to process
  @param    wavestep            Width of one win (same for obj/sky spectra)
  @param    common_start        (output) Minimum common wavelength
  @param    common_end          (output) Maximum common wavelength

  @return   The merged sky spectrum, or NULL if provided sky images had no
            good pixels

  This function subtracts sky from (only) the specified @em row of the input
  spectrum. Before subtraction, the fluxes are properly scaled according 
  to extraction slit lengths.
**/
/*----------------------------------------------------------------------------*/

static cpl_image *
subtract_sky_row(cpl_image *obj, cpl_image *obj_noise, 
         double obj_start, double obj_end, double obj_slit,
         const cpl_image *sky1, const cpl_image *sky1_noise, 
         double sky1_start, double sky1_end, double sky1_slit,
         const cpl_image *sky2, const cpl_image *sky2_noise, 
         double sky2_start, double sky2_end, double sky2_slit,
         int row, double wavestep, 
         double *common_start, double *common_end)
{
    int first_bin_obj;
    int first_bin_sky1;
    int first_bin_sky2;
    int nbins;
    cpl_image *common_obj = NULL;     /* Extract the common wavelength range ... */
    cpl_image *common_sky1 = NULL;    /* ... to these 1D images                  */
    cpl_image *common_sky2 = NULL;
    cpl_image *common_obj_noise = NULL;
    cpl_image *common_sky1_noise = NULL;
    cpl_image *common_sky2_noise = NULL;
    bool is_good1, is_good2;          /* Do the two sky images contain valid pixels? */

    cpl_image *common_sky         = NULL;   /* The combined sky spectrum normalized
                           to object slit length (returned) */
    cpl_image *common_sky_noise   = NULL;

    cpl_image *temp               = NULL;

    *common_start = uves_max_double(obj_start, uves_max_double(sky1_start, sky2_start));
    *common_end   = uves_min_double(obj_end  , uves_min_double(sky1_end  , sky2_end  ));

    if (*common_start <= *common_end)
    {
        nbins = 1 + uves_round_double((*common_end - *common_start) / wavestep);
        
        uves_msg_debug("Lower sky range: %f - %f w.l.u.", sky1_start, sky1_end);
        uves_msg_debug("Upper sky range: %f - %f w.l.u.", sky2_start, sky2_end);
        uves_msg_debug("Object sky range: %f - %f w.l.u.", obj_start, obj_end);
        uves_msg_debug("Sky/object common wavelength range in order %d: "
               "%f - %f w.l.u. (%d bins)", 
               row, *common_start, *common_end, nbins);
        
        first_bin_obj  = 1 + uves_round_double((*common_start - obj_start )/wavestep);
        first_bin_sky1 = 1 + uves_round_double((*common_start - sky1_start)/wavestep);
        first_bin_sky2 = 1 + uves_round_double((*common_start - sky2_start)/wavestep);
        
        /* Extract common bins, normalize sky windows to object slit length */
        check( common_obj       = cpl_image_extract(obj, 
                            first_bin_obj, row,
                            first_bin_obj + nbins-1, row), 
           "Error extracting common rows (object)");
        
        check( common_obj_noise = cpl_image_extract(obj_noise, 
                            first_bin_obj, row,
                            first_bin_obj + nbins-1, row), 
           "Error extracting common rows (object noise)");
        
        if (sky1 != NULL)
        {
            check( common_sky1 =
               cpl_image_extract(sky1, 
                         first_bin_sky1, row,
                         first_bin_sky1 + nbins-1, row),
               "Error extracting common rows (lower sky)");
            
            check( common_sky1_noise =
               cpl_image_extract(sky1_noise, 
                         first_bin_sky1, row,
                         first_bin_sky1 + nbins-1, row), 
               "Error extracting common rows (lower sky noise)");
            
            check(( cpl_image_multiply_scalar(common_sky1      , obj_slit / sky1_slit),
                cpl_image_multiply_scalar(common_sky1_noise, obj_slit / sky1_slit)),
               "Error normalizing sky flux");
            
            is_good1 = 
            cpl_image_count_rejected(common_sky1) <
            cpl_image_get_size_x(common_sky1)*
            cpl_image_get_size_y(common_sky1) &&  
            /* Note order of evaluation. cpl_image_get_min() would fail if
               there were no good pixels */
            cpl_image_get_min(common_sky1_noise) > 0;
        }
        else
        {
            is_good1 = false;
        }
        if (sky2 != NULL)
        {
            check( common_sky2       = cpl_image_extract(sky2, 
                                 first_bin_sky2, row,
                                 first_bin_sky2 + nbins-1, row), 
               "Error extracting common rows (upper sky)");
            
            check( common_sky2_noise = cpl_image_extract(sky2_noise, 
                                 first_bin_sky2, row,
                                 first_bin_sky2 + nbins-1, row), 
               "Error extracting common rows (upper sky noise)");
            
            check(( cpl_image_multiply_scalar(common_sky2      , obj_slit / sky2_slit),
                cpl_image_multiply_scalar(common_sky2_noise, obj_slit / sky2_slit)),
               "Error normalizing sky flux");

            is_good2 = 
            cpl_image_count_rejected(common_sky2) <
            cpl_image_get_size_x(common_sky2)*
            cpl_image_get_size_y(common_sky2) &&  
            cpl_image_get_min(common_sky2_noise) > 0;
        }
        else
        {
            is_good2 = false;
        }
               
        
        /* Optimally average the two sky windows 
           (one of which might not have been extracted) */
        if (is_good1 && is_good2)
        {
                    check( common_sky = 
                           uves_average_images(common_sky1, common_sky1_noise,
                                               common_sky2, common_sky2_noise,
                                               &common_sky_noise),
                           "Error combining sky windows");
                }
        else if (is_good1 && !is_good2)
        {
            common_sky       = cpl_image_duplicate(common_sky1);
            common_sky_noise = cpl_image_duplicate(common_sky1_noise);
        }
        else if (!is_good1 && is_good2)
        {
            common_sky       = cpl_image_duplicate(common_sky2);
            common_sky_noise = cpl_image_duplicate(common_sky2_noise);
        }
        else
        {
            common_sky = NULL;
        }
        
        if (common_sky != NULL)
        {   
            /* Do the subtraction, threshold to [0 ; oo [  */
	  /* Commented out as we should not lower threshold to 0
            check(( cpl_image_subtract (common_obj, common_sky),
                cpl_image_threshold(common_obj, 
                        0, DBL_MAX,
                        0, DBL_MAX)),
              "Error subtracting combined sky");
	  */
            check(( cpl_image_subtract (common_obj, common_sky)),
              "Error subtracting combined sky");
            
            /*  Propagate noise:  
            obj_noise := sqrt( obj_noise^2 + sky_noise^2 )  */
            check(( cpl_image_power(common_obj_noise, 2),
                cpl_image_power(common_sky_noise, 2),
                cpl_image_add  (common_obj_noise, common_sky_noise),
                cpl_image_power(common_obj_noise, 0.5)),
              "Error propagating noise during sky subtration");
            
            /* Copy results to relevant row of input spectrum */
            check(( cpl_image_copy(obj,
                       common_obj,
                       1, row),
                cpl_image_copy(obj_noise,
                       common_obj_noise,
                       1, row)),
              "Error writing subtracted flux to row %d of spectrum", row);
        }

    } /* Object and both sky windows do have an overlap in this order */
    else
    {
        int x;
        
        uves_msg_low("Extracted object and sky spectra have no overlap in order #%d. "
             "Order marked as bad", row);
        
        for (x = 1; x <= cpl_image_get_size_x(obj); x++)
        {
            check(( cpl_image_reject(obj      , x, row),
                cpl_image_reject(obj_noise, x, row)),
              "Error rejecting sky-subtracted spectrum "
              "at (x, row) = (%d, %d)", x, row);
        }
    }
    
  cleanup:
    uves_free_image(&common_obj);
    uves_free_image(&common_sky1);
    uves_free_image(&common_sky2);
    uves_free_image(&common_obj_noise);
    uves_free_image(&common_sky_noise);
    uves_free_image(&common_sky1_noise);
    uves_free_image(&common_sky2_noise);
    uves_free_image(&temp);
    if (cpl_error_get_code() != CPL_ERROR_NONE)
    {
        uves_free_image(&common_sky);
    }
    
    return common_sky;
}

/*----------------------------------------------------------------------------*/
/**
  @brief    Measure object offset w.r.t. slit center
  @param    back_subbed          The input image
  @param    ordertable           The order table. Used to get min/max order numbers
  @param    order_locations      Polynomial describing the order locations
  @param    search_range         The half search window height (in pixels)
  @param    nsamples             Number of sample points per order
  @param    doffset              (out) If non-null, the measured width (1 sigma)
                                 of the object is returned.
  @return   The average object offset in relation to the center of the order

  This function finds the flux centroid (in y-direction) of @em nsamples 
  median collapsed chunks in each order and returns the mean centroid
  position (average of all orders).

**/
/*----------------------------------------------------------------------------*/
static double get_offset(const cpl_image *back_subbed, 
             const cpl_table *ordertable, 
             const polynomial *order_locations,
             double search_range, int nsamples, double *doffset)
{
    cpl_image *chunk = NULL;          /* Chunks                          */
    cpl_image *chunk_col = NULL;      /* Chunks median collapsed along x */

    int minorder, maxorder;
    int order, x, nx, ny;
    double sum = 0, sum_o = 0, sum_oo = 0;         /* Zero'th, first and 
                              second moment of offsets */
    int s_r_int = uves_round_double(search_range); /* Search range as an integer */

    passure( back_subbed != NULL, " ");
    passure( ordertable != NULL, " ");
    passure( order_locations != NULL, " ");
    /* doffset may be NULL */
    assure( nsamples >= 1, CPL_ERROR_ILLEGAL_INPUT, 
        "Illegal number of sample points per order: %d", nsamples);

    minorder = cpl_table_get_column_min(ordertable, "Order");
    maxorder = cpl_table_get_column_max(ordertable, "Order");
    nx = cpl_image_get_size_x(back_subbed);
    ny = cpl_image_get_size_y(back_subbed);

    sum    = 0;
    sum_o  = 0;
    sum_oo = 0;
    for (order = minorder; order <= maxorder; order++)
    {
        int stepx = nx / nsamples;
        
        for (x = stepx/2; x <= nx; x += stepx)
        {
            int y = uves_round_double(
            uves_polynomial_evaluate_2d(order_locations, x, order));

            if (1 <= y - s_r_int && y + s_r_int <= ny)
            {
                double offset;

                /* Get centroid.pos. of median collapsed window */
                
                chunk = 
                cpl_image_extract(back_subbed,
                          uves_max_int(1 , x - stepx/2), 
                          y - s_r_int,
                          uves_min_int(nx, x + stepx/2),
                          y + s_r_int);
                
                chunk_col = 
                /* Result is single column image */
                cpl_image_collapse_median_create(chunk,
                                 1,
                                 0, 0); /* No filtering */
                
                
                /* Offset in world coordinates: row=1 in 'chunk_col'
                   corresponds to row=(y - s_r_int) in 'back_subbed' */
                offset = (y - s_r_int - 1) + 
                cpl_image_get_centroid_y_window(chunk_col,
                                1, 1,
                                1, 
                                cpl_image_get_size_y(chunk_col));

                /* Get offset relative to slit center */
                offset -= y;

                uves_free_image(&chunk);
                uves_free_image(&chunk_col);
                
                sum    += 1;
                sum_o  += offset;
                sum_oo += offset*offset;
            }
        }
    }
    
    /* This should never happen, but if it does 
       fail cleanly instead of dividing by zero */
    assure( sum > 0, CPL_ERROR_ILLEGAL_OUTPUT,
        "No evaluation points inside image!");

    if (doffset != NULL)
    {
        *doffset = sqrt(sum_oo/(1.0*sum) - 
                (sum_o*sum_o) / (sum*1.0*sum));
    }
    
  cleanup:
    uves_free_image(&chunk);
    uves_free_image(&chunk_col);

    return (1.0*sum_o) / sum;
}


/*----------------------------------------------------------------------------*/
/**
  @brief    Get (normalized) object to flat-field blaze function ratio
  @param    spectrum                The extracted, flat-fielded spectrum
                                    in (pixel, order)-space
  @param    spectrum_noise          Noise associated with @em spectrum

  @return   The correction factor

  Due to the different blaze function profile of flat/science spectra
  (for unknown reason) there will be order scale ripples when merging. This
  function 

  This function
  - Normalizes in each order the flux to 1:
    norm(x) = f(x) / median_x(f).
  - Fits for each x0 a low-order polynomial to norm(x0, order).
    The result is a correction factor close to 1 that corrects for
    the different shape of blaze functions (but doesn't change the
    overall normalization)

**/
/*----------------------------------------------------------------------------*/
static cpl_image *
uves_get_blaze_ratio(const cpl_image *spectrum,
             const cpl_image *spectrum_noise)
{
    int nx, ny;
    int smooth_x, smooth_y;

    cpl_image *blaze_ratio       = NULL;
    cpl_image *blaze_ratio_noise = NULL;

    cpl_table *values = NULL;
    polynomial *p = NULL;

    passure( spectrum       != NULL, " ");
    passure( spectrum_noise != NULL, " ");
    
    nx = cpl_image_get_size_x(spectrum);
    ny = cpl_image_get_size_y(spectrum);
    
    blaze_ratio       = cpl_image_duplicate(spectrum);
    blaze_ratio_noise = cpl_image_duplicate(spectrum_noise);
    assure_mem( blaze_ratio );
    assure_mem( blaze_ratio_noise );
    
    /* Normalize each row in ratio to median = 1,
       so that the overall normalization doesn't change 
    */


    {
    int x, y;
    
    for (y = 1; y <= ny; y++)
        {
        double median = cpl_image_get_median_window(blaze_ratio,
                                1, y,
                                nx, y);
        
        if (median == 0)
            {
            /* The cpl_image_get_median_window function is broken;
               it doesn't take bad pixels into account. That sometimes
               leads to a zero median */
            
            /* This mostly happens for the first and last orders */
            
            double max_noise = cpl_image_get_max(blaze_ratio_noise);

            for (x = 1; x <= nx; x++)
                {
                cpl_image_set(blaze_ratio      , x, y, 1);

                /* effectively exclude from fit: */
                cpl_image_set(blaze_ratio_noise, x, y, max_noise);
                }
            }
        else
            {
            /* Divide this row by median,
               Exclude pixels deviating more than a factor of, say, 5 */
            double exclude = 2;

            for (x = 1; x <= nx; x++)
                {
                int pis_rejected1, pis_rejected2;
                double val1, val2;
                
                val1 = cpl_image_get(blaze_ratio      , 
                             x, y, &pis_rejected1);
                val2 = cpl_image_get(blaze_ratio_noise,
                             x, y, &pis_rejected2);
                
                if (!pis_rejected1 && !pis_rejected2 && 
                    val1/median < exclude && val1/median > 1/exclude)
                    {
                    cpl_image_set(blaze_ratio      , 
                              x, y, val1 / median); 
                    cpl_image_set(blaze_ratio_noise, 
                              x, y, val2 / median); 
                    }
                else
                    {
                    /* Set to 1, then reject. This is to deal with
                       a plotter that might plot bad pixels */
                    cpl_image_set   (blaze_ratio      , x, y, 1);

                    cpl_image_reject(blaze_ratio      , x, y);
                    cpl_image_reject(blaze_ratio_noise, x, y);
                    }
                }
            }
        }

    uves_plot_image_rows(blaze_ratio, 1, ny, ny/10,
                 "x", "y", "ratio (normalized to 1)");
    }

    smooth_x = nx / 20 + 1;   /* >0 */
    smooth_y = ny / 20 + 1;   /* >0 */
    check( uves_filter_image_median(&blaze_ratio,
                    smooth_x, smooth_y,       /* x-radius, y-radius */
                    false),         /* extrapolate border pixels? */
       "Error creating smoothed ratio");

    uves_plot_image_rows(blaze_ratio, 1, ny, ny/10, "x", "y", "ratio (smoothed)");
    

    /* For each x, fit polynomial as function of y.
     * Use kappa-sigma clipping to eliminate single order
     * spectral featues. This should leave only the
     * systematics (i.e. the ratio of obj/flat blaze profiles)
     */
    {
    int x, y;
    double interpolated=0;
    for (x = 1; x <= nx; x++)
        {
        int current_row;

        /* Table rows are removed when kappa-sigma clipping,
           so we have to create a new table for each column */

        uves_free_table(&values);
        values = cpl_table_new(ny);
        cpl_table_new_column(values, "Y", CPL_TYPE_INT);
        cpl_table_new_column(values, "Ratio", CPL_TYPE_DOUBLE);
        cpl_table_new_column(values, "dRatio", CPL_TYPE_DOUBLE);
        
        assure_mem( values );

        current_row = 0;
        for (y = 1; y <= ny; y++)
            {
            double dratio=0.;
            int pis_rejected1=0, pis_rejected2=0;

            double ratio  = cpl_image_get(blaze_ratio      ,
                           x, y, &pis_rejected1);
            check_nomsg(dratio = cpl_image_get(blaze_ratio_noise, 
                                               x, y, &pis_rejected2));
/*
            uves_msg("x=%d, y=%d, rej1=%d rej2=%d",
                     x,y,pis_rejected1,pis_rejected2);
*/
            if (!pis_rejected1 && !pis_rejected2)
                {
                cpl_table_set_int   (values, "Y"     , 
                             current_row, y);
                cpl_table_set_double(values, "Ratio" , 
                             current_row, ratio);
                cpl_table_set_double(values, "dRatio", 
                             current_row, dratio);
                current_row += 1;
                }
            else
                {
                /* Ignore current order */
                }
            }
         
        cpl_table_set_size(values, current_row);
        if(current_row>UVES_MIN_LINE_ROWS_TO_MAKE_FIT)
        {
            int degree = 2;
            double kappa = 2;
            
            uves_polynomial_delete(&p);
            p = uves_polynomial_regression_1d(values,
                              "Y", "Ratio", "dRatio",
                              degree,
                              NULL, NULL, /* fit, residual^2 */
                              NULL,       /* mse             */
                              kappa);

            /* If fitting failed because there were too few points,
             * or matrix was singular */
            if (cpl_error_get_code() == CPL_ERROR_ILLEGAL_INPUT || 
            cpl_error_get_code() == CPL_ERROR_ILLEGAL_OUTPUT)
            {
                uves_error_reset();
                
                /* Then set p(x) = 1 
                 * by fitting a line through (1,1) - (2,1)
                 */


                /* The table is in a 'dirty' state (contains
                   temporary columns) if fitting routine failed
                   (that routine is not exception safe),
                   so don't try to reuse current table */
                
                uves_free_table(&values);
                values = cpl_table_new(2);
                cpl_table_new_column(values, "Y", CPL_TYPE_INT);
                cpl_table_new_column(values, "Ratio", CPL_TYPE_DOUBLE);
                cpl_table_set_int   (values, "Y"     , 0, 1);
                cpl_table_set_double(values, "Ratio" , 0, 1);
                cpl_table_set_int   (values, "Y"     , 1, 2);
                cpl_table_set_double(values, "Ratio" , 1, 1);
                
                degree = 2;
                kappa = -1;
                uves_polynomial_delete(&p);
                p = uves_polynomial_regression_1d(values,
                                  "Y", "Ratio", NULL,
                                  degree,       
                                  NULL, NULL, /* fit, residual^2 */
                                  NULL,       /* mse             */
                                  kappa);                  
            }
            assure( cpl_error_get_code() == CPL_ERROR_NONE, cpl_error_get_code(),
                "Could not fit %d. degree polynomial to column %d", degree, x);

        } else {
           interpolated=UVES_BLAZE_DUMMY_VAL;
        }

        for (y = 1; y <= ny; y++)
            {
               if(current_row>UVES_MIN_LINE_ROWS_TO_MAKE_FIT) {
                  interpolated = uves_polynomial_evaluate_1d(p, y);
               }
               cpl_image_set(blaze_ratio, x, y, fabs(interpolated));
            }

        }

    /* post smooth */
    check( uves_filter_image_median(&blaze_ratio, 
                    2*smooth_x, 2*smooth_y, /* x-radius, y-radius */
                    false),                 /* extrapolate at border? */
           "Error creating smoothed ratio");
    
    uves_plot_image_rows(blaze_ratio, 1, ny, ny/10, "x", "y", "ratio (poly. fit)");

    }

    
    /*
      printf("ratio\n");
      cpl_stats_dump(cpl_stats_new_from_image(ratio, CPL_STATS_ALL), CPL_STATS_ALL, stdout);
      printf("image\n");
      cpl_stats_dump(cpl_stats_new_from_image(image, CPL_STATS_ALL), CPL_STATS_ALL, stdout);
      printf("noise\n");
      cpl_stats_dump(cpl_stats_new_from_image(noise, CPL_STATS_ALL), CPL_STATS_ALL, stdout);
      
      passure(false, "");
      
    */

  cleanup:
    uves_free_table(&values);
    uves_polynomial_delete(&p);
    uves_free_image(&blaze_ratio_noise);

    if (cpl_error_get_code() != CPL_ERROR_NONE)
    {
        uves_free_image(&blaze_ratio);
    }
    return blaze_ratio;
}


/**@}*/