/* $Id: fors_science.c,v 1.42 2013-08-14 15:04:57 cgarcia Exp $
 *
 * This file is part of the FORS Data Reduction Pipeline
 * Copyright (C) 2002-2010 European Southern Observatory
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 */

/*
 * $Author: cgarcia $
 * $Date: 2013-08-14 15:04:57 $
 * $Revision: 1.42 $
 * $Name: not supported by cvs2svn $
 */

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

#include <cmath>
#include <stdexcept>
#include <string>
#include <sstream>
#include <iostream>
#include <algorithm>
#include <vector>
#include <utility>
#include <math.h>
#include <string.h>
#include <cpl.h>
#include <moses.h>
#include <fors_dfs.h>
#include <fors_tools.h>
#include <fors_qc.h>
#include <fors_header.h>
#include <fors_utils.h>
#include "fors_preprocess.h"
#include "fors_subtract_bias.h"
#include "fors_response.h"
#include "fors_detected_slits.h"
#include "global_distortion.h"

static int fors_science_create(cpl_plugin *);
static int fors_science_exec(cpl_plugin *);
static int fors_science_destroy(cpl_plugin *);
static int fors_science(cpl_parameterlist *, cpl_frameset *);

static int fors_science_response_fill_ignore
(const cpl_table * flux_table, const cpl_table * telluric_table,
 const std::string& grism_name,
 const std::string& resp_ignore_mode, const std::string& resp_ignore_lines,
 std::vector<double>& resp_ignore_lines_list,
 std::vector<std::pair<double, double> >& resp_ignore_ranges_list);
static cpl_image * fors_science_map_disp_profile
(const cpl_frameset * frameset, const char * flat_disp_profile_tag, 
 fors::detected_slits& det_slits, mosca::wavelength_calibration& wave_cal,
 double firstLambda, double lastLambda, double dispersion);
static void fors_science_correct_flat_sed
(cpl_image *spectra,  cpl_table *objects,
 cpl_image * mapped_flat_disp_profile,
 cpl_propertylist * flat_sed_header,
 cpl_propertylist * specphot_header);
static bool fors_science_response_apply_flat_corr
(cpl_frameset * frameset, const char * flat_disp_profile_tag,
 std::string& resp_use_flat_sed, cpl_table * grism_table, int standard);
static bool fors_science_photcal_apply_flat_corr
(cpl_frameset * frameset, const char * specphot_tag, 
 const char * master_specphot_tag, const char * flat_disp_profile_tag,
 int* fluxcal, bool response_apply_flat_corr, int standard);


static char fors_science_description[] =
"This recipe is used to reduce scientific spectra using the extraction\n"
"mask and the products created by the recipe fors_calib. The spectra are\n"
"bias subtracted, flat fielded (if a normalised flat field is specified)\n"
"and remapped eliminating the optical distortions. The wavelength calibration\n"
"can be optionally upgraded using a number of sky lines: if no sky lines\n"
"catalog of wavelengths is specified, an internal one is used instead.\n"
"If the alignment to the sky lines is performed, the input dispersion\n"
"coefficients table is upgraded and saved to disk, and a new CCD wavelengths\n"
"map is created.\n"
"This recipe accepts both FORS1 and FORS2 frames. A grism table (typically\n"
"depending on the instrument mode, and in particular on the grism used)\n"
"may also be specified: this table contains a default recipe parameter\n" 
"setting to control the way spectra are extracted for a specific instrument\n"
"mode, as it is used for automatic run of the pipeline on Paranal and in\n" 
"Garching. If this table is specified, it will modify the default recipe\n" 
"parameter setting, with the exception of those parameters which have been\n" 
"explicitly modifyed on the command line. If a grism table is not specified,\n"
"the input recipe parameters values will always be read from the command\n" 
"line, or from an esorex configuration file if present, or from their\n" 
"generic default values (that are rarely meaningful).\n" 
"In the table below the MXU acronym can be read alternatively as MOS\n"
"and LSS, depending on the instrument mode of the input data. The acronym\n"
"SCI on products should be read STD in case of standard stars observations\n"
"A CURV_COEFF table is not (yet) expected for LSS data.\n"
"Either a scientific or a standard star exposure can be specified in input.\n"
"Only in case of a standard star exposure input, the atmospheric extinction\n"
"table and a table with the physical fluxes of the observed standard star\n"
"must be specified in input, and a spectro-photometric table is created in\n"
"output. This table can then be input again to this recipe, always with an\n"
"atmospheric extinction table, and if a photometric calibration is requested\n"
"then flux calibrated spectra (in units of erg/cm/cm/s/Angstrom) are also\n" 
"written in output.\n\n"

"Input files:\n\n"
"  DO category:               Type:       Explanation:         Required:\n"
"  SCIENCE_MXU                Raw         Scientific exposure      Y\n"
"  or STANDARD_MXU            Raw         Standard star exposure   Y\n"
"  MASTER_BIAS                Calib       Master bias              Y\n"
"  GRISM_TABLE                Calib       Grism table              .\n"
"  MASTER_SKYLINECAT          Calib       Sky lines catalog        .\n"
"\n"
"  MASTER_NORM_FLAT_MXU       Calib       Normalised flat field    .\n"
"  DISP_COEFF_MXU             Calib       Inverse dispersion       Y\n"
"  CURV_COEFF_MXU             Calib       Spectral curvature       Y\n"
"  SLIT_LOCATION_MXU          Calib       Slits positions table    Y\n"
"  FLAT_DISP_PROFILE_MXU      Calib       Slits dispersion profile .\n"
"\n"
"  or, in case of LSS-like MOS/MXU data,\n"
"\n"
"  MASTER_NORM_FLAT_LONG_MXU  Calib       Normalised flat field    .\n"
"  DISP_COEFF_LONG_MXU        Calib       Inverse dispersion       Y\n"
"  SLIT_LOCATION_LONG_MXU     Calib       Slits positions table    Y\n"
"  GLOBAL_DISTORTION_TABLE    Calib       Global distortion        .\n"
"\n"
"  In case STANDARD_MXU is specified in input,\n"
"\n"
"  EXTINCT_TABLE              Calib       Atmospheric extinction   Y\n"
"  STD_FLUX_TABLE             Calib       Standard star flux       Y\n"
"  TELLURIC_CONTAMINATION     Calib       Telluric regions list    .\n"
"\n"
"  The following input files are mandatory if photometric calibrated"
"  spectra are desired:\n"
"\n"
"  EXTINCT_TABLE              Calib       Atmospheric extinction   Y\n"
"  SPECPHOT_TABLE             Calib       Response curves          Y\n"
"\n"
"  If requested for standard star data, the SPECPHOT_TABLE can be dropped:\n"
"  in this case the correction is applied using the SPECPHOT_TABLE produced\n"
"  in the same run.\n"
"\n"
"Output files:\n"
"\n"
"  DO category:               Data type:  Explanation:\n"
"  REDUCED_SCI_MXU            FITS image  Extracted scientific spectra\n"
"  REDUCED_SKY_SCI_MXU        FITS image  Extracted sky spectra\n"
"  REDUCED_ERROR_SCI_MXU      FITS image  Errors on extracted spectra\n"
"  UNMAPPED_SCI_MXU           FITS image  Sky subtracted scientific spectra\n"
"  MAPPED_SCI_MXU             FITS image  Rectified scientific spectra\n"
"  MAPPED_ALL_SCI_MXU         FITS image  Rectified science spectra with sky\n"
"  MAPPED_SKY_SCI_MXU         FITS image  Rectified sky spectra\n"
"  UNMAPPED_SKY_SCI_MXU       FITS image  Sky on CCD\n"
"  OBJECT_TABLE_SCI_MXU       FITS table  Positions of detected objects\n"
"\n"
"  Only if the global sky subtraction is requested:\n"
"  GLOBAL_SKY_SPECTRUM_MXU    FITS table  Global sky spectrum\n"
"\n"
"  Only if the sky-alignment of the wavelength solution is requested:\n"
"  SKY_SHIFTS_LONG_SCI_MXU    FITS table  Sky lines offsets (LSS-like data)\n"
"  or SKY_SHIFTS_SLIT_SCI_MXU FITS table  Sky lines offsets (MOS-like data)\n"
"  DISP_COEFF_SCI_MXU         FITS table  Upgraded dispersion coefficients\n"
"  WAVELENGTH_MAP_SCI_MXU     FITS image  Upgraded wavelength map\n"
"\n"
"  Only if a STANDARD_MXU is specified in input:\n"
"  SPECPHOT_TABLE             FITS table  Efficiency and response curves\n"
"\n"
"  Only if a photometric calibration was requested:\n"
"  REDUCED_FLUX_SCI_MXU       FITS image  Flux calibrated scientific spectra\n"
"  REDUCED_FLUX_ERROR_SCI_MXU FITS image  Errors on flux calibrated spectra\n"
"  MAPPED_FLUX_SCI_MXU        FITS image  Flux calibrated slit spectra\n\n";

#define fors_science_exit(message)            \
{                                             \
if ((const char *)message != NULL) cpl_msg_error(recipe, message);  \
cpl_free(exptime);                            \
cpl_image_delete(dummy);                      \
cpl_image_delete(mapped);                     \
cpl_image_delete(mapped_sky);                 \
cpl_image_delete(mapped_cleaned);             \
cpl_image_delete(skylocalmap);                \
cpl_image_delete(skymap);                     \
cpl_image_delete(smapped);                    \
cpl_table_delete(offsets);                    \
cpl_table_delete(photcal);                    \
cpl_table_delete(sky);                        \
fors_image_delete(&bias);                     \
cpl_image_delete(spectra);                    \
cpl_image_delete(coordinate);                 \
fors_image_delete(&norm_flat);                \
cpl_image_delete(rainbow);                    \
cpl_image_delete(rectified);                  \
cpl_image_delete(wavemap);                    \
cpl_propertylist_delete(qclist);              \
cpl_propertylist_delete(header);              \
cpl_propertylist_delete(save_header);         \
cpl_table_delete(grism_table);                \
cpl_table_delete(idscoeff);                   \
cpl_table_delete(maskslits);                  \
cpl_table_delete(overscans);                  \
cpl_table_delete(polytraces);                 \
cpl_table_delete(slits);                      \
cpl_table_delete(wavelengths);                \
cpl_vector_delete(lines);                     \
cpl_msg_indent_less();                        \
return -1;                                    \
}


#define fors_science_exit_memcheck(message)   \
{                                             \
if ((const char *)message != NULL) cpl_msg_info(recipe, message);   \
cpl_free(exptime);                            \
cpl_free(instrume);                           \
cpl_image_delete(dummy);                      \
cpl_image_delete(mapped);                     \
cpl_image_delete(mapped_cleaned);             \
cpl_image_delete(mapped_sky);                 \
cpl_image_delete(skylocalmap);                \
cpl_image_delete(skymap);                     \
cpl_image_delete(smapped);                    \
cpl_table_delete(offsets);                    \
cpl_table_delete(photcal);                    \
cpl_table_delete(sky);                        \
fors_image_delete(&bias);                     \
cpl_image_delete(spectra);                    \
cpl_image_delete(coordinate);                 \
fors_image_delete(&norm_flat);                \
cpl_image_delete(rainbow);                    \
cpl_image_delete(rectified);                  \
cpl_image_delete(wavemap);                    \
cpl_propertylist_delete(header);              \
cpl_propertylist_delete(save_header);         \
cpl_table_delete(grism_table);                \
cpl_table_delete(idscoeff);                   \
cpl_table_delete(maskslits);                  \
cpl_table_delete(overscans);                  \
cpl_table_delete(polytraces);                 \
cpl_table_delete(slits);                      \
cpl_table_delete(wavelengths);                \
cpl_vector_delete(lines);                     \
cpl_msg_indent_less();                        \
return 0;                                     \
}


/**
 * @brief    Build the list of available plugins, for this module. 
 *
 * @param    list    The plugin list
 *
 * @return   0 if everything is ok, -1 otherwise
 *
 * Create the recipe instance and make it available to the application 
 * using the interface. This function is exported.
 */

int cpl_plugin_get_info(cpl_pluginlist *plist)
{
    cpl_recipe *recipe = static_cast<cpl_recipe *>(cpl_calloc(1, sizeof *recipe ));
    cpl_plugin *plugin = &recipe->interface;

    cpl_plugin_init(plugin,
                    CPL_PLUGIN_API,
                    FORS_BINARY_VERSION,
                    CPL_PLUGIN_TYPE_RECIPE,
                    "fors_science",
                    "Extraction of scientific spectra",
                    fors_science_description,
                    "Carlo Izzo",
                    PACKAGE_BUGREPORT,
                    fors_get_license(),
                    fors_science_create,
                    fors_science_exec,
                    fors_science_destroy);

    cpl_pluginlist_append(plist, plugin);
    
    return 0;
}


/**
 * @brief    Setup the recipe options    
 *
 * @param    plugin  The plugin
 *
 * @return   0 if everything is ok
 *
 * Defining the command-line/configuration parameters for the recipe.
 */

