/* $Id: fors_img_sky_flat_impl.c,v 1.26 2013-09-10 19:16:36 cgarcia Exp $
 *
 * This file is part of the FORS Library
 * 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-10 19:16:36 $
 * $Revision: 1.26 $
 * $Name: not supported by cvs2svn $
 */

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

#include <fors_img_sky_flat_impl.h>

#include "fiera_config.h"
#include <fors_stack.h>
#include <fors_dfs.h>
#include <fors_qc.h>
#include <fors_utils.h>
#include "fors_saturation.h"
#include "fors_detmodel.h"
#include "fors_ccd_config.h"
#include "fors_overscan.h"
#include "fors_subtract_bias.h"

#include <cpl.h>

/**
 * @addtogroup fors_img_sky_flat
 */

/**@{*/

const char *const fors_img_sky_flat_name = "fors_img_sky_flat";
const char *const fors_img_sky_flat_description_short = "Compute master img_sky_flat frame";
const char *const fors_img_sky_flat_author = "Jonas M. Larsen";
const char *const fors_img_sky_flat_email = PACKAGE_BUGREPORT;
const char *const fors_img_sky_flat_description = 
"Input files:\n"
"  DO category:               Type:       Explanation:             Number:\n"
"  SKY_FLAT_IMG               Raw         Jittered sky flat fields    1+\n"
"  MASTER_BIAS                FITS image  Master bias                 1\n"
"\n"
"Output files:\n"
"  DO category:               Data type:  Explanation:\n"
"  MASTER_SKY_FLAT_IMG        FITS image  Master sky flat field\n";
/**
 * @brief    Define recipe parameters
 * @param    parameters     parameter list to fill
 */

void fors_img_sky_flat_define_parameters(cpl_parameterlist *parameters)
{
    char *context = cpl_sprintf("fors.%s", fors_img_sky_flat_name);
    
    fors_stack_define_parameters(parameters, context, "median");

    cpl_free((void *)context);

    return;
}

#undef cleanup
#define cleanup \
do { \
    cpl_frameset_delete(sflat_frames); \
    cpl_frameset_delete(master_bias_frame); \
    fors_stack_method_delete(&sm); \
    cpl_free((void *)context); \
    fors_image_delete(&master_bias); \
    fors_image_delete(&master_sky_flat); \
    fors_image_list_delete(&sflats, fors_image_delete); \
    cpl_propertylist_delete(qc); \
    fors_setting_delete(&setting); \
} while (0)

/**
 * @brief    Do the processing
 *
 * @param    frames         input frames
 * @param    parameters     The parameters list
 *
 * @return   0 if everything is ok
 */
