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

#include "vimos_spec_idp.h"
#include "vimos_science_impl.h"
#include "vimos_dfs.h"
#include "irplib_sdp_spectrum.h"
#include <moses.h>
#include "math.h"
#include <algorithm>
#include <vector>

#define PI 3.14159265358979323846

/*----------------- Private Structs Declarations -------------------------*/

const double angstrom_to_nm = 0.1;

typedef struct{

	cpl_parameterlist * parlist;
	cpl_frameset * frames;
	cpl_frame * science_frame;

	cpl_image * science_img;
	cpl_image * science_img_err;

	cpl_image * science_fluxcal_img;
	cpl_image * science_fluxcal_img_err;

	cpl_image * sky_img;

	cpl_array * wlens_angstrom;

	const char * recipe;
	const char * pipe_id;

	double spec_sye;
	double aperture;
}idp_save_info;

const char * n_combine_tag = "ESO PRO DATANCOM";

typedef std::pair<double, double> ra_dec;


/*---------------------- Forward Declarations ------------------------------*/

static
double get_property_or_nan(const cpl_propertylist * l, const char * name);

static
cpl_error_code save_spectra(idp_save_info& info, const cpl_size mplex,
		const std::vector<ra_dec>& ra_decs_for_extension);

static
bool get_ra_dec_for_extension(cpl_frameset * frames,
                const cpl_size multiplex_idx, std::vector<ra_dec>& data_out);

static
double extract_aperture_for_extension(const cpl_frame * science,
		const cpl_size index);

static
std::vector<double> extract_spec_sye(const cpl_frameset * fset,
		const cpl_size num_final_extensions);

static
cpl_array * get_wlens_angstrom(const cpl_frame * frame, const cpl_size num_wav);

static
bool check_posang_from_science(const cpl_frameset * fset);

/*---------------------------- Public Functions ------------------------------*/

/*----------------------------------------------------------------------------*/
/**
 * @brief
 *   Save the science products as IDP-compliant 1D spectra
 *
 * @param parlist  Recipe parameter list
 * @param frames   Recipe frameset containing both input and output frames
 *
 * @return Error code
 *
 * This function must be called after the main products have been generated. It
 * takes care of loading the products and generates the IDP-compliant 1D spectra.
 *
 */
/*----------------------------------------------------------------------------*/
cpl_error_code vimos_spec_idp_save(cpl_parameterlist * parlist, cpl_frameset * frames){

	const cpl_error_code start_code = cpl_error_get_code();
	if(start_code)
		return start_code;

    const int generate = dfs_get_parameter_bool(parlist,
                             "vimos.vmmosscience.generate_idp", NULL);

    if(!generate) return CPL_ERROR_NONE;

    if(!check_posang_from_science(frames))
        return CPL_ERROR_NONE;


    const char * science_frame_tag = "MOS_SCIENCE_REDUCED";
    const char * science_err_frame_tag = "MOS_SCI_ERROR_REDUCED";

	const char * science_fluxcal_frame_tag = "MOS_SCIENCE_FLUX_REDUCED";
    const char * science_fluxcal_err_frame_tag = "MOS_SCI_ERROR_FLUX_REDUCED";

    const char *sky_frame_tag = "MOS_SCI_SKY_REDUCED";

	cpl_frame * science_frame = cpl_frameset_find(frames, science_frame_tag);
	cpl_frame * science_err_frame = cpl_frameset_find(frames, science_err_frame_tag);

	cpl_frame * science_frame_fluxcal = cpl_frameset_find(frames, science_fluxcal_frame_tag);
	cpl_frame * science_err_frame_fluxcal = cpl_frameset_find(frames, science_fluxcal_err_frame_tag);

	cpl_frame * sky_frame = cpl_frameset_find(frames, sky_frame_tag);

	if(science_frame == NULL){
		cpl_msg_error(cpl_func, "%s not available, aborting IDP generation", science_frame_tag);
		return CPL_ERROR_NONE;
	}

	if(science_err_frame == NULL){
		cpl_msg_error(cpl_func, "%s not available, aborting IDP generation", science_err_frame_tag);
		return CPL_ERROR_NONE;
	}

	if(science_frame_fluxcal == NULL){
		cpl_msg_error(cpl_func, "%s not available, aborting IDP generation", science_fluxcal_frame_tag);
		return CPL_ERROR_NONE;
	}

	if(science_err_frame_fluxcal == NULL){
		cpl_msg_error(cpl_func, "%s not available, aborting IDP generation", science_fluxcal_err_frame_tag);
		return CPL_ERROR_NONE;
	}

	if(sky_frame == NULL){
		cpl_msg_error(cpl_func, "%s not available, aborting IDP generation", sky_frame_tag);
		return CPL_ERROR_NONE;
	}

	const char *recipe = "vmmosscience";
    std::string version = std::string(PACKAGE) + "-" + PACKAGE_VERSION;

	const cpl_size mplex_sz = cpl_frame_get_nextensions(science_frame);
	const cpl_size mplex_all_valid_sz = std::max(cpl_size(1), mplex_sz);
	std::vector<double> spec_syes = extract_spec_sye(frames, mplex_all_valid_sz);

	if(spec_syes.size() != (size_t)mplex_all_valid_sz){
		cpl_msg_error(cpl_func, "Unable to extract the correct number of SPEC_SYE values aborting");
		return CPL_ERROR_ILLEGAL_INPUT;
	}

	for(cpl_size mplex_idx = std::min(cpl_size(1), mplex_sz); mplex_idx <= mplex_sz; ++mplex_idx){

		std::vector<ra_dec> ra_dec_for_extension;
		if(!get_ra_dec_for_extension(frames, std::max(cpl_size(1), mplex_idx), ra_dec_for_extension)){
			cpl_msg_error(cpl_func, "Unable to extract RA/DEC for extension %lld, aborting", mplex_idx);
			return CPL_ERROR_ILLEGAL_INPUT;
		}

		cpl_image * science = cpl_image_load(cpl_frame_get_filename(science_frame), CPL_TYPE_FLOAT,
				0, mplex_idx);

		if((size_t)cpl_image_get_size_y(science) != ra_dec_for_extension.size()){
			cpl_msg_error(cpl_func, "Number of detected objects differ from rows of reduced image"
					" for extension %lld, aborting", mplex_idx);
			cpl_image_delete(science);
			return CPL_ERROR_ILLEGAL_INPUT;
		}

		double aperture = extract_aperture_for_extension(science_frame, mplex_idx);
		if(std::isnan(aperture)){
			cpl_msg_error(cpl_func, "Extraction of aperture failed for %lld extension, aborting", mplex_idx);
			cpl_image_delete(science);
			return CPL_ERROR_ILLEGAL_INPUT;
		}

		cpl_image * science_err = cpl_image_load(cpl_frame_get_filename(science_err_frame),
				CPL_TYPE_FLOAT, 0, mplex_idx);
		cpl_image * sky = cpl_image_load(cpl_frame_get_filename(sky_frame), CPL_TYPE_FLOAT, 0, mplex_idx);
		cpl_image * science_fluxcal = cpl_image_load(cpl_frame_get_filename(science_frame_fluxcal),
				CPL_TYPE_FLOAT, 0, mplex_idx);
		cpl_image * science_fluxcal_err = cpl_image_load(cpl_frame_get_filename(science_err_frame_fluxcal),
				CPL_TYPE_FLOAT, 0, mplex_idx);

		cpl_array * wlens_angstrom = get_wlens_angstrom(science_frame, cpl_image_get_size_x(science));

		cpl_error_code spectra_err_code = wlens_angstrom == NULL ?  CPL_ERROR_ILLEGAL_OUTPUT : CPL_ERROR_NONE;

		if(spectra_err_code)
					cpl_msg_error(cpl_func, "Unable to extract wavelengths for "
							"extension %lld, aborting", mplex_idx);

		if(!spectra_err_code){

			idp_save_info info{parlist, frames, science_frame, science, science_err,
				science_fluxcal, science_fluxcal_err, sky, wlens_angstrom, recipe,
				version.c_str(), spec_syes[mplex_idx - 1], aperture};

			spectra_err_code = save_spectra(info, mplex_idx, ra_dec_for_extension);
		}

		cpl_image_delete(science);
		cpl_image_delete(science_err);

		cpl_image_delete(science_fluxcal);
		cpl_image_delete(science_fluxcal_err);

		cpl_image_delete(sky);

		cpl_array_delete(wlens_angstrom);

		if(spectra_err_code) {
			cpl_msg_error(cpl_func, "Unable to save extension %lld, aborting", mplex_idx);
			cpl_error_reset();
			return spectra_err_code;
		}
	}

	return cpl_error_get_code();
}


