/* $Id: fors_trace_flat.c,v 1.8 2013-08-20 16:58:25 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-20 16:58:25 $
 * $Revision: 1.8 $
 * $Name: not supported by cvs2svn $
 */

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

#include <math.h>
#include <cpl.h>
#include <moses.h>
#include <fors_dfs.h>

static int fors_trace_flat_create(cpl_plugin *);
static int fors_trace_flat_exec(cpl_plugin *);
static int fors_trace_flat_destroy(cpl_plugin *);
static int fors_trace_flat(cpl_parameterlist *, cpl_frameset *);

static char fors_trace_flat_description[] =
"This recipe is used to trace the edges of MOS/MXU flat field slit spectra\n"
"and determine the spectral curvature solution. The input master flat field\n"
"image, product of the recipe fors_flat, is expected to be oriented with\n"
"horizontal dispersion direction and red wavelengths on the right side.\n"
"The input slits location table should be the product of the recipe\n"
"fors_detect_spectra.\n"
"\n"
"The input master flat image is shifted one pixel down and is subtracted\n"
"from the original image. The result is a vertical gradient map. Next,\n"
"the negative values are forced positive, to obtain an absolute gradient\n"
"map. The map is passed with a horizontal median filter, and after that\n"
"the gradient peaks are traced starting from the slits positions listed\n"
"in the input slits location table. The number of pixels to the left and\n"
"to the right of the reference pixel is trivially derived from the specified\n"
"spectral range and spectral dispersion.\n"
"\n"
"The output spectral curvature table contains the coefficients of the\n"
"polynomial fitting of the found traces, while the output trace table\n"
"contains the traced spectral edges positions in CCD (Y) coordinates for\n"
"each spectrum, and their comparison with their modeling. A spatial map\n"
"is also created, where to each CCD pixel is assigned the value of the\n"
"spatial coordinate along the slit (in pixel). For more details please\n"
"refer to the FORS Pipeline User's Manual.\n"
"\n"
"Note that specifying an input GRISM_TABLE will set some of the recipe\n"
"configuration parameters to default values valid for a particular grism.\n"
"Again, see the pipeline manual for more details.\n"
"\n"
"In the table below the MXU acronym can be alternatively read as MOS.\n\n"
"Input files:\n\n"
"  DO category:               Type:       Explanation:         Required:\n"
"  MASTER_SCREEN_FLAT_MXU     Calib       Master flat frame       Y\n"
"  SLIT_LOCATION_DETECT_MXU   Calib       Slits location          Y\n"
"  GRISM_TABLE                Calib       Grism table             .\n\n"
"Output files:\n\n"
"  DO category:               Data type:  Explanation:\n"
"  CURV_TRACES_MXU            FITS table  Flat field tracings\n"
"  CURV_COEFF_MXU             FITS table  Spectral curvature table\n"
"  SPATIAL_MAP_MXU            FITS image  Map of spatial coordinate\n\n";

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

