/* $Id: fors_normalise_flat.c,v 1.10 2013/09/09 12:24:55 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/09/09 12:24:55 $
 * $Revision: 1.10 $
 * $Name:  $
 */

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

#include <math.h>
#include <cpl.h>
#include <moses.h>
#include <fors_dfs.h>
#include <fors_utils.h>
#include <fors_flat_normalise.h>
#include "wavelength_calibration.h"

static int fors_normalise_flat_create(cpl_plugin *);
static int fors_normalise_flat_exec(cpl_plugin *);
static int fors_normalise_flat_destroy(cpl_plugin *);
static int fors_normalise_flat(cpl_parameterlist *, cpl_frameset *);

static char fors_normalise_flat_description[] =
"This recipe is used to normalise a master flat field frame dividing it\n"
"by its large scale illumination trend. This recipe can be applied both\n"
"to generic multi-slit (MOS/MXU) and to long slit exposures (either LSS, or\n"
"LSS-like MOS/MXU), even if different normalisation methods are applied in\n"
"such different cases. The input master flat field image is the product\n"
"of the recipe fors_flat. The input spectral curvature table, product of\n"
"the recipe fors_detect_spectra, is only required in the case of multi-slit\n"
"data.\n"
"\n"
"In the case of multi-slit data, the flat field spectra are spatially\n"
"rectified, heavily smoothed, and then mapped back on the CCD. Then the\n"
"master flat image is divided by its smoothed counterpart. The smoothing\n"
"may be obtained either by applying a running median filter of specified\n"
"sizes, or by polynomial fitting along the dispersion direction performed\n"
"independently for each row of the spatially remapped spectra.\n"
"\n"
"In the case of long-slit data, the smoothing can still be obtained either\n"
"by applying a running median filter or by polynomial fitting, but the\n"
"polynomial fitting will be performed along the spatial direction, for\n"
"each column of the spectrum.\n"
"\n"
"In the table below the MXU acronym can be alternatively read as MOS or\n"
"LSS.\n\n"
"Input files:\n\n"
"  DO category:               Type:       Explanation:         Required:\n"
"  MASTER_SCREEN_FLAT_MXU     Calib       Master flat frame       Y\n"
"  DISP_COEFF_MXU             Calib       Wavelength calibration  Y\n"        
"  CURV_COEFF_MXU             Calib       Spectral curvature      .\n"
"  SLIT_LOCATION_MXU          Calib       Slit detection          .\n"        
"  GRISM_TABLE                Calib       Grism table             Y\n\n"
"Output files:\n\n"
"  DO category:               Data type:  Explanation:\n"
"  MASTER_NORM_FLAT_MXU       FITS image  Normalised flat field\n\n";

#define fors_normalise_flat_exit(message)     \
{                                             \
if ((const char *)message != NULL) cpl_msg_error(recipe, message);  \
cpl_image_delete(spatial);                    \
cpl_image_delete(coordinate);                 \
cpl_image_delete(smo_flat);                   \
cpl_table_delete(grism_table);                \
cpl_table_delete(maskslits);                  \
cpl_table_delete(slits);                      \
cpl_table_delete(polytraces);                 \
cpl_propertylist_delete(header);              \
cpl_msg_indent_less();                        \
return -1;                                    \
}

