/*----------------------------------------------------------------------------
   
   File name    :   dark.c
   Author       :	N. Devillard
   Created on   :	January 2001
   Description  :	ISAAC dark recipe

 ---------------------------------------------------------------------------*/

/*
	$Id: dark.c,v 1.6 2002/03/08 14:58:28 yjung Exp $
	$Author: yjung $
	$Date: 2002/03/08 14:58:28 $
	$Revision: 1.6 $
 */

/*----------------------------------------------------------------------------
                                Includes
 ---------------------------------------------------------------------------*/

#include <math.h>
#include "eclipse.h"
#include "isaacp_lib.h"

/*----------------------------------------------------------------------------
                            Private functions
 ---------------------------------------------------------------------------*/

static int isaac_dark_engine(char *, char *, int, int) ;
static int isaac_dark_avg_engine(framelist *, char *) ;
static int isaac_dark_ron_engine(char *, char *, char *) ;
static int isaac_dark_ron_save(char *, char *, char *, double, double,
		double, double, double) ;
static int isaac_dark_compare(char *, char *) ;

/*----------------------------------------------------------------------------
							Main code
 ---------------------------------------------------------------------------*/

int isaac_dark_main(void * dict)
{
    dictionary  *   d ;

    char            argname[10] ;
    char        *   name_i ;
    char        *   name_o ;

    int             avg ;
    int             ron ;

    int             nfiles ;
    int             errors ;
    int             i ;

    /* Initialize */
    avg = 1 ;
    ron = 1 ;

    d = (dictionary*)dict ;
    /* Get options */
    avg = dictionary_getint(d, "arg.average", 1) ;
    ron = dictionary_getint(d, "arg.ron", 1) ;

    /* Get input/output file names */
    nfiles = dictionary_getint(d, "arg.n", -1) ;
    if (nfiles<0) {
        e_error("missing input file name(s): aborting");
        return -1 ;
    }

    /* Loop on input file names */
    errors = 0 ;
    for (i=1 ; i<nfiles ; i++) {
        sprintf(argname, "arg.%d", i);
        name_i = dictionary_get(d, argname) ;
        name_o = dictionary_get(d, "arg.output") ;
        if (name_o == NULL) name_o = strdup(get_rootname(get_basename(name_i)));
        else name_o = strdup(get_rootname(name_o)) ;

        /*
         * Once command-line options have been cleared out, call the main
         * computing function.
         */
        errors += isaac_dark_engine(name_i, name_o, avg, ron) ;
        free(name_o) ;
    }
    return errors ;
}

/*----------------------------------------------------------------------------
						Function ANSI C code 
 ---------------------------------------------------------------------------*/
static int isaac_dark_engine(
        char    *   name_i,
        char    *   name_o,
        int         avg,
        int         ron)
{
    framelist   *   lnames ;
    int             nsettings ;
    framelist   *   sublist ;
    char            outname[FILENAMESZ] ;
    int             i, j ;

    /* Test inputs */
    if (!avg && !ron) return -1 ;

    /* Read the input ASCII file */
    if ((lnames=framelist_load(name_i)) == NULL) {
        e_error("cannot read the input ASCII file") ;
        return -1 ;
    }

    /* Number of different settings */
    if ((nsettings=framelist_labelize(lnames, isaac_dark_compare)) == -1) {
        e_error("in getting the number of different settings") ;
        framelist_del(lnames) ;
        return -1 ;
    }
    e_comment(0, "there are %d different setting(s)", nsettings) ;

    /* For each setting */
    for (i=0 ; i<nsettings ; i++) {
        if ((sublist=framelist_select(lnames, i)) == NULL) {
            e_error("cannot get files for current setting") ;
            framelist_del(lnames) ;
            framelist_del(sublist) ;
            return -1 ;
        }
        e_comment(1, "concerned files:") ;
        for (j=0 ; j<sublist->n ; j++) {
            e_comment(2, "%s", sublist->name[j]) ;
        }

        /* Compute AVG if required */
        if (avg) {
            sprintf(outname, "%s_%02d.fits", name_o, i+1) ;
            isaac_dark_avg_engine(sublist, outname) ;
        }

        /* Compute RON if required */
        if (ron) {
            for (j=0 ; j<sublist->n-1 ; j++) {
				sprintf(outname, "%s_set%02d_pair%02d_ron.paf", name_o, i+1,
						j+1) ;
            	isaac_dark_ron_engine(sublist->name[j], sublist->name[j+1], 
						outname) ;
			}
        }
		framelist_del(sublist) ;
    }

    /* Free and return */
    framelist_del(lnames) ;
    return 0 ;
}