#define fors_trace_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 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 traces (%p)\n", traces);           \
cpl_table_delete(traces);                       \
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 *list)
{
    cpl_recipe *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_trace_flat",
                    "Determine spectral curvature model",
                    fors_trace_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_trace_flat_create,
                    fors_trace_flat_exec,
                    fors_trace_flat_destroy);

    cpl_pluginlist_append(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_trace_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_trace_flat.dispersion",
                                CPL_TYPE_DOUBLE,
                                "Expected spectral dispersion (Angstrom/pixel)",
                                "fors.fors_trace_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_trace_flat.startwavelength",
                                CPL_TYPE_DOUBLE,
                                "Start wavelength in spectral extraction",
                                "fors.fors_trace_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_trace_flat.endwavelength",
                                CPL_TYPE_DOUBLE,
                                "End wavelength in spectral extraction",
                                "fors.fors_trace_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);

    /*
     * Degree of spectral curvature polynomial
     */

    p = cpl_parameter_new_value("fors.fors_trace_flat.cdegree",
                                CPL_TYPE_INT,
                                "Degree of spectral curvature polynomial",
                                "fors.fors_trace_flat",
                                0);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "cdegree");
    cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
    cpl_parameterlist_append(recipe->parameters, p);

    /*
     * Curvature solution interpolation (for MOS-like data)
     */

    p = cpl_parameter_new_value("fors.fors_trace_flat.cmode",
                                CPL_TYPE_INT,
                                "Interpolation mode of curvature solution "
                                "applicable to MOS-like data (0 = no "
                                "interpolation, 1 = fill gaps, 2 = global "
                                "model)",
                                "fors.fors_trace_flat",
                                1);
    cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "cmode");
    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_trace_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;

    return fors_trace_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_trace_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_trace_flat(cpl_parameterlist *parlist, 
                               cpl_frameset *frameset)
{

    const char *recipe = "fors_trace_flat";


    /*
     * Input parameters
     */

    double      dispersion;
    double      startwavelength;
    double      endwavelength;
    int         cdegree;
    int         cmode;

    /*
     * CPL objects
     */

    cpl_image        *master_flat = NULL;
    cpl_image        *coordinate  = NULL;
    cpl_image        *spatial     = NULL;
    cpl_table        *grism_table = NULL;
    cpl_table        *maskslits   = NULL;
    cpl_table        *slits       = NULL;
    cpl_table        *traces      = NULL;
    cpl_table        *polytraces  = NULL;
    cpl_propertylist *header      = NULL;

    /*
     * Auxiliary variables
     */

    char        version[80];
    const char *master_flat_tag;
    const char *spatial_map_tag;
    const char *slit_detect_tag;
    const char *slit_location_tag;
    const char *curv_traces_tag;
    const char *curv_coeff_tag;
    int         flat_mxu;
    int         flat_mos;
    int         flat_lss;
    int         mos;
    int         nflat;
    int         rebin;
    int         nx, ny;
    int         treat_as_lss;
    double      reference;
    double      mxpos;
    int         nslits_out_det = 0;

    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_trace_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_trace_flat.dispersion", grism_table);

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

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

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

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

    cdegree = dfs_get_parameter_int(parlist,
                    "fors.fors_trace_flat.cdegree", grism_table);

    if (cdegree < 1)
        fors_trace_flat_exit("Invalid polynomial degree");

    if (cdegree > 5)
        fors_trace_flat_exit("Max allowed polynomial degree is 5");

    cmode = dfs_get_parameter_int(parlist, "fors.fors_trace_flat.cmode", NULL);

    if (cmode < 0 || cmode > 2)
        fors_trace_flat_exit("Invalid curvature solution interpolation mode");

    cpl_table_delete(grism_table); grism_table = NULL;

    if (cpl_error_get_code())
        fors_trace_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  = flat_mxu = cpl_frameset_count_tags(frameset, 
                                                "MASTER_SCREEN_FLAT_MXU");
    nflat += flat_mos = cpl_frameset_count_tags(frameset, 
                                                "MASTER_SCREEN_FLAT_MOS");
    nflat += flat_lss = cpl_frameset_count_tags(frameset, 
                                                "MASTER_SCREEN_FLAT_LSS");

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

    mos = 0;

    if (flat_mxu) {
        master_flat_tag   = "MASTER_SCREEN_FLAT_MXU";
        slit_detect_tag   = "SLIT_LOCATION_DETECT_MXU";
        slit_location_tag = "SLIT_LOCATION_MXU";
        curv_traces_tag   = "CURV_TRACES_MXU";
        curv_coeff_tag    = "CURV_COEFF_MXU";
        spatial_map_tag   = "SPATIAL_MAP_MXU";
    }
    else if (flat_mos) {
        mos = 1;
        master_flat_tag   = "MASTER_SCREEN_FLAT_MOS";
        slit_detect_tag   = "SLIT_LOCATION_DETECT_MOS";
        slit_location_tag = "SLIT_LOCATION_MOS";
        curv_traces_tag   = "CURV_TRACES_MOS";
        curv_coeff_tag    = "CURV_COEFF_MOS";
        spatial_map_tag   = "SPATIAL_MAP_MOS";
    }
    else if (flat_lss) {
        fors_trace_flat_exit("LSS spectra are not traceable: use this recipe "
                             "just for MOS/MXU data.");
    }

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

    if (cpl_frameset_count_tags(frameset, slit_detect_tag) > 1) {
        cpl_msg_error(recipe, "Too many in input: %s", slit_detect_tag);
        fors_trace_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
     */

    header = dfs_load_header(frameset, master_flat_tag, 0);

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

    instrume = (char *)cpl_propertylist_get_string(header, "INSTRUME");
    if (instrume == NULL)
        fors_trace_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_trace_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_trace_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_trace_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);
    }


    /*
     * Check if all slits have the same X offset: in such case, abort!
     */

    if (mos)
        maskslits = mos_load_slits_fors_mos(header, &nslits_out_det);
    else
        maskslits = mos_load_slits_fors_mxu(header);

    treat_as_lss = fors_mos_is_lss_like(maskslits, nslits_out_det);

    if (treat_as_lss) {
        cpl_msg_error(recipe, "All slits have the same offset: %.2f mm\n"
                      "Spectra are not traceable: the LSS data reduction\n"
                      "strategy must be applied.", mxpos);
        fors_trace_flat_exit(NULL);
    }

    cpl_table_delete(maskslits); maskslits = NULL;


    cpl_msg_indent_less();
    cpl_msg_info(recipe, "Load input frames...");
    cpl_msg_indent_more();

    master_flat = dfs_load_image(frameset, master_flat_tag, 
                                 CPL_TYPE_FLOAT, 0, 0);
    if (master_flat == NULL)
        fors_trace_flat_exit("Cannot load master flat field frame");

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


    cpl_msg_indent_less();
    cpl_msg_info(recipe, "Determining spectral curvature...");
    cpl_msg_indent_more();

    cpl_msg_info(recipe, "Tracing master flat field spectra edges...");
    traces = mos_trace_flat(master_flat, slits, reference,
                            startwavelength, endwavelength, dispersion);

    if (!traces)
        fors_trace_flat_exit("Tracing failure");

    cpl_msg_info(recipe, "Fitting flat field spectra edges...");
    polytraces = mos_poly_trace(slits, traces, cdegree);

    if (!polytraces)
        fors_trace_flat_exit("Trace fitting failure");

    if (cmode) {
        cpl_msg_info(recipe, "Computing global spectral curvature model...");
        mos_global_trace(slits, polytraces, cmode);
    }

    if (dfs_save_table(frameset, traces, curv_traces_tag, NULL, parlist,
                       recipe, version))
        fors_trace_flat_exit(NULL);

    cpl_table_delete(traces); traces = NULL;

    nx = cpl_image_get_size_x(master_flat);
    ny = cpl_image_get_size_y(master_flat);
    coordinate = cpl_image_new(nx, ny, CPL_TYPE_FLOAT);

    spatial = mos_spatial_calibration(master_flat, slits, polytraces, 
                                      reference,
                                      startwavelength, endwavelength,
                                      dispersion, 0, coordinate);

    cpl_image_delete(master_flat); master_flat = NULL;
    cpl_image_delete(spatial); spatial = NULL;

    if (dfs_save_image(frameset, coordinate, spatial_map_tag, header,
                       parlist, recipe, version))
        fors_trace_flat_exit(NULL);

    cpl_image_delete(coordinate); coordinate = NULL;
    cpl_propertylist_delete(header); header = NULL;

    if (dfs_save_table(frameset, slits, slit_location_tag, NULL,
                       parlist, recipe, version))
        fors_trace_flat_exit(NULL);

    cpl_table_delete(slits); slits = NULL;

    if (dfs_save_table(frameset, polytraces, curv_coeff_tag, NULL,
                       parlist, recipe, version))
        fors_trace_flat_exit(NULL);

    cpl_table_delete(polytraces); polytraces = NULL;

    return 0;
}