/* $Id: vimos_utils.c,v 1.18 2015/11/27 12:15:33 jim Exp $
 *
 * This file is part of the VIMOS Pipeline
 * Copyright (C) 2015 European Southern Observatory
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 */

/*
 * $Author: jim $
 * $Date: 2015/11/27 12:15:33 $
 * $Revision: 1.18 $
 * $Name:  $
 */

/* Includes */

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <sys/time.h>
#include <time.h>
#include <libgen.h>
#include <string.h>
#include <unistd.h>

#include <cpl.h>
#include <math.h>

#include "vimos_imaging_utils.h"
#include "vimos_pfits.h"
#include <casu_utils.h>
#include <casu_stats.h>
#include <casu_fits.h>
#include <casu_tfits.h>
#include <casu_mods.h>
#include <casu_wcsutils.h>

/**
    \defgroup vimos_utils vimos_utils
    \ingroup supportroutines

    \brief
    These are utility routines of various types.

    \author
    Jim Lewis, CASU
*/

/**@{*/

/*---------------------------------------------------------------------------*/
/**
    \par Name:
        vimos_create_diffimg_stats
    \par Purpose:
        Create an empty difference image stats table
    \par Description:
        Create an empty difference image stats table
    \par Language:
        C
    \param nrows
        The number of rows for the table
    \returns
        The cpl_table pointer for the new stats table
    \author
        Jim Lewis, CASU
 */
/*---------------------------------------------------------------------------*/

cpl_table *vimos_create_diffimg_stats(int nrows) {
    cpl_table *diffimstats;

    diffimstats = cpl_table_new((cpl_size)nrows);
    cpl_table_new_column(diffimstats,"xmin",CPL_TYPE_INT);
    cpl_table_set_column_unit(diffimstats,"xmin","pixels");
    cpl_table_new_column(diffimstats,"xmax",CPL_TYPE_INT);
    cpl_table_set_column_unit(diffimstats,"xmax","pixels");
    cpl_table_new_column(diffimstats,"ymin",CPL_TYPE_INT);
    cpl_table_set_column_unit(diffimstats,"ymin","pixels");
    cpl_table_new_column(diffimstats,"ymax",CPL_TYPE_INT);
    cpl_table_set_column_unit(diffimstats,"ymax","pixels");
    cpl_table_new_column(diffimstats,"mean",CPL_TYPE_FLOAT);
    cpl_table_set_column_unit(diffimstats,"mean","ADU");
    cpl_table_new_column(diffimstats,"median",CPL_TYPE_FLOAT);
    cpl_table_set_column_unit(diffimstats,"median","ADU");
    cpl_table_new_column(diffimstats,"variance",CPL_TYPE_FLOAT);
    cpl_table_set_column_unit(diffimstats,"variance","ADU**2");
    cpl_table_new_column(diffimstats,"mad",CPL_TYPE_FLOAT);
    cpl_table_set_column_unit(diffimstats,"mad","ADU");
    return(diffimstats);
}

/*---------------------------------------------------------------------------*/
/**
    \par Name:
        vimos_difference_image
    \par Purpose:
        Create a difference/ratio image and difference/ratio image stats table
    \par Description:
        A difference/ratio image is created from an input image and a master 
        image. A global difference and RMS are calculated from the difference 
        image. Then a difference image stats table will be created. This 
        breaks up the difference image into cells and calculates some basic 
        stats in each cell.
    \par Language:
        C
    \param master
        The master calibration image
    \param prog
        The new mean image
    \param bpm 
        Input bad pixel mask
    \param ncells
        The number of cells per channel
    \param oper
        The operation to be performed:
        - 1. Subtract the images
        - 2. Divide the images
    \param global_diff 
        The median difference over the whole difference image
    \param global_rms  
        The rms difference over the whole difference image
    \param diffim
        The output difference/ratio image
    \param diffimstats 
        The output difference/ratio image statistics table
    \return   
        Nothing
    \author
        Jim Lewis, CASU
 */
/*---------------------------------------------------------------------------*/

