/*-------------------------------------------------------------------------*/
/**
   @file	dark.c
   @author	Y. Jung
   @date	January 2002
   @version	$Revision: 1.5 $
   @brief	CONICA dark recipe
*/
/*--------------------------------------------------------------------------*/

/*
   $Id: dark.c,v 1.5 2002/03/08 12:27:20 yjung Exp $
   $Author: yjung $
   $Date: 2002/03/08 12:27:20 $
   $Revision: 1.5 $
 */

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

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include "eclipse.h"
#include "conicap_lib.h"

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

static int conica_dark_engine(char *, char *, int, int) ;
static int conica_dark_avg_engine(framelist *, char *) ; 
static int conica_dark_ron_engine(char *, char *, char *) ;
static int conica_dark_compare(char *, char *) ;
static int conica_dark_ron_save(char *, char *, char *, double) ; 

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

int conica_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 += conica_dark_engine(name_i, name_o, avg, ron) ;
        free(name_o) ;
    }
    return errors ;
}

static int conica_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, conica_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) ;
			conica_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) ;
                conica_dark_ron_engine(sublist->name[j], sublist->name[j+1],
                                                outname) ;
            }
        }
        framelist_del(sublist) ;
	}

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

static int conica_dark_ron_save(
		char		*	outname, 
        char        *   frame1,
        char        *   frame2,
        double          ron)
{
	FILE    *   ron_out ;
	char    *   sval ;
	
	/* Open output PAF file (formatted ASCII, see fits/pafs.c) */
    e_comment(0, "saving results to %s", outname);
    if ((ron_out = paf_print_header(outname,
                               "CONICA/dark",
                               "Readout noise computation results")) == NULL) {
        e_error("cannot open file [%s] for output", outname) ;
        return -1 ;
    }
    /* Add date */
	if ((sval = conica_get_date_obs(frame1)) != NULL) {
		fprintf(ron_out, "DATE-OBS           \"%s\" ; #Date\n", sval) ;
	}
    /* Add ARCFILE */
    if ((sval = conica_get_arcfile(frame1)) != NULL)
        fprintf(ron_out, "ARCFILE         \"%s\" ;#\n", sval) ;
    /* Add TPL ID */
	if ((sval = conica_get_templateid(frame1)) != NULL) {
		fprintf(ron_out, "TPL.ID             \"%s\"; # Template id\n", sval) ;
	}
	
    fprintf(ron_out,"#\n") ;
    fprintf(ron_out,"# Read-out noise measurements\n") ;
    fprintf(ron_out,"#\n") ;

    /* 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\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 = conica_get_dit(frame1)) != NULL) {
        fprintf(ron_out, "DET.DIT          \"%s\"\n", sval) ; 
    }
    /* Forward DET.NDIT */
    if ((sval = conica_get_ndit(frame1)) != NULL) {
        fprintf(ron_out, "DET.NDIT         \"%s\"\n", sval) ; 
    }
    /* Forward DET.NCORRS */
    if ((sval = conica_get_rom(frame1)) != NULL) {
        fprintf(ron_out, "DET.NCORRS       \"%s\"\n", sval) ; 
    }
    /* Forward DPR.TECH */
    if ((sval = conica_get_dpr_tech(frame1)) != NULL) {
        fprintf(ron_out, "DPR.TECH         \"%s\"\n", sval) ; 
    }
    /* Forward DET.NCORRS.NAME */
    if ((sval = conica_get_rom_name(frame1)) != NULL) {
        fprintf(ron_out, "DET.NCORRS.NAME   \"%s\"\n", sval) ; 
    }
    /* Try to forward DET NDSAMPLES if found in input */
    if ((sval=qfits_query_hdr(frame1, "HIERARCH ESO DET NDSAMPLES")) != NULL) {
        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");
    
    fprintf(ron_out, "QC.RON           %.4f\n", ron);
    
    if (verbose_active()) {
        fprintf(stderr, "RON: %.2f\n", ron);
    }
    fprintf(ron_out, "\n");

    fclose(ron_out) ;
    e_comment(1, "end of read-out noise computation") ;
    return 0 ;
}

static int conica_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 (conica_header_for_image(fh) != 0) {
			e_error("filtering input header: creating empty header") ;
			qfits_header_destroy(fh) ;
			fh = NULL ;
		}
	}

	conica_pro_fits(fh, outname, "REDUCED", NULL, "OK", "cal_darks", in->n, 
			in, NULL) ; 
        	
	if (conica_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 conica_dark_ron_engine(
        char    *   frame1,
        char    *   frame2,
        char    *   outname)
{
	double			ron ;	
	image_t		*	plane1 ;
	image_t		*	plane2 ;
	double			norm ;
	double			noise ;
	char		*	s ;

	int				i ;

    /* 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 readout noise */
    image_readout_noise(plane1, NULL, &noise, NULL) ;
    image_del(plane1) ;

    /* Compute norm from NDIT */
    if ((s = conica_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 RON */
    ron = noise * norm ;

    /* Write out the PAF file */
    if (conica_dark_ron_save(outname, frame1, frame2, ron) == -1) {
        e_error("cannot write PAF file") ;
        return -1 ;
    }

	return 0 ;
}

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

    double		exptime1, 
				exptime2 ;
	double		rom1,
				rom2 ;
    char        mode1[FILENAMESZ] ;
    char        mode2[FILENAMESZ] ;
	int			expno1,
				expno2 ;
	
    char    *   s ;
        
    comparison = 1 ;    
    
	/* Compare the EXPTIME  */
	if ((s = conica_get_exptime(file1)) == NULL) {
		e_error("cannot get EXPTIME from [%s]", file1) ;
		return -1 ;
	}
	exptime1 = (double)atof(s) ;
	if ((s = conica_get_exptime(file2)) == NULL) {
		e_error("cannot get EXPTIME from [%s]", file2) ;
		return -1 ;
	}
	exptime2 = (double)atof(s) ;
	if (fabs(exptime1-exptime2) > 1e-5) comparison = 0 ; 

	/* Compare the readout mode */
	if (comparison == 1) {
		if ((s = conica_get_rom(file1)) == NULL) {
			e_error("cannot get DET.NCORRS from [%s]", file1) ;
			return -1 ;
		}
		rom1 = (double)atof(s) ;
		if ((s = conica_get_rom(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 ;
	}

	/* Compare the detector mode */
    if (comparison == 1) {
        if ((s = conica_get_mode(file1)) == NULL) {
            e_error("cannot get DET.MODE.NAME from [%s]", file1) ;
            return -1 ;
        }
        strcpy(mode1, s) ;
        if ((s = conica_get_mode(file2)) == NULL) {
            e_error("cannot get DET.MODE.NAME from [%s]", file2) ;
            return -1 ;
        }
        strcpy(mode2, s) ;
        if (strcmp(mode1, mode2) != 0) comparison = 0 ;
    }

	/* Files have to be consequtive */
	if (comparison == 1) {
		if ((s = conica_get_expno(file1)) == NULL) {
			e_error("cannot get TPL.EXPNO from [%s]", file1) ;
			return -1 ;
		}
		expno1 = (int)atoi(s) ;
		if ((s = conica_get_expno(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 ;
}