static int isaac_dark_avg_engine(
		framelist	*	in,
		char		*	outname)
{
    cube_t          *   images ;
    image_t         *   avg_dark ;
    qfits_header    *   fh ;
 
    /* Load the cube */ 
    if ((images=cube_load_strings(in->name, in->n)) == NULL) {
        return -1 ;
    }

    /* Average cube */
    if (images->np > 1) {
        avg_dark = cube_avg_linear(images) ;
    } else {
        e_warning("only 1 frame used for this group") ;
        avg_dark = image_copy(images->plane[0]) ;
    }
    cube_del(images) ;

    /* Save with correct keywords */
    if ((fh = qfits_header_read(in->name[0])) == NULL) {
        e_error("cannot read header %s: creating empty header", in->name[0]) ;
    } else {
        if (isaac_header_for_image(fh) != 0) {
            e_error("filtering input header: creating empty header") ;
            qfits_header_destroy(fh) ;
            fh = NULL ;
        }
    }

    isaac_pro_fits(fh, outname, "REDUCED", NULL, isaac_dark_result, "OK", 
			"cal_darks", in->n, in, NULL) ;

    if (isaac_add_files_history(fh, in) == -1) {
        e_warning("cannot write HISTORY keywords in out file") ;
    }

    e_comment(0, "saving file [%s]", outname);
    image_save_fits_hdrdump(avg_dark, outname, fh, BPP_DEFAULT) ;
    qfits_header_destroy(fh) ;
    image_del(avg_dark) ;
    return 0 ;
}

static int isaac_dark_ron_engine(
		char		*	frame1,
		char		*	frame2,
        char        *   outname)
{
	image_t     *   plane1 ;
    image_t     *   plane2 ;
    double          norm ;
    double          noise ;
	double			ron_whole ;
	double			ron_q_ll,
					ron_q_lr,
					ron_q_ur,
					ron_q_ul ;
	int				wmode ;
	int				zone_def[4] ;
    char        *   s ;

	/* Determine instrument mode: SW or LW */
	wmode = isaac_get_wavelength_mode(frame1);

	/* Load current planes */
	if ((plane1 = image_load(frame1)) == NULL) {
		e_error("cannot load plane") ;
		return -1 ;
	}
	if ((plane2 = image_load(frame2)) == NULL) {
		e_error("cannot load plane") ;
		image_del(plane2) ;
		return -1 ;
	}
	
	/* Subtraction */
	if (image_sub_local(plane1, plane2) == -1) {
		e_error("cannot subtract planes") ;
		image_del(plane1) ;
		image_del(plane2) ;
		return -1 ;
	}
	image_del(plane2) ;

	/* Compute norm from NDIT */
	if ((s = isaac_get_ndit(frame1)) == NULL) {
		e_error("cannot get DET.NDIT from [%s]", frame1) ;
		return -1 ;
	}
	norm = atof(s) ;
	norm *= 0.5 ;
	norm = sqrt(norm) ;
	
	/* Compute readout noise according the mode */
	switch (wmode) {
		case ISAAC_SW_MODE:
		/* Get measurement for upper-left quadrant */
		zone_def[0] = 0 ;
		zone_def[1] = plane1->lx/2 ;
		zone_def[2] = plane1->ly/2 ;
		zone_def[3] = plane1->ly-1 ;
		image_readout_noise(plane1, zone_def, &noise, NULL);
		ron_q_ul = noise ;

		/* Get measurement for upper-right quadrant */
		zone_def[0] = plane1->lx/2 ;
		zone_def[1] = plane1->lx-1 ;
		zone_def[2] = plane1->ly/2 ;
		zone_def[3] = plane1->ly-1 ;
		image_readout_noise(plane1, zone_def, &noise, NULL);
		ron_q_ur = noise ;

		/* Get measurement for lower-right quadrant */
		zone_def[0] = plane1->lx/2 ;
		zone_def[1] = plane1->lx-1 ;
		zone_def[2] = 0 ;
		zone_def[3] = plane1->ly/2-1 ;
		image_readout_noise(plane1, zone_def, &noise, NULL);
		ron_q_lr = noise ;

		/* Get measurement for lower-left quadrant */
		zone_def[0] = 0 ;
		zone_def[1] = plane1->lx/2-1 ;
		zone_def[2] = 0 ;
		zone_def[3] = plane1->ly/2-1 ;
		image_readout_noise(plane1, zone_def, &noise, NULL);
		ron_q_ll = noise ;

		/* Normalize the results */
		ron_q_ul *= norm ;
		ron_q_ur *= norm ;
		ron_q_ll *= norm ;
		ron_q_lr *= norm ;
		break ;

		case ISAAC_LW_MODE:
		zone_def[0] = 0 ;
		zone_def[1] = plane1->lx-1 ;
		zone_def[2] = 0 ;
		zone_def[3] = plane1->ly-1 ;
		
		image_readout_noise(plane1, zone_def, &noise, NULL);
		ron_whole = noise ;
		ron_whole *= norm ;
		break ;
		
		default:
		e_error("Unrecognized mode - abort") ;
	    image_del(plane1) ;
		return -1 ;
		break ;
	}
	image_del(plane1) ;

    /* Write the PAF file */
    if (isaac_dark_ron_save(outname, frame1, frame2, ron_whole, ron_q_ul, 
				ron_q_ur, ron_q_lr, ron_q_ll) == -1) {
        e_error("cannot write PAF file") ;
		return -1 ;
    }

	/* Return  */
	return 0 ;
}