/*----------------------------------------------------------------------------*/
/**
 * @brief
 *   Calculate the shift from skyalign
 *
 * @param offsets  Table with measured sky lines offsets (in pixel)
 * @param slits	   Slits positions table
 *
 * @return shift from skyalign, NAN in case of failure
 *
 * The function calculates the mean shift with respect to the sky lines. See
 * See e.g. mos_wavelength_align for more information
 *
 */
/*----------------------------------------------------------------------------*/
double vimos_spec_idp_get_sky_align_shift(const cpl_table * offsets,
		const cpl_table * slits){

	if(offsets == NULL || slits == NULL) return NAN;

    const cpl_size nslits = cpl_table_get_nrow(slits);
    const cpl_size nlines = cpl_table_get_nrow(offsets);
    double off = 0;
    int num_offs = 0;
	for (cpl_size i = 0; i < nslits; i++) {
		int rej = 0;
		if (cpl_table_get_int(slits, "length", i, &rej) == 0 || rej == 1)
			continue;

		std::stringstream ss;
		ss<<"offset"<<cpl_table_get_int(slits, "slit_id", i, &rej);

		if(rej) continue;

		for(cpl_size j = 0; j < nlines; ++j){
			const double this_off = cpl_table_get_double(offsets, ss.str().c_str(), j, &rej);
			if(rej) continue;
			off += this_off;
			num_offs ++;
		}
	}

	if(cpl_error_get_code() || num_offs == 0)
	{
		cpl_error_reset();
		return NAN;
	}

    off = std::abs(off);
	return off / (double)num_offs;
}

/*----------------------------------------------------------------------------*/
/**
 * @brief
 *   Get the correct tag for the GRISM name
 *
 * @param science_header Propertylist containing informations on the quadrant used
 *
 * @return tag to use, e.g. ESO INSO GRISx NAME where x is the quadrant if 1<=x<=4.
 * NULL in case of error
 *
 */
/*----------------------------------------------------------------------------*/
const char * get_grism_name_tag(const cpl_propertylist * science_header){

	int quadrant = cpl_propertylist_get_int(science_header, "ESO OCS CON QUAD");

	if(cpl_error_get_code()){
		cpl_error_reset();
		return NULL;
	}

    switch (quadrant) {
    case 1:
        return "ESO INS GRIS1 NAME";
    case 2:
        return "ESO INS GRIS2 NAME";
    case 3:
        return "ESO INS GRIS3 NAME";
    case 4:
        return "ESO INS GRIS4 NAME";
    default:
    	return NULL;
    }
}

/*---------------------------- Private Functions ------------------------------*/

static const char * wave_units = "angstrom";
static const char * flx_units = "10**(-16)erg.cm**(-2).s**(-1).angstrom**(-1)";
static const char * flx_units_uncalib = "adu.s**(-1)";

static const char* type_qual = "spec:Data.FluxAxis.Accuracy.QualityStatus";
static const char* type_wave = "spec:Data.SpectralAxis.Value";
static const char* type_flx = "spec:Data.FluxAxis.Value";
static const char* type_flx_err = "spec:Data.FluxAxis.Accuracy.StatError";
static const char* type_flx_uncalib = "eso:Data.FluxAxis.Value";
static const char* type_flx_err_uncalib = "eso:Data.FluxAxis.Accuracy.StatError";

static const char * ucd_qual = "meta.code.qual;meta.main";
static const char * ucd_wave = "em.wl;obs.atmos";
static const char * ucd_flux = "phot.flux.density;em.wl;src.net;meta.main";
static const char * ucd_flux_err = "stat.error;phot.flux.density;meta.main";
static const char * ucd_flux_uncalib = "phot.flux.density;em.wl;src.net;stat.uncalib";
static const char * ucd_flux_err_uncalib = "stat.error;phot.flux.density;stat.uncalib";