#define fors_normalise_flat_exit_memcheck(message)  \
{                                               \
if ((const char *)message != NULL) cpl_msg_info(recipe, message);     \
printf("free master_flat (%p)\n", master_flat); \
cpl_image_delete(master_flat);                  \
printf("free spatial (%p)\n", spatial);         \
cpl_image_delete(spatial);                      \
printf("free coordinate (%p)\n", coordinate);   \
cpl_image_delete(coordinate);                   \
printf("free smo_flat (%p)\n", smo_flat);       \
cpl_image_delete(smo_flat);                     \
printf("free grism_table (%p)\n", grism_table); \
cpl_table_delete(grism_table);                  \
printf("free maskslits (%p)\n", maskslits);     \
cpl_table_delete(maskslits);                    \
printf("free slits (%p)\n", slits);             \
cpl_table_delete(slits);                        \
printf("free polytraces (%p)\n", polytraces);   \
cpl_table_delete(polytraces);                   \
printf("free header (%p)\n", header);           \
cpl_propertylist_delete(header);                \
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 *p_list)
{
    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_normalise_flat",
                    "Normalise master flat spectrum",
                    fors_normalise_flat_description,
                    "Carlo Izzo",
                    PACKAGE_BUGREPORT,
    "This file is currently part of the FORS Instrument Pipeline\n"
    "Copyright (C) 2002-2010 European Southern Observatory\n\n"
    "This program is free software; you can redistribute it and/or modify\n"
    "it under the terms of the GNU General Public License as published by\n"
    "the Free Software Foundation; either version 2 of the License, or\n"
    "(at your option) any later version.\n\n"
    "This program is distributed in the hope that it will be useful,\n"
    "but WITHOUT ANY WARRANTY; without even the implied warranty of\n"
    "MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the\n"
    "GNU General Public License for more details.\n\n"
    "You should have received a copy of the GNU General Public License\n"
    "along with this program; if not, write to the Free Software Foundation,\n"
    "Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA\n",
                    fors_normalise_flat_create,
                    fors_normalise_flat_exec,
                    fors_normalise_flat_destroy);

    cpl_pluginlist_append(p_list, 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_normalise_flat_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 (empty) parameters list in the cpl_recipe object 
     */

    recipe->parameters = cpl_parameterlist_new(); 

    /*
     * Dispersion
     */

    p = cpl_parameter_new_value("fors.fors_normalise_flat.dispersion",
                                CPL_TYPE_DOUBLE,
                                "Expected spectral dispersion (Angstrom/pixel)",
                                "fors.fors_normalise_flat",
                                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);

    /*
     * Start wavelength for spectral extraction
     */

    p = cpl_parameter_new_value("fors.fors_normalise_flat.startwavelength",
                                CPL_TYPE_DOUBLE,
                                "Start wavelength in spectral extraction",
                                "fors.fors_normalise_flat",
                                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_normalise_flat.endwavelength",
                                CPL_TYPE_DOUBLE,
                                "End wavelength in spectral extraction",
                                "fors.fors_normalise_flat",
                                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);

    /*
     * Number of knots in flat field fitting splines along spatial direction 
     */

    p = cpl_parameter_new_value("fors.fors_normalise_flat.spa_polydegree",
                                CPL_TYPE_INT,
                                "Polynomial degree for the flat field fitting "
                                "along spatial direction",
                                "fors.fors_normalise_flat",
                                -1);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "spa_polydegree");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    /*
     * Smooth box radius for flat field along spatial direction
     * (if spa_knots < 0)
     */

    p = cpl_parameter_new_value("fors.fors_normalise_flat.sradius",
                                CPL_TYPE_INT,
                                "Smooth box radius for flat field along "
                                "spatial direction (used if spa_knots < 0)",
                                "fors.fors_normalise_flat",
                                10);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "sradius");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    /*
     * Number of knots in flat field fitting splines along dispersion direction 
     */

    p = cpl_parameter_new_value("fors.fors_normalise_flat.disp_nknots",
                                CPL_TYPE_INT,
                                "Number of knots in flat field fitting "
                                "splines along dispersion direction",
                                "fors.fors_normalise_flat",
                                -1);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "disp_nknots");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    /*
     * Smooth box radius for flat field along dispersion direction
     */

    p = cpl_parameter_new_value("fors.fors_normalise_flat.dradius",
                                CPL_TYPE_INT,
                                "Smooth box radius for flat field along "
                                "dispersion direction (if disp_knots < 0)",
                                "fors.fors_normalise_flat",
                                10);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "dradius");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    /*
     * Threshold percentage for flat spline fitting with respect to the maximum
     */

    p = cpl_parameter_new_value("fors.fors_normalise_flat.splfit_threshold",
                                CPL_TYPE_DOUBLE,
                                "Threshold percentage for flat spline fitting"
                                "with respect to the maximum",
                                "fors.fors_normalise_flat",
                                0.01);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "splfit_threshold");
    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_normalise_flat_exec(cpl_plugin *plugin)
{
    cpl_recipe *recipe;
    
    if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE) 
        recipe = (cpl_recipe *)plugin;
    else 
        return -1;

    /* Issue a banner */
    fors_print_banner();

    return fors_normalise_flat(recipe->parameters, recipe->frames);
}


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