void fors_img_sky_flat(cpl_frameset *frames, const cpl_parameterlist *parameters)
{
    /* Raw */
    cpl_frameset *sflat_frames      = NULL;
    fors_image_list *sflats         = NULL;

    /* Calibration */
    cpl_frameset *master_bias_frame = NULL;
    fors_image *master_bias   = NULL; 

    /* Products */
    fors_image *master_sky_flat = NULL;
    cpl_propertylist *qc = cpl_propertylist_new();
    double saturation;

    /* Parameters */
    stack_method *sm = NULL;

    /* Other */
    char *context = cpl_sprintf("fors.%s", fors_img_sky_flat_name);
    fors_setting *setting = NULL;

    /* QCs about sky level */
    double skylevel, skylevmax, skylevmin;
    int    first_frame;

    /* Get parameters */
    sm = fors_stack_method_new(parameters, context);
    assure( !cpl_error_get_code(), return, "Could not get stacking method");

    /* Eliminate irrelevant RAW frames (header would be inherited from them) */
    cpl_frameset_erase(frames, "TEST");
    
    /* Find raw */
    sflat_frames = fors_frameset_extract(frames, SKY_FLAT_IMG);
    assure( cpl_frameset_get_size(sflat_frames) > 0, return, 
            "No %s provided", SKY_FLAT_IMG);

    /* Find calibration */
    master_bias_frame = fors_frameset_extract(frames, MASTER_BIAS);
    assure( cpl_frameset_get_size(master_bias_frame) == 1, return, 
            "One %s required. %"CPL_SIZE_FORMAT" found", 
            MASTER_BIAS, cpl_frameset_get_size(master_bias_frame));

    /* Get setting */
    setting = fors_setting_new(cpl_frameset_get_position_const(sflat_frames, 0));
    cpl_propertylist * sflat_header = 
            cpl_propertylist_load(cpl_frame_get_filename(
                    cpl_frameset_get_position(sflat_frames, 0)), 0);
    mosca::fiera_config sflat_ccd_config(sflat_header);
    cpl_propertylist_delete(sflat_header);
    assure( !cpl_error_get_code(), return, "Could not get instrument setting" );

    /* Update RON estimation from bias */
    cpl_propertylist * master_bias_header =
       cpl_propertylist_load(cpl_frame_get_filename(
               cpl_frameset_get_position(master_bias_frame, 0)), 0);
    fors_update_ccd_ron(sflat_ccd_config, master_bias_header);
    cpl_propertylist_delete(master_bias_header);
    assure( !cpl_error_get_code(), return, "Could not get RON from master bias"
            " (missing QC DET OUT? RON keywords)");

    master_bias =
            fors_image_load(cpl_frameset_get_position(master_bias_frame, 0));
    assure( !cpl_error_get_code(), return, 
            "Could not load master bias");

    /* Load raw frames  */
    fors_image_list * sflats_raw = 
            fors_image_load_list(sflat_frames);
    assure( !cpl_error_get_code(), return, "Could not load sky flat images");
    
    /* Create variances map */
    fors_image_variance_from_detmodel(sflats_raw, sflat_ccd_config);
    
    /* Subtract overscan */
    sflats = fors_subtract_prescan(sflats_raw, sflat_ccd_config);
    assure( !cpl_error_get_code(), return, "Could not subtract overscan");
    
    /* Trimm pre/overscan */
    fors_trimm_preoverscan(sflats, sflat_ccd_config);
    fors_image_list_delete(&sflats_raw, fors_image_delete);
    assure( !cpl_error_get_code(), return, "Could not trimm overscan");
    
    /* Subtract master bias */
    fors_subtract_bias_imglist(sflats, master_bias);
    saturation = fors_saturation_imglist_satper(sflats);
    assure( !cpl_error_get_code(), return, "Could not subtract master bias");

    /* Normalize to median = 1.
     *
     * (The input sky flats generally have different normalization because of
     *  the time dependent sky level. We must bring the input flats
     *  to the same normalization before stacking.)
     */
    {
        fors_image *image;

        skylevel = skylevmax = skylevmin = 0.0;

        first_frame = 1;
        
        for (image = fors_image_list_first(sflats);
             image != NULL;
             image = fors_image_list_next(sflats)) {

            skylevel = fors_image_get_median(image, NULL);

            fors_image_divide_scalar(image, skylevel, -1.0);
            
            assure( !cpl_error_get_code(), return, 
                    "Raw sky flat normalization failed");

            if (first_frame) {
                first_frame = 0;
                skylevmax = skylevmin = skylevel;
            }
            else {
                if (skylevmax < skylevel)
                    skylevmax = skylevel;
                if (skylevmin > skylevel)
                    skylevmin = skylevel;
            }
        }
    }


    /* Stack */
    master_sky_flat = fors_stack(sflats, sm);
    assure( !cpl_error_get_code(), return, "Sky flat stacking failed");

    /* QC */
    fors_qc_start_group(qc, fors_qc_dic_version, setting->instrument);
    
    fors_qc_write_group_heading(cpl_frameset_get_position(sflat_frames, 0),
                                MASTER_SKY_FLAT_IMG,
                                setting->instrument);
    assure( !cpl_error_get_code(), return, "Could not write %s QC parameters", 
            MASTER_SKY_FLAT_IMG);

    fors_qc_write_qc_double(qc,
                            saturation,
                            "QC.OVEREXPO",
                            "%",
                            "Percentage of overexposed pixels",
                            setting->instrument);

    fors_qc_write_qc_double(qc,
                            skylevmin,
                            "QC.SKYFLAT.FLUX.MIN",
                            "ADU",
                            "Median level of dimmest input flat",
                            setting->instrument);

    fors_qc_write_qc_double(qc,
                            skylevmax,
                            "QC.SKYFLAT.FLUX.MAX",
                            "ADU",
                            "Median level of brightest input flat",
                            setting->instrument);
    fors_qc_end_group();

    /* Save product */
    fors_dfs_save_image_err(frames, master_sky_flat, MASTER_SKY_FLAT_IMG,
                        qc, parameters, fors_img_sky_flat_name, 
                        cpl_frameset_get_position(sflat_frames, 0));
    assure( !cpl_error_get_code(), return, "Saving %s failed",
            MASTER_SKY_FLAT_IMG);
    
    cleanup;
    return;
}


/**@}*/