/* $Id: fors_subtract_sky_lss.c,v 1.6 2013-10-09 15:59:38 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-10-09 15:59:38 $
 * $Revision: 1.6 $
 * $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_subtract_sky_lss_create(cpl_plugin *);
static int fors_subtract_sky_lss_exec(cpl_plugin *);
static int fors_subtract_sky_lss_destroy(cpl_plugin *);
static int fors_subtract_sky_lss(cpl_parameterlist *, cpl_frameset *);

static char fors_subtract_sky_lss_description[] =
"This recipe is used to subtract the sky from wavelength calibrated\n"
"scientific spectra produced by the recipe fors_resample. A simple median\n"
"signal level is subtracted from each image column.\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 may be read STD in case of standard stars observations.\n"
"Note that only LSS or LSS-like MOS/MXU data are to be processed by this\n"
"recipe.\n\n"
"Input files:\n\n"
"  DO category:               Type:       Explanation:         Required:\n"
"  MAPPED_ALL_SCI_MXU         Raw         Scientific exposure     Y\n\n"
"Output files:\n\n"
"  DO category:               Data type:  Explanation:\n"
"  MAPPED_SCI_MXU             FITS image  Rectified scientific spectra\n"
"  MAPPED_SKY_SCI_MXU         FITS image  Rectified sky spectra\n\n";

#define fors_subtract_sky_lss_exit(message)            \
{                                             \
if ((const char *)message != NULL) cpl_msg_error(recipe, message);  \
cpl_image_delete(skymap);                     \
cpl_image_delete(sky);                        \
cpl_image_delete(spectra);                    \
cpl_propertylist_delete(header);              \
cpl_msg_indent_less();                        \
return -1;                                    \
}


#define fors_subtract_sky_lss_exit_memcheck(message)   \
{                                             \
if ((const char *)message != NULL) cpl_msg_info(recipe, message);   \
cpl_image_delete(skymap);                     \
cpl_image_delete(sky);                        \
cpl_image_delete(spectra);                    \
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_subtract_sky_lss",
                    "Subtract sky from calibrated long slit exposure",
                    fors_subtract_sky_lss_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_subtract_sky_lss_create,
                    fors_subtract_sky_lss_exec,
                    fors_subtract_sky_lss_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_subtract_sky_lss_create(cpl_plugin *plugin)
{
    cpl_recipe    *recipe;
/* Uncomment in case parameters are defined
    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(); 

    return 0;
}


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

static int fors_subtract_sky_lss_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_subtract_sky_lss(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_subtract_sky_lss_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_subtract_sky_lss(cpl_parameterlist *parlist, 
                                 cpl_frameset *frameset)
{

    const char *recipe = "fors_subtract_sky_lss";


    /*
     * Input parameters (none)
     */

    /*
     * CPL objects
     */

    cpl_image        *spectra   = NULL;
    cpl_image        *skymap    = NULL;
    cpl_image        *sky       = NULL;
    cpl_table        *maskslits = NULL;

    cpl_propertylist *header    = NULL;

    /*
     * Auxiliary variables
     */

    char        version[80];
    char       *instrume = NULL;
    const char *mapped_science_tag;
    const char *mapped_science_sky_tag;
    const char *mapped_sky_tag;
    int         mxu, mos, lss;
    int         treat_as_lss = 0;
    int         nscience;
    double      mxpos;
    int         nx, ny;
    int         standard;
    float      *data;
    float      *sdata;
    int         i, j;


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

    cpl_msg_set_indentation(2);


    /* 
     * Get configuration parameters (none)
     */

    /* 
     * 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, "MAPPED_ALL_SCI_MXU");
    mos = cpl_frameset_count_tags(frameset, "MAPPED_ALL_SCI_MOS");
    lss = cpl_frameset_count_tags(frameset, "MAPPED_ALL_SCI_LSS");
    standard = 0;

    if (mxu + mos + lss == 0) {
        mxu = cpl_frameset_count_tags(frameset, "MAPPED_ALL_STD_MXU");
        mos = cpl_frameset_count_tags(frameset, "MAPPED_ALL_STD_MOS");
        lss = cpl_frameset_count_tags(frameset, "MAPPED_ALL_STD_LSS");
        standard = 1;
    }

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

    nscience = mxu + mos + lss;

    if (nscience > 1)
        fors_subtract_sky_lss_exit("More than one scientific frame in input"); 

    if (mxu) {
        if (standard) {
            cpl_msg_info(recipe, "MXU data found");
            mapped_science_tag     = "MAPPED_STD_MXU";
            mapped_science_sky_tag = "MAPPED_ALL_STD_MXU";
            mapped_sky_tag         = "MAPPED_SKY_STD_MXU";
        }
        else {
            cpl_msg_info(recipe, "MXU data found");
            mapped_science_tag     = "MAPPED_SCI_MXU";
            mapped_science_sky_tag = "MAPPED_ALL_SCI_MXU";
            mapped_sky_tag         = "MAPPED_SKY_SCI_MXU";
        }
    }

    if (lss) {
        if (standard) {
            cpl_msg_info(recipe, "LSS data found");
            mapped_science_tag      = "MAPPED_STD_LSS";
            mapped_science_sky_tag  = "MAPPED_ALL_STD_LSS";
            mapped_sky_tag          = "MAPPED_SKY_STD_LSS";
        }
        else {
            cpl_msg_info(recipe, "LSS data found");
            mapped_science_tag      = "MAPPED_SCI_LSS";
            mapped_science_sky_tag  = "MAPPED_ALL_SCI_LSS";
            mapped_sky_tag          = "MAPPED_SKY_SCI_LSS";
        }
    }

    if (mos) {
        if (standard) {
            cpl_msg_info(recipe, "MOS data found");
            mapped_science_tag      = "MAPPED_STD_MOS";
            mapped_science_sky_tag  = "MAPPED_ALL_STD_MOS";
            mapped_sky_tag          = "MAPPED_SKY_STD_MOS";
        }
        else {
            cpl_msg_info(recipe, "MOS data found");
            mapped_science_tag      = "MAPPED_SCI_MOS";
            mapped_science_sky_tag  = "MAPPED_ALL_SCI_MOS";
            mapped_sky_tag          = "MAPPED_SKY_SCI_MOS";
        }
    }

    /*
     * Loading input data
     */

    cpl_msg_info(recipe, "Load mapped scientific exposure...");
    cpl_msg_indent_more();

    spectra = dfs_load_image(frameset, mapped_science_sky_tag, 
                             CPL_TYPE_FLOAT, 0, 0);

    if (spectra == NULL)
        fors_subtract_sky_lss_exit("Cannot load input frame");

    header = dfs_load_header(frameset, mapped_science_sky_tag, 0);

    if (header == NULL)
        fors_subtract_sky_lss_exit("Cannot load input frame header");

    instrume = (char *)cpl_propertylist_get_string(header, "INSTRUME");
    if (instrume == NULL)
        fors_subtract_sky_lss_exit("Missing keyword INSTRUME in scientific header");
    instrume = cpl_strdup(instrume);

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

    cpl_free(instrume); instrume = NULL;

    cpl_msg_indent_less();

    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_info(recipe, "All MOS slits have the same offset: %.2f\n"
                         "The LSS data reduction strategy is applied.",
                         mxpos);
        else
            fors_subtract_sky_lss_exit("This recipe can only be used "
                                       "with LSS-like data");
    }

    nx = cpl_image_get_size_x(spectra);
    ny = cpl_image_get_size_y(spectra);

    skymap = cpl_image_new(nx, ny, CPL_TYPE_FLOAT);
    sky    = cpl_image_collapse_median_create(spectra, 0, 0, 1);

    data   = cpl_image_get_data(skymap);

    for (i = 0; i < ny; i++) {
        sdata  = cpl_image_get_data(sky);
        for (j = 0; j < nx; j++) {
            *data++ = *sdata++;
        }
    }

    cpl_image_delete(sky); sky = NULL;
    cpl_image_subtract(spectra, skymap);

    if (dfs_save_image(frameset, skymap, mapped_sky_tag, header,
                       parlist, recipe, version))
        fors_subtract_sky_lss_exit(NULL);

    cpl_image_delete(skymap); skymap = NULL;

    if (dfs_save_image(frameset, spectra, mapped_science_tag, header,
                       parlist, recipe, version))
        fors_subtract_sky_lss_exit(NULL);

    cpl_image_delete(spectra); spectra = NULL;

    cpl_propertylist_delete(header); header = NULL;

    if (cpl_error_get_code()) {
        cpl_msg_error(cpl_error_get_where(), "%s", cpl_error_get_message());
        fors_subtract_sky_lss_exit(NULL);
    }
    else 
        return 0;
}