/*-------------------------------------------------------------------------*/
/**
   @file	classif.c
   @author	Y. Jung
   @date	July 2000
   @version	$Revision: 1.6 $ 
   @brief	ISAAC common functions for frames classification
*/
/*--------------------------------------------------------------------------*/

/*
	$Id: classif.c,v 1.6 2002/03/04 15:33:31 yjung Exp $
	$Author: yjung $
	$Date: 2002/03/04 15:33:31 $
	$Revision: 1.6 $
*/

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

#include <math.h>
#include <ctype.h>
#include <unistd.h>

#include "isaacp_lib.h"

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

#define INLINESZ		1024

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

static int off_comp(double, double, double) ;
static int double_sort(const void * d1, const void * d2) ;
static int * spectro_jitter_classif_engine(double *, int, double, int *) ;

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

/*-------------------------------------------------------------------------*/
/**
  @brief	Get rid of non science frames in a frame list
  @param	frame_list	frame list
  @param	nframes		nb of frames in the input frame list
  @param	nb_clean_frames	nb of frames in output list
  @return	The new list
 */
/*--------------------------------------------------------------------------*/
char ** reject_non_science_frames(
		char	**	frame_list,
		int			nframes,
		int		*	nb_clean_frames)
{
	char	**	clean_frame_list ;
	char	*	value ;
	int			i, j ;	
	
	/* Reject non-science and non-standard frames - use header */
    /* Count the valid files */
    *nb_clean_frames = 0 ;
    for (i=0 ; i<nframes ; i++) {
        if ((value=isaac_get_dpr_type(frame_list[i])) == NULL) {
			e_warning("cannot read DPR TYPE") ;
		} else if ((!strcmp(value, "OBJECT")) || (!strcmp(value, "STD"))) 
            (*nb_clean_frames) ++ ;
    }

    /* Test if some files have been selected */
    if (*nb_clean_frames == 0) {
        e_comment(1, "No SCIENCE frames found") ;
        return NULL ;
    }

    /* Allocate the new frame list */
    clean_frame_list = malloc((*nb_clean_frames)*sizeof(char*)) ;
    j = 0 ;
    for (i=0 ; i<nframes ; i++) {
        if ((value=isaac_get_dpr_type(frame_list[i])) != NULL) 
			if ((!strcmp(value, "OBJECT")) || (!strcmp(value, "STD"))) {
				clean_frame_list[j] = malloc(FILENAMESZ) ;
				strcpy(clean_frame_list[j], frame_list[i]) ;
				j++ ;
			}
	}
	return clean_frame_list ;
}


/*-------------------------------------------------------------------------*/
/**
  @brief	Get rid of half-cycle frames in a frame list
  @param	frame_list	frame list
  @param	nframes		nb of frames in the input frame list
  @param	nb_clean_frames	nb of frames in output list
  @return	The new list
 */
/*--------------------------------------------------------------------------*/
char ** reject_hc_frames(
		char	**	frame_list,
		int			nframes,
		int		*	nb_clean_frames)
{
	char	**	clean_frame_list ;
	char	*	value ;
	int			i, j ;	

	/* Reject non-science and non-standard frames - use header */
    /* Count the valid files */
    *nb_clean_frames = 0 ;
    for (i=0 ; i<nframes ; i++) 
		if ((value=isaac_get_detector_frame_type(frame_list[i])) == NULL) {
			e_warning("cannot read DET FRAME TYPE") ;
		} else if (!strcmp(qfits_pretty_string(value), "INT")) 
			(*nb_clean_frames)++ ;

    /* Test if some files have been selected */
    if (*nb_clean_frames == 0) {
        e_comment(1, "No INT frames found") ;
        return NULL ;
    }

    /* Allocate the new frame list */
    clean_frame_list = malloc((*nb_clean_frames)*sizeof(char*)) ;
    j = 0 ;
    for (i=0 ; i<nframes ; i++) {
		if ((value=isaac_get_detector_frame_type(frame_list[i])) != NULL) 
			if (!strcmp(qfits_pretty_string(value), "INT")) {
				clean_frame_list[j] = malloc(FILENAMESZ) ;
				strcpy(clean_frame_list[j], frame_list[i]) ;
				j++ ;
			}
	}
	return clean_frame_list ;
}


/*-------------------------------------------------------------------------*/
/**
  @brief    Classification engine for spectroscopic jitter (SW + LW)
  @param	frame_list	input frames names. 
  @param 	nframes		Nb of frames.
  @param	offsets_source	OFFSOURCE_HEADER or OFFSOURCE_FILE
  @param	keyword_used_for_classif	keyword used for the classification. 
  			If it begins with "isaac", use isaac_get_cumoffsety().
  @param	offsets_sourcefile	ASCII file with the offsets list
  @param	offsets_total_list	array in which a nframes offs list is written.
  			MODIFIED.
  @param	nb_classified_cubes	nb of cubes. MODIFIED.
  @param	cube_id	int array that gives for each frame the corresp cube.
  			MODIFIED
  @param	status_classification	ALGO_FAILED, ALGO_SKIPPED or ALGO_OK.
  			MODIFIED.
  @return   Classified cubes
 */
