/*---------------------------------------------------------------------------
   
   File name 	:	wavelength.c
   Author 		:	Y.Jung
   Created on	:	July 2000
   Description	:	ISAAC common function to handle with wl calibration

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

/*
	$Id: wavelength.c,v 1.8 2001/12/18 09:31:49 yjung Exp $
	$Author: yjung $
	$Date: 2001/12/18 09:31:49 $
	$Revision: 1.8 $
*/

/*---------------------------------------------------------------------------
   								Includes
 ---------------------------------------------------------------------------*/
#include <math.h>
#include <ctype.h>
#include <unistd.h>

#include "isaacp_lib.h"
#include "wavelength.h"

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

/*-------------------------------------------------------------------------*/
/**
  @name     spectro_jitter_wlcalib
  @memo     Compute the wavelength calibration for jitter recipes
  @param  	nb_coeff	nb of coeffs requested for the output polynomial
  @param	use_arc	Flag that specifie that arc table has to be used
  @param	arc_file	Arc table
  @param	framelist	Cleaned list of input frames
  @param	nbframes	nb of frames in the list
  @param	spec_length	image liength in spec direction (in pixels)
  @param	sky_frame	Frame used for sky calibration
  @param	filter_sky_flag	Flag to use or not filter
  @param	discard_lo	
  @param	discard_hi	regions to discard before calibration
  @param	discard_le	
  @param	discard_ri	
  @param	status_arc	status for arc calibration. MODIFIED
  @param	status_sky	status for sky calibration. MODIFIED
  @param	status		status for calibration. MODIFIED
  @return	dispersion polynomial coefficients 
  @doc
  The sky lines are assumed to be vertical
 */
/*--------------------------------------------------------------------------*/
double * spectro_jitter_wlcalib(
		int			nb_coeff,
		int			use_arc,
		char	*	arc_file,
		char	**	framelist,
		int			nbframes,
		int			spec_length,
		char	*	sky_frame,
		int			filter_sky_flag,
		int			discard_lo,
		int			discard_hi,
		int			discard_le,
		int			discard_ri,
		int		*	status_arc,
		int		*	status_sky,
		int		*	status)
{
    qfits_table	    *   arc_table ;
    double          *   disprel ;
    double          *   disprel2 ;
    double          *   filter ;
    image_t        	*   in ;
    image_t        	*   filtered ;
	int					coef ;
	double				wave_min,
						wave_max ;
	int					i ;
	double				null_val ;

	/* Initialize */
	*status = ALGO_FAILED ;
	disprel = NULL ;	
	null_val = 0.0 ;
	
	/* Test the required poly size	 */
	if (nb_coeff != 2) {
		e_error("Only 1st degree poly supported") ;
		return NULL ;
	}
	
    /* Try to calibrate using the arc wavelength cal. table if provided */
    if (use_arc) {
		e_comment(0, "Calibration using arc table") ;
        if ((arc_table=qfits_table_open(arc_file, 1)) == NULL) {
            e_warning("cannot open arc table") ;
            *status_arc = ALGO_FAILED ;
		} else {
            if ((disprel=(double*)qfits_query_column_data(arc_table, 
							3, (void*)&null_val)) == NULL) {
                e_warning("cannot query column of arc table") ;
                *status_arc = ALGO_FAILED ;
            } else {
                *status_arc = ALGO_OK ;
                *status = ALGO_OK ;
            }
            qfits_table_close(arc_table) ;
        }
	}

    /* Use the sky if the calibration is not done yet */
    if (*status != ALGO_OK) {
		e_comment(0, "Calibration using sky lines") ;
        /* First estimation using a physical model */
        if ((disprel=isaac_get_disprel_estimate_from_list(framelist,
                        								nbframes,
                        								spec_length,
                        								nb_coeff-1)) == NULL) {
            e_warning("cannot estimate the dispersion relation") ;
            *status_sky = ALGO_FAILED ;
            *status = ALGO_FAILED ;
			return NULL ;
        }
        /* Try to calibrate using the sky */
        if (*status_sky != ALGO_FAILED) {
			/* Find out wave_min and wave_max */
			wave_min = 0 ;
            for (i=0 ; i<nb_coeff ; i++) wave_min += disprel[i] ;
			wave_max = 0 ;
			coef     = 1 ;
			for (i=0 ; i<nb_coeff ; i++) {
				wave_max += disprel[i] * coef ;
				coef *= spec_length ;
			}

            /* Load the sky frame */
			if (sky_frame == NULL) {
				e_warning("No sky frame provided - use the physical mod.") ;
				*status_sky = ALGO_FAILED ;
				*status = ALGO_OK ;
				return disprel ;	
			}
			if ((in = image_load(sky_frame)) == NULL) {
                e_warning("cannot load the sky - use the physical model") ;
                *status_sky = ALGO_FAILED ;
                *status = ALGO_OK ;
                return disprel ;
            }
			
            /* Filter the sky frame */
           	if (filter_sky_flag) {
				filter = malloc(3*sizeof(double)) ;
            	filter[0] = 1.0 ;
            	filter[1] = 2.0 ;
            	filter[2] = 1.0 ;
            	filtered = image_filter3x1(in, filter) ;
            	free(filter) ;
			} else filtered = image_copy(in) ;
            image_del(in) ;

            if ((disprel2 = spectro_compute_disprel_from_table(filtered,
                                                   discard_lo,
                                                   discard_hi,
                                                   discard_le,
                                                   discard_ri,
                                                   0, /* remove thermal */
                                                   "oh",
                                                   wave_min,
                                                   wave_max)) == NULL) {
                e_warning("cannot compute the dispersion relation") ;
                e_warning("use the solution given by the physical model") ;
                *status_sky = ALGO_FAILED ;
                *status = ALGO_OK ;
            } else {
                free(disprel) ;
				disprel = disprel2 ;
                *status_sky = ALGO_OK ;
                *status = ALGO_OK ;
            }
            image_del(filtered) ;
        }
    }
    return disprel ;
}