static int isaac_dark_ron_save(
		char		*	name_o,
		char		*	frame1,
		char		*	frame2,
		double			ron_whole,
		double			ron_q_ul,
		double			ron_q_ur,
		double			ron_q_lr,
		double			ron_q_ll)
{
	FILE	*	ron_out ;
	char	*	sval ;
	int			wmode ;

	/* Get the wavelength mode */
	wmode = isaac_get_wavelength_mode(frame1);
	
	/* Open output PAF file (formatted ASCII, see fits/pafs.c) */
	e_comment(0, "saving results to %s", name_o);
	ron_out = paf_print_header(name_o,
							   "ISAAC/darks",
							   "Readout noise computation results");
	if (ron_out == NULL) {
		e_error("cannot open file [%s] for output: aborting RON");
		return -1 ;
	}

	/* Add PRO.CATG */
	if ((sval = isaac_get_pro_catg_value(isaac_dark_ron)) != NULL)
    	fprintf(ron_out, "PRO.CATG       \"%s\" ;# Product category\n", sval);

	/* Add date */
	if ((sval = isaac_get_date_obs(frame1)) != NULL) 
    	fprintf(ron_out, "DATE-OBS        \"%s\" ;# Date\n", sval) ;

	/* Add ARCFILE */
	if ((sval = isaac_get_arcfile(frame1)) != NULL)
		fprintf(ron_out, "ARCFILE         \"%s\" ;#\n", sval) ;

	/* Add TPL ID */
    if ((sval = isaac_get_templateid(frame1)) != NULL)
		fprintf(ron_out, "TPL.ID          \"%s\" ;# Template ID\n", sval) ;
    
	/* Add MJD-OBS for file classification */
	if ((sval = qfits_query_hdr(frame1, "MJD-OBS")) == NULL) {
		fprintf(ron_out, "MJD-OBS             0.0 ; # could not find value\n");
	} else {
		fprintf(ron_out, "MJD-OBS             %s ; # Obs start\n", sval);
	}

	/* Add input list of frames */
	fprintf(ron_out, "\n");
	fprintf(ron_out, "PRO.REC1.RAW1.NAME   \"%s\" ;#\n", get_basename(frame1));
	fprintf(ron_out, "PRO.REC1.RAW2.NAME   \"%s\" ;#\n", get_basename(frame2));
	fprintf(ron_out, "\n");

	fprintf(ron_out, "\n");
	/* Forward DET.DIT */
	if ((sval = isaac_get_dit(frame1)) != NULL) {
		fprintf(ron_out, "DET.DIT          %s\n", sval);
	}
	/* Forward DET.NDIT */
	if ((sval = isaac_get_ndit(frame1)) != NULL) {
		fprintf(ron_out, "DET.NDIT         %s\n", sval);
	}
	/* Forward DET.NCORRS */
	if ((sval = isaac_get_romode_id(frame1)) != NULL) {
		fprintf(ron_out, "DET.NCORRS       %s\n", sval);
	}
	/* Forward DPR.TECH */
	sval = isaac_get_dpr_tech(frame1);
	if (sval!=NULL) {
		fprintf(ron_out, "DPR.TECH         \"%s\"\n", sval);
	}
	/* Forward DET.NCORRS.NAME */
	sval = isaac_get_romode_name(frame1);
	if (sval!=NULL) {
		fprintf(ron_out, "DET.MODE.NAME  \"%s\"\n", sval);
	}
	/* Try to forward DET NDSAMPLES if found in input */
	sval = qfits_query_hdr(frame1, "HIERARCH ESO DET NDSAMPLES");
	if (sval!=NULL) {
		sval = qfits_pretty_string(sval);
		fprintf(ron_out, "DET.NDSAMPLES    %s\n", sval);
	}

	fprintf(ron_out,
			"\n"
			"#\n"
			"# Warning:\n"
			"# Read-out noise is measured by computing\n"
			"# pixel standard deviations over a large number\n"
			"# of randomly picked (Poisson-scattered) areas,\n"
			"# which explains why you will get different values\n"
			"# out of each recipe execution. If the method is\n"
			"# correct these values should not vary much, though.\n"
			"#\n"
		"\n");
	
	switch (wmode) {
		case ISAAC_SW_MODE:
		fprintf(ron_out, "QC.UL.RON        %.4f\n", ron_q_ul);
		fprintf(ron_out, "QC.UR.RON        %.4f\n", ron_q_ur);
		fprintf(ron_out, "QC.LR.RON        %.4f\n", ron_q_lr);
		fprintf(ron_out, "QC.LL.RON        %.4f\n", ron_q_ll);
		if (verbose_active()) {
			fprintf(stderr, "RON: %.2f %.2f %.2f %.2f\n",
					ron_q_ul,
					ron_q_ur,
					ron_q_lr,
					ron_q_ll);
		}
		break ;

		case ISAAC_LW_MODE:
		fprintf(ron_out, "QC.RON           %.4f\n", ron_whole);
		if (verbose_active()) {
			fprintf(stderr, "RON: %.2f\n", ron_whole);
		}
		break ;

		default:
		break ;
	}
	fprintf(ron_out, "\n");
	fclose(ron_out) ;
	e_comment(1, "end of read-out noise computation") ;
	return 0 ;
}