static int fors_science_create(cpl_plugin *plugin)
{
    cpl_recipe    *recipe;
    cpl_parameter *p;


    /* 
     * Check that the plugin is part of a valid recipe 
     */

    if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE) 
        recipe = (cpl_recipe *)plugin;
    else 
        return -1;

    /* 
     * Create the parameters list in the cpl_recipe object 
     */

    recipe->parameters = cpl_parameterlist_new(); 


    /*
     * Dispersion
     */

    p = cpl_parameter_new_value("fors.fors_science.dispersion",
                                CPL_TYPE_DOUBLE,
                                "Resampling step (Angstrom/pixel)",
                                "fors.fors_science",
                                0.0);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "dispersion");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    /*
     * Sky lines alignment
     */

    p = cpl_parameter_new_value("fors.fors_science.skyalign",
                                CPL_TYPE_INT,
                                "Polynomial order for sky lines alignment, "
                                "or -1 to avoid alignment",
                                "fors.fors_science",
                                -1);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "skyalign");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    /*
     * Start wavelength for spectral extraction
     */

    p = cpl_parameter_new_value("fors.fors_science.startwavelength",
                                CPL_TYPE_DOUBLE,
                                "Start wavelength in spectral extraction",
                                "fors.fors_science",
                                0.0);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "startwavelength");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    /*
     * End wavelength for spectral extraction
     */

    p = cpl_parameter_new_value("fors.fors_science.endwavelength",
                                CPL_TYPE_DOUBLE,
                                "End wavelength in spectral extraction",
                                "fors.fors_science",
                                0.0);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "endwavelength");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    /*
     * Apply flat field
     */

    p = cpl_parameter_new_value("fors.fors_science.flatfield",
                                CPL_TYPE_BOOL,
                                "Apply flat field",
                                "fors.fors_science",
                                TRUE);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "flatfield");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    /*
     * Global sky subtraction
     */

    p = cpl_parameter_new_value("fors.fors_science.skyglobal",
                                CPL_TYPE_BOOL,
                                "Subtract global sky spectrum from CCD",
                                "fors.fors_science",
                                FALSE);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "skyglobal");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    /*
     * Local sky subtraction on extracted spectra
     */

/*** New sky subtraction (search NSS)
    p = cpl_parameter_new_value("fors.fors_science.skymedian",
                                CPL_TYPE_INT,
                                "Degree of sky fitting polynomial for "
                                "sky subtraction from extracted "
                                "slit spectra (MOS/MXU only, -1 to disable it)",
                                "fors.fors_science",
                                0);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "skymedian");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);
***/

    p = cpl_parameter_new_value("fors.fors_science.skymedian",
                                CPL_TYPE_BOOL,
                                "Sky subtraction from extracted slit spectra",
                                "fors.fors_science",
                                FALSE);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "skymedian");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    /*
     * Local sky subtraction on CCD spectra
     */

    p = cpl_parameter_new_value("fors.fors_science.skylocal",
                                CPL_TYPE_BOOL,
                                "Sky subtraction from CCD slit spectra",
                                "fors.fors_science",
                                TRUE);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "skylocal");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    /*
     * Cosmic rays removal
     */

    p = cpl_parameter_new_value("fors.fors_science.cosmics",
                                CPL_TYPE_BOOL,
                                "Eliminate cosmic rays hits (only if global "
                                "sky subtraction is also requested)",
                                "fors.fors_science",
                                FALSE);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "cosmics");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    /*
     * Slit margin
     */

    p = cpl_parameter_new_value("fors.fors_science.slit_margin",
                                CPL_TYPE_INT,
                                "Number of pixels to exclude at each slit "
                                "in object detection and extraction",
                                "fors.fors_science",
                                3);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "slit_margin");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    /*
     * Extraction radius
     */

    p = cpl_parameter_new_value("fors.fors_science.ext_radius",
                                CPL_TYPE_INT,
                                "Maximum extraction radius for detected "
                                "objects (pixel)",
                                "fors.fors_science",
                                12);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "ext_radius");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    /*
     * Contamination radius
     */

    p = cpl_parameter_new_value("fors.fors_science.cont_radius",
                                CPL_TYPE_INT,
                                "Minimum distance at which two objects "
                                "of equal luminosity do not contaminate "
                                "each other (pixel)",
                                "fors.fors_science",
                                0);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "cont_radius");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    /*
     * Object extraction method
     */

    p = cpl_parameter_new_value("fors.fors_science.ext_mode",
                                CPL_TYPE_INT,
                                "Object extraction method: 0 = aperture, "
                                "1 = Horne optimal extraction",
                                "fors.fors_science",
                                1);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "ext_mode");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    /*
     * Order of polynomial modeling the instrument response.
     */

    p = cpl_parameter_new_value("fors.fors_science.resp_fit_nknots",
                                CPL_TYPE_INT,
                                "Number of knots in spline fitting of the "
                                "instrument response. "
                                "(-1: No fitting. -2: Read from grism table)",
                                "fors.fors_science",
                                -2);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "resp_fit_nknots");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    p = cpl_parameter_new_value("fors.fors_science.resp_fit_degree",
                                CPL_TYPE_INT,
                                "Degree of polynomial in fitting of the "
                                "instrument response. "
                                "(-1: No fitting. -2: Read from grism table)",
                                "fors.fors_science",
                                -2);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "resp_fit_degree");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    p = cpl_parameter_new_value("fors.fors_science.resp_ignore_mode",
                                CPL_TYPE_STRING,
                                "Types of lines/regions to ignore in response. "
                                "Valid ones are 'stellar_absorption', "
                                "'telluric' and 'command_line' (from parameter resp_ignore_lines)",
                                "fors.fors_science",
                                "stellar_absorption,telluric,command_line");
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "resp_ignore_mode");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    p = cpl_parameter_new_value("fors.fors_science.resp_ignore_points",
                                CPL_TYPE_STRING,
                                "Extra lines/regions to ignore in response. "
                                "Use a comma separated list of values. A range "
                                "can be specified like 4500.0-4600.0",
                                "fors.fors_science",
                                "");
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "resp_ignore_points");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    p = cpl_parameter_new_value("fors.fors_science.resp_use_flat_sed",
                                CPL_TYPE_STRING,
                                "Use the flat SED to normalise the observed spectra. "
                                "Value are true, false, grism_table.",
                                "fors.fors_science",
                                "grism_table");
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "resp_use_flat_sed");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    return 0;
}


/**
 * @brief    Execute the plugin instance given by the interface
 *
 * @param    plugin  the plugin
 *
 * @return   0 if everything is ok
 */

static int fors_science_exec(cpl_plugin *plugin)
{
    cpl_recipe  *   recipe ;
    int             status = 1;

    /* Get the recipe out of the plugin */
    if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE)
        recipe = (cpl_recipe *)plugin ;
    else return -1 ;

    /* Issue a banner */
    fors_print_banner();

    try
    {
        status = fors_science(recipe->parameters, recipe->frames);
    }
    catch(std::exception& ex)
    {
        cpl_msg_error(cpl_func, "Recipe error: %s", ex.what());
    }
    catch(...)
    {
        cpl_msg_error(cpl_func, "An uncaught error during recipe execution");
    }

    return status;
}


/**
 * @brief    Destroy what has been created by the 'create' function
 *
 * @param    plugin  The plugin
 *
 * @return   0 if everything is ok
 */

static int fors_science_destroy(cpl_plugin *plugin)
{
    cpl_recipe *recipe;
    
    if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE) 
        recipe = (cpl_recipe *)plugin;
    else 
        return -1;

    cpl_parameterlist_delete(recipe->parameters); 

    return 0;
}


/**
 * @brief    Interpret the command line options and execute the data processing
 *
 * @param    parlist     The parameters list
 * @param    frameset    The set-of-frames
 *
 * @return   0 if everything is ok
 */