static const char * wave_clmn_name = "WAVE";
static const char * flux_clmn_name = "FLUX";
static const char * flux_clmn_err_name = "ERR";

static cpl_array *
extract_flux(const cpl_image * img, const cpl_size idx){

	const cpl_size sz_x = cpl_image_get_size_x(img);
	cpl_array * to_ret = cpl_array_new(sz_x, CPL_TYPE_DOUBLE);
	for(cpl_size x = 0; x < sz_x; ++x){
		int rej = 0;
		const double val = cpl_image_get(img, x + 1, idx + 1, &rej);
		if(!rej)
			cpl_array_set(to_ret, x, val);
		else
			cpl_array_set_invalid(to_ret, x);
	}

	return to_ret;
}

static cpl_array *
make_qual(const cpl_array * flx_data) {
    const cpl_size flux_sz = cpl_array_get_size(flx_data);
    cpl_array * qual_data = cpl_array_new(flux_sz, CPL_TYPE_INT);
    for (cpl_size x = 0; x < flux_sz; ++x) {
        int reject = 0;
        const double val = cpl_array_get(flx_data, x, &reject);
        if (reject)  // 1 datum is invalid, -1 there was an error
            cpl_array_set_invalid(qual_data, x);
        else
            cpl_array_set(qual_data, x, val == -1.0 ? 1 : 0);
    }
    return qual_data;
}

static
irplib_sdp_spectrum * extract_1D_spectrum(idp_save_info& info, const cpl_size idx){

	irplib_sdp_spectrum * spectrum = irplib_sdp_spectrum_new();
	const cpl_size sz = cpl_image_get_size_x(info.science_img);

	cpl_error_code fail = irplib_sdp_spectrum_set_nelem(spectrum, sz);
	if(fail){
		irplib_sdp_spectrum_delete(spectrum);
		return NULL;
	}

	cpl_array * flx_data = extract_flux(info.science_fluxcal_img, idx);
	cpl_array * flx_err_data = extract_flux(info.science_fluxcal_img_err, idx);

	cpl_array * bkg_data = extract_flux(info.sky_img, idx);

	cpl_array * flx_red_data = extract_flux(info.science_img, idx);
	cpl_array * flx_red_data_err = extract_flux(info.science_img_err, idx);

    cpl_array * qual_data = make_qual(flx_data);

	fail = irplib_sdp_spectrum_add_column(spectrum, wave_clmn_name, CPL_TYPE_DOUBLE, wave_units,
			NULL, type_wave, ucd_wave, info.wlens_angstrom);

    if(!fail)
    {
        fail = irplib_sdp_spectrum_replace_column_comment(spectrum, wave_clmn_name,
                      "TUCD", "Air wavelength");
    }

	if(fail){
		irplib_sdp_spectrum_delete(spectrum);
		return NULL;
	}

	fail = irplib_sdp_spectrum_add_column(spectrum, flux_clmn_name, CPL_TYPE_DOUBLE, flx_units,
			NULL, type_flx, ucd_flux, flx_data);

	if(fail){
		irplib_sdp_spectrum_delete(spectrum);
		return NULL;
	}

	fail = irplib_sdp_spectrum_add_column(spectrum, flux_clmn_err_name, CPL_TYPE_DOUBLE, flx_units,
				NULL, type_flx_err, ucd_flux_err, flx_err_data);

	if(fail){
		irplib_sdp_spectrum_delete(spectrum);
		return NULL;
	}

	fail = irplib_sdp_spectrum_add_column(spectrum, "BGFLUX", CPL_TYPE_DOUBLE, flx_units_uncalib,
			NULL, type_flx_uncalib, ucd_flux_uncalib, bkg_data);

	if(fail){
		irplib_sdp_spectrum_delete(spectrum);
		return NULL;
	}

	fail = irplib_sdp_spectrum_add_column(spectrum, "FLUX_REDUCED", CPL_TYPE_DOUBLE, flx_units_uncalib,
			NULL, type_flx_uncalib, ucd_flux_uncalib, flx_red_data);

	if(fail){
		irplib_sdp_spectrum_delete(spectrum);
		return NULL;
	}

	fail = irplib_sdp_spectrum_add_column(spectrum, "ERR_REDUCED", CPL_TYPE_DOUBLE, flx_units_uncalib,
			NULL, type_flx_err_uncalib, ucd_flux_err_uncalib, flx_red_data_err);

	if(fail){
		irplib_sdp_spectrum_delete(spectrum);
		return NULL;
	}

	fail = irplib_sdp_spectrum_add_column(spectrum, "QUAL", CPL_TYPE_INT, "",
			NULL, type_qual, ucd_qual, qual_data);

	if(fail){
		irplib_sdp_spectrum_delete(spectrum);
		return NULL;
	}

	cpl_array_delete(flx_data);
	cpl_array_delete(flx_err_data);

	cpl_array_delete(flx_red_data);
	cpl_array_delete(flx_red_data_err);

	cpl_array_delete(bkg_data);

	cpl_array_delete(qual_data);

	return spectrum;
}

static
double get_int_or_nan(const cpl_propertylist * p, const char * name){

	double r = cpl_propertylist_get_int(p, name);

	if(cpl_error_get_code()){
		cpl_error_reset();
		return NAN;
	}

	return r;
}

static
double get_double_or_nan(const cpl_propertylist * p, const char * name){

	double r = cpl_propertylist_get_double(p, name);

	if(cpl_error_get_code()){
		cpl_error_reset();
		return NAN;
	}

	return r;
}