static int fors_normalise_flat_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_normalise_flat(cpl_parameterlist *parlist, 
                               cpl_frameset *frameset)
{

    const char *recipe = "fors_normalise_flat";


    /*
     * Input parameters
     */

    double      dispersion;
    double      startwavelength;
    double      endwavelength;
    int         spa_polyorder;
    int         disp_nknots;
    int         sradius;
    int         dradius;
    float       splfit_threshold;

    /*
     * CPL objects
     */

    fors_image       *master_flat_f= NULL;
    cpl_image        *smo_flat     = NULL;
    cpl_image        *coordinate   = NULL;
    cpl_image        *spatial      = NULL;
    cpl_table        *grism_table  = NULL;
    cpl_table        *slits        = NULL;
    cpl_table        *polytraces   = NULL;
    cpl_table        *maskslits    = NULL;
    cpl_propertylist *header       = NULL;

    /*
     * Auxiliary variables
     */

    char        version[80];
    const char *master_flat_tag = NULL;
    const char *master_norm_flat_tag = NULL;
    const char *slit_location_tag = NULL;
    const char *curv_coeff_tag = NULL;
    const char *disp_coeff_tag = NULL;
    int         mxu, mos, lss;
    int         nflat;
    int         rebin;
    int         nx, ny;
    int         treat_as_lss = 0;
    double      reference;

    const char *instrume = NULL;


    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_normalise_flat_exit("Too many in input: GRISM_TABLE"); 

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

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

    if (dispersion <= 0.0)
        fors_normalise_flat_exit("Invalid spectral dispersion value");

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

    endwavelength = dfs_get_parameter_double(parlist,
                    "fors.fors_normalise_flat.endwavelength", grism_table);
    if (endwavelength > 1.0) {
        if (endwavelength < 3000.0 || endwavelength > 13000.0)
            fors_normalise_flat_exit("Invalid wavelength");
        if (startwavelength < 1.0)
            fors_normalise_flat_exit("Invalid wavelength interval");
    }

    if (startwavelength > 1.0)
        if (endwavelength - startwavelength <= 0.0)
            fors_normalise_flat_exit("Invalid wavelength interval");

    spa_polyorder = dfs_get_parameter_int(parlist,
            "fors.fors_normalise_flat.spa_polydegree", NULL);
    disp_nknots = dfs_get_parameter_int(parlist,
            "fors.fors_normalise_flat.disp_nknots", NULL);
    sradius = dfs_get_parameter_int(parlist,
            "fors.fors_normalise_flat.sradius", NULL);
    dradius = dfs_get_parameter_int(parlist, 
            "fors.fors_normalise_flat.dradius", NULL);
    splfit_threshold = dfs_get_parameter_double(parlist, 
            "fors.fors_normalise_flat.splfit_threshold", NULL);

    cpl_table_delete(grism_table); grism_table = NULL;

    if (cpl_error_get_code())
        fors_normalise_flat_exit("Failure reading the configuration "
                                 "parameters");


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

    nflat  = mxu = cpl_frameset_count_tags(frameset, "MASTER_SCREEN_FLAT_MXU");
    nflat += mos = cpl_frameset_count_tags(frameset, "MASTER_SCREEN_FLAT_MOS");
    nflat += lss = cpl_frameset_count_tags(frameset, "MASTER_SCREEN_FLAT_LSS");

    if (nflat == 0) {
        fors_normalise_flat_exit("Missing input master flat field frame");
    }
    if (nflat > 1) {
        cpl_msg_error(recipe, "Too many input flat frames (%d > 1)", nflat);
        fors_normalise_flat_exit(NULL);
    }

    if (mxu) {
        master_flat_tag      = "MASTER_SCREEN_FLAT_MXU";
        master_norm_flat_tag = "MASTER_NORM_FLAT_MXU";
        slit_location_tag    = "SLIT_LOCATION_MXU";
        curv_coeff_tag       = "CURV_COEFF_MXU";
        disp_coeff_tag       = "DISP_COEFF_MXU";
    }
    else if (mos) {
        master_flat_tag      = "MASTER_SCREEN_FLAT_MOS";
        master_norm_flat_tag = "MASTER_NORM_FLAT_MOS";
        slit_location_tag    = "SLIT_LOCATION_MOS";
        curv_coeff_tag       = "CURV_COEFF_MOS";
        disp_coeff_tag       = "DISP_COEFF_MOS";
    }
    else if (lss) {
        master_flat_tag      = "MASTER_SCREEN_FLAT_LSS";
        master_norm_flat_tag = "MASTER_NORM_FLAT_LSS";
        disp_coeff_tag       = "DISP_COEFF_LSS";
    }

    header = dfs_load_header(frameset, master_flat_tag, 0);

    if (mos || mxu) {
        int nslits_out_det = 0;
        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!
         */

        treat_as_lss = fors_mos_is_lss_like(maskslits, nslits_out_det); 

        cpl_table_delete(maskslits); maskslits = NULL;

        if (treat_as_lss)
            cpl_msg_warning(recipe, "All MOS slits have the same offset. "
                            "The LSS data reduction strategy is applied!");
    }

    /* Check for required inputs */
    if (cpl_frameset_count_tags(frameset, disp_coeff_tag) != 0) {
        cpl_msg_error(recipe, "One (and only one) %s must be in inputs", 
                      slit_location_tag);
        fors_normalise_flat_exit(NULL);
    }
    if (!(lss || treat_as_lss)) {
        if (cpl_frameset_count_tags(frameset, curv_coeff_tag) != 1) {
            cpl_msg_error(recipe, "One (and only one) %s must be in inputs", 
                          curv_coeff_tag);
            fors_normalise_flat_exit(NULL);
        }


        if (cpl_frameset_count_tags(frameset, slit_location_tag) != 0) {
            cpl_msg_error(recipe, "One (and only one) %s must be in inputs", 
                          slit_location_tag);
            fors_normalise_flat_exit(NULL);
        }
    }

    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 (!dfs_equal_keyword(frameset, "ESO DET CHIP1 ID"))
        cpl_msg_warning(cpl_func,"Input frames are not from the same chip");


    /*
     * Get the reference wavelength and the rebin factor along the
     * dispersion direction from the master flat frame
     */

    if (header == NULL)
        fors_normalise_flat_exit("Cannot load master flat frame header");

    instrume = cpl_propertylist_get_string(header, "INSTRUME");
    if (instrume == NULL)
        fors_normalise_flat_exit("Missing keyword INSTRUME in master "
                                 "flat header");

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

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

    if (cpl_error_get_code() != CPL_ERROR_NONE)
        fors_normalise_flat_exit("Missing keyword ESO INS GRIS1 WLEN "
                                 "in master flat frame header");

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

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

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

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

    if (cpl_error_get_code() != CPL_ERROR_NONE)
        fors_normalise_flat_exit("Missing keyword ESO DET WIN1 BINX "
                                 "in master flat header");

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

    /* Load flat frame */
    cpl_msg_indent_less();
    cpl_msg_info(recipe, "Load input frames...");
    cpl_msg_indent_more();

    cpl_frame * master_flat_frame = cpl_frameset_find(frameset, master_flat_tag);
    master_flat_f = fors_image_load(master_flat_frame);
    if (master_flat_f == NULL)
        fors_normalise_flat_exit("Cannot load master flat field frame");
    mosca::image master_flat(cpl_image_duplicate(master_flat_f->data),
                             cpl_image_power_create(master_flat_f->variance, 0.5),
                             true, mosca::X_AXIS);

    /* Load wavelength calibration */
    cpl_table * disp_coeff = dfs_load_table(frameset, disp_coeff_tag, 1);
    if (disp_coeff == NULL)
        fors_normalise_flat_exit("Cannot load wave calib table");
    mosca::wavelength_calibration wave_cal(disp_coeff, reference);

    cpl_msg_indent_less();
    cpl_msg_info(recipe, "Perform flat field normalisation...");
    cpl_msg_indent_more();

    /* The normaliser object */
    fors::flat_normaliser normaliser;

    if (lss || treat_as_lss) {

        /* FIXME:
         * The LSS data calibration is still dirty: it doesn't apply
         * any spatial rectification, and only in future an external
         * spectral curvature model would be provided in input. Here
         * and there temporary solutions are adpted, such as accepting
         * the preliminary wavelength calibration.
         */

        /*
         * Flat field normalisation is done directly on the master flat
         * field (without spatial rectification first). The spectral
         * curvature model may be provided in input, in future releases.
         */

        normaliser.lss_normalise(master_flat, wave_cal, sradius, dradius, 
                                 spa_polyorder, disp_nknots, splfit_threshold);
    }
    else
    {

        /*
         * This is the generic MOS/MXU handling
         */

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

        polytraces = dfs_load_table(frameset, curv_coeff_tag, 1);
        if (slits == NULL)
            fors_normalise_flat_exit("Cannot load spectral curvature table");

        nx = master_flat.size_dispersion();
        ny = master_flat.size_spatial();

        coordinate = cpl_image_new(nx, ny, CPL_TYPE_FLOAT);
        spatial = mos_spatial_calibration(master_flat.get_cpl_image(),
                slits, polytraces, 
                reference,
                startwavelength, endwavelength,
                dispersion, 0, coordinate);

        cpl_image_delete(spatial); spatial = NULL;

        normaliser.mos_normalise(master_flat, wave_cal, coordinate, slits, polytraces,
                startwavelength, endwavelength,
                dispersion, sradius, dradius,
                spa_polyorder, disp_nknots, splfit_threshold);

        cpl_image_delete(coordinate); coordinate = NULL;
        cpl_table_delete(polytraces); polytraces = NULL;
        cpl_table_delete(slits); slits = NULL;
    }


    /* Save the normalised image */
    fors_image * norm_flat = fors_image_new(
            cpl_image_duplicate(master_flat.get_cpl_image()),
            cpl_image_power_create(master_flat.get_cpl_image_err(), 2.));
    fors_dfs_save_image_err(frameset, norm_flat,
                            master_norm_flat_tag, header, parlist,
                            recipe, master_flat_frame);
    if(cpl_error_get_code() != CPL_ERROR_NONE)
        fors_normalise_flat_exit(NULL);
    cpl_propertylist_delete(header); header = NULL;

    /* Get the spectral shape of the slits and save it */
    cpl_image * wave_profiles = normaliser.get_wave_profiles_im();
    header = cpl_propertylist_new();
    cpl_propertylist_append_string(header, CPL_DFS_PRO_TYPE, "REDUCED");
    cpl_propertylist_append_string(header, CPL_DFS_PRO_CATG, "FLAT_DISP_PROFILE") ;
    cpl_dfs_save_image(frameset, NULL, parlist, frameset, 
                       master_flat_frame, wave_profiles, CPL_BPP_IEEE_FLOAT,
                       recipe, header, NULL, PACKAGE "/" PACKAGE_VERSION,  
                       "flat_disp_profile.fits");
    if(cpl_error_get_code() != CPL_ERROR_NONE)
        fors_normalise_flat_exit(NULL);
    cpl_propertylist_delete(header); header = NULL;
    


    return 0;
}