/*----------------------------------------------------------------------------
   
   File name    :   lw_detlin.c
   Author       :	N. Devillard
   Created on   :	April 2001
   Description  :	ISAAC LW detector linearity test.

Inputs:
	- A list of frames to process, with various DITs like:
	0.1384
	0.2
	0.3
	0.4
	0.1384
	0.5
	0.6
	0.7
	0.8
	0.1384
	0.9
	1.0
	1.1
	1.2
	0.1384
	- A list of corresponding dark frames (same DITs as the ones used
	above).

Process:
	- Subtract darks from input frames.
	- Check the stability of the level in the DIT=0.1384 frames.
	  exit if changes too much (1% level).
	- Use all frames but 0.1384 frames. Fit to each pixel the function
	  DIT = a*flux + b*flux^2 + c*flux^3
	- Determine a, b, c, fit error and chi-square estimate of the goodness
	  of fit.
	- Construct 4 images: a, b, c, goodness of fit.

Outputs:
	- Image of a, b, c coefficients
	- Image of the goodness of fit.


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

/*
   $Id: lw_detlin.c,v 1.13 2001/10/22 12:23:34 ndevilla Exp $
   $Author: ndevilla $
   $Date: 2001/10/22 12:23:34 $
   $Revision: 1.13 $
 */

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

#include <stdio.h>
#include <string.h>
#include "eclipse.h"
#include "isaacp_lib.h"

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

#define FRAME_DARK		1
#define FRAME_LAMP		2

#define NPARTS			6	

/*----------------------------------------------------------------------------
                            Private functions
 ---------------------------------------------------------------------------*/
static int isaac_lw_detlin_engine(char * name_i, char * name_o);
static cube_t * isaac_lw_detlin_load(char * listname, double ** ditval);
static int isaac_lw_detlin_save(cube_t*, char*, char*, int);


static int part=0 ;

/*----------------------------------------------------------------------------
                            Main code
 ---------------------------------------------------------------------------*/
int isaac_lw_detlin_main(void * dict)
{
	dictionary	*	d ;
	char		*	name_i ;
	char		*	name_o ;

	d = (dictionary*)dict ;

	/* No options */

	/* Get input/output file names */
	name_i = dictionary_get(d, "arg.1");
	if (name_i==NULL) {
		e_error("missing input file name: aborting");
		return -1 ;
	}
	name_o = dictionary_get(d, "arg.output");
	if (name_o==NULL) {
		name_o = "detlin" ;
	}

	return isaac_lw_detlin_engine(name_i, name_o);
}

static int isaac_lw_detlin_engine(char * name_i, char * name_o)
{
	cube_t	*	detlin ;
	double	*	ditval ;
	cube_t	*	fitres ;
	int			datancom ;
	int			sta ;

	/* Load inputs */
	detlin = isaac_lw_detlin_load(name_i, &ditval);
	if (detlin==NULL) {
		return -1 ;
	}
	datancom = detlin->np ;
	part ++ ;
	e_comment(0, "-> part %d of %d: fitting polynomials (long)", part, NPARTS);
	fitres = detector_linearity_fit(detlin, ditval);
	/* Discard data */
	cube_del(detlin);
	free(ditval);

	if (fitres==NULL) {
		e_error("fitting function to planes: aborting");
		return -1 ;
	}

	/* Save results */
	sta = isaac_lw_detlin_save(fitres, name_i, name_o, datancom);
	cube_del(fitres);

	e_comment(0, "done");
	return sta ;
}


/*
 * Load all input frames, check that there are as many darks as linearity
 * images and that they have the corresponding integration times.
 * The files with identical DITS are also checked for intensity variations.
 */
