/*----------------------------------------------------------------------------
   
   File name    :   detlin.c
   Author       :	Y. Jung
   Created on   :	June 2001
   Description  :	CONICA detector linearity test.

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

/*
   $Id: detlin.c,v 1.5 2001/10/22 12:23:26 ndevilla Exp $
   $Author: ndevilla $
   $Date: 2001/10/22 12:23:26 $
   $Revision: 1.5 $
 */

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

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

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

/*----------------------------------------------------------------------------
                                Defines
 ---------------------------------------------------------------------------*/

#define	FRAME_INT		2001
#define	FRAME_STDEV		2002
#define	FRAME_DARK		2003

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

static int conica_detlin_engine(char *, char *, int);
static double3 * conica_detlin_get_DITs(framelist *, framelist *, framelist *) ;

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

int conica_detlin_main(void * dict)	
{
	dictionary  *   d ;
	int				display ;
	
	char            argname[10] ;
	char        *   name_i ;
	char        *   name_o ;
	int				nfiles ;

	int				errors ;
	int				i ;
	
	d = (dictionary*)dict ;
	/* Get options */
	display = dictionary_getint(d, "arg.display", 0) ;
	
	/* 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_detlin_engine(name_i, name_o, display) ;
        free(name_o) ;
    }
    return errors ;
}

static int conica_detlin_engine(
		char 	* 	name_i, 
		char 	* 	name_o, 
		int			display)
{
	framelist		*	input_list ;
	framelist		*	int_list ;
	framelist		*	dark_list ;
	framelist		*	stdev_list ;
	double3			*	dits_means_var ;
	double3			*	local_double3 ;
	double				stdev_val ;
	cube_t			*	cube_int ;
	cube_t			*	cube_dark ;
	cube_t			*	cube_stdev ;
	cube_t 			*	linearity_cu ;
	double			*	fitted ;
	int					dark ;
	int					stdev ;
	char			*	sval ;
	double				gain,
						ron ;
	gnuplot_ctrl    *   handle ;
	char                cmd[FILENAMESZ] ;
	char                outname[FILENAMESZ] ;
	int				i ;

	/* Initialize */
	dark = 1 ;
	stdev = 1 ;
	
	/* Load input framelist */
	if ((input_list = framelist_load(name_i)) == NULL) {
		e_error("cannot load input frame list") ;
		return -1 ;
	}
	
	/* Assign labels to frames */
	for (i=0 ; i<input_list->n ; i++) {
		if ((sval = conica_get_frame_type(input_list->name[i])) == NULL) {
			e_error("cannot read frame type in header - abort") ;
			framelist_del(input_list) ;
			return -1 ;
		} else {
			if (!strcmp(sval, "STDEV")) input_list->label[i] = FRAME_STDEV ;
			else if (!strcmp(sval, "INT")) input_list->label[i] = FRAME_INT ;
			else e_warning("invalid frame type !") ;
		}
		if ((sval=qfits_query_hdr(input_list->name[i], "DPR.TYPE")) != NULL) {
			if (!strcmp(sval, "DARK")) input_list->label[i] = FRAME_DARK ;
		} else {
			e_warning("cannot read DPR TYPE keyword") ;
			dark = 0 ;
		}
	}
	/* Separate frames INT STDEV and DARK */
	if ((int_list = framelist_select(input_list, FRAME_INT)) == NULL) {
		e_error("no INT frames") ;
		framelist_del(input_list) ;
		return -1 ;
	}
	if ((stdev_list = framelist_select(input_list, FRAME_STDEV)) == NULL) {
		e_warning("no STDEV frames") ;
		stdev = 0 ;
	}
	if (dark == 1) {
		if ((dark_list = framelist_select(input_list, FRAME_DARK)) == NULL) {
			e_warning("no DARK frames") ;
			dark = 0 ;
		}
	}
	framelist_del(input_list) ;

	if (stdev == 1) {
		if (int_list->n != stdev_list->n) {
			e_error("nb of ST_DEV frames <> nb of INT frames !") ;
			stdev = 0 ;
			framelist_del(stdev_list) ;	
			stdev_list = NULL ;
		}
	}
	if (dark == 1) {
		if (int_list->n != dark_list->n) {
			e_error("nb of INT frames <> nb of DARK frames !") ;
			dark = 0 ;
			framelist_del(dark_list) ;
			dark_list = NULL ;
		}
	}

	/* Check if the DITS are the same in both lists and get them	 */
	if ((dits_means_var = conica_detlin_get_DITs(int_list, stdev_list, 
					dark_list)) == NULL) {
		e_error("cannot get DIT values") ;
		framelist_del(int_list) ;
		if (stdev == 1) framelist_del(stdev_list) ;
		if (dark == 1) framelist_del(dark_list) ;
		return -1 ;
	}

	/* Load INT frames */
	e_comment(0, "Load INT frames") ;
	if ((cube_int = cube_load_strings(int_list->name, int_list->n)) == NULL) {
		e_error("cannot load INT frames") ;
		framelist_del(int_list) ;
		if (stdev == 1) framelist_del(stdev_list) ;
		if (dark == 1) framelist_del(dark_list) ;
		double3_del(dits_means_var) ;
		return -1 ;
	}
	framelist_del(int_list) ;

	/* subtract INT by dark if dark is present */
	if (dark == 1) {
		e_comment(0, "Subtract dark") ;
		cube_dark = cube_load_strings(dark_list->name, dark_list->n) ;
		cube_sub(cube_int, cube_dark) ;
		cube_del(cube_dark) ;
	}
	if (dark == 1) framelist_del(dark_list) ;
	
	/* Compute linearity */
	e_comment(0, "Compute linearity") ;
	if ((linearity_cu = detector_linearity_fit(cube_int, 
									dits_means_var->x)) == NULL) {
		e_error("computing linearity") ;
		if (stdev == 1) framelist_del(stdev_list) ;
		double3_del(dits_means_var) ;
		cube_del(cube_int) ;
		return -1 ;
	}

	/* Write out linearity coefficients frames */
	e_comment(0, "Write out output_lin table") ;
	sprintf(outname, "%s_a.fits", name_o) ;
	image_save_fits(linearity_cu->plane[0], outname, BPP_DEFAULT) ;
	sprintf(outname, "%s_b.fits", name_o) ;
	image_save_fits(linearity_cu->plane[1], outname, BPP_DEFAULT) ;
	sprintf(outname, "%s_c.fits", name_o) ;
	image_save_fits(linearity_cu->plane[2], outname, BPP_DEFAULT) ;
	sprintf(outname, "%s_rms.fits", name_o) ;
	image_save_fits(linearity_cu->plane[3], outname, BPP_DEFAULT) ;
	cube_del(linearity_cu) ;

	if (stdev == 0) {
		e_warning("No stdev frames - no Gain and RON computation") ;
		double3_del(dits_means_var) ;
		cube_del(cube_int) ;
		return -1 ;
	}
	
	/* Load ST_DEV frames */
	e_comment(0, "Load ST_DEV frames") ;
	if ((cube_stdev = cube_load_strings(stdev_list->name, 
										stdev_list->n)) == NULL) {
		e_error("cannot load STDEV frames") ;
		framelist_del(stdev_list) ;
		double3_del(dits_means_var) ;
		cube_del(cube_int) ; 
		return -1 ;
	}
	framelist_del(stdev_list) ;

	/* Fill dits_means_var */
	for (i=0 ; i<dits_means_var->n ; i++) {
		dits_means_var->y[i] = image_getmean(cube_int->plane[i]) ;
		stdev_val = image_getmean(cube_stdev->plane[i]) ;
		dits_means_var->z[i] = stdev_val * stdev_val ;
	}
	cube_del(cube_int) ; 
	cube_del(cube_stdev) ; 

	/* Comput Gain and RON */
	e_comment(0, "Compute Gain and RON") ;

	local_double3 = double3_new(dits_means_var->n) ;
	for (i=0 ; i<dits_means_var->n ; i++) {
		local_double3->x[i] = dits_means_var->y[i] ;
		local_double3->y[i] = dits_means_var->z[i] ;
	}
	fitted = fit_slope_robust(local_double3) ;
	double3_del(local_double3) ;
	gain = 1/fitted[1] ;
	ron = sqrt(fitted[0]) * gain ;
	e_comment(1, "COMPUTED GAIN: %g\n", gain) ;
	e_comment(1, "COMPUTED RON : %g\n", ron) ;

	if (display) {
		handle = gnuplot_init() ;
		gnuplot_setstyle(handle, "points") ;
		gnuplot_set_xlabel(handle, "counts[ADU]") ;
		gnuplot_set_ylabel(handle, "var[ADU^2]") ;
		gnuplot_plot_xy(handle, dits_means_var->y, dits_means_var->z, 
				dits_means_var->n, "") ;
		printf("press enter to continue\n") ;
		while (getchar() != '\n') {}
		sprintf(cmd, "replot %g+%g*x\n", fitted[0], fitted[1]) ;
		gnuplot_cmd(handle, cmd) ;
		printf("press enter to continue\n") ;
		while (getchar() != '\n') {}
		gnuplot_close(handle) ;
	}
	free(fitted) ;
	double3_del(dits_means_var) ;

	return 0 ;
}