void vimos_difference_image(cpl_image *master, cpl_image *prog, 
                            unsigned char *bpm, int ncells, int oper, 
                            float *global_diff, float *global_rms, 
                            cpl_image **diffim, 
                            cpl_table **diffimstats) {
    float *ddata,*work,mean,sig,med,mad;
    long nx,ny,npts;
    int nc1,nc2,nr,idx,idy;
    int icx,icy,indy1,indy2,indx1,indx2,jp,jcx,jj,jcy,ii;
    const char *fctid = "vimos_difference_image";

    /* Initialise the output */

    *diffim = NULL;
    *diffimstats = NULL;
    *global_diff = 0.0;
    *global_rms = 0.0;

    /* Is there a programme frame or a master? */

    if (prog == NULL || master == NULL)
        return;

    /* Start by subtracting the master image from the programme image */

    switch (oper) {
    case 1:
        *diffim = cpl_image_subtract_create(prog,master);
        break;
    case 2:
        *diffim = cpl_image_divide_create(prog,master);
        break;
    default:
        *diffim = NULL;
        cpl_msg_error(fctid,"Invalid operation requested %" CPL_SIZE_FORMAT,
                      (cpl_size)oper);
        break;
    }      
    if (*diffim == NULL)
        return;

    /* Work out a median difference */

    ddata = cpl_image_get_data_float(*diffim);
    nx = (int)cpl_image_get_size_x(*diffim);
    ny = (int)cpl_image_get_size_y(*diffim);
    npts = nx*ny;
    casu_medmad(ddata,bpm,npts,global_diff,global_rms);
    *global_rms *= 1.48;

    /* Work out how to divide the image */

    switch (ncells) {
    case 1:
        nc1 = 1;
        nc2 = 1;
        break;
    case 2:
        nc1 = 2;
        nc2 = 1;
        break;
    case 4:
        nc1 = 2;
        nc2 = 2;
        break;
    case 8:
        nc1 = 4;
        nc2 = 2;
        break;
    case 16:
        nc1 = 4;
        nc2 = 4;
        break;
    case 32:
        nc1 = 8;
        nc2 = 4;
        break;
    case 64:
        nc1 = 8;
        nc2 = 8;
        break;
    default:
        nc1 = 8;
        nc2 = 8;
        break;
    }

    /* Create a difference image stats table */

    *diffimstats = vimos_create_diffimg_stats(nc1*nc2);

    /* How big is each cell? */

    idx = nx/nc1;
    idy = ny/nc2;
    work = cpl_malloc(idx*idy*sizeof(*work));

    /* Now loop for each cell */

    nr = 0;
    for (icy = 0; icy < nc2; icy++) {
        indy1 = idy*icy;
        indy2 = min(ny,indy1+idy-1);
        for (icx = 0; icx < nc1; icx++) {
            indx1 = idx*icx;
            indx2 = min(nx,indx1+idx-1);
            jp = 0;
            for (jcy = indy1; jcy < indy2; jcy++) {
                jj = jcy*nx;
                for (jcx = indx1; jcx < indx2; jcx++) {
                    ii = jj + jcx;
                    if (bpm != NULL && bpm[ii] == 0)
                        work[jp++] = ddata[ii];
                    else if (bpm == NULL)
                        work[jp++] = ddata[ii];
                }
            }
            (void)casu_meansig(work,NULL,(long)jp,&mean,&sig);
            (void)casu_medmad(work,NULL,(long)jp,&med,&mad);
            cpl_table_set_int(*diffimstats,"xmin",(cpl_size)nr,indx1+1);
            cpl_table_set_int(*diffimstats,"xmax",(cpl_size)nr,indx2+1);
            cpl_table_set_int(*diffimstats,"ymin",(cpl_size)nr,indy1+1);
            cpl_table_set_int(*diffimstats,"ymax",(cpl_size)nr,indy2+1);
            cpl_table_set_float(*diffimstats,"mean",(cpl_size)nr,mean);
            cpl_table_set_float(*diffimstats,"median",(cpl_size)nr,med);
            cpl_table_set_float(*diffimstats,"variance",(cpl_size)nr,
                                (sig*sig));
            cpl_table_set_float(*diffimstats,"mad",(cpl_size)(nr++),mad);
        }
    }
    cpl_free(work);
}       

/*---------------------------------------------------------------------------*/
/**
    \par Name:
        vimos_get_det_frameset
    \par Purpose:
        Extract a frameset for a particular detector
    \par Description:
        VIMOS input files contain the data for just one of the four detectors
        in the focal plane. This routine takes a frameset of a whole list
        of input files that have been 'labelised' by detector name and
        extracts the subset that match a given detector name.
    \par Language:
        C
    \param in
        The input frameset. This must have been labelised by detector name
        before calling this routine.
    \param labels
        The 'by detector name' labels
    \param nlab
        The number of labels in the list
    \param testname
        The name of the detector whose frames you want.
    \return   
        The frameset containing the subset of frames requested. If none
        of the frames in the input frameset match, then NULL is returned.
    \author
        Jim Lewis, CASU
 */
/*---------------------------------------------------------------------------*/

cpl_frameset *vimos_get_det_frameset(cpl_frameset *in, cpl_size *labels, 
                                     cpl_size nlab, char *testname) {
    cpl_frameset *res;
    cpl_propertylist *p;
    cpl_frame *fr;
    int i,match;
    char name[16];

    for (i = 0; i < nlab; i++) {
        res = cpl_frameset_extract(in,labels,(cpl_size)i);
        fr = cpl_frameset_get_position(res,0);
        p = cpl_propertylist_load(cpl_frame_get_filename(fr),
                                  cpl_frame_get_nextensions(fr));
        vimos_pfits_get_chipname(p,name);
        match = (! strcmp(name,testname));
        cpl_propertylist_delete(p);
        if (match) 
            return(res);
        else
            cpl_frameset_delete(res);
    }
    return(NULL);
}

/*---------------------------------------------------------------------------*/
/**
    \par Name:
        vimos_get_obs_frameset
    \par Purpose:
        Extract a frameset for a particular observation
    \par Description:
        VIMOS input files contain the data for just one of the four detectors
        in the focal plane. This routine takes a frameset of a whole list
        of input files that have been 'labelised' by observation time and
        extracts the subset for a given label. The subset is reordered
        into a standard output detector order.
    \par Language:
        C
    \param in
        The input frameset. This must have been labelised by MJD
        before calling this routine.
    \param labels
        The 'by MJD' labels
    \param nlab
        The number of labels in the list
    \param ilab
        The label number you want to extract.
    \param vimos_names
        The names of the detectors in use.
    \return   
        The frameset containing the subset of frames requested in
        the correct order. If a frame doesn't exist for a particular
        detector, then a dummy frame with the tag DUMMY is put 
        as a place holder in the output frameset
    \author
        Jim Lewis, CASU
 */
/*---------------------------------------------------------------------------*/