static int fors_science(cpl_parameterlist *parlist, cpl_frameset *frameset)
{

    const char *recipe = "fors_science";


    /*
     * Input parameters
     */

    double      dispersion;
    int         skyalign;
    const char *wcolumn = "WLEN";
    const char *resp_ignore_mode;
    const char *resp_ignore_lines;
    double      startwavelength;
    double      endwavelength;
    int         flatfield;
    int         skyglobal;
    int         skylocal;
    int         skymedian;
    int         cosmics;
    int         slit_margin;
    int         ext_radius;
    int         cont_radius;
    int         ext_mode;
    int         resp_fit_nknots;
    int         resp_fit_degree;
    int         photometry = 0;


    /*
     * CPL objects
     */

    fors_image_list  *all_science;
    cpl_image       **images;

    fors_image       *bias           = NULL;
    fors_image       *norm_flat      = NULL;
    fors_image       *science_ima    = NULL;
    fors_image       *smapped_ima    = NULL;
    cpl_image        *science_ima_nosky =NULL;
    cpl_image        *science_ima_nosky_var =NULL;
    cpl_image        *spectra        = NULL;
    cpl_image        *rectified      = NULL;
    cpl_image        *coordinate     = NULL;
    cpl_image        *rainbow        = NULL;
    cpl_image        *mapped         = NULL;
    cpl_image        *mapped_var     = NULL;
    cpl_image        *mapped_sky     = NULL;
    cpl_image        *mapped_sky_var = NULL;
    cpl_image        *mapped_cleaned = NULL;
    cpl_image        *smapped        = NULL;
    cpl_image        *smapped_var    = NULL;
    cpl_image        *smapped_nosky  = NULL;
    cpl_image        *smapped_nosky_var  = NULL;
    cpl_image        *wavemap        = NULL;
    cpl_image        *skymap         = NULL;
    cpl_image        *skylocalmap    = NULL;
    cpl_image        *dummy          = NULL;

    cpl_table        *grism_table    = NULL;
    cpl_table        *overscans      = NULL;
    cpl_table        *wavelengths    = NULL;
    cpl_table        *idscoeff       = NULL;
    cpl_table        *slits          = NULL;
    cpl_table        *maskslits      = NULL;
    cpl_table        *polytraces     = NULL;
    cpl_table        *offsets        = NULL;
    cpl_table        *sky            = NULL;
    cpl_table        *photcal        = NULL;
    cpl_table        *response_interp = NULL;

    cpl_vector       *lines          = NULL;

    cpl_propertylist *header         = NULL;
    cpl_propertylist *save_header    = NULL;
    cpl_propertylist *qclist         = NULL;
    const cpl_frame  *ref_sci_frame  = NULL;



    /*
     * Auxiliary variables
     */

    char        version[80];
    const char *instrume = NULL;
    const char *wheel4;
    const char *science_tag = NULL;
    const char *master_norm_flat_tag = NULL;
    const char *disp_coeff_tag = NULL;
    const char *disp_coeff_sky_tag = NULL;
    const char *wavelength_map_sky_tag = NULL;
    const char *curv_coeff_tag = NULL;
    const char *slit_location_tag = NULL;
    const char *reduced_science_tag = NULL;
    const char *reduced_flux_science_tag = NULL;
    const char *reduced_sky_tag = NULL;
    const char *reduced_error_tag = NULL;
    const char *reduced_flux_error_tag = NULL;
    const char *mapped_science_tag = NULL;
    const char *mapped_flux_science_tag = NULL;
    const char *unmapped_science_tag = NULL;
    const char *mapped_science_sky_tag = NULL;
    const char *mapped_sky_tag = NULL;
    const char *unmapped_sky_tag = NULL;
    const char *global_sky_spectrum_tag = NULL;
    const char *flat_disp_profile_tag = NULL;
    const char *object_table_tag = NULL;
    const char *skylines_offsets_tag = NULL;
    const char *global_distortion_table_tag = NULL;
    const char *specphot_tag;
    const char *master_specphot_tag = "MASTER_SPECPHOT_TABLE";
    const char *telluric_contamination_tag = "TELLURIC_CONTAMINATION";
    int         mxu, mos, lss;
    int         treat_as_lss = 0;
    int         have_phot = 0;
    int         nscience;
    double     *exptime = NULL;
    double      alltime;
    double      airmass = -1;
    double      mxpos;
    double      mean_rms;
    int         nlines;
    int         rebin;
    double     *line;
    int         nx, ny;
    int         ccd_xsize, ccd_ysize;
    double      ref_wave;
    double      gain;
    double      ron;
    int         standard;
    int         highres;
    int         i;
    double      wstart;
    double      wstep;
    int         wcount;
    int         global_distortion = 0;


    snprintf(version, 80, "%s-%s", PACKAGE, PACKAGE_VERSION);

    cpl_msg_set_indentation(2);


    /* 
     * Get configuration parameters
     */

    cpl_msg_info(recipe, "Recipe %s configuration parameters:", recipe);
    cpl_msg_indent_more();

    if (cpl_frameset_count_tags(frameset, "GRISM_TABLE") > 1)
        fors_science_exit("Too many in input: GRISM_TABLE");

    grism_table = dfs_load_table(frameset, "GRISM_TABLE", 1);

    dispersion = dfs_get_parameter_double(parlist, 
                    "fors.fors_science.dispersion", grism_table);

    if (dispersion <= 0.0)
        fors_science_exit("Invalid resampling step");

    skyalign = dfs_get_parameter_int(parlist, 
                    "fors.fors_science.skyalign", NULL);

    if (skyalign > 2)
        fors_science_exit("Max polynomial degree for sky alignment is 2");

    startwavelength = dfs_get_parameter_double(parlist, 
                    "fors.fors_science.startwavelength", grism_table);
    if (startwavelength < 3000.0 || startwavelength > 13000.0)
        fors_science_exit("Invalid wavelength");

    endwavelength = dfs_get_parameter_double(parlist, 
                    "fors.fors_science.endwavelength", grism_table);
    if (endwavelength < 3000.0 || endwavelength > 13000.0)
        fors_science_exit("Invalid wavelength");

    if (endwavelength - startwavelength <= 0.0)
        fors_science_exit("Invalid wavelength interval");

    flatfield = dfs_get_parameter_bool(parlist, "fors.fors_science.flatfield", 
                                       NULL);

    skyglobal = dfs_get_parameter_bool(parlist, "fors.fors_science.skyglobal", 
                                       NULL);
    skylocal  = dfs_get_parameter_bool(parlist, "fors.fors_science.skylocal", 
                                       NULL);
    skymedian = dfs_get_parameter_bool(parlist, "fors.fors_science.skymedian", 
                                       NULL);
/* NSS
    skymedian = dfs_get_parameter_int(parlist, "fors.fors_science.skymedian", 
                                       NULL);
*/

    if (skylocal && skyglobal)
        fors_science_exit("Cannot apply both local and global sky subtraction");

    if (skylocal && skymedian)
        fors_science_exit("Cannot apply sky subtraction both on extracted "
                          "and non-extracted spectra");

    cosmics = dfs_get_parameter_bool(parlist, 
                                     "fors.fors_science.cosmics", NULL);

    if (cosmics)
        if (!(skyglobal || skylocal))
            fors_science_exit("Cosmic rays correction requires "
                              "either skylocal=true or skyglobal=true");

    slit_margin = dfs_get_parameter_int(parlist, 
                                        "fors.fors_science.slit_margin",
                                        NULL);
    if (slit_margin < 0)
        fors_science_exit("Value must be zero or positive");

    ext_radius = dfs_get_parameter_int(parlist, "fors.fors_science.ext_radius",
                                       NULL);
    if (ext_radius < 0)
        fors_science_exit("Value must be zero or positive");

    cont_radius = dfs_get_parameter_int(parlist, 
                                        "fors.fors_science.cont_radius",
                                        NULL);
    if (cont_radius < 0)
        fors_science_exit("Value must be zero or positive");

    ext_mode = dfs_get_parameter_int(parlist, "fors.fors_science.ext_mode",
                                       NULL);
    if (ext_mode < 0 || ext_mode > 1)
        fors_science_exit("Invalid object extraction mode");

    resp_fit_nknots = dfs_get_parameter_int(parlist, "fors.fors_science.resp_fit_nknots",
                                            grism_table);
    if (resp_fit_nknots >= 0 &&  resp_fit_nknots < 2)
        fors_science_exit("Invalid instrument response spline knots");

    resp_fit_degree = dfs_get_parameter_int(parlist, "fors.fors_science.resp_fit_degree",
                                            grism_table);
    
    if (resp_fit_degree >= 0 && resp_fit_degree < 1)
        fors_science_exit("Invalid instrument response polynomial degree");
    
    if (resp_fit_degree > 0 && resp_fit_nknots > 0 )
        fors_science_exit("Only spline or polynomial fitting of response allowed, but not both");

    if (resp_fit_degree < 0 && resp_fit_nknots < 0 )
        fors_science_exit("Either spline or polynomial fitting of response has to be specified");

    resp_ignore_mode = dfs_get_parameter_string(parlist, 
                    "fors.fors_science.resp_ignore_mode", NULL);

    resp_ignore_lines = dfs_get_parameter_string(parlist, 
                    "fors.fors_science.resp_ignore_points", NULL);

    std::string resp_use_flat_sed = dfs_get_parameter_string(parlist, 
                    "fors.fors_science.resp_use_flat_sed", NULL);

    if (cpl_error_get_code())
        fors_science_exit("Failure getting the configuration parameters");

    
    /* 
     * Check input set-of-frames
     */

    cpl_msg_indent_less();
    cpl_msg_info(recipe, "Check input set-of-frames:");
    cpl_msg_indent_more();

    mxu = cpl_frameset_count_tags(frameset, "SCIENCE_MXU");
    mos = cpl_frameset_count_tags(frameset, "SCIENCE_MOS");
    lss = cpl_frameset_count_tags(frameset, "SCIENCE_LSS");
    standard = 0;

    if (mxu + mos + lss == 0) {
        mxu = cpl_frameset_count_tags(frameset, "STANDARD_MXU");
        mos = cpl_frameset_count_tags(frameset, "STANDARD_MOS");
        lss = cpl_frameset_count_tags(frameset, "STANDARD_LSS");
        standard = 1;
    }

    if (mxu + mos + lss == 0)
        fors_science_exit("Missing input scientific frame");

    nscience = mxu + mos + lss;

    if (mxu && mxu < nscience)
        fors_science_exit("Input scientific frames must be of the same type"); 

    if (mos && mos < nscience)
        fors_science_exit("Input scientific frames must be of the same type"); 

    if (lss && lss < nscience)
        fors_science_exit("Input scientific frames must be of the same type"); 

    if (mxu) {
        cpl_msg_info(recipe, "MXU data found");
        if (standard) {
            science_tag              = "STANDARD_MXU";
            reduced_science_tag      = "REDUCED_STD_MXU";
            reduced_flux_science_tag = "REDUCED_FLUX_STD_MXU";
            unmapped_science_tag     = "UNMAPPED_STD_MXU";
            mapped_science_tag       = "MAPPED_STD_MXU";
            mapped_flux_science_tag  = "MAPPED_FLUX_STD_MXU";
            mapped_science_sky_tag   = "MAPPED_ALL_STD_MXU";
            skylines_offsets_tag     = "SKY_SHIFTS_SLIT_STD_MXU";
            wavelength_map_sky_tag   = "WAVELENGTH_MAP_STD_MXU";
            disp_coeff_sky_tag       = "DISP_COEFF_STD_MXU";
            mapped_sky_tag           = "MAPPED_SKY_STD_MXU";
            unmapped_sky_tag         = "UNMAPPED_SKY_STD_MXU";
            object_table_tag         = "OBJECT_TABLE_STD_MXU";
            reduced_sky_tag          = "REDUCED_SKY_STD_MXU";
            reduced_error_tag        = "REDUCED_ERROR_STD_MXU";
            reduced_flux_error_tag   = "REDUCED_FLUX_ERROR_STD_MXU";
            specphot_tag             = "SPECPHOT_TABLE";
        }
        else {
            science_tag              = "SCIENCE_MXU";
            reduced_science_tag      = "REDUCED_SCI_MXU";
            reduced_flux_science_tag = "REDUCED_FLUX_SCI_MXU";
            unmapped_science_tag     = "UNMAPPED_SCI_MXU";
            mapped_science_tag       = "MAPPED_SCI_MXU";
            mapped_flux_science_tag  = "MAPPED_FLUX_SCI_MXU";
            mapped_science_sky_tag   = "MAPPED_ALL_SCI_MXU";
            skylines_offsets_tag     = "SKY_SHIFTS_SLIT_SCI_MXU";
            wavelength_map_sky_tag   = "WAVELENGTH_MAP_SCI_MXU";
            disp_coeff_sky_tag       = "DISP_COEFF_SCI_MXU";
            mapped_sky_tag           = "MAPPED_SKY_SCI_MXU";
            unmapped_sky_tag         = "UNMAPPED_SKY_SCI_MXU";
            object_table_tag         = "OBJECT_TABLE_SCI_MXU";
            reduced_sky_tag          = "REDUCED_SKY_SCI_MXU";
            reduced_error_tag        = "REDUCED_ERROR_SCI_MXU";
            reduced_flux_error_tag   = "REDUCED_FLUX_ERROR_SCI_MXU";
            specphot_tag             = "SPECPHOT_TABLE";
        }

        master_norm_flat_tag    = "MASTER_NORM_FLAT_MXU";
        disp_coeff_tag          = "DISP_COEFF_MXU";
        curv_coeff_tag          = "CURV_COEFF_MXU";
        slit_location_tag       = "SLIT_LOCATION_MXU";
        global_sky_spectrum_tag = "GLOBAL_SKY_SPECTRUM_MXU";
        flat_disp_profile_tag   = "FLAT_SED_MXU";

        if (!cpl_frameset_count_tags(frameset, master_norm_flat_tag)) {
            master_norm_flat_tag   = "MASTER_NORM_FLAT_LONG_MXU";
            disp_coeff_tag         = "DISP_COEFF_LONG_MXU";
            slit_location_tag      = "SLIT_LOCATION_LONG_MXU";
            flat_disp_profile_tag   = "FLAT_SED_LONG_MXU";
        }
    }

    if (lss) {
        cpl_msg_info(recipe, "LSS data found");

        if (cosmics && !skyglobal)
            fors_science_exit("Cosmic rays correction for LSS "
                              "data requires --skyglobal=true");

        if (standard) {
            science_tag              = "STANDARD_LSS";
            reduced_science_tag      = "REDUCED_STD_LSS";
            reduced_flux_science_tag = "REDUCED_FLUX_STD_LSS";
            unmapped_science_tag     = "UNMAPPED_STD_LSS";
            mapped_science_tag       = "MAPPED_STD_LSS";
            mapped_flux_science_tag  = "MAPPED_FLUX_STD_LSS";
            mapped_science_sky_tag   = "MAPPED_ALL_STD_LSS";
            skylines_offsets_tag     = "SKY_SHIFTS_SLIT_STD_LSS";
            wavelength_map_sky_tag   = "WAVELENGTH_MAP_STD_LSS";
            disp_coeff_sky_tag       = "DISP_COEFF_STD_LSS";
            mapped_sky_tag           = "MAPPED_SKY_STD_LSS";
            unmapped_sky_tag         = "UNMAPPED_SKY_STD_LSS";
            object_table_tag         = "OBJECT_TABLE_STD_LSS";
            reduced_sky_tag          = "REDUCED_SKY_STD_LSS";
            reduced_error_tag        = "REDUCED_ERROR_STD_LSS";
            reduced_flux_error_tag   = "REDUCED_FLUX_ERROR_STD_LSS";
            specphot_tag             = "SPECPHOT_TABLE";
        }
        else {
            science_tag              = "SCIENCE_LSS";
            reduced_science_tag      = "REDUCED_SCI_LSS";
            reduced_flux_science_tag = "REDUCED_FLUX_SCI_LSS";
            unmapped_science_tag     = "UNMAPPED_SCI_LSS";
            mapped_science_tag       = "MAPPED_SCI_LSS";
            mapped_flux_science_tag  = "MAPPED_FLUX_SCI_LSS";
            mapped_science_sky_tag   = "MAPPED_ALL_SCI_LSS";
            skylines_offsets_tag     = "SKY_SHIFTS_SLIT_SCI_LSS";
            wavelength_map_sky_tag   = "WAVELENGTH_MAP_SCI_LSS";
            disp_coeff_sky_tag       = "DISP_COEFF_SCI_LSS";
            mapped_sky_tag           = "MAPPED_SKY_SCI_LSS";
            unmapped_sky_tag         = "UNMAPPED_SKY_SCI_LSS";
            object_table_tag         = "OBJECT_TABLE_SCI_LSS";
            reduced_sky_tag          = "REDUCED_SKY_SCI_LSS";
            reduced_error_tag        = "REDUCED_ERROR_SCI_LSS";
            reduced_flux_error_tag   = "REDUCED_FLUX_ERROR_SCI_LSS";
            specphot_tag             = "SPECPHOT_TABLE";
        }

        master_norm_flat_tag        = "MASTER_NORM_FLAT_LSS";
        disp_coeff_tag              = "DISP_COEFF_LSS";
        slit_location_tag           = "SLIT_LOCATION_LSS";
        global_sky_spectrum_tag     = "GLOBAL_SKY_SPECTRUM_LSS";
        flat_disp_profile_tag       = "FLAT_SED_LSS";
    }

    if (mos) {
        cpl_msg_info(recipe, "MOS data found");
        if (standard) {
            science_tag              = "STANDARD_MOS";
            reduced_science_tag      = "REDUCED_STD_MOS";
            reduced_flux_science_tag = "REDUCED_FLUX_STD_MOS";
            unmapped_science_tag     = "UNMAPPED_STD_MOS";
            mapped_science_tag       = "MAPPED_STD_MOS";
            mapped_flux_science_tag  = "MAPPED_FLUX_STD_MOS";
            mapped_science_sky_tag   = "MAPPED_ALL_STD_MOS";
            skylines_offsets_tag     = "SKY_SHIFTS_SLIT_STD_MOS";
            wavelength_map_sky_tag   = "WAVELENGTH_MAP_STD_MOS";
            disp_coeff_sky_tag       = "DISP_COEFF_STD_MOS";
            mapped_sky_tag           = "MAPPED_SKY_STD_MOS";
            unmapped_sky_tag         = "UNMAPPED_SKY_STD_MOS";
            object_table_tag         = "OBJECT_TABLE_STD_MOS";
            reduced_sky_tag          = "REDUCED_SKY_STD_MOS";
            reduced_error_tag        = "REDUCED_ERROR_STD_MOS";
            reduced_flux_error_tag   = "REDUCED_FLUX_ERROR_STD_MOS";
            specphot_tag             = "SPECPHOT_TABLE";
        }
        else {
            science_tag              = "SCIENCE_MOS";
            reduced_science_tag      = "REDUCED_SCI_MOS";
            reduced_flux_science_tag = "REDUCED_FLUX_SCI_MOS";
            unmapped_science_tag     = "UNMAPPED_SCI_MOS";
            mapped_science_tag       = "MAPPED_SCI_MOS";
            mapped_flux_science_tag  = "MAPPED_FLUX_SCI_MOS";
            mapped_science_sky_tag   = "MAPPED_ALL_SCI_MOS";
            skylines_offsets_tag     = "SKY_SHIFTS_SLIT_SCI_MOS";
            wavelength_map_sky_tag   = "WAVELENGTH_MAP_SCI_MOS";
            disp_coeff_sky_tag       = "DISP_COEFF_SCI_MOS";
            mapped_sky_tag           = "MAPPED_SKY_SCI_MOS";
            unmapped_sky_tag         = "UNMAPPED_SKY_SCI_MOS";
            object_table_tag         = "OBJECT_TABLE_SCI_MOS";
            reduced_sky_tag          = "REDUCED_SKY_SCI_MOS";
            reduced_error_tag        = "REDUCED_ERROR_SCI_MOS";
            reduced_flux_error_tag   = "REDUCED_FLUX_ERROR_SCI_MOS";
            specphot_tag             = "SPECPHOT_TABLE";
        }

        master_norm_flat_tag    = "MASTER_NORM_FLAT_MOS";
        disp_coeff_tag          = "DISP_COEFF_MOS";
        curv_coeff_tag          = "CURV_COEFF_MOS";
        slit_location_tag       = "SLIT_LOCATION_MOS";
        global_sky_spectrum_tag = "GLOBAL_SKY_SPECTRUM_MOS";
        flat_disp_profile_tag   = "FLAT_SED_MOS";

        if (!cpl_frameset_count_tags(frameset, master_norm_flat_tag)) {
            master_norm_flat_tag   = "MASTER_NORM_FLAT_LONG_MOS";
            disp_coeff_tag         = "DISP_COEFF_LONG_MOS";
            slit_location_tag      = "SLIT_LOCATION_LONG_MOS";
            flat_disp_profile_tag   = "FLAT_SED_LONG_MOS";
        }
    }

    if (cpl_frameset_count_tags(frameset, "MASTER_BIAS") == 0)
        fors_science_exit("Missing required input: MASTER_BIAS");

    if (cpl_frameset_count_tags(frameset, "MASTER_BIAS") > 1)
        fors_science_exit("Too many in input: MASTER_BIAS");

    if (skyalign >= 0)
        if (cpl_frameset_count_tags(frameset, "MASTER_SKYLINECAT") > 1)
            fors_science_exit("Too many in input: MASTER_SKYLINECAT");

    if (cpl_frameset_count_tags(frameset, telluric_contamination_tag) > 1)
        fors_science_exit("Too many in input: TELLURIC_CONTAMINATION");
    
    if (cpl_frameset_count_tags(frameset, flat_disp_profile_tag) > 1)
        fors_science_exit("Too many in input: FLAT_DISP_PROFILE_*");

    if (cpl_frameset_count_tags(frameset, disp_coeff_tag) == 0) {
        cpl_msg_error(recipe, "Missing required input: %s", disp_coeff_tag);
        fors_science_exit(NULL);
    }

    if (cpl_frameset_count_tags(frameset, disp_coeff_tag) > 1) {
        cpl_msg_error(recipe, "Too many in input: %s", disp_coeff_tag);
        fors_science_exit(NULL);
    }

    if (cpl_frameset_count_tags(frameset, slit_location_tag) == 0) {
        cpl_msg_error(recipe, "Missing required input: %s",
                      slit_location_tag);
        fors_science_exit(NULL);
    }

    if (cpl_frameset_count_tags(frameset, slit_location_tag) > 1) {
        cpl_msg_error(recipe, "Too many in input: %s", slit_location_tag);
        fors_science_exit(NULL);
    }

    if (cpl_frameset_count_tags(frameset, master_norm_flat_tag) > 1) {
        if (flatfield) {
            cpl_msg_error(recipe, "Too many in input: %s", 
                          master_norm_flat_tag);
            fors_science_exit(NULL);
        }
        else {
            cpl_msg_warning(recipe, "%s in input are ignored, "
                            "since flat field correction was not requested", 
                            master_norm_flat_tag);
        }
    }

    if (cpl_frameset_count_tags(frameset, master_norm_flat_tag) == 1) {
        if (!flatfield) {
            cpl_msg_warning(recipe, "%s in input is ignored, "
                            "since flat field correction was not requested", 
                            master_norm_flat_tag);
        }
    }

    if (cpl_frameset_count_tags(frameset, master_norm_flat_tag) == 0) {
        if (flatfield) {
            cpl_msg_error(recipe, "Flat field correction was requested, "
                          "but no %s are found in input",
                          master_norm_flat_tag);
            fors_science_exit(NULL);
        }
    }

    if (standard) {
    
        if (cpl_frameset_count_tags(frameset, "EXTINCT_TABLE") == 0) {
            cpl_msg_warning(recipe, "An EXTINCT_TABLE was not found in input: "
                            "instrument response curve will not be produced.");
            standard = 0;
        }

        if (cpl_frameset_count_tags(frameset, "EXTINCT_TABLE") > 1)
            fors_science_exit("Too many in input: EXTINCT_TABLE");
    
        if (cpl_frameset_count_tags(frameset, "STD_FLUX_TABLE") == 0) {
            cpl_msg_warning(recipe, "A STD_FLUX_TABLE was not found in input: "
                            "instrument response curve will not be produced.");
            standard = 0;
        }

        if (cpl_frameset_count_tags(frameset, "STD_FLUX_TABLE") > 1)
            fors_science_exit("Too many in input: STD_FLUX_TABLE");

        if (!dfs_equal_keyword(frameset, "ESO OBS TARG NAME")) {
            cpl_msg_warning(recipe, "The target name of observation does not "
                            "match the standard star catalog: "
                            "instrument response curve will not be produced.");
            standard = 0;
        }
        //If we still "have" a standard, perform also photometry
        if(standard)
           photometry = 1;
    }

   
    have_phot  = cpl_frameset_count_tags(frameset, specphot_tag);
    have_phot += cpl_frameset_count_tags(frameset, master_specphot_tag);

    if (!standard) {
        if (have_phot == 0) {
            cpl_msg_info(recipe, 
                    "A SPECPHOT_TABLE was not found in input: "
                    "no photometric calibrated "
                    "spectra will be produced.");
            photometry = 0;
        }
        else
        {
            cpl_msg_info(recipe, 
                    "A SPECPHOT_TABLE was found in input: "
                    "photometric calibrated "
                    "spectra will be produced.");            
            photometry = 1;
        }

        if (photometry) {
            if (cpl_frameset_count_tags(frameset, "EXTINCT_TABLE") != 1)
                fors_science_exit("One and only one EXTINCT_TABLE is needed "
                        "to calibrate in photometry");


            if (have_phot > 1)
                fors_science_exit("Too many in input: SPECPHOT_TABLE");
        }
    }

    if (!dfs_equal_keyword(frameset, "ESO INS GRIS1 ID"))
        cpl_msg_warning(cpl_func,"Input frames are not from the same grism");

    if (!dfs_equal_keyword(frameset, "ESO INS FILT1 ID"))
        cpl_msg_warning(cpl_func,"Input frames are not from the same filter");

    if (cpl_frameset_count_tags(frameset, "SPECPHOT_TABLE"))
    {
        cpl_frameset * frameset_nostd = cpl_frameset_duplicate(frameset);
        cpl_frameset_erase(frameset_nostd, "SPECPHOT_TABLE");
        cpl_frameset_erase(frameset_nostd, "GLOBAL_DISTORTION_TABLE");
        if (!dfs_equal_keyword(frameset_nostd, "ESO DET CHIP1 ID"))
            cpl_msg_warning(cpl_func,"Input frames are not from the same chip."
                               " This does not apply to specphot table.");
        cpl_frameset_delete(frameset_nostd);
    }
    else
    {
        if (!dfs_equal_keyword(frameset, "ESO DET CHIP1 ID"))
            cpl_msg_warning(cpl_func,"Input frames are not from the same chip");
    }

    {
        cpl_frameset * frameset_detector = cpl_frameset_duplicate(frameset);
        cpl_frameset_erase(frameset_detector, "GRISM_TABLE");
        if (!dfs_equal_keyword(frameset_detector, "ESO DET CHIP1 NAME"))
            cpl_msg_warning(cpl_func,"Input frames are not from the same chip mosaic");
        cpl_frameset_delete(frameset_detector);
    }
 
    cpl_msg_indent_less();


    /*
     * Loading input data
     */

    /* Classify frames */
    fors_dfs_set_groups(frameset);
    const cpl_frame * bias_frame = 
            cpl_frameset_find_const(frameset, "MASTER_BIAS");
    
    exptime = (double * )cpl_calloc(nscience, sizeof(double));


    //Get the reference frame to inherit science frames
    ref_sci_frame = cpl_frameset_find_const(frameset, science_tag);

    if (nscience > 1) {

        cpl_msg_info(recipe, "Load %d scientific frames and median them...",
                     nscience);
        cpl_msg_indent_more();


        header = dfs_load_header(frameset, science_tag, 0);

        if (header == NULL)
            fors_science_exit("Cannot load scientific frame header");

        alltime = exptime[0] = cpl_propertylist_get_double(header, "EXPTIME");

        if (cpl_error_get_code() != CPL_ERROR_NONE)
            fors_science_exit("Missing keyword EXPTIME in scientific "
                              "frame header");

        if (standard || photometry) {
            airmass = fors_get_airmass(header);
            if (airmass < 0.0) 
                fors_science_exit("Missing airmass information in "
                                  "scientific frame header");
        }

        cpl_propertylist_delete(header); header = NULL;

        cpl_msg_info(recipe, "Scientific frame 1 exposure time: %.2f s", 
                     exptime[0]);

        for (i = 1; i < nscience; i++) {

            header = dfs_load_header(frameset, NULL, 0);

            if (header == NULL)
                fors_science_exit("Cannot load scientific frame header");
    
            exptime[i] = cpl_propertylist_get_double(header, "EXPTIME");

            alltime += exptime[i];
    
            if (cpl_error_get_code() != CPL_ERROR_NONE)
                fors_science_exit("Missing keyword EXPTIME in scientific "
                                  "frame header");
    
            cpl_propertylist_delete(header); header = NULL;

            cpl_msg_info(recipe, "Scientific frame %d exposure time: %.2f s", 
                         i+1, exptime[i]);
        }

        all_science = fors_image_list_new();
        cpl_frameset * science_frames = fors_frameset_extract(frameset, science_tag);

        for (i = 0; i < nscience; i++) 
        {
            cpl_frame * this_science_frame = 
                    cpl_frameset_get_position(science_frames, i);
            fors_image * this_science_ima = 
                    fors_image_load_preprocess(this_science_frame, bias_frame);

            if (this_science_ima) {
                fors_image_divide_scalar(this_science_ima, exptime[i], 0);
                fors_image_list_insert(all_science, this_science_ima); 
            }
            else
                fors_science_exit("Cannot load scientific frame");
        }

        hdrl_image *combined;
        cpl_image *contrib;
        hdrl_imagelist * all_science_hdrl = fors_image_list_to_hdrl(all_science);
        hdrl_parameter * combine_par = hdrl_collapse_median_parameter_create();
        hdrl_imagelist_collapse(all_science_hdrl, combine_par, 
                                &combined, &contrib);
        
        science_ima = fors_image_from_hdrl(combined);
        
        fors_image_multiply_scalar(science_ima, alltime, 0);

        hdrl_imagelist_delete(all_science_hdrl);
        fors_image_list_delete(&all_science, fors_image_delete);
        hdrl_image_delete(combined);
        hdrl_parameter_delete(combine_par);
        cpl_image_delete(contrib);
        cpl_frameset_delete(science_frames);
    }
    else {
        cpl_msg_info(recipe, "Load scientific exposure...");
        cpl_msg_indent_more();

        header = dfs_load_header(frameset, science_tag, 0);

        if (header == NULL)
            fors_science_exit("Cannot load scientific frame header");

        if (standard || photometry) {
            airmass = fors_get_airmass(header);
            if (airmass < 0.0) 
                fors_science_exit("Missing airmass information in "
                                  "scientific frame header");
        }

        /*
         * Insert here a check on supported filters:
         */

        wheel4 = cpl_propertylist_get_string(header, 
                                                     "ESO INS OPTI9 TYPE");
        if (cpl_error_get_code() != CPL_ERROR_NONE) {
            fors_science_exit("Missing ESO INS OPTI9 TYPE in flat header");
        }

        if (strcmp("FILT", wheel4) == 0) {
            wheel4 = cpl_propertylist_get_string(header,
                                                         "ESO INS OPTI9 NAME");
            cpl_msg_error(recipe, "Unsupported filter: %s", wheel4);
            fors_science_exit(NULL);
        }

        alltime = exptime[0] = cpl_propertylist_get_double(header, "EXPTIME");

        if (cpl_error_get_code() != CPL_ERROR_NONE)
            fors_science_exit("Missing keyword EXPTIME in scientific "
                              "frame header");

        cpl_propertylist_delete(header); header = NULL;

        cpl_msg_info(recipe, "Scientific frame exposure time: %.2f s", 
                     exptime[0]);

        const cpl_frame * science_frame = 
                cpl_frameset_find_const(frameset, science_tag);
        science_ima = fors_image_load_preprocess(science_frame, bias_frame);
        if(science_ima == NULL)
            fors_science_exit("Could not preprocess science image");
    }

    if (science_ima->data == NULL)
        fors_science_exit("Cannot load scientific frame");

    cpl_free(exptime); exptime = NULL;

    cpl_msg_indent_less();


    /*
     * Get the reference wavelength and the rebin factor along the
     * dispersion direction from a scientific exposure
     */

    header = dfs_load_header(frameset, science_tag, 0);

    if (header == NULL)
        fors_science_exit("Cannot load scientific frame header");

    instrume = cpl_propertylist_get_string(header, "INSTRUME");
    if (instrume == NULL)
        fors_science_exit("Missing keyword INSTRUME in scientific header");

    if (instrume[4] == '1')
        snprintf(version, 80, "%s/%s", "fors1", VERSION);
    if (instrume[4] == '2')
        snprintf(version, 80, "%s/%s", "fors2", VERSION);

    ref_wave = cpl_propertylist_get_double(header, "ESO INS GRIS1 WLEN");

    if (cpl_error_get_code() != CPL_ERROR_NONE)
        fors_science_exit("Missing keyword ESO INS GRIS1 WLEN in scientific "
                        "frame header");

    if (ref_wave < 3000.0)   /* Perhaps in nanometers... */
        ref_wave *= 10;

    if (ref_wave < 3000.0 || ref_wave > 13000.0) {
        cpl_msg_error(recipe, "Invalid central wavelength %.2f read from "
                      "keyword ESO INS GRIS1 WLEN in scientific frame header",
                      ref_wave);
        fors_science_exit(NULL);
    }

    cpl_msg_info(recipe, "The central wavelength is: %.2f", ref_wave);

    rebin = cpl_propertylist_get_int(header, "ESO DET WIN1 BINX");

    if (cpl_error_get_code() != CPL_ERROR_NONE)
        fors_science_exit("Missing keyword ESO DET WIN1 BINX in scientific "
                        "frame header");

    if (rebin != 1) {
        dispersion *= rebin;
        cpl_msg_warning(recipe, "The rebin factor is %d, and therefore the "
                        "resampling step used is %f A/pixel", rebin, 
                        dispersion);
        ext_radius /= rebin;
        cpl_msg_warning(recipe, "The rebin factor is %d, and therefore the "
                        "extraction radius used is %d pixel", rebin, 
                        ext_radius);
    }

    gain = cpl_propertylist_get_double(header, "ESO DET OUT1 CONAD");

    if (cpl_error_get_code() != CPL_ERROR_NONE)
        fors_science_exit("Missing keyword ESO DET OUT1 CONAD in scientific "
                          "frame header");

    cpl_msg_info(recipe, "The gain factor is: %.2f e-/ADU", gain);

    ron = cpl_propertylist_get_double(header, "ESO DET OUT1 RON");

    if (cpl_error_get_code() != CPL_ERROR_NONE)
        fors_science_exit("Missing keyword ESO DET OUT1 RON in scientific "
                          "frame header");

    ron /= gain;     /* Convert from electrons to ADU */

    cpl_msg_info(recipe, "The read-out-noise is: %.2f ADU", ron);

    if (mos || mxu) {
        int nslits_out_det = 0;
/* goto skip; */
        if (mos)
            maskslits = mos_load_slits_fors_mos(header, &nslits_out_det);
        else
            maskslits = mos_load_slits_fors_mxu(header);

        /*
         * Check if all slits have the same X offset: in such case,
         * treat the observation as a long-slit one!
         */

        mxpos = cpl_table_get_column_median(maskslits, "xtop");
     
        treat_as_lss = fors_mos_is_lss_like(maskslits, nslits_out_det);

        cpl_table_delete(maskslits); maskslits = NULL;
/* skip: */

        if (treat_as_lss) {
            cpl_msg_warning(recipe, "All MOS slits have the same offset: %.2f\n"
                            "The LSS data reduction strategy is applied!",
                            mxpos);
        }
    }

    if (lss || treat_as_lss) {
        if (skylocal) {
            if (cosmics)
                fors_science_exit("Cosmic rays correction for LSS or LSS-like "
                                  "data requires --skyglobal=true");
        }
        global_distortion_table_tag = "GLOBAL_DISTORTION_TABLE";
        if (cpl_frameset_count_tags(frameset, global_distortion_table_tag) == 0)
            global_distortion = 0;
        else if (cpl_frameset_count_tags(frameset, global_distortion_table_tag) > 1)
        {
            cpl_msg_error(recipe, "Too many in input: %s", curv_coeff_tag);
            fors_science_exit(NULL);
        }
        else
            global_distortion = 1;

    }
    else {
        if (cpl_frameset_count_tags(frameset, curv_coeff_tag) == 0) {
            cpl_msg_error(recipe, "Missing required input: %s", curv_coeff_tag);
            fors_science_exit(NULL);
        }

        if (cpl_frameset_count_tags(frameset, curv_coeff_tag) > 1) {
            cpl_msg_error(recipe, "Too many in input: %s", curv_coeff_tag);
            fors_science_exit(NULL);
        }
    }
    
    std::string grism_name = 
            cpl_propertylist_get_string(header, "ESO INS GRIS1 NAME");

    cpl_propertylist_delete(header); header = NULL;


    /*
     * Remove the master bias
     */

    cpl_msg_info(recipe, "Remove the master bias...");

    bias = fors_image_load(bias_frame);
    if (bias == NULL)
        fors_science_exit("Cannot load master bias");

    fors_subtract_bias(science_ima, bias);
    if (cpl_error_get_code() != CPL_ERROR_NONE)
        fors_science_exit("Cannot remove bias from scientific frame");

    fors_image_delete(&bias); bias = NULL;

    /* Flat-field the science */
    cpl_msg_indent_less();
    cpl_msg_info(recipe, "Load normalised flat field (if present)...");
    cpl_msg_indent_more();

    if (flatfield) {

        const cpl_frame * flat_frame = 
                cpl_frameset_find_const(frameset, master_norm_flat_tag);
        norm_flat = fors_image_load(flat_frame);
        
        //Substitute the zeros with ones, to avoid inf in division
        cpl_size npix = fors_image_get_size_x(norm_flat) *
                        fors_image_get_size_y(norm_flat);
        float * data_p = cpl_image_get_data_float(norm_flat->data);
        float * var_p  = cpl_image_get_data_float(norm_flat->variance);
        for(int pix = 0; pix < npix; ++pix)
        {
            if(data_p[pix] <= 0)
            {
                data_p[pix] = FLT_MAX;
                var_p[pix] = FLT_MAX;
            }
        }

        if (norm_flat) {
            cpl_msg_info(recipe, "Apply flat field correction...");
            fors_image_divide(science_ima, norm_flat);
            if (cpl_error_get_code() != CPL_ERROR_NONE) {
                cpl_msg_error(recipe, "Failure of flat field correction: %s",
                              cpl_error_get_message());
                fors_science_exit(NULL);
            }
            fors_image_delete(&norm_flat); norm_flat = NULL;
        }
        else {
            cpl_msg_error(recipe, "Cannot load input %s for flat field "
                          "correction", master_norm_flat_tag);
            fors_science_exit(NULL);
        }

    }
    spectra = science_ima->data;
    ccd_xsize = nx = cpl_image_get_size_x(science_ima->data);
    ccd_ysize = ny = cpl_image_get_size_y(science_ima->data);


    if (skyalign >= 0) {
        cpl_msg_indent_less();
        cpl_msg_info(recipe, "Load input sky line catalog...");
        cpl_msg_indent_more();

        wavelengths = dfs_load_table(frameset, "MASTER_SKYLINECAT", 1);

        if (wavelengths) {

            /*
             * Cast the wavelengths into a (double precision) CPL vector
             */

            nlines = cpl_table_get_nrow(wavelengths);

            if (nlines == 0)
                fors_science_exit("Empty input sky line catalog");

            if (cpl_table_has_column(wavelengths, wcolumn) != 1) {
                cpl_msg_error(recipe, "Missing column %s in input line "
                              "catalog table", wcolumn);
                fors_science_exit(NULL);
            }

            line = (double* )cpl_malloc(nlines * sizeof(double));
    
            for (i = 0; i < nlines; i++)
                line[i] = cpl_table_get(wavelengths, wcolumn, i, NULL);

            cpl_table_delete(wavelengths); wavelengths = NULL;

            lines = cpl_vector_wrap(nlines, line);
        }
        else {
            cpl_msg_info(recipe, "No sky line catalog found in input - fine!");
        }
    }


    /*
     * Load the slit location table
     */

    slits = dfs_load_table(frameset, slit_location_tag, 1);
    if (slits == NULL)
        fors_science_exit("Cannot load slits location table");


    /*
     * Load the wavelength calibration table
     */

    idscoeff = dfs_load_table(frameset, disp_coeff_tag, 1);
    if (idscoeff == NULL)
        fors_science_exit("Cannot load wavelength calibration table");

    cpl_msg_indent_less();
    cpl_msg_info(recipe, "Processing scientific spectra...");
    cpl_msg_indent_more();

    /*
     * Load the spectral curvature table, or provide a dummy one
     * in case of LSS or LSS-like data (single slit with flat curvature)
     */

    if (!(lss || treat_as_lss)) {
        polytraces = dfs_load_table(frameset, curv_coeff_tag, 1);
        if (polytraces == NULL)
            fors_science_exit("Cannot load spectral curvature table");
    }
    else
    {
        //Provide a dummy curvature table
        int * slit_id = cpl_table_get_data_int(slits, "slit_id");
        double * ytop    = cpl_table_get_data_double(slits, "ytop");
        double * ybottom = cpl_table_get_data_double(slits, "ybottom");
        polytraces = cpl_table_new(2);
        cpl_table_new_column(polytraces, "slit_id", CPL_TYPE_INT);
        cpl_table_new_column(polytraces, "c0", CPL_TYPE_DOUBLE);
        cpl_table_set_int(polytraces, "slit_id", 0, slit_id[0]);
        cpl_table_set_int(polytraces, "slit_id", 1, slit_id[0]);
        cpl_table_set_double(polytraces, "c0", 0, ytop[0]);
        cpl_table_set_double(polytraces, "c0", 1, ybottom[0]);
    }

    /*
     * This one will also generate the spatial map from the spectral 
     * curvature table (in the case of multislit data)
     */

    if (lss || treat_as_lss) {
        if(global_distortion) 
        {
            cpl_msg_info(cpl_func, "Applying global distortion to LSS-like data");
            cpl_table * distortion_table = dfs_load_table(frameset, 
                                              global_distortion_table_tag, 1);
            mosca::global_distortion global_dist(distortion_table); 
            smapped = global_dist.calibrate_spatial(science_ima->data, slits, ref_wave,
                    startwavelength, endwavelength, dispersion);
            smapped_var = global_dist.calibrate_spatial(science_ima->variance, slits, ref_wave,
                    startwavelength, endwavelength, dispersion);
            if(smapped == NULL)
                fors_science_exit("Could not apply spatial calibration");
            cpl_table_delete(distortion_table);
            smapped_ima = fors_image_new(smapped, smapped_var);
        }
        else 
        {
            int first_row = cpl_table_get_double(slits, "ybottom", 0, NULL);
            int last_row = cpl_table_get_double(slits, "ytop", 0, NULL);
            smapped = cpl_image_extract(science_ima->data, 1, first_row, nx, last_row);
            smapped_var = cpl_image_extract(science_ima->variance, 1, first_row, nx, last_row);
            smapped_ima = fors_image_new(smapped, smapped_var);
        }
    }
    else {
        cpl_msg_info(cpl_func, "Correcting spatial distortion");
        coordinate = cpl_image_new(nx, ny, CPL_TYPE_FLOAT);

        smapped = mos_spatial_calibration(science_ima->data, slits, polytraces, ref_wave,
                                          startwavelength, endwavelength,
                                          dispersion, 1, coordinate);
        smapped_var = mos_spatial_calibration(science_ima->variance, slits, polytraces, ref_wave,
                                          startwavelength, endwavelength,
                                          dispersion, 1, coordinate);
        smapped_ima = fors_image_new(smapped, smapped_var);
    }


    /*
     * Generate a rectified wavelength map from the wavelength calibration 
     * table
     */

    rainbow = mos_map_idscoeff(idscoeff, nx, ref_wave, startwavelength, 
                               endwavelength);

    if (dispersion > 1.0)
        highres = 0;
    else
        highres = 1;

    if (skyalign >= 0) {
        if (skyalign) {
            cpl_msg_info(recipe, "Align wavelength solution to reference "
            "skylines applying %d order residual fit...", skyalign);
        }
        else {
            cpl_msg_info(recipe, "Align wavelength solution to reference "
            "skylines applying median offset...");
        }

        offsets = mos_wavelength_align(smapped_ima->data, slits, ref_wave, 
                                       startwavelength, endwavelength, 
                                       idscoeff, lines, highres, skyalign, 
                                       rainbow, 4);

        if (offsets) {
            if (standard)
                cpl_msg_warning(recipe, "Alignment of the wavelength solution "
                                "to reference sky lines may be unreliable in "
                                "this case!");

            fors_dfs_save_table(frameset, offsets, skylines_offsets_tag, NULL, 
                               parlist, recipe, ref_sci_frame);
            if(cpl_error_get_code() != CPL_ERROR_NONE)
                fors_science_exit(NULL);

            cpl_table_delete(offsets); offsets = NULL;
        }
        else {
            if (cpl_error_get_code()) {
                if (cpl_error_get_code() == CPL_ERROR_INCOMPATIBLE_INPUT) {
                    cpl_msg_error(recipe, "The IDS coeff table is "
                    "incompatible with the input slit position table.");
                }
                cpl_msg_error(cpl_func, "Error found in %s: %s",
                              cpl_error_get_where(), cpl_error_get_message());
                fors_science_exit(NULL);
            }
            cpl_msg_warning(recipe, "Alignment of the wavelength solution "
                            "to reference sky lines could not be done!");
            skyalign = -1;
        }

    }

    if (lss || treat_as_lss) {
        int first_row = cpl_table_get_double(slits, "ybottom", 0, NULL);
        int last_row = cpl_table_get_double(slits, "ytop", 0, NULL);
        int ylow, yhig;

        ylow = first_row;
        yhig = last_row;

        wavemap = cpl_image_new(ccd_xsize, ccd_ysize, CPL_TYPE_FLOAT);
        cpl_image_copy(wavemap, rainbow, 1, ylow);

    }
    else {
        wavemap = mos_map_wavelengths(coordinate, rainbow, slits, 
                                      polytraces, ref_wave, 
                                      startwavelength, endwavelength,
                                      dispersion);
    }

    cpl_image_delete(rainbow); rainbow = NULL;
    cpl_image_delete(coordinate); coordinate = NULL;

    /*
     * Here the wavelength calibrated slit spectra are created. This frame
     * contains sky_science.
     */

    mapped_sky = mos_wavelength_calibration(smapped_ima->data, ref_wave,
                                            startwavelength, endwavelength,
                                            dispersion, idscoeff, 1);

    mapped_sky_var = mos_wavelength_calibration(smapped_ima->variance, ref_wave,
                                                startwavelength, endwavelength,
                                                dispersion, idscoeff, 1);

    cpl_msg_indent_less();
    cpl_msg_info(recipe, "Check applied wavelength against skylines...");
    cpl_msg_indent_more();

    mean_rms = mos_distortions_rms(mapped_sky, lines, startwavelength,
                                   dispersion, 6, highres);

    cpl_vector_delete(lines); lines = NULL;

    cpl_msg_info(recipe, "Mean residual: %f", mean_rms);

    mean_rms = cpl_table_get_column_mean(idscoeff, "error");

    cpl_msg_info(recipe, "Mean model accuracy: %f pixel (%f A)",
                 mean_rms, mean_rms * dispersion);

    header = cpl_propertylist_new();
    cpl_propertylist_update_double(header, "CRPIX1", 1.0);
    cpl_propertylist_update_double(header, "CRPIX2", 1.0);
    cpl_propertylist_update_double(header, "CRVAL1", 
                                   startwavelength + dispersion/2);
    cpl_propertylist_update_double(header, "CRVAL2", 1.0);
    /* cpl_propertylist_update_double(header, "CDELT1", dispersion);
    cpl_propertylist_update_double(header, "CDELT2", 1.0); */
    cpl_propertylist_update_double(header, "CD1_1", dispersion);
    cpl_propertylist_update_double(header, "CD1_2", 0.0);
    cpl_propertylist_update_double(header, "CD2_1", 0.0);
    cpl_propertylist_update_double(header, "CD2_2", 1.0);
    cpl_propertylist_update_string(header, "CTYPE1", "LINEAR");
    cpl_propertylist_update_string(header, "CTYPE2", "PIXEL");

    /* Perform time normalisation */
    cpl_propertylist_update_string(header, "BUNIT", "ADU/s");
    dummy = cpl_image_divide_scalar_create(mapped_sky, alltime);
    cpl_image * dummy_var = 
        cpl_image_divide_scalar_create(mapped_sky_var, alltime * alltime);
    fors_image * fors_mapped_sky = fors_image_new(
                   cpl_image_duplicate(dummy),
                   cpl_image_duplicate(dummy_var));
    fors_dfs_save_image_err(frameset, fors_mapped_sky, mapped_science_sky_tag,
             header, parlist, recipe, ref_sci_frame);
    if(cpl_error_get_code() != CPL_ERROR_NONE)
        fors_science_exit(NULL);
    cpl_image_delete(dummy); dummy = NULL;
    cpl_image_delete(dummy_var); dummy_var = NULL;
    fors_image_delete(&fors_mapped_sky);

/*    if (skyglobal == 0 && skymedian < 0) {    NSS */
    if (skyglobal == 0 && skymedian == 0 && skylocal == 0) {
        cpl_image_delete(mapped_sky); mapped_sky = NULL;
    }

    if (skyglobal || skylocal) {

        cpl_msg_indent_less();

        if (skyglobal) {
            cpl_msg_info(recipe, "Global sky determination...");
            cpl_msg_indent_more();
            skymap = cpl_image_new(nx, ny, CPL_TYPE_FLOAT);
            if (lss || treat_as_lss) {
                sky = mos_sky_map_super(science_ima->data, wavemap, dispersion, 
                                        2.0, 50, skymap);
            }
            else {
                sky = mos_sky_map_super(science_ima->data, wavemap, dispersion, 
                                        2.0, 50, skymap);
            }
            if (sky) {
                //TODO: This assumes that the sky computation has no error
                science_ima_nosky = cpl_image_subtract_create(science_ima->data, skymap);
                science_ima_nosky_var = cpl_image_duplicate(science_ima->variance);
            }
            else {
                cpl_image_delete(skymap); skymap = NULL;
            }
        }
        else {
            cpl_msg_info(recipe, "Local sky determination...");
            cpl_msg_indent_more();
            science_ima_nosky = cpl_image_duplicate(science_ima->data);
            skymap = mos_subtract_sky(science_ima_nosky, slits, polytraces, ref_wave,
                           startwavelength, endwavelength, dispersion);
            //TODO: This assumes that the sky computation has no error
            science_ima_nosky_var = cpl_image_duplicate(science_ima->variance);
        }

        if (skymap) {
            if (skyglobal) {
                /* Perform time normalisation */
                cpl_table_divide_scalar(sky, "sky", alltime);

                fors_dfs_save_table(frameset, sky, global_sky_spectrum_tag, 
                                    NULL, parlist, recipe, ref_sci_frame);
                if(cpl_error_get_code() != CPL_ERROR_NONE)
                    fors_science_exit(NULL);
    
                cpl_table_delete(sky); sky = NULL;
            }

            save_header = dfs_load_header(frameset, science_tag, 0);

            /* Perform time normalisation */
            cpl_image_divide_scalar(skymap, alltime);
            
            fors_dfs_save_image(frameset, skymap, unmapped_sky_tag,
                                save_header, parlist, recipe, ref_sci_frame);
            if(cpl_error_get_code() != CPL_ERROR_NONE)
                fors_science_exit(NULL);

            cpl_image_delete(skymap); skymap = NULL;

            /* Saving unmapped science */
            fors_image * fors_science_ima_nosky = fors_image_new(
                cpl_image_duplicate(science_ima_nosky),
                cpl_image_duplicate(science_ima_nosky_var));
            fors_dfs_save_image_err(frameset, fors_science_ima_nosky, 
                unmapped_science_tag, save_header, parlist, recipe, ref_sci_frame);
            if(cpl_error_get_code() != CPL_ERROR_NONE)
                fors_science_exit(NULL);
            fors_image_delete(&fors_science_ima_nosky);

            cpl_propertylist_delete(save_header); save_header = NULL;

            if (cosmics) {
                cpl_msg_info(recipe, "Removing cosmic rays...");
                mos_clean_cosmics(science_ima_nosky, gain, -1., -1.);
            }

            if (lss || treat_as_lss) {
                if(global_distortion) 
                {
                    cpl_table * distortion_table = dfs_load_table(frameset, 
                                                      global_distortion_table_tag, 1);
                    mosca::global_distortion global_dist(distortion_table); 
                    smapped_nosky = global_dist.calibrate_spatial(science_ima_nosky, slits, ref_wave,
                            startwavelength, endwavelength, dispersion);
                    smapped_nosky_var = global_dist.calibrate_spatial(science_ima_nosky_var, slits, ref_wave,
                            startwavelength, endwavelength, dispersion);
                    if(smapped_nosky == NULL || smapped_nosky_var == NULL)
                        fors_science_exit("Could not apply spatial calibration");
                    cpl_table_delete(distortion_table);
                }
                else 
                {
                    int first_row = cpl_table_get_double(slits, "ybottom", 0, NULL);
                    int last_row = cpl_table_get_double(slits, "ytop", 0, NULL);
                    smapped_nosky = cpl_image_extract(science_ima_nosky, 1, first_row, nx, last_row);
                    smapped_nosky_var = cpl_image_extract(science_ima_nosky_var, 1, first_row, nx, last_row);
                }
            }
            else {
                smapped_nosky = mos_spatial_calibration(science_ima_nosky, slits, polytraces, 
                                                  ref_wave, startwavelength, 
                                                  endwavelength, dispersion, 
                                                  1, NULL);
                smapped_nosky_var = mos_spatial_calibration(science_ima_nosky_var, slits, polytraces, 
                                                  ref_wave, startwavelength, 
                                                  endwavelength, dispersion, 1, NULL);
            }
        }
        else {
            cpl_msg_warning(recipe, "Sky subtraction failure");
            if (cosmics)
                cpl_msg_warning(recipe, "Cosmic rays removal not performed!");
            cosmics = skylocal = skyglobal = 0;
        }
    }
    else
    {
        smapped_nosky = cpl_image_duplicate(smapped_ima->data);
        smapped_nosky_var = cpl_image_duplicate(smapped_ima->variance);
    }

    //cpl_image_delete(spectra); spectra = NULL;

    if (skyalign >= 0) {
        save_header = dfs_load_header(frameset, science_tag, 0);
        cpl_propertylist_update_string(save_header, "BUNIT", "Angstrom");
        fors_dfs_save_image(frameset, wavemap, wavelength_map_sky_tag,
                            save_header, parlist, recipe, ref_sci_frame);
        if(cpl_error_get_code() != CPL_ERROR_NONE)
            fors_science_exit(NULL);
        cpl_propertylist_delete(save_header); save_header = NULL;
    }

    cpl_image_delete(wavemap); wavemap = NULL;

    mapped = mos_wavelength_calibration(smapped_nosky, ref_wave,
                                        startwavelength, endwavelength,
                                        dispersion, idscoeff, 1);

    mapped_var = mos_wavelength_calibration(smapped_nosky_var, ref_wave,
                                        startwavelength, endwavelength,
                                        dispersion, idscoeff, 1);

    //cpl_image_delete(smapped); smapped = NULL;

/*    if (skymedian >= 0) {    NSS */
    if (skymedian) {
            cpl_msg_indent_less();
            cpl_msg_info(recipe, "Local sky median determination...");
            cpl_msg_indent_more();
       
/*   NSS      skylocalmap = mos_sky_local(mapped, slits, skymedian); */
/*            skylocalmap = mos_sky_local(mapped, slits, 0);        */
            skylocalmap = mos_sky_local_old(mapped, slits);
            //TODO: The variance remains the same because the sky is assumed to have no error
            cpl_image_subtract(mapped, skylocalmap);
/*
            if (dfs_save_image(frameset, skylocalmap, mapped_sky_tag, header, 
                               parlist, recipe, version))
                fors_science_exit(NULL);
*/
            cpl_image_delete(skylocalmap); skylocalmap = NULL;
    }

    //Determine whether the response will contain the flat correction
    bool response_apply_flat_corr = 
            fors_science_response_apply_flat_corr
                (frameset, flat_disp_profile_tag, resp_use_flat_sed, 
                        grism_table, standard);

    //Determine wheter the flux calibration will use the flat correction
    bool photcal_apply_flat_corr = fors_science_photcal_apply_flat_corr
            (frameset, specphot_tag, master_specphot_tag,
             flat_disp_profile_tag, &photometry, 
             response_apply_flat_corr, standard);

    //TODO: Place this in a better place. The whole recipe needs refactoring...
    cpl_image        * mapped_flat_disp_profile = NULL;
    cpl_propertylist * flat_disp_header = NULL;
    if(photcal_apply_flat_corr || response_apply_flat_corr)
    {
        fors::detected_slits det_slits = 
                fors::detected_slits_from_tables(slits, polytraces, nx);
        mosca::wavelength_calibration wave_cal(idscoeff, ref_wave);
        mapped_flat_disp_profile = 
                fors_science_map_disp_profile(frameset, flat_disp_profile_tag, 
                          det_slits, wave_cal,
                          startwavelength, endwavelength, dispersion);
        flat_disp_header = dfs_load_header(frameset, flat_disp_profile_tag, 0);
    }
    
/*    if (skyglobal || skymedian >= 0 || skylocal) {   NSS */
    if (skyglobal || skymedian || skylocal) {

        skylocalmap = cpl_image_subtract_create(mapped_sky, mapped);

        cpl_image_delete(mapped_sky); mapped_sky = NULL;

        /* Perform time normalisation */
        cpl_propertylist_update_string(header, "BUNIT", "ADU/s");
        dummy = cpl_image_divide_scalar_create(skylocalmap, alltime);
        fors_dfs_save_image(frameset, dummy, mapped_sky_tag, header,
                parlist, recipe, ref_sci_frame);
        if(cpl_error_get_code() != CPL_ERROR_NONE)
            fors_science_exit(NULL);
        cpl_image_delete(dummy); dummy = NULL;

        cpl_msg_indent_less();
        cpl_msg_info(recipe, "Object detection...");
        cpl_msg_indent_more();

        if (cosmics || nscience > 1) {
            dummy = mos_detect_objects(mapped, slits, slit_margin, ext_radius, 
                                       cont_radius);
        }
        else {
            mapped_cleaned = cpl_image_duplicate(mapped);
            mos_clean_cosmics(mapped_cleaned, gain, -1., -1.);
            dummy = mos_detect_objects(mapped_cleaned, slits, slit_margin, 
                                       ext_radius, cont_radius);

            cpl_image_delete(mapped_cleaned); mapped_cleaned = NULL;
        }

        cpl_image_delete(dummy); dummy = NULL;

        fors_dfs_save_table(frameset, slits, object_table_tag, NULL, parlist, 
                           recipe, ref_sci_frame);
        if(cpl_error_get_code() != CPL_ERROR_NONE)
            fors_science_exit(NULL);

        cpl_msg_indent_less();
        cpl_msg_info(recipe, "Object extraction...");
        cpl_msg_indent_more();

        
        images = mos_extract_objects(mapped, mapped_var, skylocalmap, slits, 
                                     ext_mode, ron, gain, 1);

        cpl_image_delete(skylocalmap); skylocalmap = NULL;

        if (images) {
            if (standard) {
                cpl_msg_info(cpl_func, "Computing response...");

                if(response_apply_flat_corr)
                    cpl_msg_info(cpl_func, "  An additional flat-sed-corrected"
                            " response will be computed");
                cpl_table *ext_table  = NULL;
                cpl_table *flux_table = NULL;
                cpl_table *telluric_table = NULL;
                double flat_sed_norm_factor;

                ext_table = dfs_load_table(frameset, "EXTINCT_TABLE", 1);
                flux_table = dfs_load_table(frameset, "STD_FLUX_TABLE", 1);
                
                if(cpl_frameset_count_tags(frameset, telluric_contamination_tag))
                    telluric_table = dfs_load_table(frameset, "TELLURIC_CONTAMINATION", 1);

                std::vector<double> resp_ignore_lines_list;
                std::vector<std::pair<double, double> > resp_ignore_ranges_list;
                
                if(fors_science_response_fill_ignore(flux_table,
                                                     telluric_table,
                                                     grism_name,
                                                     resp_ignore_mode,
                                                     resp_ignore_lines,
                                                     resp_ignore_lines_list,
                                                     resp_ignore_ranges_list))
                    fors_science_exit("Cannot parse the response ignored lines");

                photcal = fors_compute_response(images[0], 
                                                mapped_flat_disp_profile,
                                                flat_disp_header,
                                                slits,
                                                startwavelength,
                                                dispersion, gain,
                                                alltime, ext_table,
                                                airmass, flux_table,
                                                resp_ignore_lines_list,
                                                resp_ignore_ranges_list,
                                                resp_fit_nknots,
                                                resp_fit_degree,
                                                response_interp,
                                                flat_sed_norm_factor);

                cpl_table_delete(ext_table);
                cpl_table_delete(flux_table);
                if(cpl_frameset_count_tags(frameset, telluric_contamination_tag))
                    cpl_table_delete(telluric_table);

                if (photcal) {

                    float *data;
                    char   keyname[30];

                    qclist = dfs_load_header(frameset, science_tag, 0);

                    if (qclist == NULL)
                        fors_science_exit("Cannot reload scientific "
                                "frame header");

                    /*
                     * QC1 parameters
                     */

                    wstart = 3700.;
                    wstep  = 400.;
                    wcount = 15;

                    dummy = cpl_image_new(wcount, 1, CPL_TYPE_FLOAT);
                    data = cpl_image_get_data_float(dummy);
                    map_table(dummy, wstart, wstep, photcal, 
                              "WAVE", "EFFICIENCY");

                    for (i = 0; i < wcount; i++) {
                        sprintf(keyname, "QC.SPEC.EFFICIENCY%d.LAMBDA", 
                                i + 1);
                        if (fors_header_write_double(qclist, 
                                wstart + wstep * i,
                                keyname, "Angstrom",
                                "Wavelength of "
                                "efficiency evaluation")) {
                            fors_science_exit("Cannot write qc param");
                        }

                        sprintf(keyname, "QC.SPEC.EFFICIENCY%d", i + 1);
                        if (fors_header_write_double(qclist,
                                data[i],
                                keyname, "e-/photon",
                                "Efficiency")) {
                            fors_science_exit("Cannot write qc param");
                        }
                    }

                    if(mapped_flat_disp_profile != NULL)
                    {
                        if (fors_header_write_int(qclist,
                                1,
                                "QC.RESP.FLAT_SED_CORR", "",
                                "Response corrected from flat dispersion profile")) {
                            fors_science_exit("Cannot write qc param");
                        }
                        if (fors_header_write_double(qclist,
                                flat_sed_norm_factor,
                                "QC.RESP.FLAT_SED_NORM", "",
                                "Normalisation factor applied to flat sed")) {
                            fors_science_exit("Cannot write qc param");
                        }
                    }
                        
                    cpl_image_delete(dummy); dummy = NULL;

                    fors_dfs_save_table(frameset, photcal, specphot_tag, qclist,
                                       parlist, recipe, ref_sci_frame);
                    //TODO: Create a fors_dfs_save_tables to save more than one
                    cpl_table_save(response_interp, NULL, NULL,
                                   "specphot_table.fits", CPL_IO_EXTEND);

                    if(cpl_error_get_code() != CPL_ERROR_NONE)
                        fors_science_exit(NULL);

                    cpl_propertylist_delete(qclist); qclist = NULL;

                    if (have_phot) {
                        cpl_table_delete(photcal); photcal = NULL;
                        cpl_table_delete(response_interp); response_interp = NULL;
                    }
                }
                else
                {
                    fors_science_exit("Couldn't compute the response");
                }
            }

            if (photometry) {
                cpl_image *calibrated;
                cpl_table *ext_table;
                cpl_propertylist * specphot_header;

                ext_table = dfs_load_table(frameset, "EXTINCT_TABLE", 1);

                if (cpl_frameset_count_tags(frameset, specphot_tag) == 0) 
                {
                    response_interp = dfs_load_table(frameset, 
                            master_specphot_tag, 2);
                    specphot_header = dfs_load_header(frameset, master_specphot_tag, 0);
                }
                else 
                {
                    response_interp = dfs_load_table(frameset, specphot_tag, 2);
                    specphot_header = dfs_load_header(frameset, specphot_tag, 0);
                }

                cpl_image * science_images = cpl_image_duplicate(images[0]);

                if(photcal_apply_flat_corr)
                    fors_science_correct_flat_sed(science_images, slits, 
                            mapped_flat_disp_profile, flat_disp_header, specphot_header);
                
                calibrated = mos_apply_photometry(science_images, response_interp, 
                                                  ext_table, startwavelength, 
                                                  dispersion, gain, alltime, 
                                                  airmass);
                cpl_propertylist_update_string(header, "BUNIT", 
                                   "10^(-16) erg/(cm^2 s Angstrom)");

                fors_dfs_save_image(frameset, calibrated,
                                    reduced_flux_science_tag, header,
                                    parlist, recipe, ref_sci_frame);
                if(cpl_error_get_code() != CPL_ERROR_NONE)
                {
                    cpl_image_delete(calibrated);
                    fors_science_exit(NULL);
                }

                cpl_image_delete(science_images);
                cpl_image_delete(calibrated);

                /* Apply photometry calibration to error */

 
                science_images = cpl_image_duplicate(images[0]);
                cpl_image * science_images_err = cpl_image_duplicate(images[2]);
                if(photcal_apply_flat_corr)
                {
                    //We assume that the mapped flat profile has no associated error
                    fors_science_correct_flat_sed(science_images, slits, 
                            mapped_flat_disp_profile, flat_disp_header, specphot_header);
                    fors_science_correct_flat_sed(science_images_err, slits, 
                            mapped_flat_disp_profile, flat_disp_header, specphot_header);
                }

                calibrated = mos_propagate_photometry_error(science_images, 
                                                  science_images_err, response_interp,
                                                  ext_table, startwavelength,
                                                  dispersion, gain, alltime,
                                                  airmass);

                cpl_propertylist_update_string(header, "BUNIT", 
                                   "10^(-16) erg/(cm^2 s Angstrom)");

                fors_dfs_save_image(frameset, calibrated,
                                    reduced_flux_error_tag, header,
                                    parlist, recipe, ref_sci_frame);
                if(cpl_error_get_code() != CPL_ERROR_NONE)
                {
                    cpl_image_delete(calibrated);
                    fors_science_exit(NULL);
                }

                cpl_image_delete(science_images);
                cpl_image_delete(science_images_err);
                cpl_table_delete(ext_table);
                cpl_image_delete(calibrated);
            }

    
            /* Perform time normalisation */
            cpl_propertylist_update_string(header, "BUNIT", "ADU/s");
            cpl_image_divide_scalar(images[0], alltime);

            fors_dfs_save_image(frameset, images[0], reduced_science_tag, 
                                header, parlist, recipe, ref_sci_frame);
            if(cpl_error_get_code() != CPL_ERROR_NONE)
                fors_science_exit(NULL);

            /* Perform time normalisation */
            cpl_image_divide_scalar(images[1], alltime);

            fors_dfs_save_image(frameset, images[1], reduced_sky_tag, header,
                                parlist, recipe, ref_sci_frame);
            if(cpl_error_get_code() != CPL_ERROR_NONE)
                fors_science_exit(NULL);

            cpl_image_delete(images[1]);

            /* Perform time normalisation */
            cpl_propertylist_update_string(header, "BUNIT", "ADU/s");
            cpl_image_divide_scalar(images[2], alltime);

            fors_dfs_save_image(frameset, images[2], reduced_error_tag, header,
                               parlist, recipe, ref_sci_frame);
            if(cpl_error_get_code() != CPL_ERROR_NONE)
                fors_science_exit(NULL);

            cpl_image_delete(images[0]);
            cpl_image_delete(images[2]);

            cpl_free(images);
        }
        else {
            cpl_msg_warning(recipe, "No objects found: the products "
                            "%s, %s, and %s are not created", 
                            reduced_science_tag, reduced_sky_tag, 
                            reduced_error_tag);
        }

    }

    if(mapped_flat_disp_profile != NULL)
    {
        cpl_image_delete(mapped_flat_disp_profile);
        cpl_propertylist_delete(flat_disp_header);
    }
    cpl_table_delete(slits); slits = NULL;

    if (skyalign >= 0) {
        fors_dfs_save_table(frameset, idscoeff, disp_coeff_sky_tag, NULL, 
                           parlist, recipe, ref_sci_frame);
        if(cpl_error_get_code() != CPL_ERROR_NONE)
            fors_science_exit(NULL);
    }

    cpl_table_delete(idscoeff); idscoeff = NULL;

    if (photometry && photcal) {
        cpl_image *calibrated;
        cpl_table *ext_table; 


        ext_table = dfs_load_table(frameset, "EXTINCT_TABLE", 1);

        calibrated = mos_apply_photometry(mapped, response_interp,
                ext_table, startwavelength,
                dispersion, gain, alltime,
                airmass);

        cpl_propertylist_update_string(header, "BUNIT", 
                                       "10^(-16) erg/(cm^2 s Angstrom)");

        fors_dfs_save_image(frameset, calibrated,
                mapped_flux_science_tag, header,
                parlist, recipe, ref_sci_frame);
        if(cpl_error_get_code() != CPL_ERROR_NONE)
        {
            cpl_image_delete(calibrated);
            fors_science_exit(NULL);
        }

        cpl_table_delete(ext_table);
        cpl_image_delete(calibrated);
    }

    /* Perform time normalisation */
    cpl_propertylist_update_string(header, "BUNIT", "ADU/s");
    cpl_image_divide_scalar(mapped, alltime);
    cpl_image_divide_scalar(mapped_var, alltime * alltime);

    /* Saving mapped science */
    fors_image * fors_mapped = fors_image_new(
            cpl_image_duplicate(mapped), cpl_image_duplicate(mapped_var));
    fors_dfs_save_image_err(frameset, fors_mapped, mapped_science_tag, header,
            parlist, recipe, ref_sci_frame);
    if(cpl_error_get_code() != CPL_ERROR_NONE)
        fors_science_exit(NULL);
    fors_image_delete(&fors_mapped);

    cpl_table_delete(photcal); photcal = NULL;
    cpl_table_delete(response_interp); response_interp = NULL;
    cpl_image_delete(mapped); mapped = NULL;
    cpl_propertylist_delete(header); header = NULL;
    cpl_table_delete(polytraces); polytraces = NULL;
    cpl_table_delete(grism_table); grism_table = NULL;


    if (cpl_error_get_code()) {
        cpl_msg_error(cpl_func, "Error found in %s: %s",
                      cpl_error_get_where(), cpl_error_get_message());
        fors_science_exit(NULL);
    }
    else 
        return 0;
    return 0;
}