static
bool add_exposures(irplib_sdp_spectrum * dest, cpl_frameset * allframes, const cpl_propertylist * source){

	const char * total_exptime_tag = "ESO PRO EXPTTOT";
    const char * mjd = "MJD-OBS";
    const char * exptime_s = "EXPTIME";

	const double mjd_obs = get_double_or_nan(source, mjd);
	if(std::isnan(mjd_obs)) return false;

	bool success = irplib_sdp_spectrum_copy_exptime(dest, source, total_exptime_tag) == CPL_ERROR_NONE;
	success = success && irplib_sdp_spectrum_copy_texptime(dest, source, total_exptime_tag) == CPL_ERROR_NONE;

	const double sec_in_day = 86400.0;
	cpl_frameset_iterator * f_iter = cpl_frameset_iterator_new(allframes);
	const cpl_frame * frame = cpl_frameset_iterator_get_const(f_iter);
    double max_mjd_obs = 0;
    double exptime_last = 0;;
	while (frame != NULL)
	{
	    if(std::string(cpl_frame_get_tag(frame)) == "MOS_SCIENCE")
	    {
	        cpl_propertylist * raw_science_plist = cpl_propertylist_load(cpl_frame_get_filename(frame), 0);
	        double this_mjd_obs = get_double_or_nan(raw_science_plist, mjd);
	        if(std::isnan(this_mjd_obs)) return false;
	        if(this_mjd_obs > max_mjd_obs)
	        {
	            max_mjd_obs = this_mjd_obs;
	            exptime_last = get_double_or_nan(raw_science_plist, exptime_s);
	        }
	    }
        cpl_frameset_iterator_advance(f_iter, 1);
        frame = cpl_frameset_iterator_get_const(f_iter);
	}
	cpl_frameset_iterator_delete(f_iter);
    const double mjd_end = max_mjd_obs + exptime_last / sec_in_day;
	success = success && irplib_sdp_spectrum_set_mjdobs(dest, mjd_obs) == CPL_ERROR_NONE;
	success = success && irplib_sdp_spectrum_set_mjdend(dest, mjd_end) == CPL_ERROR_NONE;
	success = success && irplib_sdp_spectrum_set_telapse(dest, (mjd_end - mjd_obs) * sec_in_day) == CPL_ERROR_NONE;
	success = success && irplib_sdp_spectrum_set_tmid(dest, (mjd_end + mjd_obs) / 2.0) == CPL_ERROR_NONE;

	return success;
}

static bool
add_pro(irplib_sdp_spectrum * spec, const cpl_propertylist * plist){

	const char * pro_soft_tag = "ESO PRO REC1 PIPE ID";

	bool success = irplib_sdp_spectrum_copy_procsoft(spec, plist, pro_soft_tag) == CPL_ERROR_NONE;
	success = success && irplib_sdp_spectrum_set_referenc(spec, "") == CPL_ERROR_NONE;
	success = success && irplib_sdp_spectrum_set_prodcatg(spec, "SCIENCE.SPECTRUM") == CPL_ERROR_NONE;

    int isci = 1;
    while(1)
    {
        std::stringstream ss;
        ss<<"ESO PRO REC1 RAW"<<isci<<" NAME";
        if(cpl_propertylist_has(plist, ss.str().c_str()))
        {
            success = success && irplib_sdp_spectrum_copy_prov(spec, isci, plist, ss.str().c_str()) == CPL_ERROR_NONE;
            isci++;
        }
        else
            break;
    }

	return success;
}

static bool
add_obs(irplib_sdp_spectrum * spec, const cpl_propertylist * plist){

	const char * prog_id_tag = "ESO OBS PROG ID";
	const char * ob_id_tag = "ESO OBS ID";
	const char * obs_tech = "ESO PRO TECH";

	bool success = irplib_sdp_spectrum_copy_progid(spec, plist, prog_id_tag) == CPL_ERROR_NONE;
	success = success && irplib_sdp_spectrum_copy_obid(spec, 1, plist, ob_id_tag) == CPL_ERROR_NONE;
	success = success && irplib_sdp_spectrum_copy_obstech(spec, plist, obs_tech) == CPL_ERROR_NONE;
	success = success && irplib_sdp_spectrum_set_fluxcal(spec, "ABSOLUTE/UNCALIBRATED") == CPL_ERROR_NONE;
	return success;
}

static bool valid_flx(const cpl_array * wave, const cpl_array * flux,
		const cpl_size i, double * wave_out_a){

	int rej_w = 0;
	const double w = cpl_array_get(wave, i, &rej_w);

	if(rej_w) return false;

	int rej_f = 0;
	const double f = cpl_array_get(flux, i, &rej_f);

	if(f == -1.0 || rej_f) return false;

	*wave_out_a = w;

	return true;
}

static
const char * get_filt_name_tag(const cpl_propertylist * plist){

	int quadrant = cpl_propertylist_get_int(plist, "ESO OCS CON QUAD");

	if(cpl_error_get_code()){
		cpl_error_reset();
		return NULL;
	}

    switch (quadrant) {
    case 1:
        return "ESO INS FILT1 NAME";
    case 2:
        return "ESO INS FILT2 NAME";
    case 3:
        return "ESO INS FILT3 NAME";
    case 4:
        return "ESO INS FILT4 NAME";
    default:
    	return NULL;
    }
}


static bool add_specbin(irplib_sdp_spectrum * dest, const cpl_propertylist * plist){

	const double cd11 = get_double_or_nan(plist, "CD1_1") * angstrom_to_nm;
	if(std::isnan(cd11)) return false;

	return !irplib_sdp_spectrum_set_specbin(dest, cd11);
}

static bool
add_flux_data(irplib_sdp_spectrum * dest, const cpl_propertylist * plist){

	const char * grism_name_tag = get_grism_name_tag(plist);

	if(!grism_name_tag){
		cpl_msg_warning(cpl_func, "GRISM name extraction failed, aborting IDP generation");
		return false;
	}

	bool success = !irplib_sdp_spectrum_set_specsys(dest, "TOPOCENT");
	success = success && !irplib_sdp_spectrum_set_contnorm(dest, CPL_FALSE);
	success = success && !irplib_sdp_spectrum_set_totflux(dest, CPL_FALSE);
	success = success && !irplib_sdp_spectrum_set_fluxcal(dest, "ABSOLUTE");
	success = success && !irplib_sdp_spectrum_set_fluxerr(dest, -2.0);
	success = success && !irplib_sdp_spectrum_copy_dispelem(dest, plist, grism_name_tag);
	success = success && add_specbin(dest, plist);
	success = success && !irplib_sdp_spectrum_set_voclass(dest, "SPECTRUM V2.0");
	success = success && !irplib_sdp_spectrum_set_vopub(dest, "ESO/SAF");
	success = success && !irplib_sdp_spectrum_set_vopub(dest, "ESO/SAF");
	success = success && !irplib_sdp_spectrum_copy_object(dest, plist, "OBJECT");
	success = success && !irplib_sdp_spectrum_copy_title(dest, plist, "OBJECT");

	return success;
}