cpl_frameset *vimos_get_obs_frameset(cpl_frameset *in, cpl_size *labels, 
                                     cpl_size nlab, cpl_size ilab,
                                     char *vimos_names[]) {
    cpl_frameset *res,*test;
    cpl_propertylist *p;
    cpl_frame *fr,*fr2;
    int i,j,match;
    char name[16];

    test = cpl_frameset_extract(in,labels,ilab);
    res = cpl_frameset_new();
    for (i = 0; i < VIMOS_NEXTN; i++) {
        match = 0;
        for (j = 0; j < cpl_frameset_get_size(test); j++) {
            fr = cpl_frameset_get_position(test,(cpl_size)j);
            p = cpl_propertylist_load(cpl_frame_get_filename(fr),
                                      cpl_frame_get_nextensions(fr));
            vimos_pfits_get_chipname(p,name);
            cpl_propertylist_delete(p);
            if (strncmp(name,vimos_names[i],16) == 0) {
                match = 1;
                break;
            }
        }
        if (match == 1) {
            fr2 = cpl_frame_duplicate(fr);
            cpl_frameset_insert(res,fr2);
        } else {
            fr2 = cpl_frame_new();
            cpl_frame_set_tag(fr2,"DUMMY");
            cpl_frameset_insert(res,fr2);
            cpl_msg_warning("vimos_get_obs_frameset",
                            "Observation for detector %s missing in set %d",
                            vimos_names[i],(int)ilab);
        }
    }
    cpl_frameset_delete(test);
    return(res);
}


/*---------------------------------------------------------------------------*/
/**
    \par Name:
        vimos_compare_names
    \par Purpose:
        Comparison routine used to labelise a frameset by detector name
    \par Description:
        Compare the detector names for two frames. Return the string
        comparison.
    \par Language:
        C
    \param fr1
        The first frame.
    \param fr2
        The second frame
    \return   
        The result of the string comparison of the detectors names for
        the two input frames
    \author
        Jim Lewis, CASU
 */
/*---------------------------------------------------------------------------*/

int vimos_compare_names(const cpl_frame *fr1, const cpl_frame *fr2) {
    cpl_propertylist *pp;
    char name1[16],name2[16];

    /* Test entries */

    if (fr1 == NULL || fr2 == NULL) 
        return(-1);

    /* Load the propertylist for each frame and extract the chip name
       NB: we have to test the number of extensions here. If 
       we are dealing with tile compressed images, the FITS header will
       be in the first extension. If it's not compressed then the info will
       be in the primary */

    pp = cpl_propertylist_load(cpl_frame_get_filename(fr1),
                               cpl_frame_get_nextensions(fr1));
    (void)vimos_pfits_get_chipname(pp,name1);
    cpl_propertylist_delete(pp);
    pp = cpl_propertylist_load(cpl_frame_get_filename(fr2),
                               cpl_frame_get_nextensions(fr2));
    (void)vimos_pfits_get_chipname(pp,name2);
    cpl_propertylist_delete(pp);

    /* Compare the names */

    if (strcmp(name1,name2)) 
        return(0);
    else
        return(1);
}
/*---------------------------------------------------------------------------*/
/**
    \par Name:
        vimos_compare_mjds
    \par Purpose:
        Comparison routine used to labelise a frameset by mjd
    \par Description:
        Compare the mjds for two frames. These have to be within
        0.5 a second because of inconsistencies in the reported 
        MJD for the single exposure raw frames. Return the comparison.
    \par Language:
        C
    \param fr1
        The first frame.
    \param fr2
        The second frame
    \return   
        The result of the comparison of the MJDs for the two input frames
    \author
        Jim Lewis, CASU
 */
/*---------------------------------------------------------------------------*/

int vimos_compare_mjds(const cpl_frame *fr1, const cpl_frame *fr2) {
    cpl_propertylist *pp;
    double mjd1,mjd2,adiff;

    /* Test entries */

    if (fr1 == NULL || fr2 == NULL) 
        return(-1);

    /* Load the propertylist for each frame and extract the Modified Julian
       Date. NB: we have to test the number of extensions here. If 
       we are dealing with tile compressed images, the FITS header will
       be in the first extension. If it's not compressed then the info will
       be in the primary */

    pp = cpl_propertylist_load(cpl_frame_get_filename(fr1),
                               cpl_frame_get_nextensions(fr1));
    (void)vimos_pfits_get_mjd(pp,&mjd1);
    cpl_propertylist_delete(pp);
    pp = cpl_propertylist_load(cpl_frame_get_filename(fr2),
                               cpl_frame_get_nextensions(fr2));
    (void)vimos_pfits_get_mjd(pp,&mjd2);
    cpl_propertylist_delete(pp);

    /* Round the dates to 1.0e-6 */

    adiff = fabs(mjd1 - mjd2);

    /* Compare the times and see if they are less than 0.5 seconds apart */

    return(adiff < 5.79e-6);
}

/*---------------------------------------------------------------------------*/
/**
    \par Name:
        vimos_compare_tpl_start
    \par Purpose:
        Comparison routine used to labelise a frameset by TPL START
    \par Description:
        Compare the TPL START for two frames.
    \par Language:
        C
    \param fr1
        The first frame.
    \param fr2
        The second frame
    \return   
        The result of the comparison of the TPL START for the two input frames
    \author
        Matteo Salmistraro, ESO
 */
/*---------------------------------------------------------------------------*/

int vimos_compare_tpl_start(const cpl_frame *fr1, const cpl_frame *fr2) {

    const char * lst1;
    const char * lst2;

    /* Test entries */

    if (fr1 == NULL || fr2 == NULL) 
        return(-1);

    /* Load the propertylist for each frame and extract the Modified Julian
       Date. NB: we have to test the number of extensions here. If 
       we are dealing with tile compressed images, the FITS header will
       be in the first extension. If it's not compressed then the info will
       be in the primary */

    const char * fname1 = cpl_frame_get_filename(fr1);
    cpl_propertylist * pp1 = cpl_propertylist_load(fname1,
                               cpl_frame_get_nextensions(fr1));

    (void)vimos_pfits_get_tpl_start(pp1, &lst1);


    const char * fname2 = cpl_frame_get_filename(fr2);
    cpl_propertylist * pp2 = cpl_propertylist_load(fname2,
                               cpl_frame_get_nextensions(fr2));

    (void)vimos_pfits_get_tpl_start(pp2, &lst2);

    int diff = strcmp(lst1, lst2);

    cpl_propertylist_delete(pp1);
    cpl_propertylist_delete(pp2);

    return diff == 0;
}

