/*----------------------------------------------------------------------------
   
   File name    :   sp_extract.c
   Author       :	Y. Jung
   Created on   :	Feb. 2001
   Description  :	ISAAC spectrum extraction 

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

/*
	$Id: sp_extract.c,v 1.24 2002/03/04 14:13:19 yjung Exp $
	$Author: yjung $
	$Date: 2002/03/04 14:13:19 $
	$Revision: 1.24 $
 */

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

#include "eclipse.h"
#include "isaacp_lib.h"

/*---------------------------------------------------------------------------
                            Function prototypes
 ---------------------------------------------------------------------------*/

static int sp_extract_engine(char *, char *, int, int, int, int, int) ;
static int sp_extract_write_table(char *, int, pixelvalue **, char *) ;

/*----------------------------------------------------------------------------
                            Main code
 ---------------------------------------------------------------------------*/
int isaac_sp_extract_main(void * dict)
{
	dictionary	*	d ;
		
	int				spec_width ;
	int				sky_dist ;
	int				sky_width ;
	int				spec_pos ;
	int				no_graph ;

	char			argname[10] ;
	char    	*	name_i ;
    char    	*	name_o ;
    int     		nfiles ;
	
	int				errors ;
	int				i ;
	 
	d = (dictionary*)dict ;
   	/* Get options */
	spec_width = dictionary_getint(d, "arg.width", 10) ;
	sky_dist   = dictionary_getint(d, "arg.sky_dist", 10) ;
	sky_width  = dictionary_getint(d, "arg.sky_width", 30) ;
	spec_pos   = dictionary_getint(d, "arg.spec_pos", 512) ;
	no_graph   = dictionary_getint(d, "arg.nograph", 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") ;
		printf("%s\n", get_rootname(get_basename(name_i))) ;
		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 += sp_extract_engine(name_i,
                            		name_o,
									spec_width,
									sky_dist,
									sky_width,
									spec_pos,
									no_graph) ;
		free(name_o) ;
	}
	return errors ;
}


/*---------------------------------------------------------------------------
  							Function codes
 ---------------------------------------------------------------------------*/

/*-------------------------------------------------------------------------*/
/**
  @name		sp_extract_engine
  @memo		extract a spectrum at a given position
  @param	imagename	image name
  @param	tablename	output table name
  @param	spec_width	spectrum width
  @param	sky_dist	sky distance
  @param	sky_width	sky width
  @param	spec_pos	spectrum position
  @param	no_graph	flag to disable the gnuplot output
  @return	int -1 in error case
  @doc		the spectrum in the input image is an horizontal line.
  			spectrum_position is the distance in pix. from the bottom.
 */