static int fors_science_response_fill_ignore
(const cpl_table * flux_table, const cpl_table * telluric_table,
 const std::string& grism_name,
 const std::string& resp_ignore_mode, const std::string& resp_ignore_lines,
 std::vector<double>& resp_ignore_lines_list,
 std::vector<std::pair<double, double> >& resp_ignore_ranges_list)
{
    //Reading response mode
    bool mask_stellar_absorption = false;
    bool mask_telluric = false;
    bool mask_commandline = false;
    std::string mode(resp_ignore_mode);
    while(mode.length() > 0)
    {
        //Parsing ignore_lines (values are separated by comma)
        std::string::size_type found = mode.find(',');
        std::string mode_str;
        if(found != std::string::npos)
        {
            mode_str = mode.substr(0, found);
            mode = mode.substr(found+1);
        }
        else
        {
            mode_str = mode;
            mode = "";
        }
        if(mode_str == "stellar_absorption")
            mask_stellar_absorption = true;
        if(mode_str == "telluric")
            mask_telluric = true;
        if(mode_str == "command_line")
            mask_commandline = true;
    }
    
    //Adding lines from the standard star table
    if(mask_stellar_absorption)
    {
        if(cpl_table_has_column(flux_table, "STLLR_ABSORP"))
        {
            cpl_size stdtable_size = cpl_table_get_nrow(flux_table);
            for(cpl_size irow = 0; irow < stdtable_size; ++irow)
                if(cpl_table_get_int(flux_table, "STLLR_ABSORP", irow, NULL))
                    resp_ignore_lines_list.push_back
                    (cpl_table_get_float(flux_table, "WAVE", irow, NULL));
        }
        else
            cpl_msg_warning(cpl_func, " Column STLLR_ABSORP not found in std "
                    "star table. Value 'stellar_absorption' in 'resp_ignore_mode' is ignored.");
    }
    
    //Adding regions from the telluric contamination table
    if(mask_telluric && telluric_table != NULL)
    {
        if(cpl_table_has_column(telluric_table, grism_name.c_str()))
        {
            cpl_size telltable_size = cpl_table_get_nrow(telluric_table);
            for(cpl_size irow = 0; irow < telltable_size; ++irow)
                if(cpl_table_get_int(telluric_table, grism_name.c_str(), irow, NULL))
                {
                    double wave_start = 
                            cpl_table_get_float(telluric_table, "START_WAVE", irow, NULL);
                    double wave_end = 
                            cpl_table_get_float(telluric_table, "END_WAVE", irow, NULL);
                    resp_ignore_ranges_list.push_back(std::make_pair(wave_start, wave_end));
                }
        }
    }
    
    //Adding lines and ranges from the command line
    if(mask_commandline)
    {
        std::string lines(resp_ignore_lines);
        while(lines.length() > 0)
        {
            std::string::size_type found = lines.find(',');
            std::string line_str;
            if(found != std::string::npos)
            {
                line_str = lines.substr(0, found);
                lines = lines.substr(found+1);
            }
            else
            {
                line_str = lines;
                lines = "";
            }
            //We have a line or an interval. Let's check which of them
            std::string::size_type found_interval = line_str.find('-');
            if(found_interval != std::string::npos)     //It is an interval
            {
                double wave_start;
                std::istringstream iss1(line_str.substr(0, found_interval));
                //We cannot use simpy iss1 >> std::ws && iss1.eof() because in
                //some STL implementations iss >> std::ws gives a failure
                //if there are no more characters (the standard is not clear)
                //See http://llvm.org/bugs/show_bug.cgi?id=19497
                if ( !(iss1 >> wave_start) || !(iss1.eof() || (iss1 >> std::ws 
&& iss1.eof())) )
                {
                    cpl_msg_error(cpl_func, "Cannot interpret number in resp_ignore_lines");
                    return 1;
                }
                double wave_end;
                std::istringstream iss2(line_str.substr(found_interval + 1));
                if ( !(iss2 >> wave_end) || !(iss2.eof() || (iss2 >> std::ws && iss2.eof())) )
                {
                    cpl_msg_error(cpl_func, "Cannot interpret number in resp_ignore_lines");
                    return 1;
                }
                resp_ignore_ranges_list.push_back(std::make_pair(wave_start, wave_end));
            
            }
            else   //It is a single line
            {
                double wave;
                std::istringstream iss(line_str);
                if ( !(iss >> wave) || !(iss.eof() || (iss >> std::ws && iss.eof())) )
                {
                    cpl_msg_error(cpl_func, "Cannot interpret number in resp_ignore_lines");
                    return 1;
                }
                resp_ignore_lines_list.push_back(wave);
            }
        }
    }
    
    std::sort(resp_ignore_lines_list.begin(), resp_ignore_lines_list.end());  
    std::sort(resp_ignore_ranges_list.begin(), resp_ignore_ranges_list.end());  
    
    cpl_msg_indent_more();
    std::string all_lines;
    if(resp_ignore_lines_list.size() != 0)
    {
        std::ostringstream oss;
        for(size_t i = 0 ; i < resp_ignore_lines_list.size(); i++)
            oss<<resp_ignore_lines_list[i]<<" ";
        cpl_msg_info(cpl_func, "Total list of lines to ignore in reponse: %s", 
                     oss.str().c_str());
    }
    if(resp_ignore_ranges_list.size() != 0)
    {
        std::ostringstream oss;
        for(size_t i = 0 ; i < resp_ignore_ranges_list.size(); i++)
            oss<<resp_ignore_ranges_list[i].first<<"-"<<resp_ignore_ranges_list[i].second<<" ";
        cpl_msg_info(cpl_func, "Total list of ranges to ignore in reponse: %s", 
                     oss.str().c_str());
    }
    
    return 0;
}