/*---------------------------------------------------------------------------*/
/**
    \par Name:
        vimos_compare_lsts
    \par Purpose:
        Comparison routine used to labelise a frameset by LST
    \par Description:
        Compare the lst for two frames. These _should_ be identical
        but we'll let them differ by 0.5 a second just in case.
        Return the comparison.
    \par Language:
        C
    \param fr1
        The first frame.
    \param fr2
        The second frame
    \return
        The result of the comparison of the LSTs for the two input frames
    \author
        Jim Lewis, CASU
 */
/*---------------------------------------------------------------------------*/

int vimos_compare_lsts(const cpl_frame *fr1, const cpl_frame *fr2) {
    cpl_propertylist *pp;
    double lst1,lst2,adiff;

    /* Test entries */

    if (fr1 == NULL || fr2 == NULL)
        return(-1);

    /* Load the propertylist for each frame and extract the Modified Julian
       Date. NB: we have to test the number of extensions here. If
       we are dealing with tile compressed images, the FITS header will
       be in the first extension. If it's not compressed then the info will
       be in the primary */

    const char * fname1 = cpl_frame_get_filename(fr1);
    pp = cpl_propertylist_load(fname1,
                               cpl_frame_get_nextensions(fr1));
    (void)vimos_pfits_get_lst(pp,&lst1);
    cpl_propertylist_delete(pp);

    const char * fname2 = cpl_frame_get_filename(fr2);
    pp = cpl_propertylist_load(fname2,
                               cpl_frame_get_nextensions(fr2));
    (void)vimos_pfits_get_lst(pp,&lst2);
    cpl_propertylist_delete(pp);

    /* Look at the difference ... */

    adiff = fabs(lst1 - lst2);

    /* Compare the times and see if they are less than 2 seconds apart.
       This might fail if you have two exposures on either side of midnight. */

    return(adiff < 2.0);
}


/*---------------------------------------------------------------------------*/
/**
    \par Name:
        vimos_load_names
    \par Purpose:
        Load the standard list of detectors
    \par Description:
        Look at the chip name in one frame. From the name of the chip,
        decide which group of chip names is appropriate.
    \par Language:
        C
    \param frm
        The input frame
    \param vimos_names
        The output string array with the names to be used
    \return   
        CASU_OK if all went well. CASU_FATAL if no name was found in the
        input frame header.
    \author
        Jim Lewis, CASU
 */
/*---------------------------------------------------------------------------*/

int vimos_load_names(cpl_frame *frm, char *vimos_names[]) {
    const char *fctid = "vimos_load_names";
    char chipname[16];
    cpl_propertylist *p;
    int i;
    
    /* Check to make sure you have a valid propertylist */

    if (frm == NULL) {
        cpl_msg_error(fctid,"Unable to get name set from supplied frame");
        return(CASU_FATAL);
    }
    p = cpl_propertylist_load(cpl_frame_get_filename(frm),
                              cpl_frame_get_nextensions(frm));

    /* Get the name of the detector. If it starts with the string "CCD" then
       it's the old set. */

    vimos_pfits_get_chipname(p,chipname);
    if (! strncmp(chipname,"CCD",3)) {
        for (i = 0; i < 4; i++)
            vimos_names[i] = (char *)vimos_oldnames[i];
    } else {
        for (i = 0; i < 4; i++)
            vimos_names[i] = (char *)vimos_newnames[i];
    }
    cpl_propertylist_delete(p);
    return(CASU_OK);
}

/*---------------------------------------------------------------------------*/
/**
    \par Name:
        vimos_load_trimreg
    \par Purpose:
        Load the standard trim regions for a detector set
    \par Description:
        Look at the chip name in one frame. From the name of the chip,
        decide which set of trim regions is appropriate
    \par Language:
        C
    \param p
        The input propertylist
    \param which
        If set to 1 then the old chips are being used. If 2 then it's the
        new chips
    \return   
        CASU_OK if all went well.
    \author
        Jim Lewis, CASU
 */
/*---------------------------------------------------------------------------*/

int vimos_load_trimreg(cpl_propertylist *p, int *which) {
    char chipname[16];
    
    /* Get the name of the detector. If it starts with the string "CCD" then
       it's the old set. */

    vimos_pfits_get_chipname(p,chipname);
    if (! strncmp(chipname,"CCD",3)) 
        *which = 1;
    else
        *which = 2;

    return(CASU_OK);
}

/*---------------------------------------------------------------------------*/
/**
    \par Name:
        vimos_wcsfit
    \par Purpose:
        Do a WCS fit to an individual image
    \par Description:
        For each image in a list, create a source catalogue if one isn't 
        given in the argument list. Fit a WCS to the source positions.
    \par Language:
        C
    \param in
        A list of casu_fits structures for the input images
    \param incat
        A list of casu_tfits structures for the input source catalogues.
    \param nf
        The number of images in the lists
    \param catname
        The name of the standard catalogue
    \param catpath
        The full path to the standard index file
    \param cdssearch
        The CDS catalogue to be searched for standards. 0 if using local
        catalogues.
    \param cacheloc
        The location of the standard star cache
    \param keepms
        If set, then the matched standards files will be kept
    \param trim
        If set, then we will ignore objects in the source catalogue bands 
	where the data is likely to be rubbish. Value is the number of the
	extension that we're currently processing. Only set this for standard
	star exposures as science exposures will have already been trimmed
    \param mstds
        The matched standards catalogue
    \returns
        Nothing
    \author
        Jim Lewis, CASU
 */