/*-------------------------------------------------------------------------*/
/**
  @name     isaac_get_disprel_estimate
  @memo     from a physical model of the instrument, find out the
  wavelength range associated to a given instrument configuration.
  Information are fetched from the FITS header.
  @param    filename    input file name
  @param    nbpix       nb of pixels to calibrate
  @param    poly_deg    polynomial degree
  @return   newly allocated array of 3 polynomical coefficients, such
  as  wave = c[0] + c[1] * pix + c[2] * pix * pix
  @doc  returns NULL in case of error
 */
/*--------------------------------------------------------------------------*/
double *
isaac_get_disprel_estimate(char * filename, int nbpix, int poly_deg)
{
    double  *   c ;
    double      central_wavelength ;
    char    *   s ;
    char        objective[32];
    char        resolution[32];

    double3 *   plist ;
    double  *   coefs ;

    int         i ;

    /* Get various information from the FITS header */
    central_wavelength = isaac_get_central_wavelength(filename);
    if (central_wavelength < 0.0) {
        e_error("cannot get central wavelength from [%s]", filename);
        return NULL ;
    }

    s = isaac_get_objective(filename) ;
    if (s==NULL) {
        e_error("cannot get objective used from [%s]", filename);
        return NULL ;
    }
    strcpy(objective, strlwc(s));

    s = isaac_get_resolution(filename);
    if (s==NULL) {
        e_error("cannot get resolution used from [%s]", filename);
        return NULL ;
    }
    strcpy(resolution, s);

    /* c is an array of nbpix doubles with the wavelength for each pix. */
    if ((c=isaac_physical_model(central_wavelength,
                                objective,
                                resolution,
                                nbpix)) == NULL) {
            e_error("cannot compute the physical model calibration") ;
            return NULL ;
    }

    /* A polynomial fit is computed  */
    if (poly_deg<0) {
        e_error("invalid polynomial degree specified") ;
        free(c) ;
        return NULL ;
    }

    plist = double3_new(nbpix);
    for (i=0 ; i<nbpix ; i++) {
        plist->x[i] = (double)i ;
        plist->y[i] = c[i] ;
    }
    free(c) ;
    coefs=fit_1d_poly(poly_deg, plist, NULL);
    double3_del(plist) ;
    return coefs ;
}


/*-------------------------------------------------------------------------*/
/**
  @name     isaac_get_disprel_estimate_from_list
  @memo     the same as isaac_get_disprel_estimate, but on several
  frames (in case the keywords are not present in the first file, try
  with the second, and so on...)
  @param    lnames  list of files names
  @param    nbfiles number of files
  @param    nbpix   number of pixels to calibrate
  @param    poly_deg    polynomial degree
  @return   the same as isaac_get_disprel_estimate
 */