static cpl_image * fors_science_map_disp_profile
(const cpl_frameset * frameset, const char * flat_disp_profile_tag, 
 fors::detected_slits& det_slits, mosca::wavelength_calibration& wave_cal,
 double firstLambda, double lastLambda, double dispersion)
{
    const cpl_frame * flat_disp_profile_frame = 
            cpl_frameset_find_const(frameset, flat_disp_profile_tag);
    const char * filename = cpl_frame_get_filename(flat_disp_profile_frame);

    cpl_image * flat_disp_profile =  
            cpl_image_load(filename, CPL_TYPE_FLOAT, 0, 0); 

    int nl = (lastLambda - firstLambda) / dispersion;
    
    cpl_image * mapped_flat_disp_profile = 
            cpl_image_new(nl, cpl_image_get_size_y(flat_disp_profile), CPL_TYPE_FLOAT);
    
    for(size_t i_slit = 0; i_slit < det_slits.size(); ++i_slit)
    {
        int slit_begin_pos = 
                det_slits[i_slit].get_position_spatial_corrected();
        int slit_end_pos = slit_begin_pos +  
                det_slits[i_slit].get_length_spatial_corrected();
        int slit_mean_pos = (slit_end_pos - slit_begin_pos) / 2;

        for(cpl_size i_wave = 0; i_wave < nl; ++i_wave)
        {
            int null;
            double wave = firstLambda + i_wave * dispersion;
            double pixel = wave_cal.get_pixel(slit_mean_pos, wave);
            int i_pix = std::ceil(pixel+0.5);
            if(i_pix>= 0 && i_pix <cpl_image_get_size_x(flat_disp_profile))
                cpl_image_set(mapped_flat_disp_profile, i_wave+1, i_slit+1, 
                        cpl_image_get(flat_disp_profile, i_pix+1, i_slit+1, &null));
        }
    }
    return mapped_flat_disp_profile;
}