/*---------------------------------------------------------------------------*/
    
void vimos_wcsfit(casu_fits **in, casu_tfits **incat, int nf, 
                  char *catname, char *catpath, int cdssearch,
                  char *cacheloc, int keepms, int trim,
                  casu_tfits **mstds) {
    int status,nstd,ncat,slevel,n,i,j,oldnew,imin,imax,jmin,jmax;
    float *x,*y,*ra,*dec;
    double r,d,*rac,*dra,*decc,*ddec,*ra_d,*dec_d;
    casu_tfits *tcat;
    cpl_table *stdscat,*cat,*tmp,*tmp2,*matchstds;
    cpl_propertylist *p;
    cpl_wcs *wcs;
    const char *fctid = "vimos_wcsfit";

    /* Loop for all the images */

    for (j = 0; j < nf; j++) {

        /* Define a few things */

        status = CASU_OK;
        tcat = incat[j];

        /* Get some standard stars */

        (void)casu_getstds(casu_fits_get_ehu(in[j]),1,catpath,catname,cdssearch,
                           cacheloc,&stdscat,&status);
        if (status != CASU_OK) {
            freetable(stdscat);
            cpl_msg_error(fctid,
                          "Failed to find any standards for %s[%" CPL_SIZE_FORMAT "]",
                          casu_fits_get_filename(in[j]),
                          (cpl_size)casu_fits_get_nexten(in[j]));
            continue;
        }
        nstd = (int)cpl_table_get_nrow(stdscat);
	
	/* If we want to ignore the objects in the trim regions, then
	   do that now. (NB: standard observations only as science images 
	   will have already been trimmed) */

        cat = casu_tfits_get_table(tcat);
        ncat = (int)cpl_table_get_nrow(cat);
	cpl_table_select_all(cat);
	cpl_table_and_selected_float(cat,"Classification",CPL_NOT_EQUAL_TO,0.0);
	tmp = NULL;
	if (trim) {
	    (void)vimos_load_trimreg(casu_fits_get_ehu(in[j]),&oldnew);
	    if (oldnew == 1) {
		imin = trimreg_old[trim-1][0];
		imax = trimreg_old[trim-1][1];
		jmin = trimreg_old[trim-1][2];
		jmax = trimreg_old[trim-1][3];
	    } else {
		imin = trimreg_new[trim-1][0];
		imax = trimreg_new[trim-1][1];
		jmin = trimreg_new[trim-1][2];
		jmax = trimreg_new[trim-1][3];
	    }
	    (void)cpl_table_and_selected_float(cat,"X_coordinate",
					       CPL_GREATER_THAN,(float)imin);
	    (void)cpl_table_and_selected_float(cat,"X_coordinate",
					       CPL_LESS_THAN,(float)imax);
	    (void)cpl_table_and_selected_float(cat,"Y_coordinate",
					       CPL_GREATER_THAN,(float)jmin);
	    (void)cpl_table_and_selected_float(cat,"Y_coordinate",
					       CPL_LESS_THAN,(float)jmax);
	    tmp = cpl_table_extract_selected(cat);
	    cat = tmp;
	}

        /* If there are too many objects in the catalogue then first restrict
           ourselves by ellipticity. Cut so that there are similar numbers of
           objects in the standards and the object catalogues by retaining the
           brighter objects */

        if (ncat > 500 && ncat > 2.0*nstd) {
            tmp = cpl_table_duplicate(cat);
	    cpl_table_unselect_all(tmp);
            (void)cpl_table_or_selected_float(tmp,"Ellipticity",CPL_LESS_THAN,0.5);
            tmp2 = cpl_table_extract_selected(tmp);
            ncat = (int)cpl_table_get_nrow(tmp2);
            freetable(tmp);
            p = cpl_propertylist_new();
            cpl_propertylist_append_bool(p,"Isophotal_flux",TRUE);
            cpl_table_sort(tmp2,(const cpl_propertylist *)p);
            cpl_propertylist_delete(p);
            slevel = min(ncat,max(1,min(5000,max(500,2*nstd))));
            tmp = cpl_table_extract(tmp2,1,(cpl_size)slevel);
            freetable(tmp2);
            ncat = (int)cpl_table_get_nrow(tmp);
            cat = tmp;
        }

        /* Now match this against the catalogue */

        (void)casu_matchstds(cat,stdscat,300.0,&matchstds,&status);
        freetable(stdscat);
        freetable(tmp);
        if (status != CASU_OK) {
            freetable(matchstds);
            cpl_msg_error(fctid,"Failed to match standards to catalogue");
            return;
        }

        /* Fit the plate solution */

        (void)casu_platesol(casu_fits_get_ehu(in[j]),
                            casu_tfits_get_ehu(tcat),matchstds,6,1,&status);
        if (keepms) {
            cpl_table_duplicate_column(matchstds,"RA_calc",matchstds,
                                       "RA");
            cpl_table_duplicate_column(matchstds,"diffRA",matchstds,"RA");
            cpl_table_duplicate_column(matchstds,"Dec_calc",matchstds,
                                       "Dec");
            cpl_table_duplicate_column(matchstds,"diffDec",matchstds,"Dec");
            mstds[j] = casu_tfits_wrap(matchstds,tcat,NULL,NULL);
        } else {
            freetable(matchstds);
        }
        if (status != CASU_OK) {
            cpl_msg_error(fctid,"Failed to fit WCS");
            continue;
        }

        /* Update the RA and DEC of the objects in the object catalogue */

        cat = casu_tfits_get_table(tcat);
        n = (int)cpl_table_get_nrow(cat);
        wcs = cpl_wcs_new_from_propertylist(casu_fits_get_ehu(in[j]));
        if (wcs == NULL) {
            cpl_msg_error(fctid,"Failed to fill RA and Dec in catalogue");
            return;
        }
        x = cpl_table_get_data_float(cat,"X_coordinate"); 
        y = cpl_table_get_data_float(cat,"Y_coordinate");
        ra = cpl_table_get_data_float(cat,"RA");
        dec = cpl_table_get_data_float(cat,"DEC");
        for (i = 0; i < n; i++) {
            casu_xytoradec(wcs,(double)x[i],(double)y[i],&r,&d);
            ra[i] = (float)r;
            dec[i] = (float)d;
        }
        cpl_wcs_delete(wcs);

        /* If keeping the matchstds catalogue, then add some stuff to it */

        if (keepms) {
            rac = cpl_table_get_data_double(matchstds,"RA_calc");
            dra = cpl_table_get_data_double(matchstds,"diffRA");
            decc = cpl_table_get_data_double(matchstds,"Dec_calc");
            ddec = cpl_table_get_data_double(matchstds,"diffDec");

            /* Now compute the equatorial coordinates and compare with the 
               standard values */

            n = (int)cpl_table_get_nrow(matchstds);
            x = cpl_table_get_data_float(matchstds,"X_coordinate");
            y = cpl_table_get_data_float(matchstds,"Y_coordinate");
            ra_d = cpl_table_get_data_double(matchstds,"RA");
            dec_d = cpl_table_get_data_double(matchstds,"Dec");
            wcs = cpl_wcs_new_from_propertylist(casu_fits_get_ehu(in[j]));
            for (i = 0; i < n; i++) {
                casu_xytoradec(wcs,(double)x[i],(double)y[i],&r,&d);
                rac[i] = (float)r;
                decc[i] = (float)d;
                dra[i] = rac[i] - ra_d[i];
                ddec[i] = decc[i] - dec_d[i];
            }
            cpl_wcs_delete(wcs);
        }           
    }
}