static double calc_snr(const irplib_sdp_spectrum * spectrum){
    const cpl_array * flx = irplib_sdp_spectrum_get_column_data(spectrum, flux_clmn_name);
    const cpl_array * err = irplib_sdp_spectrum_get_column_data(spectrum, flux_clmn_err_name);

    cpl_array * div = cpl_array_duplicate(flx);
    cpl_array_divide(div, err);
    int valid;
    for(int i=0; i<cpl_array_get_size(flx); ++i)
        if(cpl_array_get_double(flx, i, &valid)==-1 &&
           cpl_array_get_double(err, i, &valid)==-1)
            cpl_array_set_invalid(div, i);

    const double ret = cpl_array_get_median(div);
    cpl_array_delete(div);

    if(cpl_error_get_code()) 
    {
        cpl_error_reset();
        return NAN;
    }

    return ret;
}

static bool add_snr(irplib_sdp_spectrum * spectrum){

	const double snr = calc_snr(spectrum);
	if(std::isnan(snr)) return false;
	return irplib_sdp_spectrum_set_snr(spectrum, snr) == CPL_ERROR_NONE;
}

static bool
add_wave_data(irplib_sdp_spectrum * dest){

	const cpl_array * wave = irplib_sdp_spectrum_get_column_data(dest, wave_clmn_name);
	const cpl_array * flux = irplib_sdp_spectrum_get_column_data(dest, flux_clmn_name);


	double wave_min_valid_a = NAN;
	double wave_max_valid_a = NAN;

	for(cpl_size i = 0; i < cpl_array_get_size(wave); ++i){
		if(valid_flx(wave, flux, i, &wave_min_valid_a)){
			break;
		}
	}

	for(cpl_size i = cpl_array_get_size(wave) - 1; i >= 0; --i){
		if(valid_flx(wave, flux, i, &wave_max_valid_a)){
			break;
		}
	}

    int null_val_min = 0, null_val_max = 0;
    const double tdmin = cpl_array_get_double(wave, 0, &null_val_min);
    const double tdmax = cpl_array_get_double(wave,
                                              cpl_array_get_size(wave)-1, &null_val_max);

	if(std::isnan(wave_min_valid_a)) return false;
	if(std::isnan(wave_max_valid_a)) return false;
	if(wave_min_valid_a >= wave_max_valid_a) return false;
    if(null_val_min || null_val_max) return false;

	const double wave_max_valid_nm = wave_max_valid_a * angstrom_to_nm;
	const double wave_min_valid_nm = wave_min_valid_a * angstrom_to_nm;

	const double spec_val = (wave_max_valid_nm + wave_min_valid_nm) / 2.0;
	const double spec_bw = wave_max_valid_nm - wave_min_valid_nm;

	bool success = !irplib_sdp_spectrum_set_wavelmin(dest, wave_min_valid_nm);
	success = success && !irplib_sdp_spectrum_set_wavelmax(dest, wave_max_valid_nm);

	success = success && !irplib_sdp_spectrum_set_specbw(dest, spec_bw);
	success = success && !irplib_sdp_spectrum_set_nelem(dest, cpl_array_get_size(wave));
	success = success && !irplib_sdp_spectrum_set_specval(dest, spec_val);

	success = success && !irplib_sdp_spectrum_set_tdmin(dest, tdmin);
	success = success && !irplib_sdp_spectrum_set_tdmax(dest, tdmax);

	return success;
}


static std::vector<ra_dec> get_ra_dec(const cpl_table * objects)
{
    cpl_size maxobjects = 1;
    std::stringstream name_stream;
    name_stream<<"object_"<<maxobjects;

    /*
     * Count objects to extract
     */

    while (cpl_table_has_column(objects, name_stream.str().c_str())) {
        maxobjects++;
        name_stream.str(std::string());
        name_stream<<"object_"<<maxobjects;
    }


    cpl_size nslits = cpl_table_get_nrow(objects);
    std::vector<ra_dec> to_ret;

    for (cpl_size i = 0; i < nslits; i++) {
        for (cpl_size j = 1; j < maxobjects; j++) {
        	std::stringstream obj_name_stream;
			obj_name_stream<<"object_"<<j;
            if (cpl_table_is_valid(objects, obj_name_stream.str().c_str(), i)){

            	std::stringstream ra_stream; ra_stream << "ra_"<<j;
            	std::stringstream dec_stream; dec_stream << "dec_"<<j;

            	std::string ra_clmn_name = ra_stream.str();
            	std::string dec_clmn_name = dec_stream.str();

            	bool missing_clmn = !cpl_table_has_column(objects, ra_clmn_name.c_str()) ||
            			!cpl_table_has_column(objects, dec_clmn_name.c_str());

            	if(missing_clmn){
            		cpl_msg_warning(cpl_func, "A Ra/Dec column is missing, aborting IDP generation");
            		return {};
            	}

            	const bool is_ra_valid = cpl_table_is_valid(objects, ra_clmn_name.c_str(), i);
            	const bool is_dec_valid = cpl_table_is_valid(objects, dec_clmn_name.c_str(), i);

            	if(is_ra_valid && is_dec_valid){
            		const double ra = cpl_table_get_double(objects, ra_stream.str().c_str(), i, NULL);
            		const double dec = cpl_table_get_double(objects, dec_stream.str().c_str(), i, NULL);
            		to_ret.emplace_back(ra, dec);
            	}
            	else
            	{
            		cpl_msg_warning(cpl_func, "A Ra/Dec column is invalid while object is valid, "
            				"aborting IDP generation");
            		return {};
            	}
            }
        }
    }

    if(cpl_error_get_code()){
    	cpl_error_reset();
    	to_ret.clear();
    }
    return to_ret;
}