#define MAX_COLNAME      (80)

static void fors_science_correct_flat_sed
(cpl_image *spectra,  cpl_table *objects,
 cpl_image * mapped_flat_disp_profile,
 cpl_propertylist * flat_sed_header,
 cpl_propertylist * specphot_header)
{

    cpl_size nx = cpl_image_get_size_x(spectra);
    cpl_size nslits = cpl_table_get_nrow(objects);
    int   maxobjects = 1;
    char  name[MAX_COLNAME];
    int   null;
    snprintf(name, MAX_COLNAME, "object_%d", maxobjects);
    while (cpl_table_has_column(objects, name)) {
        maxobjects++;
        snprintf(name, MAX_COLNAME, "object_%d", maxobjects);
    }
    maxobjects--;
    int idx_obj = 0;
    //TODO: Move it away from here. It is also repeated in the science photometric correction
    for (cpl_size i_slit = 0; i_slit < nslits; i_slit++) {
        /* Get the proper normalisation factor. The flat sed has to be
         * normalises by the same normalisation factor applied in the
         * creation of specphot.
         * So the original factor applied to the flat sed has to be removed 
         * and then aaply the one in the specphot */
        std::ostringstream norm_key;
        norm_key<< "ESO QC FLAT SED"<<i_slit+1<<" NORM";
        double flat_sed_norm_orig =
            cpl_propertylist_get_double(flat_sed_header, norm_key.str().c_str());
        double specphot_flat_sed_norm = 
            cpl_propertylist_get_double(specphot_header, "ESO QC RESP FLAT_SED_NORM");
        double flat_sed_final_norm_factor = flat_sed_norm_orig / specphot_flat_sed_norm;

        for (int i_obj = 1; i_obj <= maxobjects; i_obj++) {
            snprintf(name, MAX_COLNAME, "object_%d", i_obj);
            if (cpl_table_is_valid(objects, name, i_slit))
            {
                /* Divide the target spectrum by the corresponding slit profile */
                for (cpl_size i_pix = 0; i_pix < nx; i_pix++)
                {
                    double profile_val  = cpl_image_get(mapped_flat_disp_profile, i_pix+1, i_slit+1, &null);
                    if(profile_val != 0)
                        cpl_image_set(spectra, i_pix+1, idx_obj+1, 
                                cpl_image_get(spectra, i_pix+1, idx_obj+1, &null) / profile_val *
                                flat_sed_final_norm_factor);
                }
                idx_obj++;
            }
        }
    }
}