/*--------------------------------------------------------------------------*/
int spectro_jitter_classification(
		char	**	frame_list,
		int			nframes,
		int			offsets_source,
		char	*	keyword_used_for_classif,	
		char	*	offsets_sourcefile,
		double	**	offsets_total_list,	
		int		*	nb_classified_cubes,
		int		**	cube_id,
		int		*	status_classification)
{
    double      *   sorted ;
    double          offset_min,
                    offset_max,
                    offset_thresh ;
    int             nb_offsets ;
    double          currentoff ;

    int             i ;

  	/* Get the offsets */
	if (offsets_source == OFFSOURCE_HEADER) {
        /* Read the offsets in the input files headers */
        *offsets_total_list=get_cumoffsety_from_fits_header(nframes,
                                        frame_list,
                                        keyword_used_for_classif) ;
        if (*offsets_total_list == NULL) {
            e_error("cannot get the offsets") ;
            *status_classification = ALGO_FAILED ;
            return -1 ;
        }
    } else if (offsets_source == OFFSOURCE_FILE) {
        e_comment(1, "Read offsets file") ;
        *offsets_total_list=get_offsets_from_file(offsets_sourcefile,
                                            nframes,
                                            frame_list) ;
        if (*offsets_total_list == NULL) {
            e_error("cannot get the offsets from file [%s]", offsets_sourcefile) ;
            *status_classification = ALGO_FAILED ;
            return -1 ;
        }
    } else {
        e_warning("bad specified offsets source") ;
        *status_classification = ALGO_FAILED ;
        return -1 ;
    }
	
    /* Find the number of different values taken by the offset */
    /* and the max and min values of the offsets */
    sorted = malloc(nframes*sizeof(double)) ;
    for (i=0 ; i<nframes ; i++) {
        sorted[i] = (*offsets_total_list)[i] ;
    }
    qsort(sorted, nframes, sizeof(double), double_sort) ;
    nb_offsets = 1 ;
    currentoff = sorted[0] ;
    for (i=1 ; i<nframes ; i++) {
        if (currentoff != sorted[i]) {
            nb_offsets++ ;
            currentoff = sorted[i] ;
        }
	}
    offset_min = sorted[0] ;
    offset_max = sorted[nframes-1] ;
    offset_thresh = (offset_min + offset_max) / 2.0 ;
    free(sorted) ;

    /* Separate the offset total list in offset lists */
    /* Write the parameters in the blackboard */
    if (nb_offsets < 2) {
        e_error("not enough different offsets: [%d]", nb_offsets) ;
        *status_classification = ALGO_FAILED ;
        return -1 ;
    } else {
        /* Classify the frames according their position in the list */

        /* Find number of cubes  */
        if ((*cube_id = spectro_jitter_classif_engine(*offsets_total_list,
											nframes,
											offset_thresh,
											nb_classified_cubes)) == NULL) {
            e_error("cannot classify the frames") ;
            *status_classification = ALGO_FAILED ;
            return -1 ;
		}
        e_comment(1, "Nb of combined images : %d\n", *nb_classified_cubes/2) ;
    }

    /* Return */
	*status_classification = ALGO_OK ;
	return 0 ;
}


/*-------------------------------------------------------------------------*/
/**
  @brief    Compare two doubles for qsort
  @param    d1 first double
  @param    d2 second double
  @return   int
 */
/*--------------------------------------------------------------------------*/
static int double_sort(const void * d1, const void * d2)
{   
    if (*(double*)d1 > *(double*)d2) return 1 ;
    else if (*(double*)d1 == *(double*)d2) return 0 ;
    else return -1 ;
}


/*-------------------------------------------------------------------------*/
/**
  @brief	match an offsets sequence with a defined pattern	
  @param	list	Offsets list
  @param	noffsets	nb of offsets in offset list
  @param	thresh	Offset threshold
  @param	nb_classified_cubes	Nb of cubes.
  @return	array of noffsets int whose values go from 0 to nb_classified_cubes-1
  Compares the offsets to a threshold and return a number of
  batches according to the offsets positions. A sequence <<>> gives 2
  batches (<< and >>).  A sequence >>><<<>><< gives 4 batches (>>>,
  <<<, >>, <<). X '>' (resp '<') has to be followed by X (at least)
  '<' (resp. '>').
 */
/*--------------------------------------------------------------------------*/
static int * spectro_jitter_classif_engine(
		double	*	list,
		int			noffsets,
		double		thresh,
		int		*	nb_classified_cubes)
{
	int	*	cube_id ;
	int     last_cube ;
	
    int     i,    /* Current frame id */
			j,    /* nb of frames in the first cube */
			k,    /* nb of frames in the second cube */
			l ;

	/* Initialize */
    *nb_classified_cubes = 0 ;
    i = 0 ;

    /* Allocate output array */
   	cube_id = malloc(noffsets*sizeof(int)) ;

	/* Initialize cube_id */
	for (i=0 ; i<noffsets ; i++) cube_id[i] = -1 ;	
	
	i = 0 ;
    while (i < noffsets) {
        j = 0 ;
        /* Count the number of successive '+' or '-' (j) */
        while ((i+j<noffsets) && (!off_comp(list[i], list[i+j], thresh))) j++ ;
        
		if (i+j >= noffsets) i = noffsets ;
        else {
            k = 0 ;
            /* Check if there are j '-' or '+' (k) */
            while ((i+j+k<noffsets) 
					&& (!off_comp(list[i+j], list[i+j+k], thresh)) 
                    && (k<j)) k++ ;
			last_cube = 1 ;
			if (i+j+k < noffsets) {
                for (l=i+j+k ; l<noffsets ; l++) {
                    if (off_comp(list[i+j], list[l], thresh)) {
                        last_cube = 0 ;
                        break ;
                    }
                }
            }

            if (last_cube == 0) {
				for (l=0 ; l<j ; l++) cube_id[i+l] = *nb_classified_cubes ;
				for (l=0 ; l<k ; l++) cube_id[i+j+l] = *nb_classified_cubes + 1 ;
                (*nb_classified_cubes) += 2 ;
                i += j+k ;
            } else {
                for (l=0 ; l<j ; l++) cube_id[i+l] = *nb_classified_cubes ;
				for (l=0 ; l<noffsets - (i+j) ; l++) 
					cube_id[i+j+l] = *nb_classified_cubes + 1 ;
                (*nb_classified_cubes) +=2 ;
                i = noffsets ;
            }
        }
    }
	
	/* Nb of cubes found should be even */
    if (*nb_classified_cubes % 2) {
        free(cube_id) ;
        *nb_classified_cubes = 0 ;
        return NULL ;
    }
    return cube_id ;
}