static bool get_ra_dec_for_extension(cpl_frameset * frames, const cpl_size multiplex_idx, std::vector<ra_dec>& data_out){

    const char *object_table_tag = "OBJECT_SCI_TABLE";

    const cpl_frame * objects_table_frame = cpl_frameset_find_const(frames, object_table_tag);

    if(objects_table_frame == NULL){
    	cpl_msg_error(cpl_func, "%s unavailable, aborting IDP generation", object_table_tag);
    	return false;
    }

    cpl_table * obs = cpl_table_load(cpl_frame_get_filename(objects_table_frame), multiplex_idx, 1);

    if(obs == NULL){
    	cpl_msg_error(cpl_func, "Cannot load %s, aborting IDP generation", object_table_tag);
    	return false;
    }

    data_out = get_ra_dec(obs);

    cpl_table_delete(obs);

    return !data_out.empty();
}

static double get_stacked_effron(const double detron,
        const cpl_propertylist * science_plist){

    const double nstacked = get_int_or_nan(science_plist, n_combine_tag);
    return detron / sqrt(nstacked);
}

static bool add_ron_from_bias(irplib_sdp_spectrum * dest, idp_save_info& info,
        const cpl_propertylist * science_plist){

    const char * qc_ron_tag = "ESO QC RON";

	cpl_frame * bias = cpl_frameset_find(info.frames, "MASTER_BIAS");
	if(bias == NULL) return false;


	cpl_propertylist * l = cpl_propertylist_load(cpl_frame_get_filename(bias), 0);
	bool success = irplib_sdp_spectrum_copy_detron(dest, l, qc_ron_tag) == CPL_ERROR_NONE;
	if(success){
	    const double effron = get_stacked_effron(irplib_sdp_spectrum_get_detron(dest), science_plist);
	    success = success && !std::isnan(effron);
	    success = success && irplib_sdp_spectrum_set_effron(dest, effron) == CPL_ERROR_NONE;
	}
	cpl_propertylist_delete(l);
	return success;
}

static bool add_stacking_pars(irplib_sdp_spectrum * dest, const cpl_propertylist * plist){

	bool success = !irplib_sdp_spectrum_copy_ncombine(dest, plist, n_combine_tag);
	success = success && !irplib_sdp_spectrum_set_mepoch(dest, CPL_FALSE);
	return success;
}

/* makes use of crder, so must run after add_drs_err() */
static
bool add_radec_err(cpl_propertylist * dest, cpl_propertylist * source) {
    // VLT positioning rms error [3 arcsec] projected equally along RA and DEC
    // rms positioning error: 3 arcsec
    const double err_vlt_ra = 3. / 3600. / sqrt(2.);
    const double err_vlt_dec = 3. / 3600. / sqrt(2.);

    const double posang = get_property_or_nan(source, "ESO ADA POSANG");
    const double crder = get_property_or_nan(dest, idp_output_crder_tag);
    if (std::isnan(posang) || std::isnan(crder)) return false;

    const double cos_posang = cos(posang * PI / 180.0);
    const double sin_posang = sin(posang * PI / 180.0);

    // the error due to the wavelength calibration projected along the RA and
    // DEC axes
    const double err_disp_ra = std::abs(crder / 3600. * sin_posang);
    const double err_disp_dec = std::abs(crder / 3600. * cos_posang);

    cpl_table * slits = mos_load_slits_vimos(source, 0);
    if (!slits) return false;
    const double med_slit_width = cpl_table_get_column_median(slits, "ywidth");
    cpl_table_delete(slits);

    // the possible offset of a target from the center of the slit along the
    // dispersion axis (0.3 x slit width [arcsec]), projected along the RA and
    // DEC axes
    const double err_slit_ra = abs((med_slit_width / 3.) / 3600. * sin_posang);
    const double err_slit_dec = abs((med_slit_width / 3.) / 3600. * cos_posang);

    const double ra_err = sqrt(pow(err_vlt_ra, 2.) + pow(err_disp_ra, 2.) +
            pow(err_slit_ra, 2.));
    const double dec_err = sqrt(pow(err_vlt_dec, 2.) + pow(err_disp_dec, 2.) +
            pow(err_slit_dec, 2.));

    bool success = !cpl_propertylist_update_double(dest, "RA_ERR", ra_err);
    success = success && !cpl_propertylist_set_comment(dest, "RA_ERR",
            "Error on spectroscopic target position [deg]");
    success = success && !cpl_propertylist_update_double(dest, "DEC_ERR",
            dec_err);
    success = success && !cpl_propertylist_set_comment(dest, "DEC_ERR",
            "Error on spectroscopic target position [deg]");

    return success;
}

static bool add_filter(cpl_propertylist * dest, const cpl_propertylist * source){

	const char * filt_name_tag = get_filt_name_tag(source);

	if(!filt_name_tag){
		cpl_msg_warning(cpl_func, "FILTER name extraction failed, aborting IDP generation");
		return false;
	}

	const bool has_filter_name =
			cpl_propertylist_has(source, filt_name_tag) &&
			cpl_propertylist_get_type(source, filt_name_tag) == CPL_TYPE_STRING;

	const char * filt_name = "Free";
	if(has_filter_name){
		filt_name = cpl_propertylist_get_string(source, filt_name_tag);
	}
	return cpl_propertylist_update_string(dest, "FILTER", filt_name) == CPL_ERROR_NONE;
}

static bool
add_from_master_screen_flat(irplib_sdp_spectrum * dest,
		idp_save_info& info,
		const cpl_size extension_idx){

	const char * tag = "MOS_MASTER_SCREEN_FLAT";
	const cpl_frame * frame = cpl_frameset_find_const(info.frames, tag);
	if(frame == NULL){
		cpl_msg_warning(cpl_func, "Unable to extract %s, abort IDP generation", tag);
		return false;
	}

	const cpl_size ext = std::max(cpl_size(0), (extension_idx - 1)* 2 + 1);
	cpl_propertylist * l = cpl_propertylist_load(cpl_frame_get_filename(frame), ext);

	if(l == NULL){
		cpl_msg_warning(cpl_func, "Unable to extract propertylist for %s, abort IDP generation", tag);
		return false;
	}

	bool success = !irplib_sdp_spectrum_copy_lamnlin(dest, l, idp_input_screen_flat_nwave);
	success = success && !irplib_sdp_spectrum_copy_specres(dest, l, idp_input_screen_flat_resolution);
	success = success && !irplib_sdp_spectrum_copy_lamrms(dest, l, idp_input_screen_flat_accuracy);
	cpl_propertylist_delete(l);

	if(!success) return success;

	const double spec_err = irplib_sdp_spectrum_get_lamrms(dest) / sqrt(irplib_sdp_spectrum_get_lamnlin(dest));

	return !irplib_sdp_spectrum_set_specerr(dest, spec_err);
}