/*--------------------------------------------------------------------------*/
static int sp_extract_engine(
        char    *   imagename,
        char    *   tablename,
        int         spec_width,
        int         sky_dist,
        int         sky_width,
        int         spec_pos,
		int			no_graph)
{
    image_t        	*   inimage;
    pixelvalue      *   array[2] ;
	double			*	plot_x,
					*	plot_y ;
	
	double			*	wavecal ;
    image_t        	*   extr_line ;
    int                 low_side, up_side ;
    int                 spec_length, slit_length ;
    int                 sky_pos[4] ;
    pixelvalue          sky_estim ;
	double				median_1 ;
	double				median_2 ;
	int					wave_flag ;
	char				outtablename[1024] ;
	
	int					i ;

	wave_flag = 1 ;
	
    /* Load the input image */
    if ((inimage = image_load(imagename)) == NULL) {
        e_error("in loading image [%s] - aborting", imagename) ;
		return -1 ;
    }
	
	/* Wavelength calibration */
	if ((wavecal=isaac_get_disprel_estimate(imagename,	
											inimage->lx,
											2)) == NULL) {
		e_warning("cannot estimate the dispersion model") ;
		wave_flag = 0 ;	
	}

    /* Verify the validity of spec_pos */
    spec_length = inimage->lx ;
    slit_length = inimage->ly ;
    if (spec_pos <= 1 || spec_pos >= inimage->ly) {
		e_error("invalid spectrum position: [%d] - aborting",spec_pos) ;
		image_del(inimage) ;
		free(wavecal) ;
		return -1 ;
	}

    /* Allocate the ouput arrays */
	array[0] = malloc(spec_length*sizeof(pixelvalue)) ;
	array[1] = malloc(spec_length*sizeof(pixelvalue)) ;

	if (sky_width == 0) {
		e_comment(1, "No sky background subtraction") ;
	}

    /* Extract the spectrum and get rid of the residual sky */
    for (i=0 ; i<spec_length ; i++) {
        /* Spectrum position */
        low_side = (int)(spec_pos - (spec_width-1)/2) ;
		up_side = low_side + spec_width ;
        if (low_side < 1) low_side = 1 ;
		if (up_side > slit_length) up_side = slit_length ;
		spec_width = up_side - low_side ;
		
        /* Positions of the residual sky */
        sky_pos[1] = (int)(spec_pos - sky_dist) ;
        sky_pos[0] = sky_pos[1] - sky_width ;
        sky_pos[2] = (int)(spec_pos + sky_dist) ;
        sky_pos[3] = sky_pos[2] + sky_width ;

		if ((extr_line = image_getvig(inimage,
												i+1,
												low_side,
												i+1,
												up_side)) == NULL) {
			e_error("error in line extraction - aborting") ;
			free(array[0]) ;
			free(array[1]) ;
			image_del(inimage) ;
			free(wavecal) ;
			return -1 ;
		}
		
		if ((sky_pos[0] < 1) && (sky_width !=0)) {
			median_1 = image_getmedian_vig(inimage,
												i+1,
												sky_pos[2],
												i+1,
												sky_pos[3]) ;
			sky_estim = median_1 ;
		} else if ((sky_pos[3] > slit_length) && (sky_width !=0)) {
			median_1 = image_getmedian_vig(inimage,
												i+1,
												sky_pos[0],
												i+1,
												sky_pos[1]) ;
			sky_estim = median_1 ;
		} else if (sky_width !=0) {
			median_1 = image_getmedian_vig(inimage,
												i+1,
												sky_pos[0],
												i+1,
												sky_pos[1]) ;
			median_2 = image_getmedian_vig(inimage,
												i+1,
												sky_pos[2],
												i+1,
												sky_pos[3]) ;
			sky_estim=(median_1 + median_2)/2 ;
		} else {
			sky_estim = 0 ;
		}

        array[1][i]=(pixelvalue)image_getsumpix(extr_line) ;
        array[1][i] = (pixelvalue)array[1][i] - spec_width*sky_estim ;
    	image_del(extr_line) ;
	}
	image_del(inimage) ;

	/* Fill array[0] */
	if (wave_flag) {
		for (i=0 ; i<spec_length ; i++) {
			array[0][i] = (pixelvalue)(wavecal[0] +
				wavecal[1]*(i+1) +
				wavecal[2]*(i+1)*(i+1)) ;
		}
		free(wavecal) ;
	} else {
		for (i=0 ; i<spec_length ; i++) {
			array[0][i] = (pixelvalue)(i+1) ;
		}
	}

	/* Write the output table */
	sprintf(outtablename, "%s.tfits", tablename) ;
	if (sp_extract_write_table(outtablename,
							spec_length,
							(pixelvalue**)array,
							imagename) == -1) {
		e_error("cannot write output FITS table") ;
	}
	
	if (no_graph == 0) {
		/* Display the plot  */
		plot_x = malloc(spec_length * sizeof(double));
		plot_y = malloc(spec_length * sizeof(double));
	
		for (i=0 ; i<spec_length ; i++) {
			plot_x[i] = array[0][i];
			plot_y[i] = array[1][i];
		}
		gnuplot_plot_once("Extracted spectrum",
						  "lines",
						  "wavelength",
						  "spectrum",
						  plot_x,
						  plot_y,
						  spec_length);
		free(plot_x) ;
		free(plot_y) ;
	}
	
    /* Free and return */
	free(array[0]) ;
	free(array[1]) ;
	return 0 ;
}