/*---------------------------------------------------------------------------*/
/**
    \par Name:
        vimos_testfrms
    \par Purpose:
        Test the frames in a frameset to make sure they are what you 
        expect
    \par Description:
        Each of the frames in a frameset are checked to make sure they
        have the right number of extensions and that each image/table
        is loadable.
    \par Language:
        C
    \param frms
        The input frameset
    \param nextn_expected
        The number of extensions expected in each image/table file
    \param isimg
        1 if these are image frames. 0 we check for tables. -1 could be either
    \param checkwcs
        TRUE if we need to check for a valid WCS
    \return
        The number of errors
    \author
        Jim Lewis, CASU
 */
/*---------------------------------------------------------------------------*/

int vimos_testfrms(cpl_frameset *frms, int nextn_expected, int isimg, 
                   int checkwcs) {
    int i,nf,nerr;
    cpl_frame *fr;

    /* Return immediately if given nonsense */

    if (frms == NULL) 
        return(0);

    /* Loop for each frame in the frameset */

    nf = cpl_frameset_get_size(frms);
    nerr = 0;
    for (i = 0; i < nf; i++) {
        fr = cpl_frameset_get_position(frms,i);
        nerr += vimos_testfrm_1(fr,nextn_expected,isimg,checkwcs);
    }
    
    /* Return value */

    return(nerr);
}

/*---------------------------------------------------------------------------*/
/**
    \par Name:
        vimos_testfrm_1
    \par Purpose:
        Test a frame to make sure it is what you expect
    \par Description:
        The frame is checked to make sure it has the right number of 
        extensions and that each image/table is loadable.
    \par Language:
        C
    \param fr
        The input frame
    \param nextn_expected
        The number of extensions expected in each image/table file
    \param isimg
        1 if this is an image frame. 0 we check for tables. -1 we work
        out what it is by looking at the first extension.
    \param checkwcs
        TRUE if we need to check for a valid WCS
    \return
        The number of errors
    \author
        Jim Lewis, CASU
 */
/*---------------------------------------------------------------------------*/