static int isaac_dark_compare(
        char    *   file1,
        char    *   file2)
{
    int         comparison ;

    double      dit1,
                dit2 ;
	double		ndit1,
				ndit2 ;
    double      rom1,
                rom2 ;
    int         expno1,
                expno2 ;

    char    *   s ;

    comparison = 1 ;

    /* Compare the DIT  */
    if ((s = isaac_get_dit(file1)) == NULL) {
        e_error("cannot get DET.DIT from [%s]", file1) ;
        return -1 ;
    }
    dit1 = (double)atof(s) ;
    if ((s = isaac_get_dit(file2)) == NULL) {
        e_error("cannot get DET.DIT from [%s]", file2) ;
        return -1 ;
    }
    dit2 = (double)atof(s) ;
    if (fabs(dit1-dit2) > 1e-5) comparison = 0 ;

    /* Compare the NDIT  */
    if ((s = isaac_get_ndit(file1)) == NULL) {
        e_error("cannot get DET.DIT from [%s]", file1) ;
        return -1 ;
    }
    ndit1 = (double)atof(s) ;
    if ((s = isaac_get_ndit(file2)) == NULL) {
        e_error("cannot get DET.DIT from [%s]", file2) ;
        return -1 ;
    }
    ndit2 = (double)atof(s) ;
    if (fabs(ndit1-ndit2) > 1e-5) comparison = 0 ;

    /* Compare the readout mode */
    if (comparison == 1) {
        if ((s = isaac_get_romode_id(file1)) == NULL) {
            e_error("cannot get DET.NCORRS from [%s]", file1) ;
            return -1 ;
        }
        rom1 = (double)atof(s) ;
        if ((s = isaac_get_romode_id(file2)) == NULL) {
            e_error("cannot get DET.NCORRS from [%s]", file2) ;
            return -1 ;
        }
        rom2 = (double)atof(s) ;
        if (fabs(rom1-rom2) > 1e-5) comparison = 0 ;
    }

    /* Files have to be consequtive */
    if (comparison == 1) {
        if ((s = isaac_get_current_exp_nb(file1)) == NULL) {
            e_error("cannot get TPL.EXPNO from [%s]", file1) ;
            return -1 ;
        }
        expno1 = (int)atoi(s) ;
        if ((s = isaac_get_current_exp_nb(file2)) == NULL) {
            e_error("cannot get TPL.EXPNO from [%s]", file2) ;
            return -1 ;
        }
        expno2 = (int)atoi(s) ;
        if (fabs(expno1 - expno2) > 1) comparison = 0 ;
    }
    return comparison ;
}