/*-------------------------------------------------------------------------*/
/**
  @brief    compares two doubles in regard to a given threshold
  @param    off1 first double
  @param    off2 second double
  @param    thresh  threshold
  @return   1 if threshold is between offsets, 0 otherwise
 */
/*--------------------------------------------------------------------------*/
static int off_comp(double off1, double off2, double thresh)
{
    if (((off1>thresh) && (off2<thresh)) || ((off1<thresh) && (off2>thresh))) 
		return 1 ;
    else return 0 ;
}


/*-------------------------------------------------------------------------*/
/**
  @brief    separate a frame list into nb_cubes new cubes
  @param    in  frame list 
  @param	nbframes number of frames in the list
  @param    t_switch    switch table. values from 0 to nbcubes-1
  @param    nb_cubes    nb of cubes (ie nb of values taken by t_switch)
  @return   newly allocated list of nb_cubes 
  destroys the input cube after use
 */
/*--------------------------------------------------------------------------*/
cube_t ** spectro_jitter_framelist_to_cubes(
		char	**	in,
		int			nframes,
		int		*	t_switch,
		int			nb_cubes)
{
    cube_t **  classified ;
	cube_t	*	incube ;
    int		*	nb_offsets ;
	
	int         i, j ;
    int     *   ind ;
 
	/* Allocate nb_offsets */
	nb_offsets = calloc(nb_cubes, sizeof(int)) ;
	
	/* Find for each cube how many planes it has */
	for (i=0 ; i<nb_cubes ; i++) 
		for (j=0 ; j<nframes ; j++) 
			if (t_switch[j] == i) nb_offsets[i] ++ ;	
	
	/* Load framelist in incube */
	if ((incube = cube_load_strings(in, nframes)) == NULL) {
		e_error("cannot load input frames") ;
		return NULL ;
	}
	
	/* Allocate output array of cubes */
    classified = malloc(nb_cubes * sizeof(cube_t*)) ;
    for (i=0 ; i<nb_cubes ; i++) {
        classified[i] = cube_new(incube->lx, incube->ly, nb_offsets[i]);
    }
	free(nb_offsets) ;

    /*
     * Distribute planes in classified output cubes according to
     * t_switch table.
     */

    /* Allocate ind */
    ind = calloc(nb_cubes, sizeof(int)) ;
    for (i=0 ; i<nframes ; i++) {
        if (t_switch[i] == -1) {
            image_del(incube->plane[i]) ;
            incube->plane[i] = NULL ;
        } else {
            classified[t_switch[i]]->plane[ind[t_switch[i]]] = incube->plane[i] ;
            incube->plane[i] = NULL ;
            ind[t_switch[i]]++ ;
        }
    }
    free(ind) ;
    cube_del_shallow(incube);
    return classified ;
}


/*-------------------------------------------------------------------------*/
/**
  @brief	Get an offset from an offsets list
  @param	list	offsets list
  @param	nboff	nb of offsets in the list
  @param	cube_id	array that identifies cubes
  @param	cube	specifies the cube
  @param	iter	specifies which offset of the cube
  @return	the found offset
 */
/*--------------------------------------------------------------------------*/
double get_offset_in_list(
		double	*	list,
		int			noff,
		int		*	cube_id,
		int			cube,
		int			iter)
{
	int		curr_frame ;
	double	offset ;
	int		i ;

	offset = 0.0 ;
	curr_frame = 0 ;
	
	for (i=0 ; i<noff ; i++) {
		if (cube_id[i] == cube) {
			if (curr_frame == iter) return list[i] ;
			else curr_frame++ ;
		}
	}
	e_error("cube or image specified is out of range") ;
	return 0.0 ;
}


/*-------------------------------------------------------------------------*/
/**
  @brief	create a list of HC frames from a list of frames
  @param	list	list of frames
  @param	nb_frames	nb of frames in list
  @param	nb_hc		nb of hc frames 
  @param	type		1 for HC1, 2 for HC2
  @return	list of hc frames
 */
/*--------------------------------------------------------------------------*/
char ** get_hc_frames(
		char	**	list,
		int			nb_frames,
		int		*	nb_hc,
		int			type)
{
	char	**	out_list ;
	int			frame_type ;
	int			cur_frame_type ;
	
	int			i, j ;
	
	/* Initialize */
	*nb_hc = 0 ;

	/* HC1 or HC2 ? */
	if (type == 1) frame_type = ISAAC_HALF_CYCLE_1 ;
	else if (type == 2) frame_type = ISAAC_HALF_CYCLE_2 ;
	else return NULL ;

	/* Count the valid files */
	for (i=0 ; i<nb_frames ; i++) {
		cur_frame_type = isaac_get_frame_type(list[i]) ;
		if (cur_frame_type == frame_type) (*nb_hc) ++ ;
	}

    /* Test if some files have been selected */
    if (*nb_hc == 0) {  
        e_comment(1, "No HC frames found") ;
        return NULL ;
    }
	
    /* Allocate the new frame list */
    out_list = malloc((*nb_hc)*sizeof(char*)) ;
	j = 0 ;
    for (i=0 ; i<nb_frames ; i++) {
		cur_frame_type = isaac_get_frame_type(list[i]) ;
		if (cur_frame_type == frame_type) {
			out_list[j] = malloc(FILENAMESZ) ;
			strcpy(out_list[j], list[i]) ;
			j++ ;
		}
	}
	return out_list ;
}


/*-------------------------------------------------------------------------*/
/**
  @brief    get an offset list from a list of FITS headers
  @param    nframes number of frames
  @param    filenames   list of frames
  @param    keyword keyword to look for ("isaac-default" for ISAAC std)
  @return   a list of offsets
  ISAAC specific
 */