/*-------------------------------------------------------------------------*/
/**
  @name		conica_detlin_get_DITs
  @memo		Read DIT values 
  @param	list1 	first frame list
  @param	list2	second frame list
  @param	list3	third frame list (may be NULL)
  @return	list of DITs (NULL if not the same DITs in both lists)
 */
/*--------------------------------------------------------------------------*/
static double3 * conica_detlin_get_DITs(
		framelist	*	list1,
		framelist	*	list2,
		framelist	*	list3)
{
	double3	*	dits ;
	char	*	sval ;
	int			i ;
	
	dits = double3_new(list1->n) ;
	/* For each input frame */
	for (i=0 ; i<list1->n ; i++) {
		/* Get the DIT in first list */
		if ((sval = conica_get_dit(list1->name[i])) == NULL) {
			e_error("cannot read DIT in input file") ;
			double3_del(dits) ;
			return NULL ;
		} else dits->x[i] = (double)atof(sval) ;
	}
	
	/* If list 2 (stdev) is provided, verify the DITs */
	if (list2 != NULL) {
		for (i=0 ; i<list2->n ; i++) {
			if ((sval = conica_get_dit(list2->name[i])) == NULL) {
				e_error("cannot read DIT in input file") ;
				double3_del(dits) ;
				return NULL ;
			} else {
				if (fabs(dits->x[i]-(double)atof(sval)) > 1e-3) {
					e_error("not the same DITs in stdev frames") ;
					double3_del(dits) ;
					return NULL ;
				}
			}
		}
	}
	
	/* If list 3 (dark) is provided, verify the DITs */
	if (list3 != NULL) {
		for (i=0 ; i<list3->n ; i++) {
			if ((sval = conica_get_dit(list3->name[i])) == NULL) {
				e_error("cannot read DIT in input file") ;
				double3_del(dits) ;
				return NULL ;
			} else {
				if (fabs(dits->x[i]-(double)atof(sval)) > 1e-3) {
					e_error("not the same DITs in dark frames") ;
					double3_del(dits) ;
					return NULL ;
				}
			}
		}
	}
	return dits ;
}