int vimos_testfrm_1(cpl_frame *fr, int nextn_expected, int isimg,
                    int checkwcs) {
    int nextn,nerr,j,simple,jj,isimg2;
    casu_fits *test;
    casu_tfits *testt;
    cpl_propertylist *p;
    cpl_wcs *wcs;
    const char *fctid="vimos_testfrm_1";

    /* Return immediately if given nonsense */

    if (fr == NULL)
        return(0);

    /* Test to see how many extensions there are and compare to see
       if it matches the number expected */

    nextn = cpl_frame_get_nextensions(fr);
    simple = 0;
    if (isimg == 1 && nextn == 0) {
        simple = 1;
    } else if (nextn != nextn_expected) {
        cpl_msg_error(fctid,"Frame %s has %" CPL_SIZE_FORMAT " extensions, expected %" CPL_SIZE_FORMAT "\n",
                      cpl_frame_get_filename(fr),(cpl_size)nextn,
                      (cpl_size)nextn_expected);
        return(1);
    }

    /* If we don't know what to expect, then work it out here */

    if (isimg == -1) {
        jj = (nextn == 0 ? 0 : 1);
        test = casu_fits_load(fr,CPL_TYPE_FLOAT,jj);
        isimg2 = (test == NULL ? 0 : 1);
    } else {
        isimg2 = isimg;
    }

    /* Test to see if you can load each of the extensions */

    nerr = 0;
    for (j = 1; j <= nextn; j++) {
        if (isimg2 == 1) {
            jj = (simple ? 0 : j);
            test = casu_fits_load(fr,CPL_TYPE_FLOAT,jj);
            if (test == NULL) {
                cpl_msg_error(fctid,
                              "Frame image %s[%" CPL_SIZE_FORMAT "] won't load",
                              cpl_frame_get_filename(fr),(cpl_size)jj);
                nerr++;
                continue;
            } else if (casu_is_dummy(casu_fits_get_ehu(test))) {
                cpl_msg_error(fctid,
                              "Frame image %s[%" CPL_SIZE_FORMAT "] is a dummy",
                              cpl_frame_get_filename(fr),(cpl_size)jj);
                nerr++;
                continue;
            }
            if (checkwcs) {
                p = casu_fits_get_ehu(test);
                wcs = cpl_wcs_new_from_propertylist(p);
                if (wcs == NULL) {
                    cpl_msg_error(fctid,
                                  "Frame image %s[%" CPL_SIZE_FORMAT "] WCS invalid",
                                  cpl_frame_get_filename(fr),(cpl_size)jj);
                    nerr++;
                    continue;
                } 
                freewcs(wcs);
            }
            freefits(test);
        } else {
            testt = casu_tfits_load(fr,j);
            if (testt == NULL) {
                cpl_msg_error(fctid,
                              "Frame table %s[%" CPL_SIZE_FORMAT "] won't load\n",
                              cpl_frame_get_filename(fr),(cpl_size)j);
                nerr++;
                continue;
            } else if (casu_is_dummy(casu_tfits_get_ehu(testt))) {
                cpl_msg_error(fctid,
                              "Frame table %s[%" CPL_SIZE_FORMAT "] is a dummy",
                              cpl_frame_get_filename(fr),(cpl_size)j);
                nerr++;
                continue;
            }
            freetfits(testt);
        }
    }
    return(nerr);
}

/*---------------------------------------------------------------------------*/
/**
    \par Name:
        vimos_testrdgn
    \par Purpose:
        Test the read/gain table
    \par Description:
        The read/gain table is tested to make sure it's readable and to
        make sure there is one row for each detector. No testing is
        done to the actual values in the rows.
    \par Language:
        C
    \param fr
        An example input frame that can be used to find the extension names
    \param readgain
        The read/gain table
    \return
        The number of errors
    \author
        Jim Lewis, CASU
 */
/*---------------------------------------------------------------------------*/

int vimos_testrdgn(cpl_frame *fr, cpl_frame *readgain) {
    int i,nerr;
    cpl_table *rdgn;
    casu_fits *f;
    const char *fctid = "vimos_testrdgn";

    /* Can we load the table? */

    nerr = 0;
    rdgn = cpl_table_load(cpl_frame_get_filename(readgain),1,0);
    if (cpl_error_get_code() != 0) {
        cpl_msg_error(fctid,"Read/gain table %s[1] won't load\n",
                      cpl_frame_get_filename(fr));
        nerr++;
        return(nerr);
    }

    /* Loop for each extension and make sure there is an entry for each one */

    for (i = 1; i <= VIMOS_NEXTN; i++) {
        cpl_table_unselect_all(rdgn);
        f = casu_fits_load(fr,CPL_TYPE_UNSPECIFIED,i);
        cpl_table_or_selected_string(rdgn,"EXTNAME",CPL_EQUAL_TO,
                                     casu_fits_get_extname(f));
        if (cpl_table_count_selected(rdgn) != 1) {
            cpl_msg_error(fctid,"No read/gain entry for %s",
                          casu_fits_get_extname(f));
            nerr++;
        }
        freefits(f);
    }
    cpl_table_delete(rdgn);
    return(nerr);
}

/*---------------------------------------------------------------------------*/
/**
    \par Name:
        vimos_getrdgn
    \par Purpose:
        Get the readnoise and gain estimate for the current detector
    \par Description:
        The read/gain table is searched for the current detector name.
        The readnoise and gain estimates are returned for that detector.
    \par Language:
        C
    \param readgain
        Input frame with the master read/gain table
    \param extname
        The extension name for which we want read/gain info
    \param readnoise
        The output readnoise in ADUs
    \param gain
        The output gain in e-/ADUs
    \author
        Jim Lewis, CASU
 */
/*---------------------------------------------------------------------------*/

void vimos_getrdgn(cpl_frame *readgain, char *extname, float *readnoise,
                   float *gain) {
    cpl_table *rgtab,*subset;
    int isnull;

    /* Read the input table */

    rgtab = cpl_table_load(cpl_frame_get_filename(readgain),1,0);

    /* Search the table for the correct row. NB: the read/gain table has
       already been tested so there should be no error here */

    cpl_table_unselect_all(rgtab);
    cpl_table_or_selected_string(rgtab,"EXTNAME",CPL_EQUAL_TO,extname);
    subset = cpl_table_extract_selected(rgtab);
    *readnoise = cpl_table_get_float(subset,"READNOISE",0,&isnull);
    *gain = cpl_table_get_float(subset,"GAIN",0,&isnull);
    cpl_table_delete(subset);
    cpl_table_delete(rgtab);
}

/*---------------------------------------------------------------------------*/
/**
    \par Name:
        vimos_get_groupsof4
    \par Purpose:
        Divide a frameset of VIMOS images into a 'group of 4' structure.
    \par Description:
        The input frames in a frameset are divided by their observation
        times into groups of 4 image that belong to the same exposure.
    \par Language:
        C
    \param fset
        Input frameset
    \param gr
        The output 'group of 4' array 
    \param ng
        The output number of groups
    \param vimos_names
        The set of names for the detectors
    \author
        Jim Lewis, CASU
 */