/*--------------------------------------------------------------------------*/
double * get_cumoffsety_from_fits_header(
        int         nframes,
        char    **  filenames,
        char    *   keyword)
{
    int			i ;
    double	*	ret_offs ;
    char	*   frame_name ;
    char	*   val ;
    double		dval ;
    int			use_isaac_dy ;
    
    use_isaac_dy = 0 ;
    if (!strncmp(keyword, "isaac", 5)) use_isaac_dy = 1 ;

    ret_offs = (double*)malloc(nframes*sizeof(double)) ;
    
    for (i=0 ; i<nframes ; i++) {
        frame_name = filenames[i] ;
        if (use_isaac_dy) {
            val = isaac_get_cumoffsety(frame_name);
        } else {
            val = qfits_query_hdr(frame_name, keyword);
        }
        if (val==NULL) {
            e_error("getting Y offset from file [%s]", frame_name);
            free(ret_offs) ;
            return NULL ;
        }
        sscanf(val, "%lg", &dval);
        ret_offs[i] = dval ;
    }
    return ret_offs ;
}


/*-------------------------------------------------------------------------*/
/**
  @memo     get an offset list from a ascii file 
  @param    offsets_sourcefile  ASCII file with the offsets list
  @param    nframes number of input frames
  @param    filelist    list of input frame names
  @return   a list of offsets
  ISAAC specific
 */
/*--------------------------------------------------------------------------*/
double * get_offsets_from_file(
		char    *   offsets_sourcefile,
        int         nframes,
        char    **  filelist)
{
    double          *   ret_offs ;
    double              dval ;
    FILE            *   file ;

    int                 i ;

    /* Test if ASCII file exists */
    if (!file_exists(offsets_sourcefile)) {
        e_error("offset file [%s] not found", offsets_sourcefile) ;
        return NULL ;
    }

    /* Open the ascii file */
    if ((file = fopen(offsets_sourcefile, "r")) == NULL) {
        e_error("cannot open offsets file") ;
        return NULL ;
    }

    /* Allocate the output array */
    ret_offs = (double*)malloc(nframes*sizeof(double)) ;

    for (i=0 ; i<nframes ; i++) {
        fscanf(file, "%lg\n", &dval) ;
        e_comment(1, "File [%s] -> offset [%lg]", get_basename(filelist[i]), 
				dval) ;
        ret_offs[i] = dval ;
    }
    return ret_offs ;
}

/*-------------------------------------------------------------------------*/
/**
  @brief    compare the grating, the slit and the wl in input files
            headers 
  @param    file1   first file
  @param    file2   second file
  @return   -1 in error case, 0 for different settings, 1 for equal
            settings
 */
/*--------------------------------------------------------------------------*/
int compare_settings(
        char    *   file1,
        char    *   file2)
{
    int         comparison ;

    double      wl1,
                wl2 ;
    char        grat_name1[FILENAMESZ] ;
    char        grat_name2[FILENAMESZ] ;
    char        opti_id1[FILENAMESZ] ;
    char        opti_id2[FILENAMESZ] ;
        
    char    *   s ;
        
    comparison = 1 ;    
        
    /* Get the slit used */
    s = isaac_get_optical_id(file1) ;
    if (s==NULL) {  
        e_error("cannot get optical id from [%s]", file1) ;
        return -1 ; 
    }               
    strcpy(opti_id1, s) ;
    s = isaac_get_optical_id(file2) ;
    if (s==NULL) {  
        e_error("cannot get optical id from [%s]", file2) ;
        return -1 ;
    }
    strcpy(opti_id2, s) ;

    if (strcmp(opti_id1, opti_id2) != 0) comparison = 0 ;

    /* Get the grating name */
    if (comparison == 1) {
        s = isaac_get_resolution(file1) ;
        if (s==NULL) {
            e_error("cannot get resolution from [%s]", file1) ;
            return -1 ;
        }
        strcpy(grat_name1, s) ;
        s = isaac_get_resolution(file2) ;
        if (s==NULL) {
            e_error("cannot get resolution from [%s]", file2) ;
            return -1 ;
        }
        strcpy(grat_name2, s) ;

        if (strcmp(grat_name1, grat_name2) != 0) comparison = 0 ;
    }

    /* Compare the central wavelength */
    if (comparison == 1) {
        wl1 = isaac_get_central_wavelength(file1) ;
        if (wl1 == -1.00) {
            e_error("cannot get central wavelength from [%s]", file1) ;
            return -1 ;
        }
        wl2 = isaac_get_central_wavelength(file2) ;
        if (wl2 == -1.00) {
            e_error("cannot get entral wavelength from [%s]", file2) ;
            return -1 ;
        }
		if (fabs(wl1 - wl2) > 1e-4) comparison = 0 ;
    }

    return comparison ;
}


/*-------------------------------------------------------------------------*/
/**
  @brief    identify the activated lamps in a frames list   
  @param    lnames  list of frames
  @return   array: array[i]: 0 if no activated lamps in lnames[i], 1
  if xenon lamp is activated, 2 if argon lamp is activated, 3 if both
  lamps are activated
 */
/*--------------------------------------------------------------------------*/
int * find_activated_lamps(framelist * lnames)
{
    int     *   lamps ;

    int         i ;

    /* Allocate the output array */
    lamps = malloc(lnames->n*sizeof(int)) ;

    for (i=0 ; i<lnames->n ; i++) {
        if (isaac_is_xenon_lamp_active(lnames->name[i]) == 1) {
            if (isaac_is_argon_lamp_active(lnames->name[i]) == 1) {
                lamps[i] = 3 ;
            } else if (isaac_is_argon_lamp_active(lnames->name[i]) == 0) {
                lamps[i] = 1 ;
            } else {
                e_error("cannot check if argon lamp is activated") ;
                free(lamps) ;
                return NULL ;
            }
        } else {
            if (isaac_is_argon_lamp_active(lnames->name[i]) == 1) {
                lamps[i] = 2 ;
            } else if (isaac_is_argon_lamp_active(lnames->name[i]) == 0) {
                lamps[i] = 0 ;
            } else {
                e_error("cannot check if xenon lamp is activated") ;
                free(lamps) ;
                return NULL ;
            }
        }
    }

    /* Free and return */
    return lamps ;
}