static bool add_drs_err(cpl_propertylist * spec_plist,
		const cpl_propertylist * science_plist, const double lam_rms){

	if(std::isnan(lam_rms)) return false;

	const double cd11 = get_double_or_nan(science_plist, "CD1_1");
	const double pszx = get_double_or_nan(science_plist, "ESO DET CHIP1 PSZX");
	const double binx = get_int_or_nan(science_plist, "ESO DET WIN1 BINX");

	const double crder = lam_rms / cd11 * pszx * binx;

	if(std::isnan(crder)) return false;

	bool success = !cpl_propertylist_update_double(spec_plist, idp_output_csyer_tag, 3.0);
	success = success && !cpl_propertylist_set_comment(spec_plist, idp_output_csyer_tag, "Systematic error in position [arcsec]");
	success = success && !cpl_propertylist_update_double(spec_plist, idp_output_crder_tag, crder);
	success = success && !cpl_propertylist_set_comment(spec_plist, idp_output_crder_tag, "Random error in position [arcsec]");

	return success;
}

static cpl_error_code append_metadata(irplib_sdp_spectrum * dest, idp_save_info& info,
		const cpl_size multiplex_idx, cpl_propertylist * spec_plist, const double ra,
		const double dec){

	cpl_propertylist * science_plist = cpl_propertylist_load(cpl_frame_get_filename(info.science_frame), multiplex_idx);

	if(science_plist == NULL)
		return CPL_ERROR_ILLEGAL_INPUT;

	bool success = add_exposures(dest, info.frames, science_plist);
	success = success && add_obs(dest, science_plist);
	success = success && add_pro(dest, science_plist);
	success = success && add_wave_data(dest);
	success = success && add_ron_from_bias(dest, info, science_plist);
	success = success && add_flux_data(dest, science_plist);
	success = success && add_snr(dest);
	success = success && add_stacking_pars(dest, science_plist);
	success = success && add_filter(spec_plist, science_plist);
	success = success && add_from_master_screen_flat(dest, info, multiplex_idx);
	success = success && !irplib_sdp_spectrum_set_specsye(dest, info.spec_sye);
	success = success && !irplib_sdp_spectrum_set_aperture(dest, info.aperture);
	success = success && !irplib_sdp_spectrum_set_ra(dest, ra);
	success = success && !irplib_sdp_spectrum_set_dec(dest, dec);

	if(!success) {
		cpl_propertylist_delete(science_plist);
		return CPL_ERROR_ILLEGAL_OUTPUT;
	}

	const double lam_rms = irplib_sdp_spectrum_get_lamrms(dest);

	success = add_drs_err(spec_plist, science_plist, lam_rms);
	success = success && add_radec_err(spec_plist, science_plist);

	cpl_propertylist_delete(science_plist);

	if(success) return CPL_ERROR_NONE;
	return CPL_ERROR_ILLEGAL_OUTPUT;
}

static cpl_error_code
save_spectrum(idp_save_info& info, const cpl_size spec_idx, const cpl_size multiplex_ex,
		const double ra, const double dec){

	if(std::isnan(ra) || std::isnan(dec)){
		cpl_msg_warning(cpl_func, "RA/DEC unable to extract for %lld, extension %lld", spec_idx, multiplex_ex);
		return CPL_ERROR_ILLEGAL_OUTPUT;
	}

	if(std::isnan(info.spec_sye)){
			cpl_msg_warning(cpl_func, "SPEC_SYE unable to extract for %lld, extension %lld", spec_idx, multiplex_ex);
			return CPL_ERROR_ILLEGAL_OUTPUT;
	}

	if(std::isnan(info.aperture)){
			cpl_msg_warning(cpl_func, "APERTURE unable to extract for %lld, extension %lld", spec_idx, multiplex_ex);
			return CPL_ERROR_ILLEGAL_OUTPUT;
	}

	irplib_sdp_spectrum * spec = extract_1D_spectrum(info, spec_idx);

	if(spec == NULL)
		return CPL_ERROR_ILLEGAL_OUTPUT;

	std::string multiplex_bit = multiplex_ex == 0 ? "_" : "_" + std::to_string(multiplex_ex) + "_" ;

	std::string science_tag =  std::string(cpl_frame_get_tag(info.science_frame)) + "_IDP"+
			multiplex_bit + std::to_string(spec_idx + 1);

	std::string fname = science_tag + ".fits";
	std::transform(fname.begin(), fname.end(), fname.begin(), [](const unsigned char c){return std::tolower(c);});

	cpl_propertylist * plist = cpl_propertylist_new();
	cpl_propertylist_update_string(plist, "ESO PRO CATG", science_tag.c_str());

	cpl_error_code err = append_metadata(spec, info, multiplex_ex, plist, ra, dec);

	if(!err)
		err = irplib_dfs_save_spectrum(info.frames, NULL, info.parlist, info.frames, info.science_frame,
			spec, info.recipe, plist, NULL, NULL, info.pipe_id, "PRO-1.15", fname.c_str());
	else
		cpl_msg_error(cpl_func, "unable to extract metadata for spectrum %lld, extension %lld", spec_idx, multiplex_ex);

	irplib_sdp_spectrum_delete(spec);
	cpl_propertylist_delete(plist);
	return err;
}

static cpl_error_code
save_spectra(idp_save_info& info, const cpl_size mplex, const std::vector<ra_dec>& ra_decs_for_extension){

	const cpl_size sz_y = cpl_image_get_size_y(info.science_img);
	for(cpl_size y = 0; y < sz_y; ++y){
		const ra_dec& ra_dec_coords = ra_decs_for_extension[y];
		const cpl_error_code err = save_spectrum(info, y, mplex, ra_dec_coords.first, ra_dec_coords.second);
		if(err) {
			cpl_msg_error(cpl_func, "Unable to save spectrum %lld, aborting", y);
			return err;
		}
	}
	return CPL_ERROR_NONE;
}