static cube_t * isaac_lw_detlin_load(char * listname, double ** ditval)
{
	framelist	*	in_list ;
	framelist	*	dark_list ;
	framelist	*	lamp_list ;
	int				i, j ;
	char		*	sval ;
	int				err ;
	int				n_dark, n_lamp ;
	cube_t		*	lampcube ;
	char		*	dark_integ,
				*	lamp_integ ;
	char		*	init_dit ;
	int			*	same_dit ;
	int				n_same_dit ;
	double		*	level_same_dit ;
	double		*	ditval_load ;
	image_t		*	lamp_1,
				*	dark_1 ;
	int				lx, ly, np ;

	/* Load framelist */
	part++;
	e_comment(0, "-> part %d of %d: frame identification", part, NPARTS);
	in_list = framelist_load(listname);
	if (in_list==NULL) {
		e_error("cannot load %s", listname);
		return NULL ;
	}
	e_comment(1, "framelist [%s] parsed Ok", listname);
	/* Assign labels to frames */
	err=0 ;
	n_dark=n_lamp=0 ;
	for (i=0 ; i<in_list->n ; i++) {
		sval = qfits_query_hdr(in_list->name[i], "DPR.TYPE");
		if (sval==NULL) {
			e_error("no DPR TYPE for frame %s", in_list->name[i]);
			err++ ;
		} else {
			if (!strcmp(sval, "DARK")) {
				/* Frame is a dark */
				in_list->label[i] = FRAME_DARK ;
				n_dark++ ;
			} else if (!strcmp(sval, "LAMP")) {
				/* Frame is a linearity image */
				in_list->label[i] = FRAME_LAMP ;
				n_lamp++ ;
			} else {
				/* Invalid DPR TYPE for this frame */
				e_error("invalid DPR TYPE for frame %s: [%s]",
						in_list->name[i],
						sval);
				err++;
			}
		}
	}
	/* Check that there are as many darks as input images */
	if (n_dark!=n_lamp) {
		e_error("inconsistent data set: %d darks for %d images",
				n_dark,
				n_lamp);
		err++;
	}
	if (err) {
		e_error("%d error(s) parsing list %s", err, listname);
		framelist_del(in_list);
		return NULL ;
	}
	e_comment(1, "all frames correctly labelled");
	/* Create new framelists for linearity and dark frames */
	lamp_list = framelist_select(in_list, FRAME_LAMP);
	dark_list = framelist_select(in_list, FRAME_DARK);
	/* Discard initial list */
	framelist_del(in_list);

	/*
	 * Check out that they have consistent integration times
	 * Remember which frames have the same integration time.
	 */
	part++;
	e_comment(0, "-> part %d of %d: checking DIT consistency", part, NPARTS);
	err=0 ;
	same_dit = malloc(lamp_list->n * sizeof(int));
	ditval_load = malloc(lamp_list->n * sizeof(double));
	n_same_dit=0 ;
	for (i=0 ; i<lamp_list->n ; i++) {
		/* Get integration time for lamps */
		lamp_integ = qfits_query_hdr(lamp_list->name[i], "DET.DIT");
		if (lamp_integ==NULL) {
			e_error("frame %s has no DET.DIT", lamp_list->name[i]);
			err++ ;
			break ;
		}
		e_comment(1, "LAMP %s DIT %s",
				  get_basename(lamp_list->name[i]),
				  lamp_integ);
		/* If first loaded frame, record integration time */
		if (i==0) {
			init_dit = strdup(lamp_integ);
			same_dit[0] = 1 ;
			n_same_dit++ ;
		} else {
			/* Store integration time */
			if (!strcmp(init_dit, lamp_integ)) {
				same_dit[i] = 1 ;
				n_same_dit ++ ;
			} else {
				same_dit[i] = 0 ;
			}
		}
		ditval_load[i] = (double)atof(lamp_integ);
		/* Get integration time for dark */
		dark_integ = qfits_query_hdr(dark_list->name[i], "DET.DIT");
		if (dark_integ==NULL) {
			e_error("frame %s has no DET.DIT: aborting", dark_list->name[i]);
			err++ ;
			break ;
		}
		e_comment(1, "DARK %s DIT %s",
				  get_basename(dark_list->name[i]),
				  lamp_integ);
		/* Compare DIT for lamp and dark */
		if (strcmp(lamp_integ, dark_integ)) {
			e_error("DIT inconsistency");
			e_error("file %s has DIT=%s", lamp_list->name[i], lamp_integ);
			e_error("file %s has DIT=%s", dark_list->name[i], dark_integ);
			err++ ;
		}
	}
	if (init_dit!=NULL)
		free(init_dit);

	/* Check that there are frames with identical DITs */
	if (n_same_dit<1) {
		e_error("no two frames with identical DIT");
		err++ ;
	}
	/* Report errors */
	if (err) {
		e_error("%d error(s) in data set", err);
		free(ditval_load);
		free(same_dit);
		framelist_del(lamp_list);
		framelist_del(dark_list);
		return NULL ;
	}
	e_comment(1, "DIT consistency Ok");

	/* Compute level in frames of identical DIT */
	part++;
	e_comment(0, "-> part %d of %d: checking lamp stability", part, NPARTS);
	level_same_dit = malloc(n_same_dit * sizeof(double));
	j=0 ;
	lx=ly=0 ;
	for (i=0 ; i<n_lamp ; i++) {
		if (same_dit[i]) {
			/* Load lamp frame */
			lamp_1 = image_load(lamp_list->name[i]);
			if (lamp_1==NULL) {
				e_error("loading frame %s: aborting", lamp_list->name[i]);
				free(ditval_load);
				free(level_same_dit) ;
				free(same_dit);
				framelist_del(lamp_list);
				framelist_del(dark_list);
				return NULL ;
			}
			/* Load dark frame */
			dark_1 = image_load(dark_list->name[i]);
			if (dark_1==NULL) {
				e_error("loading frame %s: aborting", dark_list->name[i]);
				free(ditval_load);
				free(level_same_dit);
				free(same_dit);
				framelist_del(lamp_list);
				framelist_del(dark_list);
				return NULL ;
			}
			if (lx==0 || ly==0) {
				lx = lamp_1->lx ;
				ly = lamp_1->ly ;
			}
			/* Subtract dark from lamp */
			image_sub_local(lamp_1, dark_1);
			/* Discard dark frame */
			image_del(dark_1);
			/* Record level in subtracted frame */
			level_same_dit[j] = image_getmean(lamp_1);
			/* Discard lamp frame */
			image_del(lamp_1);
			e_comment(1, "level for LAMP %02d: %g",
						i+1,
						level_same_dit[j]);
			j++ ;
		}
	}
	/* Check level in frames of identical DIT */
	e_comment(1, "checking level in frames");
	err=0 ;
	for (i=1 ; i<n_same_dit ; i++) {
		if (((level_same_dit[i] - level_same_dit[0]) / level_same_dit[0])>0.01) {
			e_error("level difference # %d too high", i+1);
			err++ ;
		}
	}
	/* Discard levels */
	free(level_same_dit);
	if (err) {
		e_error("too much difference in frames of identical DIT: aborting");
		free(ditval_load);
		free(same_dit);
		framelist_del(lamp_list);
		framelist_del(dark_list);
		return NULL ;
	} else {
		e_comment(1, "lamp level check Ok");
	}

	/* Load frames and subtract them as they load */
	part++;
	e_comment(0, "-> part %d of %d: load dark-subtracted frames", part, NPARTS);
	np = lamp_list->n - n_same_dit + 1;
	same_dit[0] = 0 ;
	lampcube = cube_new(lx, ly, np);
	j=0 ;
	err=0 ;
	for (i=0 ; i<lamp_list->n ; i++) {
		if (same_dit[i]==0) {
			e_comment(1, "loading/subtracting DIT %g", ditval_load[i]);
			lamp_1 = image_load(lamp_list->name[i]);
			if (lamp_1==NULL) {
				e_error("loading frame %s", lamp_list->name[i]);
				err++ ;
				break ;
			}
			dark_1 = image_load(dark_list->name[i]);
			if (dark_1==NULL) {
				e_error("loading frame %s", dark_list->name[i]);
				err++ ;
				break ;
			}
			image_sub_local(lamp_1, dark_1);
			image_del(dark_1);
			lampcube->plane[j]=lamp_1 ;
			j++ ;
			if (j>np) {
				e_error("internal: counting planes");
				err++ ;
				break ;
			}
		}
	}
	framelist_del(lamp_list);
	framelist_del(dark_list);
	free(same_dit);
	if (err) {
		e_error("loading data: aborting");
		free(ditval_load);
		cube_del(lampcube);
		return NULL ;
	}
	e_comment(1, "frame loading Ok");

	/* Assign list of loaded DITs */
	(*ditval) = ditval_load ;
	/* Returned loaded cube */
	return lampcube ;
}