/*-------------------------------------------------------------------------*/
/**
  @brief    for each setting, look if it is a medium or low resol
            setting
  @param    lnames  list of file names
  @param	nsettings	number of settings			
  @return   array of nsettings chars : l or m for low or medium resolution
 */
/*--------------------------------------------------------------------------*/
char * get_resolutions(
        framelist	*  	lnames,
        int         	nsettings)
{
    char    *   resolution ;

    char        reso[81] ;
    char    *   s ;
    int         i,
                j ;

    /* Allocate resolution table */
    resolution = malloc(nsettings*sizeof(char)) ;

    for (i=0 ; i<nsettings ; i++) {
        j=0 ;
        while (j<lnames->n) {
            if (lnames->label[j] == i) {
                s = isaac_get_resolution(lnames->name[j]) ;
                if (s==NULL) {
                    e_error("cannot get resolution from [%s]", lnames->name[j]);
                    free(resolution) ;
                    return NULL ;
                }
                strcpy(reso, strlwc(s)) ;
                if (reso[0] == 'l') {
                    resolution[i] = 'l' ;
                    j = lnames->n ;
                } else if (reso[0] == 'm') {
                    resolution[i] = 'm' ;
                    j = lnames->n ;
                } else {
                    e_warning("unexpected resolution read in [%s]: %c",
                            lnames->name[j], reso[0]) ;
                    resolution[i] = (char)reso[0] ;
                    j = lnames->n ;
                }
            } else j++ ;
        }
    }

    /* Free and return */
    return resolution ;
}


/*-------------------------------------------------------------------------*/
/**
  @brief    find out in an ASCII file which file name corresponds to a
            master dark frame. The order number specifies, in case
            there are several files corresponding to the definition,
            which one should be returned. First file is order 1,
            second file is order 2, etc.
  @param    filename    name of an ASCII file
  @param    id          id
  @param    order       order number
  @return   pointer to statically allocated char string
  This is highly DFS dependant
	Example  :   If the input ASCII file contains:
	--- begin 'input' ---
    /data/flat1.fits    MASTER_FLAT
    /data/flat2.fits    MASTER_FLAT
    /data/dark1.fits    MASTER_DARK
    /data/dark2.fits    MASTER_DARK
   --- end 'input' -----
    
   Calling:
    
   isaac_get_calibframe_name("input", "MASTER_FLAT", 1)
   returns: "/data/flat1.fits"
   isaac_get_calibframe_name("input", "MASTER_FLAT", 2)
   returns: "/data/flat2.fits"
   etc.
 */
/*--------------------------------------------------------------------------*/
char * isaac_get_calibframe_name(char * filename, char * id, int order)
{
    FILE    *   in ;
    char        line[INLINESZ];
    char        do_id[INLINESZ];
    static char match[INLINESZ];
    int         found ;

    if (order<=0) order=1 ;

    /* Open the file for reading */
    if ((in=fopen(filename, "r"))==NULL) {
        e_error("cannot open file: %s", filename);
        return NULL ;
    }

    found=0 ;
    memset(match, 0, INLINESZ);
    while (fgets(line, INLINESZ-1, in)!=NULL) {
        sscanf(line, "%s %s", match, do_id);
        if (!strcmp(do_id, id)) {
            found++;
            if (found==order) break ;
        }
    }
    fclose(in);
    if (found==order) return match ;
    return NULL ;
}


/*-------------------------------------------------------------------------*/
/**
  @brief    for a list of FITS frames, extract the keyword value for
            each frame and return as a list of doubles.
  @param    fnames  list of file names  
  @param    nfiles  number of files
  @param    keyword a FITS keyword
  @return   list of doubles
 */
/*--------------------------------------------------------------------------*/
double * get_detconfig(char ** fnames, int nfiles, char * keyword)
{
    double      *   list ;
    char        *   val ;
    int             i ;

    list = malloc(nfiles * sizeof(double)) ;
    for (i=0 ; i<nfiles ; i++) {
        compute_status("getting FITS info from headers", i, nfiles, 1) ;
        val = qfits_query_hdr(fnames[i], keyword);
        if (val==NULL) {
            e_error("in file [%s] to get value for [%s]", fnames[i], keyword);
            free(list) ;
            return NULL ;
        }
        val = qfits_pretty_string(val);
        if (sscanf(val, "%lg", &list[i])!=1) {
            e_error("getting [%s] from file [%s]", keyword, fnames[i]) ;
            free(list);
            return NULL;
        }
    }
    return list ;
}


/*-------------------------------------------------------------------------*/
/**
  @brief	Filter half-cycle frames out of a frame list.
  @param	flist	Framelist object to purge.
  @return	int 0 if Ok, -1 otherwise.
  This function processes a framelist object to remove any half-cycle
  frame. These frames are only expected in long-wavelength mode and they
  can be detected in two ways. Either the frame type in the framelist
  contains the word 'half' (case-insensitive), or the frame name
  corresponds to a FITS file that is referenced as a half-cycle frame.
 */