typedef struct{
	double crpix1;
	double crval1;
	double cdelt1;
}wave_meta;

static bool
wave_meta_valid(const wave_meta m){
	return !std::isnan(m.cdelt1) && !std::isnan(m.crpix1) && !std::isnan(m.crval1);
}

static
wave_meta get_wave_meta(const cpl_frame * frame){

	const char * crpix_tag = "CRPIX1";
	const char * crval_tag = "CRVAL1";
	const char * cdelt_tag = "CD1_1";

	cpl_propertylist * l = cpl_propertylist_load(cpl_frame_get_filename(frame), 0);

	wave_meta to_ret{NAN, NAN, NAN};

	const bool has_values = cpl_propertylist_has(l, crpix_tag)
					&& cpl_propertylist_has(l, crval_tag)
					&& cpl_propertylist_has(l, cdelt_tag);

	if(!has_values){
		cpl_propertylist_delete(l);
		return to_ret;
	}

	to_ret.crpix1 = cpl_propertylist_get_double(l, crpix_tag);
	to_ret.crval1 = cpl_propertylist_get_double(l, crval_tag);
	to_ret.cdelt1 = cpl_propertylist_get_double(l, cdelt_tag);

	cpl_propertylist_delete(l);

	return to_ret;
}

cpl_array * get_wlens_angstrom(const cpl_frame * frame, const cpl_size num_wav){

	wave_meta cr_data = get_wave_meta(frame);

	if(!wave_meta_valid(cr_data))
		return NULL;

	cpl_array * wlens = cpl_array_new(num_wav, CPL_TYPE_DOUBLE);

	for(cpl_size w = 1; w <= num_wav; ++w){
		const double val = (w - cr_data.crpix1) * cr_data.cdelt1 + cr_data.crval1;
		cpl_array_set(wlens, w - 1, val);
	}

	return wlens;
}

static double get_property_or_nan(const cpl_propertylist * l, const char * name){
	if(l == NULL) return NAN;

	double v = cpl_propertylist_get_double(l, name);

	if(cpl_error_get_code()){
		cpl_error_reset();
		return NAN;
	}

	return v;
}

static double get_avg_spec_eye(const cpl_frame * offset_frame, const cpl_size start,
		const cpl_size end){

	double ret = 0.0;

	for(cpl_size i = start; i < end; ++i){

		cpl_propertylist * plist = cpl_propertylist_load(cpl_frame_get_filename(offset_frame), i);
		const double spec_eye = get_property_or_nan(plist, idp_input_screen_skyaling_shift);
		cpl_propertylist_delete(plist);

		if(std::isnan(spec_eye)) return NAN;

		ret += spec_eye;
	}

	return ret / (double)(end - start);
}

std::vector<double> extract_spec_sye(const cpl_frameset * fset,
		const cpl_size num_final_extensions){

	const char * sky_tag = "MOS_SCI_DISP_COEFF_SKY";
	const cpl_frame * offsets = cpl_frameset_find_const(fset, sky_tag);

	if(offsets == NULL){
	    cpl_msg_warning(cpl_func, "Unable to extract %s to calculate "
	            "SPEC_SYE, IDP generation aborted", sky_tag);
	    return {};
	}
	const cpl_size n_ext = cpl_frame_get_nextensions(offsets);
	const cpl_size num_mos = cpl_frameset_count_tags(fset, "MOS_SCIENCE");

	const cpl_size expected_extensions = num_mos * num_final_extensions;

	if(expected_extensions != n_ext)
	{
		cpl_msg_warning(cpl_func, "Expecting %lld extensions, "
				"found %lld ins %s", expected_extensions, n_ext, sky_tag);
		return {};
	}

	std::vector<double> to_ret;
	for(cpl_size i = 1; i <= n_ext; i += num_mos){
		const double spec_eye = get_avg_spec_eye(offsets, i, i + num_mos);
		if(std::isnan(spec_eye)){
			cpl_msg_warning(cpl_func, "SPEC_EYE extraction failed in extension range [%lld, %lld[", i, i + num_mos);
			return {};
		}
		to_ret.push_back(spec_eye);
	}

	return to_ret;
}

static double extract_aperture_for_extension(const cpl_frame * science, const cpl_size index){

	cpl_propertylist * plist = cpl_propertylist_load(cpl_frame_get_filename(science), index);
	if(plist == NULL || cpl_error_get_code()){
		cpl_propertylist_delete(plist);
		return NAN;
	}

	cpl_size idx = 1;
	std::vector<double> slits_dim_y;
	while(true){

		std::stringstream tag_name_stream;
		tag_name_stream<<"ESO INS SLIT"<<idx<<" DIMY";
		std::string tag_name = tag_name_stream.str();
		idx++;

		if(!cpl_propertylist_has(plist, tag_name.c_str()))
			break;

		slits_dim_y.push_back(cpl_propertylist_get_double(plist, tag_name.c_str()));
		if(cpl_error_get_code()){
			cpl_error_reset();
			slits_dim_y.clear();
			break;
		}
	}

	cpl_propertylist_delete(plist);

	if(slits_dim_y.empty()) return NAN;

	std::size_t target_idx = slits_dim_y.size() / 2;
	std::vector<double>::iterator target = slits_dim_y.begin() + target_idx;
	std::nth_element(slits_dim_y.begin(), target, slits_dim_y.end());
	const double aperture = slits_dim_y[target_idx] / 3600.0;

	return aperture;
}

static
bool check_posang_from_science(const cpl_frameset * fset){

    cpl_propertylist * science_header = dfs_load_header(fset, "MOS_SCIENCE", 0);
    if(science_header == NULL) return false;

    double posang = -1;
    const bool is_valid = is_posang_valid(science_header, posang);
    if(!is_valid){
        cpl_msg_warning(cpl_func, "ADA.POSANG %f is not allowed in VIMOS MOS, "
                "IDP generation is disabled", posang);
    }
    cpl_propertylist_delete(science_header);
    return is_valid;
}