/*---------------------------------------------------------------------------*/

void vimos_get_groupsof4(cpl_frameset *fset, char *vimos_names[],
                         groupof4 **gr, int *ng) {
    int i,nr,j,k,nn;
    cpl_size *labels,nlab;
    cpl_frameset *res;
    cpl_propertylist *p;
    cpl_frame *fr;
    char name[16];

    /* Set some things up in case of silly input */

    if (fset == NULL) {
        *gr = NULL;
        *ng = 0;
        return;
    }

    /* Label the frameset by mjd and get memory for the groups of 4*/

    labels = cpl_frameset_labelise(fset,vimos_compare_lsts,&nlab);
    *ng = (int)nlab;
    *gr = cpl_malloc((*ng)*sizeof(groupof4));

    /* Loop for each label and extract the frameset */

    for (i = 0; i < *ng; i++) {
        res = cpl_frameset_extract(fset,labels,(cpl_size)i);
        nr = (int)cpl_frameset_get_size(res);

        /* Initialise a few things for this group */

        for (j = 0; j < 4; j++) 
            (*gr)[i].inf[j] = NULL;

        /* Loop for each member of the group. Get the INS part of the
           header and store it away. */

        for (j = 0; j < nr; j++) {
            fr = cpl_frameset_get_position(res,j);
            nn = cpl_frame_get_nextensions(fr);
            p = cpl_propertylist_load(cpl_frame_get_filename(fr),nn);

            /* Extract the date-obs of the first one for a name */

            if (j == 0) 
                strcpy((*gr)[i].name,cpl_propertylist_get_string(p,"DATE-OBS"));

            /* See which detector this belongs to */

            vimos_pfits_get_chipname(p,name);
            for (k = 0; k < VIMOS_NEXTN; k++) {
                if (! strcmp(vimos_names[k],name)) {
                    (*gr)[i].inf[k] = cpl_frame_duplicate(fr);
                    break;
                }
            }
            cpl_propertylist_delete(p);
        }                    
        
        /* Tidy up */

        cpl_frameset_delete(res);
    }

    /* And tidy up a little more */
    
    cpl_free(labels);
}

/*---------------------------------------------------------------------------*/
/**
    \par Name:
        vimos_free_groupsof4
    \par Purpose:
        Free the workspace in a groups of 4 structure
    \par Description:
        Free the workspace in a groups of 4 structure
    \par Language:
        C
    \param gr
        The 'group of 4' array 
    \param ng
        The number of groups
    \author
        Jim Lewis, CASU
 */
/*---------------------------------------------------------------------------*/

void vimos_free_groupsof4(int *ngr4, groupof4 **gr) {
    int i,j;

    for (i = 0; i < *ngr4; i++) {
        for (j = 0; j < 4; j++) 
            freeframe((*gr)[i].inf[j]);
    }
    freespace(*gr);
    *ngr4 = 0;
}

/*---------------------------------------------------------------------------*/
/**
    \par Name:
        vimos_copywcs
    \par Purpose:
        Copy a WCS from one propertylist to another
    \par Description:
        Copy a full WCS from the first propertylist into the second one
    \par Language:
        C
    \param p1
        The first input propertylist. This is the source of the WCS
    \param p2
        The second input propertylist. The WCS will be copied to this
    \return
        Nothing
    \author
        Jim Lewis, CASU
 */
/*---------------------------------------------------------------------------*/

void vimos_copywcs(cpl_propertylist *p1, cpl_propertylist *p2) {
    (void)cpl_propertylist_erase_regexp(p2,CPL_WCS_REGEXP,0);
    cpl_propertylist_copy_property_regexp(p2,p1,CPL_WCS_REGEXP,0);
}

/**@}*/

/*

$Log: vimos_utils.c,v $
Revision 1.18  2015/11/27 12:15:33  jim
detabbed some lines

Revision 1.17  2015/09/15 10:51:10  jim
Fixed superficial issue

Revision 1.16  2015/09/15 10:39:05  jim
added vimos_wcscopy

Revision 1.15  2015/08/07 13:07:15  jim
Fixed copyright to ESO

Revision 1.14  2015/06/25 11:54:07  jim
Fixed includes to casu

Revision 1.13  2015/05/13 12:01:18  jim
Fixed _wcsfit routine so that even if the WCS fit fails, the matched standards
catalogue has the correct number of columns

Revision 1.12  2015/05/01 09:02:22  jim
Modified to pass detector name information differently so that Mac compiler
doesn't gripe

Revision 1.11  2015/03/13 10:22:59  jim
Fixed some compiler warnings

Revision 1.10  2015/03/10 12:45:53  jim
More leaks

Revision 1.9  2015/02/14 12:37:22  jim
Expanded and modified vimos_wcsfit

Revision 1.8  2015/01/29 12:07:20  jim
Removed some support routines from recipes to here. Modified comments

Revision 1.7  2014/12/11 12:27:28  jim
new version

Revision 1.6  2014/05/06 07:37:41  jim
superficial fixes

Revision 1.5  2013/11/21 12:49:18  jim
Fixed to have separate trim regions for old and new chipsets

Revision 1.4  2013/11/21 09:39:08  jim
detabbed

Revision 1.3  2013/11/04 06:03:27  jim
Fixed some comments

Revision 1.2  2013/10/25 09:32:29  jim
fixed some docs

Revision 1.1.1.1  2013-10-24 08:36:12  jim
Initial import


*/