/*--------------------------------------------------------------------------*/
int isaac_lw_filter_halfcycle(framelist ** flist)
{
	framelist	*	purged ;
	int				i, j ;
	int				nval ;
	int			*	frame_ok ;
	char		*	str ;

	if (flist==NULL || (*flist)==NULL) return -1 ;

	frame_ok = malloc((*flist)->n * sizeof(int));
	nval = (*flist)->n ;
	/* Try rejecting frames based on frame type */
	if ((*flist)->type!=NULL) {
		for (i=0 ; i<(*flist)->n ; i++) {
			frame_ok[i] = 1 ;
			if ((*flist)->type[i]!=NULL) {
				if (strstr((*flist)->type[i], "half")!=NULL) {
					frame_ok[i]=0 ;
					nval -- ;
				}
			}
		}
	} else {
		/* Try rejecting frames based on detector ID keyword */
		for (i=0 ; i<(*flist)->n ; i++) {
			frame_ok[i] = 1 ;
			str = isaac_get_detector_frame_type((*flist)->name[i]);
			if (str!=NULL) {
				if (!strncmp(str, "HALF", 4)) {
					frame_ok[i]=0 ;
					nval -- ;
				}
			}
		}
	}

	/* If all frames Ok, return now */
	if (nval==(*flist)->n) {
		free(frame_ok);
		return 0 ;
	}

	/* Purge framelist if needed */
	purged = calloc(1, sizeof(framelist));
	purged->filename = strdup((*flist)->filename);
	purged->n = nval ;
	purged->name = calloc(nval, sizeof(char*));
	if ((*flist)->type!=NULL) {
		purged->type = calloc(nval, sizeof(char*));
	} else {
		purged->type = NULL ;
	}

	j=0 ;
	for (i=0 ; i<(*flist)->n ; i++) {
		if (frame_ok[i]) {
			purged->name[j] = strdup((*flist)->name[i]);
			if (purged->type!=NULL) {
				if ((*flist)->type[i]!=NULL) {
					purged->type[j] = strdup((*flist)->type[i]);
				}
			}
			j++ ;
		}
	}
	free(frame_ok);
	framelist_del(*flist);
	(*flist) = purged ;
	return 0 ;
}


/*-------------------------------------------------------------------------*/
/**
  @brief	Classifies a charmatrix of files according to a token retrieval 
  @param  in		files to classify.
  @param  colno	the column number to use as an alternative to the token retrieval function
  @param	nb_of_classes	number of classes found
  @param	funk_get_tok	the token retrieval function
  @return array of charmatrixes classified by token, NULL if error
  This function is typically used to classify a group of frames in a
  charmatrix. 
  Either the token has been provided in the framelist text
  file as column colno (!=0) or it has to be retrieved by funk_get_tok which
  takes a filename of the framelist as a parameter. A typical example
  is one of the functions isaac_get_XXX in dicbkeys.h 
  The function returns an array of charmatrixes each having the unique
  token in colno (if supplied) or appended as the last column (i.e. in->lx).
 */
/*--------------------------------------------------------------------------*/
charmatrix ** classify_charmatrix_by_token(
		charmatrix 	* 	in, 
		int				colno,
		int			*	nb_of_classes,
		char		*	(*funk_get_tok)(char*))
{
	charmatrix 	**	out=NULL,
				**	c=NULL;
	char		*	currtok=NULL,
				**	toktab=NULL;	
	int				use_input_col,
					found_class;
	int			*	classcount=NULL;				
	int				i,j,k,index_newcm;
	
	/* Initialize */
	i = 0 ;

	if (!in){
		e_error("classify_charmatrix_by_token: in NULL");
		return NULL;
	}
	
	if (colno <= 0)
		use_input_col = 0;
	else 
		use_input_col = 1;
	
	if (colno >= in->lx){
		e_warning( "classify_charmatrix_by_token: colno=%d lx=%d",colno,in->lx);
		use_input_col = 0;
	}
	if (!in->ly){
		e_error("classify_charmatrix_by_token: no files");
		return NULL;
	}
	e_comment(0,"classify_charmatrix_by_token: lx=%d ly=%d colno=%d",in->lx,
			in->ly,colno);
	toktab = calloc(in->ly,sizeof(char*));
	classcount = calloc(in->ly,sizeof(int));
/* Count nb of tokens to classify */	
	*nb_of_classes = 0;
	for (j=0; j<in->ly; j++){
		if (use_input_col)
			currtok = charmatrix_elem(in,colno,j);
		else 
			currtok=(*funk_get_tok)(charmatrix_elem(in,0,j));	
		if (currtok){
			found_class = 0;
			for (i=0;i<*nb_of_classes;i++){
				if (!strcmp(currtok,toktab[i])){
					found_class=1;
					break;
				}
			}
		}
		else{
			e_warning("classify_charmatrix_by_token: Token not found in file %s",
					charmatrix_elem(in,0,j));
			found_class = 1;
		}
		if (!found_class){
			if (!currtok)
				currtok = "UNKNOWN";
			toktab[(*nb_of_classes)++] = strdup(currtok);
		}
		else{
			classcount[i]++; /* one more occurrence for this class */
		}
	}

/* create new charmatrix */		
	out = calloc(*nb_of_classes,sizeof(charmatrix*));
	if (debug_active()>1)
		e_comment(0,"classify_charmatrix_by_token: found %d different tokens",*nb_of_classes);
	for (j=0;j<*nb_of_classes;j++){
		e_comment(0,"%s %d occurences",toktab[j],classcount[j]);
		c = &out[j];
		*c = calloc(1,sizeof(charmatrix));
		(*c)->lx = in->lx;
		if (!use_input_col)
			(*c)->lx++;
		(*c)->ly = classcount[j];
		(*c)->c = calloc((*c)->lx*(*c)->ly,sizeof(char*));
		classcount[j]=0;
	}
/* fill in new charmatrix */	
	for (j=0; j<in->ly; j++){
		if (use_input_col)
			currtok = charmatrix_elem(in,colno,j);
		else 
			currtok=(*funk_get_tok)(charmatrix_elem(in,0,j));	
		if (!currtok)
			currtok = "UNKNOWN";
		found_class = 0;
		for (i=0;i<*nb_of_classes;i++){
			if (!strcmp(currtok,toktab[i])){
				found_class=1;
				break;
			}
		}
		if (found_class){
			index_newcm = classcount[i]++; 
			c = &out[i];
			for (k=0;k<in->lx;k++){
				(*c)->c[index_newcm*(*c)->lx+k] = strdup(charmatrix_elem(in,k,j));	
			}
			if (!use_input_col){
				k = (*c)->lx-1;
				(*c)->c[index_newcm*(*c)->lx+k] = strdup(currtok);
			}
		}
		else{
			e_error("classify_charmatrix_by_token: could not find token %s in table",currtok);
		}
	}
	if (debug_active()>2){
		for (i=0;i<*nb_of_classes;i++){
			e_comment(0,"class %s",out[i]->c[out[i]->lx-1]);
			for (j=0;j<out[i]->ly;j++){
				e_comment(0,"file %s",out[i]->c[j*out[i]->lx]);
			}
		}
	}
	free(classcount);
	for (j=0;j<*nb_of_classes;j++)
		free(toktab[j]);
	free(toktab);
	return out;
}