static int isaac_lw_detlin_save(
	cube_t 	* fitres,
	char 	* name_i,
	char 	* name_o,
	int		  datancom)
{
	double			med ;
	image_t		*	div ;
	char			outname[FILENAMESZ];
	char		*	refname ;
	qfits_header*	fh,
				*	fh_spec ;
	framelist	*	raw ;
				
	part++;
	e_comment(0, "-> part %d of %d: saving results", part, NPARTS);

	/* Compute B/A and find its median */
	div = image_div(fitres->plane[1], fitres->plane[0]);
	med = image_getmedian(div);
	image_del(div);
	printf("median B/A: %g\n", med);

	/* Compute C/A and find its median */
	div = image_div(fitres->plane[2], fitres->plane[0]);
	med = image_getmedian(div);
	image_del(div);
	printf("median C/A: %g\n", med);

	/* Load header from first input file */
	if (is_ascii_list(name_i)==1) {
		refname = framelist_firstname(name_i) ;
	} else {
		refname = name_i ;
	}
	if ((fh=qfits_header_read(refname))==NULL) {
		e_error("getting header from reference frame [%s]", refname);
		return -1 ;
	}
	/* Prepare header for image output */
	isaac_header_for_image(fh);

	/* Save coeff A image */
	fh_spec = qfits_header_copy(fh);
	sprintf(outname, "%s_A.fits", get_basename(name_o));
	e_comment(1, "saving image [%s]", outname);
	raw = framelist_load(name_i) ;
	isaac_pro_fits(	fh_spec,
					outname,
					NULL,	/* Not a reduced file */
					NULL,	/* No reduction level applicable */
					isaac_imag_lw_detlin_coeff_A,
					"OK",	/* Status */
					"lw_detlin",	/* Recipe ID */
					datancom,
					raw,
					NULL);
	image_save_fits_hdrdump(fitres->plane[0], outname, fh_spec, BPP_DEFAULT);
	qfits_header_destroy(fh_spec);

	/* Save coeff B image */
	fh_spec = qfits_header_copy(fh);
	sprintf(outname, "%s_B.fits", get_basename(name_o));
	e_comment(1, "saving image [%s]", outname);
	isaac_pro_fits(	fh_spec,
					outname,
					NULL,	/* Not a reduced file */
					NULL,	/* No reduction level applicable */
					isaac_imag_lw_detlin_coeff_B,
					"OK",	/* Status */
					"lw_detlin",	/* Recipe ID */
					datancom,
					raw,
					NULL);
	image_save_fits_hdrdump(fitres->plane[1], outname, fh_spec, BPP_DEFAULT);
	qfits_header_destroy(fh_spec);

	/* Save coeff C image */
	fh_spec = qfits_header_copy(fh);
	sprintf(outname, "%s_C.fits", get_basename(name_o));
	e_comment(1, "saving image [%s]", outname);
	isaac_pro_fits(	fh_spec,
					outname,
					NULL,	/* Not a reduced file */
					NULL,	/* No reduction level applicable */
					isaac_imag_lw_detlin_coeff_C,
					"OK",	/* Status */
					"lw_detlin",	/* Recipe ID */
					datancom,
					raw,
					NULL);
	image_save_fits_hdrdump(fitres->plane[2], outname, fh_spec, BPP_DEFAULT);
	qfits_header_destroy(fh_spec);


	/* Save goodness of fit image */
	fh_spec = qfits_header_copy(fh);
	sprintf(outname, "%s_Q.fits", get_basename(name_o));
	e_comment(1, "saving image [%s]", outname);
	isaac_pro_fits(	fh_spec,
					outname,
					NULL,	/* Not a reduced file */
					NULL,	/* No reduction level applicable */
					isaac_imag_lw_detlin_coeff_Q,
					"OK",	/* Status */
					"lw_detlin",	/* Recipe ID */
					datancom,
					raw,
					NULL);
	framelist_del(raw) ;
	image_save_fits_hdrdump(fitres->plane[3], outname, fh_spec, BPP_DEFAULT);
	qfits_header_destroy(fh_spec);

	/* Delete reference header */
	qfits_header_destroy(fh);
	return 0 ;
}