static bool fors_science_response_apply_flat_corr
(cpl_frameset * frameset, const char * flat_disp_profile_tag,
 std::string& resp_use_flat_sed, cpl_table * grism_table, int standard)
{
    if(standard)
    {
        bool have_flat_sed = false;
        if(cpl_frameset_count_tags(frameset, flat_disp_profile_tag) > 0)
            have_flat_sed= true;

        bool requested_in_grism_table = false;
        int null;
        if(cpl_table_get_int(grism_table, "RESP_USE_FLAT_SED", 0, &null))
            requested_in_grism_table = true;

        if(have_flat_sed)
        {
            if(resp_use_flat_sed == "false" || (resp_use_flat_sed == "grism_table" &&
                    !requested_in_grism_table))
            {
                cpl_msg_warning(cpl_func, "Flat SED is part of the input but "
                        "no correction has been requested");
                return false;
            }
            return true;
        }
        else
        {
            if(resp_use_flat_sed == "true" || (resp_use_flat_sed == "grism_table" &&
                    requested_in_grism_table))
                throw std::invalid_argument("Flat SED correction requested "
                        "but FLAT_SED it is not part of input.");
            return false;
        }
    }
    else
        return false;
    return false;

}

static bool fors_science_photcal_apply_flat_corr
(cpl_frameset * frameset, const char * specphot_tag, 
 const char * master_specphot_tag, const char * flat_disp_profile_tag,
 int* fluxcal, bool response_apply_flat_corr, int standard)
{
    
    //We are correcting a stdstar. Use what it was used to compute the response
    if(standard)
        return response_apply_flat_corr;
    else //We have a science
    {
        if(*fluxcal) //Logic only if flux calibration is going to be performed
        {
            bool have_flat_sed = false;
            if(cpl_frameset_count_tags(frameset, flat_disp_profile_tag) > 0)
                have_flat_sed= true;
            
            //Check if input response has been flat-sed-corrected 
            bool resp_has_flat_corr = false;
            cpl_table * response_interp;
            if (cpl_frameset_count_tags(frameset, master_specphot_tag) != 0) {
                response_interp = dfs_load_table(frameset, 
                        master_specphot_tag, 2);
                if(cpl_table_has_column(response_interp, "RESPONSE_FFSED"))
                    resp_has_flat_corr = true;
            }
            else if (cpl_frameset_count_tags(frameset, specphot_tag) != 0){
                response_interp = dfs_load_table(frameset, specphot_tag, 2);
                if(cpl_table_has_column(response_interp, "RESPONSE_FFSED"))
                    resp_has_flat_corr = true;
            }

            if(resp_has_flat_corr && !have_flat_sed)
            {
                cpl_msg_warning(cpl_func, "The response is corrected with the"
                   " FLAT_SED but there is no FLAT_SED provided for the "
                   "science frame. Therefore the data are not flux-calibrated");
                *fluxcal = 0;
                return false;
            }
            
            if(!resp_has_flat_corr && have_flat_sed)
            {
                cpl_msg_warning(cpl_func, "FLAT_SED in science sof ignored as"
                        " no FLAT_SED correction was applied to the response");
                return false;
            }
            
            if(!resp_has_flat_corr && !have_flat_sed)
                return false;
            
            if(resp_has_flat_corr && have_flat_sed)
                return true;
        }
        else        //No flux calibration 
            return false;
    }
    return false;
}