/*-------------------------------------------------------------------------*/
/**
  @brief  Classifies a list of files according to a token retrieval function. 
  @param  in				files to classify.
  @param	nfiles			number of files in input
  @param	funk_get_tok	the token retrieval function
  @param	toktab			token table
  @param	classcount		the number of files for each token
  @param	nb_of_classes	the number of different tokens found.
  @return array of file lists classified by token, NULL if error
  This function is typically used to classify a group of frames in a
  filelist (char**). 
  The token has to be retrieved by funk_get_tok which
  takes a filename of the framelist as a parameter. A typical example
  is one of the functions isaac_get_XXX in dicbkeys.h 
  The function returns an array of file lists char**.
  The token table is similar to the file list containing the found tokens. 
 */
/*--------------------------------------------------------------------------*/
char *** classify_charlist_by_token(
		char	** 	in, 
		int			nfiles,
		char	***	toktab,
		int		*	nb_of_classes,
		int		**	classcount,				
		char	*	(*funk_get_tok)(char*))
{
	char	 	***	out=NULL,
				***	c=NULL;
	char		*	currtok=NULL;
	int				found_class;
	int				i,j,index_newcm;

	/* Initialize */
	i = 0 ;

	if (!in){
		e_error("classify_charlist_by_token: in NULL");
		return NULL;
	}
	
	if (!nfiles){
		e_error("classify_charlist_by_token: no files");
		return NULL;
	}
	(*toktab) = calloc(nfiles,sizeof(char*));
	(*classcount) = calloc(nfiles,sizeof(int));
/* Count nb of tokens to classify */	
	*nb_of_classes = 0;
	for (j=0; j<nfiles; j++){
		currtok=(*funk_get_tok)(in[j]);	
		if (currtok){
			found_class = 0;
			for (i=0;i<*nb_of_classes;i++){
				if (!strcmp(currtok,(*toktab)[i])){
					found_class=1;
					break;
				}
			}
		}
		else{
			e_warning("classify_charlist_by_token: Token not found in file %s", in[j]);
			found_class = 1;
		}
		if (!found_class){
			if (!currtok)
				currtok = "UNKNOWN";
			(*toktab)[(*nb_of_classes)++] = strdup(currtok);
		}
		else{
			(*classcount)[i]++; /* one more occurrence for this class */
		}
	}

/* create new list of chars*/		
	out = calloc(*nb_of_classes,sizeof(char**));
	if (debug_active()>1)
		e_comment(0,"classify_charlist_by_token: found %d different tokens in %d files",
			*nb_of_classes,nfiles);
	for (j=0;j<*nb_of_classes;j++){
		out[j] = calloc((*classcount)[j]+1,sizeof(char*));
		(*classcount)[j]=0;
	}
/* fill in new lists of chars */	
	for (j=0; j<nfiles; j++){
		currtok=(*funk_get_tok)(in[j]);	
		if (!currtok)
			currtok = "UNKNOWN";
		found_class = 0;
		for (i=0;i<*nb_of_classes;i++){
			if (!strcmp(currtok,(*toktab)[i])){
				found_class=1;
				break;
			}
		}
		if (found_class){
			index_newcm = (*classcount)[i]++; 
			c = &(out[i]);
			(*c)[index_newcm] = strdup(in[j]);	
		}
		else{
			e_error("classify_charlist_by_token: could not find token %s in table",currtok);
		}
	}
	if (debug_active()>2){
		for (i=0;i<*nb_of_classes;i++){
			e_comment(0,"class %s",(*toktab)[i]);
		}
	}
	return out;
}

	
/*-------------------------------------------------------------------------*/
/*
   @brief  Classifies a framelist according to a token retrieval function. 
   @param  in				framelist to classify.
   @param	funk_get_tok	the token retrieval function
   @param	toktab			token table
   @param	nb_of_classes	the number of different tokens found.
   @return array of framelists classified by token, NULL if error
   This function is typically used to classify a group of frames in a frameist. 
   The token has to be retrieved by funk_get_tok which
   takes a filename of the framelist as a parameter. A typical example
   is one of the functions isaac_get_XXX in dicbkeys.h 
   The function returns an array of framelists.
   The token table is a list containing the found tokens. 
 */