/*--------------------------------------------------------------------------*/
double * isaac_get_disprel_estimate_from_list(
        char    **  lnames,
        int         nbfiles,
        int         nbpix,
        int         poly_deg)
{
    double  *   c ;
    double      central_wavelength ;
    char    *   s ;
    char        objective[32];
    char        resolution[32];
    double3 *   plist ;
    double  *   coefs ;
    int         found ;
    int         i ;

	/* Initialize */
    found = 0 ;
    i = 0 ;
	central_wavelength = 0.0 ;
	
	if (nbfiles == 0) return NULL ;
	
    /* Get various information from the FITS header */
    while ((i<nbfiles) && (found == 0)) {
        found = 1 ;
        central_wavelength = isaac_get_central_wavelength(lnames[i]);
        if (central_wavelength < 0.0) {
            e_warning("cannot get central wavelength from [%s]", lnames[i]);
            found = 0 ;
        }

        s = isaac_get_objective(lnames[i]) ;
        if (s==NULL) {
            e_warning("cannot get objective used from [%s]", lnames[i]);
            found = 0 ;
        } else {
            strcpy(objective, strlwc(s));
        }

        s = isaac_get_resolution(lnames[i]);
        if (s==NULL) {
            e_warning("cannot get resolution used from [%s]", lnames[i]);
            found = 0 ;
        } else {
            strcpy(resolution, s);
        }

        i++ ;
    }

    if (found == 0) {
        e_error("cannot find grating info in header files") ;
        return NULL ;
    }

    /* c is an array of nbpix doubles with the wavelength for each pix. */
    if ((c=isaac_physical_model(central_wavelength,
                                objective,
                                resolution,
                                nbpix)) == NULL) {
            e_error("cannot compute the physical model calibration") ;
            return NULL ;
    }

    /* A polynomial fit is computed  */
    if (poly_deg<0) {
        e_error("invalid polynomial degree specified") ;
        free(c) ;
        return NULL ;
    }

    plist = double3_new(nbpix);
    for (i=0 ; i<nbpix ; i++) {
        plist->x[i] = (double)i ;
        plist->y[i] = c[i] ;
    }
    free(c) ;
    coefs=fit_1d_poly(poly_deg, plist, NULL);
    double3_del(plist) ;
    return coefs ;
}


/*-------------------------------------------------------------------------*/
/**
  @name     isaac_physical_model    
  @memo     computes a wavelength calibration thanks to a physical
  model and infos from the header   
  @param    lambda_c    central wavelength
  @param    objective   slit used
  @param    resolution  resolution (medium or low)
  @param    nbpix       number of pixels to calibrate
  @return   wavelength coefficients
 */