/*-------------------------------------------------------------------------*/
/**
  @name		sp_extract_write_table
  @memo		write the FITS output table
  @param	name		name of the table to write
  @param	nb_lines	nb of lines
  @param	array		pixelvalue array (2 col.)
  @param	imagename	name of the ref image used to transfer the
  						header
  @return	int (-1 in error case)	
 */
/*--------------------------------------------------------------------------*/
static int sp_extract_write_table(
		char		*	name,
		int				nb_lines,
		pixelvalue	**	array,
		char		*	imagename)
{
	qfits_header	*	fh ;
    qfits_table	    *   table ;
	qfits_col       *   col ;
	double			*	array_double[2] ;
	framelist		*	lnames ;
	
	int					i ;

	/* Write the output qfits_table table (informations) */
	table = qfits_table_new(name, QFITS_BINTABLE, -1, 2, nb_lines) ;
	col = table->col ;
	for (i=0 ; i<table->nc ; i++) {
		qfits_col_fill(col, 1, sizeof(double), TFITS_BIN_TYPE_D, " ", " ", " ",
                " ", 0, 0.0, 0, 1.0, i*sizeof(double)) ;
		col++ ;
	}
	col = table->col ;
	sprintf(col->tlabel, "Wavelength") ;
	col++ ;
	sprintf(col->tlabel, "Extracted") ;

	/* Convert pixelvalue array to double array */
	array_double[0] = malloc(nb_lines*sizeof(double)) ;
	array_double[1] = malloc(nb_lines*sizeof(double)) ;
	for (i=0 ; i<nb_lines ; i++) {
		array_double[0][i] = (double)array[0][i] ;
		array_double[1][i] = (double)array[1][i] ;
	}

	/* WRITE THE OUTPUT FILE */
	/* Read the input header */
	if ((fh = qfits_header_read(imagename)) == NULL) {
        e_error("in writing the output fits file") ;
        qfits_table_close(table) ;
        free(array_double[0]) ;
        free(array_double[1]) ;
		return -1 ;
    }

    /* Prepare it for table output */
    if (isaac_header_for_table(fh) == -1) {
        e_error("in writing the output fits file") ;
        qfits_header_destroy(fh) ;
        qfits_table_close(table) ;
        free(array_double[0]) ;
        free(array_double[1]) ;
        return -1 ;
    }

	/* Write the PRO keywords in the header */
    if (isaac_pro_fits(fh,
                        name,
                        "REDUCED",
                        NULL,
                        isaac_spec_extract,
                        "OK",
                        "isaacp_sp_extract",
                        1,
						NULL,
						NULL) == -1) {
        e_error("in writing PRO keywords in output file") ;
        qfits_header_destroy(fh) ;
        qfits_table_close(table) ;
        free(array_double[0]) ;
        free(array_double[1]) ;
        return -1 ;
    }
	
	/* Create lnames */
	lnames = framelist_new(1) ;
	lnames->name[0] = strdup(imagename) ;
	
	/* Write the HISTORY keywords with the input file names */
    if (isaac_add_files_history(fh, lnames) == -1) {
        e_warning("cannot write HISTORY keywords in out file") ;
    }
	framelist_del(lnames) ;

	/* Write the file on disk */
    if (qfits_save_table_hdrdump((void**)array_double, table, fh) == -1) {
        e_error("cannot write file: %s", name) ;
        qfits_header_destroy(fh) ;
        qfits_table_close(table) ;
        free(array_double[0]) ;
        free(array_double[1]) ;
        return -1 ;
    }
	qfits_table_close(table) ;
    qfits_header_destroy(fh) ;
    free(array_double[0]) ;
    free(array_double[1]) ;

	e_comment(0, "File [%s] produced", name) ;
	return 0 ;
}