/*--------------------------------------------------------------------------*/
framelist * classify_framelist_by_token(
		framelist* 	in, 
		char	***	toktab,
		int		*	nb_of_classes,
		char	*	(*funk_get_tok)(char*))
{
	framelist	*	out=NULL,
				*	c=NULL;
	char		*	currtok=NULL;
	int				found_class;
	int				i,j,index_newcm;
	int			*	classcount;				

	/* Initialize */
	i = 0 ;	

	if (!in){
		e_error("classify_framelist_by_token: in NULL");
		return NULL;
	}
	
	if (!in->n){
		e_error("classify_framelist_by_token: no files");
		return NULL;
	}
	(*toktab) = calloc(in->n,sizeof(char*));
	classcount = calloc(in->n,sizeof(int));
/* Count nb of tokens to classify */	
	*nb_of_classes = 0;
	for (j=0; j<in->n; j++){
		currtok=(*funk_get_tok)(in->name[j]);	
		if (currtok){
			found_class = 0;
			for (i=0;i<*nb_of_classes;i++){
				if (!strcmp(currtok,(*toktab)[i])){
					found_class=1;
					break;
				}
			}
		}
		else{
			e_warning("classify_framelist_by_token: Token not found in file %s", 
					in->name[j]);
			found_class = 1;
		}
		if (!found_class){
			if (!currtok)
				currtok = "UNKNOWN";
			(*toktab)[(*nb_of_classes)++] = strdup(currtok);
		}
		else{
			classcount[i]++; /* one more occurrence for this class */
		}
	}

/* create new framelist */		
	out = calloc(*nb_of_classes,sizeof(framelist));
	if (debug_active()>1)
		e_comment(0,"classify_framelist_by_token: found %d different tokens in %d files",
				*nb_of_classes,in->n);
	if (!in->type)
		e_warning("classify_framelist_by_tokenin->type NULL");
	for (j=0;j<*nb_of_classes;j++){
		if (debug_active()>1)
			e_comment(0,"%s %d occurences",(*toktab)[j],classcount[j]);
		out[j].n = classcount[j];
		out[j].filename = calloc(strlen(in->filename)+1,sizeof(char));
		out[j].name = calloc(classcount[j],sizeof(char*));
		if (in->type)
			out[j].type = calloc(classcount[j],sizeof(char*));
		else
			out[j].type = NULL;
		classcount[j]=0;
	}
/* fill in new framelist */	
	for (j=0; j<in->n; j++){
		currtok=(*funk_get_tok)(in->name[j]);	
		if (!currtok)
			currtok = "UNKNOWN";
		found_class = 0;
		for (i=0;i<*nb_of_classes;i++){
			if (!strcmp(currtok,(*toktab)[i])){
				found_class=1;
				break;
			}
		}
		if (found_class){
			index_newcm = classcount[i]++; 
			c = &(out[i]);
			c->name[index_newcm] = strdup(in->name[j]);	
			if (in->type){
				if (in->type[j]){
					c->type[index_newcm]= strdup(in->type[j]);	
				}
				else
					c->type[index_newcm]=NULL;
			}
		}
		else{
			e_error("classify_framelist_by_token: could not find token %s in table",
					currtok);
		}
	}
	if (debug_active()>2){
		for (i=0;i<*nb_of_classes;i++){
			e_comment(0,"class %s",(*toktab)[i]);
			for (j=0;j<classcount[i];j++){
				e_comment(0,"file %s",out[i].name[j]);
			}
		}
	}
	free(classcount);
	return out;
}

	
/*-------------------------------------------------------------------------*/
/**
  @brief  selects files in a list of files according to a token function. 
  @param  in				files to classify (freed).
  @param	nfiles			number of files in (modified)
  @param	funk_get_tok	the token retrieval function
  @param	token			the token to match by funk_get_tok()
  @return list of files matching token, NULL if error
  This function is typically used to selct frames in a
  filelist (char**). 
  The token has to match what is retrieved by funk_get_tok which
  takes a filename of the framelist as a parameter. A typical example
  is one of the functions isaac_get_XXX in dicbkeys.h 
  The function returns a file list char**.
  The input list is deallocated and the nfiles parameter is updated
  for convenience, i.e. this function is called repeatedly whith
  different token/funk_get_tok's.
 */
/*--------------------------------------------------------------------------*/
char ** select_charlist_by_token(
		char	**	in, 
		int		*	nfiles,
		char	*	token,
		char	*	(*funk_get_tok)(char*))
{
	char	 	**	out=NULL,
				***	class=NULL;
	char		**	toktab=NULL;	
	int				i,j;
	int				nb_of_classes;
	int			*	classcount=NULL;
	int				found = 0;
	int				err = 0;
	int				in_files = 0;
	
	if (!in){
		e_error("select_charlist_by_token: in NULL");
		return NULL;
	}
	
	if (!nfiles){
		e_error("select_charlist_by_token: nfiles NULL");
		return NULL;
	}
	if (!(*nfiles)){
		e_error("select_charlist_by_token: no files");
		return NULL;
	}
	if (!token){
		e_error("select_charlist_by_token: token NULL");
		return NULL;
	}
	in_files = *nfiles;
	class = classify_charlist_by_token( 
				in,
				in_files,
				&toktab,
				&nb_of_classes,
				&classcount,
				funk_get_tok);
	if (!class) /* error occured in classify_charlist_by_token */
		return NULL;
	if (!toktab)
		return NULL;
	if (debug_active()>1)
		e_comment(0,"Selecting on token %s in %d classes",token,
			nb_of_classes);
/* find relevent class */	
	for (j=0; j<nb_of_classes; j++){
		if (!toktab[j]){
			e_error("toktab[%d] NULL",j);
			err=-1;
			break;
		}
		if (!strcmp(token,toktab[j]) && !err){
			
			*nfiles = classcount[j];
			if (!(*nfiles)){
				e_error("select_charlist_by_token: no files by token %s",
						token);
				err=-1;
			}
			out=calloc(*nfiles,sizeof(char*));
			for (i=0;i<*nfiles;i++)
				out[i]= strdup(class[j][i]);
			found = 1;
		}
	}
	if (!found)
		e_warning("token not found");
	for (i=0;i<nb_of_classes;i++){
		for (j=0;j<classcount[i];j++){
			free(class[i][j]);
		}
		free(class[i]);
		free(toktab[i]);
	}
	free(class);
	free(classcount);
	free(toktab);	
	for (i=0;i<in_files;i++)
		free(in[i]);
	free(in);
	return out;
}