/*--------------------------------------------------------------------------*/
double *
isaac_physical_model(
        double      lambda_c,
        char    *   objective,
        char    *   resolution,
        int         nbpix)
{
    double      beam_diff ;
    double      focal_length,
                pixel_size,
                pupil;
    double      gr;
    double      gr_dir;
    double      isaacLRdir = ISAAC_LR_DIR ;
    double      isaacMRdir = ISAAC_MR_DIR ;

    double      isaacLRGrating = ISAAC_LR_GRATING ;
    double      isaacMRGrating = ISAAC_MR_GRATING ;
    int         order;
    double  *   disp_rel;
    double      x;
    double      angle_in;
    double      angle_out;
    double      wsl,
                wsu,
                wsm;
    int         i;

    double      a,
                b,
                det;


    /* This module determines the dispersion relation of ISAAC for the */
    /* different configurations of objectives, gratings, detectors, and  */
    /* central wavelengths */

    /* Now reading configuration, low- or medium-resolution,  */
    /* and central wavelength, to determine, order number, grating setting, */
    /* wavelength range. */


    /* assumed Optical configuration: */

    /* Focal lens of objectives for short- and long-wavelength objective: */

    /* (Values at 77 Kelvin). */
    /* S1=f/1.75 S2=f/3.25 L1=f/1.56 L2=f/4.77 L3=f/9.98 */
    /* Pupil size: 100 mm */
    /* Pixel size: 18.5 microns (SW), 27 microns (LW) */
    /* Gratings: low resolutions, 40 gr/mm, entering at about 5 degrees */
    /* medium resol.,  210 gr/mm, entering at about 23 degrees */
    /* Beam difference: 2.72 degrees */

    /* Test: L2, medium res., wave=2120 nm, pixel 27 microns,   */
    /* The module computes: 2057.42,  2182.45,  1.22103, order 2 */
    /* Lab measure: 1.28 A/pixel measured with pixels 30 microns. */

    /* Test: S1, low res., wave=1250 nm, pixel 18.5 microns */
    /* The module computes: 913.132,  1586.49,  6.57576, order 4 */

    /* Syntax: test <central_wave_nm> <low/medium> <s1/s2/l1/l2/l3> */

    focal_length = ISAAC_FOCAL_LENGTH_MM ;
    pixel_size   = ISAAC_PIXEL_SIZE_S ;
    beam_diff    = ISAAC_BEAM_DIFF ;
    pupil        = ISAAC_PUPIL_SIZE_MM ;

    disp_rel = malloc(nbpix * sizeof(double)) ;
    if (disp_rel == NULL) {
        e_error("memory allocation failure: aborting") ;
        return NULL ;
    }

    if (objective[0] == 's') {
        if (objective[1] == '1') focal_length = pupil*ISAAC_FLGTH_S1 ;
        if (objective[1] == '2') focal_length = pupil*ISAAC_FLGTH_S2 ;
        pixel_size = ISAAC_PIXEL_SIZE_S ;
    }

    if (objective[0] == 'l') {
        if (objective[1] == '1') focal_length = pupil*ISAAC_FLGTH_L1;
        if (objective[1] == '2') focal_length = pupil*ISAAC_FLGTH_L2;
        if (objective[1] == '3') focal_length = pupil*ISAAC_FLGTH_L3;
        pixel_size = ISAAC_PIXEL_SIZE_M ;
    }

    focal_length *= 1e-3; /* Convert mm to m. */

    /* Resolution to lower case letters */
	strcpy(resolution, strlwc(resolution)) ;
	
    /* Display configuration */

    e_comment(1, "configuration: ") ;
    if (resolution[0]=='l') e_comment(2, "configuration: low") ;
    if (resolution[0]=='m') e_comment(2, "configuration: medium") ;
    e_comment(2, "lambda_c     : %g", lambda_c) ;
    e_comment(2, "objective    : %s", objective) ;
    e_comment(2, "focal length : %g", focal_length) ;

    /* Will be eventually loaded from the IDF */
    if (resolution[0] == 'l') { gr = isaacLRGrating; gr_dir = isaacLRdir;}
	else if (resolution[0] == 'm') { gr = isaacMRGrating; gr_dir = isaacMRdir;}
	else {    
		e_error("wrong grating! %c", resolution[0]) ;
        free(disp_rel) ;
        return NULL ;
    }

    gr *= 1e-6; /* To have gr in groove/nm. */
    beam_diff *= M_PI/180.0; /* Convert beam_diff in radians. */


    /* 40 gr/mm --> 5. degrees input */
    /* 210 gr/mm -->  23. degrees  */

    /* convert lambda_c : A->nm */
    lambda_c /= 10 ;

    if ( lambda_c >= 890 && lambda_c < 8000) {
        if ( lambda_c >= 890 && lambda_c < 990 ) {
            order = 6;
        } else if ( lambda_c >= 990 && lambda_c < 1100 ) {
            order = 5;
        } else if ( lambda_c >= 1100 && lambda_c < 1400 ) {
            order = 4;
        } else if ( lambda_c >= 1400 && lambda_c < 1850 ) {
            order = 3;
        } else if ( lambda_c >= 1850 && lambda_c < 2500 ) {
            order = 2;
        } else if ( lambda_c >= 2500 && lambda_c < 5500 ) {
            order = 1;
        } else {
            order = 1;
        }
    } else {
        angle_in = ANGLE_IN_DEFAULT ;
        angle_out = ANGLE_OUT_DEFAULT ;
        order = (int)(0.5+(sin(angle_in)+sin(angle_out))/(lambda_c*gr));
    }

    e_comment(2, "order    : %d", order) ;


    /* The following is the solution of the set of equations : */
    /* (1)  sin(angle_in) + sin(angle_out) = order*gr*lambda_c  */
    /* (2)  angle_out - angle_in = beam_diff */
    /* for the two angles angle_in and angle_out. The system has two */
    /* solutions, only one is retained. */

    det = 2*sin(beam_diff);
    a   = order*gr*lambda_c*sin(beam_diff);
    b   = sqrt(order*gr*lambda_c*order*gr*lambda_c*
        (sin(beam_diff)*sin(beam_diff)+2*cos(beam_diff)-2)
        + 2*sin(beam_diff)*sin(beam_diff)*(1-cos(beam_diff)));

    angle_in = asin((a-b)/det);
    angle_out = angle_in + beam_diff;

    for (i=0; i<nbpix;i++) {
        x = (double)(i-(nbpix/2.)) * pixel_size * 1e-6; /* x is in meters */
        disp_rel[i] =
            (sin(angle_in)+sin(angle_out+atan(x/focal_length)))/(order*gr);

        /* Convert disprel from nm to angstroms */
        disp_rel[i] *= 10 ;
    }

    wsl = disp_rel[0];
    wsu = disp_rel[nbpix-1];
    wsm = (wsl+wsu)/2.;

    e_comment(2, "central wavelength computed: %g A", wsm) ;
    e_comment(2, "wave band: [%g %g] (angstroms)", wsl, wsu) ;
    e_comment(2, "average dispersion: %g A/pixel", (wsu-wsl)/(double)nbpix) ;

    return disp_rel ;
}