/* $Id: vimos_ima_science.c,v 1.45 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.45 $
 * $Name:  $
 */

/* Includes */

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

#include <stdio.h>
#include <libgen.h>
#include <string.h>

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

#include "vmutils.h"
#include "vimos_imaging_utils.h"
#include "vimos_pfits.h"
#include "vimos_dfs.h"
#include "vimos_mods.h"
#include "vimos_var.h"
#include "casu_utils.h"
#include "casu_stats.h"
#include "casu_fits.h"
#include "casu_mask.h"
#include "casu_wcsutils.h"
#include "casu_mods.h"

/* Structure definitions */

typedef struct {

    /* Input */

    int         preview_only;
    int         minphotom;
    int         savecat;
    int         savemstd;
    int         prettynames;
    int         chop_crud;
    int         cdssearch_astrom;
    int         cdssearch_photom;
    int         ignore_fringe;
    char        *cacheloc;
    float       magerrcut;

    /* Stacking parameters */

    float       stk_lthr;
    float       stk_hthr;
    int         stk_method;
    int         stk_seeing;
    int         stk_fast;
    int         stk_nfst;

    /* Source catalogue extraction parameters for stacks*/

    int         stk_cat_ipix;
    float       stk_cat_thresh;
    int         stk_cat_icrowd;
    float       stk_cat_rcore;
    int         stk_cat_nbsize;

} configstruct;

typedef struct {
    cpl_size         *labels;
    groupof4         *gr4;
    int              ngr4;
    cpl_frame        *master_bias;
    cpl_frame        *master_dark;
    cpl_frame        *master_twilight_flat;
    cpl_frame        *master_conf;
    cpl_frame        *master_mstd_phot;
    cpl_frame        *master_fringe;
    cpl_frame        *master_fringe_var;
    cpl_frame        *master_stdstar;
    casu_mask        *mask;
    cpl_frame        *phottab;
    cpl_table        *tphottab;
    cpl_frameset     *science_frames;
    cpl_frame        **product_frames_simple;
    cpl_frame        **product_frames_simple_var;
    cpl_frame        **product_frames_simple_cat;
    cpl_frame        **product_frames_simple_mstd_a;
    cpl_frame        *product_frame_stack;
    cpl_frame        *product_frame_stack_conf;
    cpl_frame        *product_frame_stack_var;
    cpl_frame        *product_frame_stack_cat;
    cpl_frame        *product_frame_stack_mstd_a;
    cpl_frame        *product_frame_stack_mstd_p;
    float            *gaincors;
    char             *catpath_a;
    char             *catname_a;
    char             *catpath_p;
    char             *catname_p;
    cpl_frame        *schlf_n;
    cpl_frame        *schlf_s;
    cpl_frame        *readgain;
    cpl_frameset     *prod_im;
    cpl_frameset     *prod_tab;

    /* Level 1 stuff */

    casu_fits        *fbias;
    casu_fits        *fdark;
    casu_fits        *fflat;
    casu_fits        *fconf;
    casu_fits        *fconf2;
    casu_fits        *ffringe;
    casu_fits        *ffringe_var;

    casu_fits        **sci_fits;
    casu_fits        **sci_vars;
    casu_tfits       **sci_cats;
    casu_tfits       **sci_mstd_a;

} memstruct;

/* List of data products types */

enum {SIMPLE_FILE,
      SIMPLE_VAR,
      SIMPLE_CAT,
      STACK_FILE,
      STACK_CONF,
      STACK_VAR,
      STACK_CAT,
      SIMPLE_MSTD_A,
      STACK_MSTD_A,
      STACK_MSTD_P
};

/* Recipe name */

#define VIMOS_RECIPENAME "vimos_ima_science"

/* Function prototypes */

static int vimos_ima_science_create(cpl_plugin *);
static int vimos_ima_science_exec(cpl_plugin *);
static int vimos_ima_science_destroy(cpl_plugin *);
static int vimos_ima_science(cpl_parameterlist *, cpl_frameset *);

static int vimos_sci_phot_stdstar(cpl_frame *stdstar, casu_fits **stack,
                                   casu_tfits **stackcat,
                                   const char *airmasskey);

static int vimos_sci_save_simple(casu_fits *obj, cpl_frameset *framelist,
                                 cpl_parameterlist *parlist, 
                                 cpl_frame *template, int isfirst,
                                 const char *tag, char *fname, char *cname, 
                                 char *assoc, cpl_frame **product_frame);
static int vimos_sci_save_stack(casu_fits *stack, cpl_frameset *framelist,
                                cpl_parameterlist *parlist, 
                                cpl_frame *template, int fnametype, char *cname,
                                int prodtype, char *assoc[], char *photosys,
                                cpl_frame **product_frame);
static int vimos_sci_save_stack_conf(casu_fits *stack, cpl_frameset *framelist,
                                     cpl_parameterlist *parlist, 
                                     cpl_frame *template, int fnametype,
                                     char *cname, cpl_frame **product_frame);
static int vimos_sci_save_cat(casu_tfits *scat, cpl_frameset *framelist,
                              cpl_parameterlist *parlist, 
                              cpl_frame *template, int fnametype, 
                              int ptype, int fnumber, char *cname, 
                              char *photosys, cpl_frame **product_frame);
static void vimos_sci_product_name(char *template, int producttype, 
                                   int nametype, int fnumber, char *outfname);
static void vimos_sci_update_hdr(cpl_frameset *simple, cpl_frame *stack, 
                                  int typef);
static void vimos_sci_update_hdr_stack(casu_fits *stack[], casu_fits *conf,
                                       casu_fits *var, casu_tfits *cat);
static void vimos_sci_init(memstruct *ps);
static void vimos_sci_tidy(memstruct *ps, int level);

static char vimos_ima_science_description[] =
"vimos_ima_science -- VIMOS science recipe.\n\n"
"Take a list of science frames and correct them to remove bias, dark\n"
"and flat field signatures. Optionally defringe the reddest filter\n"
"images. Stack jitter sequences. Astrometrically and photometrically\n"
"calibrate the stacks\n\n"
"The program accepts the following files in the SOF:\n\n"
"    Tag                   Description\n"
"    -----------------------------------------------------------------------\n"
"    %-21s A list of raw science images\n"
"    %-21s A master bias frame\n"
"    %-21s A master dark frame\n"
"    %-21s A master flat frame\n"
"    %-21s A master confidence map\n"
"    %-21s A master fringe map (optional)\n"
"    %-21s A master fringe variance map (optional)\n"
"    %-21s A master readgain table\n"
"    %-21s A master 2MASS index for astrometry or\n"
"    %-21s A master PPMXL index for astrometry or\n"
"    %-21s A master APASS index for astrometry or\n"
"    %-21s A master local index for astrometry\n"
"    %-21s A master APASS index for photometry or\n"
"    %-21s A master local index for photometry or\n"
"    %-21s A master processed standard star image or table\n"
"    %-21s A master photometric matched standards table (optional)\n"
"    %-21s Northern Schlegel Map\n"
"    %-21s Southern Schlegel Map\n"
"    %-21s A photometric calibration table\n"
"All of these are required unless specified otherwise."
"\n";

/**@{*/

/**
    \ingroup recipelist
    \defgroup vimos_ima_science vimos_ima_science
    \brief VIMOS science recipe

    \par Name: 
        vimos_ima_science
    \par Purpose: 
        Reduce a full pawprint of VIMOS data
    \par Description: 
        A single pawprint of VIMOS data in a single filter is corrected 
        for instrumental signature. The images are stacked and the stacked 
        image is astrometrically calibrated using a standard astrometric
        catalogue and photometrically calibrated using either the APASS 
        catalogue or a local standard file.
    \par Language:
        C
    \par Parameters:
        - General reduction parameters:
            - \b savecat (int): If set, then catalogues for the simple images are saved
            - \b savemstd (int): If set, then matched standard catalogues are saved
            - \b preview_only (bool): If set we only get a preview of the reduction
            - \b minphotom (int): The minimum number of standards for photometry
            - \b prettynames (bool): True if we're using nice file names       
            - \b chop_crud (string): Method for removing cruddy edges
                - none: Don't chop anything
                - lowconf_block: Chop blocks with low confidence
                - hardcoded: Chop hardcoded regions
                - lowconf_pix: Ignore all pixels with low confidence
                - hardconf_pix: hardcoded + lowconf_pix
            - \b cdssearch_astrom (string): Use CDS for astrometric standards?
                - none: Use no CDS catalogues. Use local catalogues instead
                - 2mass: 2MASS PSC
                - usnob: USNOB catalogue
                - ppmxl: PPMXL catalogue
                - wise: WISE catalogue
            - \b cdssearch_photom (string): Use CDS for photometric standards?
                - none: Use no CDS catalogues. Use local catalogues instead
                - apass: APASS catalogue
            - \b ignore_fringe (bool): Flag to ignore a provided fringe frame
            - \b cacheloc (string): A directory where we can put the standard
                 star cache
            - \b magerrcut (float): A cut in the magnitude errors for 
                 photometric standards
        - Parameters for source detection on stacks:
            - \b stk_cat_ipix (int): The minimum allowable size of an object
            - \b stk_cat_thresh (float): The detection threshold in sigma above sky
            - \b stk_cat_icrowd (int): If set then the deblending software will be used
            - \b stk_cat_rcore (float): The core radius in pixels
            - \b stk_cat_nbsize (int): The smoothing box size for background map estimation        
        - Parameters for stacking algorithm:
            - \b stk_lthr (float): Low rejection threshold
            - \b stk_hthr (float): Upper rejection threshold
            - \b stk_method (int): 0 -> nearest neighbour, 1 -> bilinear interpolation
            - \b stk_seeing (bool): TRUE if we want to weight images by seeing
            - \b stk_fast (string): flag for fast or slow stacking algorithm.
                 The fast algorithm is much greedier in memory.
                - auto: Let the recipe decide which we use
                - slow: Use the slow algorithm
                - fast: Use the fast algorithm
            - \b stk_nfst (int): The number of images above which we revert
                 to the slow stacking algorithm when in auto mode
    \par Input File Types:
        The following list gives the file types that can appear in the SOF
        file. The word in bold is the DO category value.
        - \b OBJECT (required): A list of raw science images for processing
        - \b MASTER_BIAS (required): A master bias frame. This is needed 
             for the bias correction of the input images.
        - \b MASTER_DARK (required): A master dark frame. This is needed 
             for the dark correction of the input images.
        - \b MASTER_TWILIGHT_FLAT (required): A master twilight flat frame
             appropriate for the filter used in the science observations.
             This is used to remove the flat field effects from the science
             frames and also for the detector level gain correction.
        - \b MASTER_FRINGE (optional): A master fringe map
        - \b MASTER_FRINGE_VAR (optional): A master fringe variance map
        - \b MASTER_READGAIN (required): A master read/gain table
        - \b MASTER_CONF (required): A master confidence map 
             appropriate for the filter used in the science observations.
             This is used for weighting and to mark bad pixels.
        - \b MASTER_2MASS_CATLAOGUE_ASTROM or \b MASTER_PPMXL_CATALOGUE_ASTROM
             or \b MASTER_APASS_CATALOGUE_ASTROM or \b MASTER_LOCCAT_ASTROM
             (required): A master astrometric standard star catalogue. This
             is not required if the --cdssearch_astrom switch is used.
        - \b MASTER_APASS_CATALOGUE_PHOTOM or \b MASTER_LOCCAT_PHOTOM
             (required): A master photometric standard star catalogue. This
             is not required if the --cdssearch_photom switch is used.
        - \b MASTER_STDSTAR (optional): A standard star processed  
             source catalogue of the same filter as the object images. If this
             exists in the sof, then it will be used as the basis for the
             photometric solution. If not, then a solution will be found from
             standard star catalogues as per above. 
        - \b MATCHSTD_PHOTOM (optional): If this is specified, then all the
             photometric calibration is derived from it rather than from a
             matching exercise to an external catalogue.
        - \b SCHLEGEL_MAP_NORTH (required): The Northern Schlegel dust map
        - \b SCHLEGEL_MAP_SOUTH (required): The Southern Schlegel dust map
    \par Output Products:
        The following list gives the output data products that are generated
        by this recipe. The word in bold gives the DPR CATG keyword value for
        each product:
        - The instrumentally corrected and astrometerically calibrated 
          individual science exposures, (\b BASIC_CALIBRATED_SCI), their
          corresponding variance maps (\b BASIC_VAR_MAP_SCI), if 
          requested, their source catalogues (\b OBJECT_CATALOGUE_SCI) and,
          if requested, their astrometric matched standards catalogues
          (MATCHSTD_ASTROM).
        - The stacked science pawprint (\b JITTERED_IMAGE_SCI).
        - The stacked science pawprint confidence map  (\b CONFIDENCE_MAP_SCI).
        - The stacked science pawprint variance map (\b JITTERED_VAR_IMAGE)
        - The stacked science pawprint source catalogue  
          (\b OBJECT_CATALOGUE_SCI_JITTERED).
        - If requested, the astrometric and photometric matched standards 
          catalogues for the stack (\b MATCHSTD_ASTROM and \b MATCHSTD_PHOTOM).
    \par Output QC Parameters:
        - \b SATURATION
             The saturation level in ADUs.
        - \b MEAN_SKY
             The mean level of the background sky over the image in ADUs.
        - \b SKY_NOISE
             The RMS of the background sky over the image in ADUs
        - \b IMAGE_SIZE
             The average size of stellar images on the image in pixels
        - \b ELLIPTICITY
             The average ellipticity of stellar images on the image
        - \b POSANG
             The average position angle in degrees from North towards East.
             NB: this value only makes sense if the ellipticity is significant
        - \b APERTURE_CORR
             The aperture correction for an aperture of radius rcore.
        - \b NOISE_OBJ
             The number of noise objects found in the image
        - \b MAGZPT
             The photometric zero point
        - \b MAGZERR
             The internal error in the photometric zero point
        - \b MAGNZPT
             The number of stars used to determine the magnitude zero point
        - \b MAGNCUT
             The number of stars cut from magnitude zero point calculation
        - \b SKYBRIGHT
             The sky brightness in mag/arcsec**2
        - \b LIMITING_MAG
             The limiting magnitude for this image for a 5 sigma detection.
        - \b WCS_DCRVAL1
             The offset of the equatorial coordinate represented by CRVAL1 
             from the raw frame to the reduced one (degrees).
        - \b WCS_DCRVAL2
             The offset of the equatorial coordinate represented by CRVAL2
             from the raw frame to the reduced one (degrees).
        - \b WCS_DTHETA
             The change in the WCS coordinate system rotation from the raw
             to the reduced frame (degrees)
        - \b WCS_SCALE
             The scale of the pixels in the reduced frames in arcseconds per
             pixel
        - \b WCS_SHEAR
             The shear of the astrometric solution in the form of the 
             difference between the rotation of the x solution and the rotation
             of the y solution (abs(xrot) - abs(yrot) in degrees)
        - \b WCS_RMS
             The average error in the WCS fit (arcsec)
    \par Notes
        None.
    \par Fatal Error Conditions:
        - NULL input frameset
        - Input frameset headers incorrect meaning that RAW and CALIB frame
          can't be distinguished
        - No science frames in the input frameset
        - No master bias, dark, flat or confidence map frame in input frameset
        - Unable to save data products
        - Unable to load images
    \par Non-Fatal Error Conditions:
        - None
    \par Conditions Leading To Dummy Products:
        - Master calibration images flagged as dummys
        - Detector for the current image is flagged as dead in all 
          science images
        - Failure of one of the processing steps
    \par Author:
        Jim Lewis, CASU
    \par Code Reference: 
        vimos_ima_science.c
*/

/* Function code */

/*---------------------------------------------------------------------------*/
/**
  @brief    Build the list of available plugins, for this module.
  @param    list    the plugin list
  @return   0 if everything is ok

  This function is exported.
 */
/*---------------------------------------------------------------------------*/

int cpl_plugin_get_info(cpl_pluginlist *list) {
    cpl_recipe  *recipe = cpl_calloc(1,sizeof(*recipe));
    cpl_plugin  *plugin = &recipe->interface;
    char alldesc[SZ_ALLDESC];
    (void)snprintf(alldesc,SZ_ALLDESC,vimos_ima_science_description,
                   VIMOS_SCI_OBJECT_RAW,VIMOS_CAL_BIAS,VIMOS_CAL_DARK,
                   VIMOS_CAL_TWILIGHT_FLAT,VIMOS_CAL_CONF,VIMOS_CAL_FRINGE,
                   VIMOS_CAL_FRINGE_VAR,VIMOS_CAL_READGAIN,VIMOS_CAL_2MASS_A,
                   VIMOS_CAL_PPMXL_A,VIMOS_CAL_APASS_A,VIMOS_CAL_LOCCAT_A,
                   VIMOS_CAL_APASS_P,VIMOS_CAL_LOCCAT_P,VIMOS_CAL_STDSTAR,
                   VIMOS_CAL_MSTD_PHOT,VIMOS_CAL_SCHL_N,VIMOS_CAL_SCHL_S,
                   VIMOS_CAL_PHOTTAB);
    cpl_plugin_init(plugin,
                    CPL_PLUGIN_API,
                    VIMOS_BINARY_VERSION,
                    CPL_PLUGIN_TYPE_RECIPE,
                    "vimos_ima_science",
                    "Science processing for imaging",
                    alldesc,
                    "Jim Lewis",
                    "jrl@ast.cam.ac.uk",
                    vimos_get_license(),
                    vimos_ima_science_create,
                    vimos_ima_science_exec,
                    vimos_ima_science_destroy);

    cpl_pluginlist_append(list,plugin);

    return(0);
}

/*---------------------------------------------------------------------------*/
/**
  @brief    Setup the recipe options
  @param    plugin  the plugin
  @return   0 if everything is ok

  Create the recipe instance and make it available to the application using the
  interface.
 */
/*---------------------------------------------------------------------------*/

static int vimos_ima_science_create(cpl_plugin *plugin) {
    cpl_recipe      *recipe;
    cpl_parameter   *p;

    /* Get the recipe out of the plugin */

    if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE)
        recipe = (cpl_recipe *)plugin;
    else 
        return(-1);

    /* Create the parameters list in the cpl_recipe object */

    recipe->parameters = cpl_parameterlist_new();

    /* Fill in flag to use save the output catalogue */

    p = cpl_parameter_new_value("vimos.vimos_ima_science.savecat",
                                CPL_TYPE_BOOL,"Save catalogue?",
                                "vimos.vimos_ima_science",FALSE);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"savecat");
    cpl_parameterlist_append(recipe->parameters,p);

    /* Fill in flag to use save the matched standard catalogues */

    p = cpl_parameter_new_value("vimos.vimos_ima_science.savemstd",
                                CPL_TYPE_BOOL,
                                "Save matched standard catalogues?",
                                "vimos.vimos_ima_science",FALSE);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"savemstd");
    cpl_parameterlist_append(recipe->parameters,p);

    /* Fill in flag to just print out how the data will be processed and 
       then exit */

    p = cpl_parameter_new_value("vimos.vimos_ima_science.preview_only",
                                CPL_TYPE_BOOL,"Preview only?",
                                "vimos.vimos_ima_science",FALSE);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"preview_only");
    cpl_parameterlist_append(recipe->parameters,p);

    /* Fill in flag to just print out how the data will be processed and 
       then exit */

    p = cpl_parameter_new_range("vimos.vimos_ima_science.minphotom",
                                CPL_TYPE_INT,
                                "Minimum stars for photometry solution",
                                "vimos.vimos_ima_science",1,1,100000);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"minphotom");
    cpl_parameterlist_append(recipe->parameters,p);

    /* Fill in flag to decide whether use predictable names or nice
       names based on input file names */

    p = cpl_parameter_new_value("vimos.vimos_ima_science.prettynames",
                                CPL_TYPE_BOOL,"Use pretty product names?",
                                "vimos.vimos_ima_science",FALSE);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"prettynames");
    cpl_parameterlist_append(recipe->parameters,p);

    /* Fill in flag to decide how or whether we chop off cruddy bits 
       around the outside of each detector */

    p = cpl_parameter_new_enum("vimos.vimos_ima_science.chop_crud",
                               CPL_TYPE_STRING,"Chop crud method",
                               "vimos.vimos_ima_science",
                               "hardconf_pix",5,"none","lowconf_block",
                               "hardcoded","lowconf_pix","hardconf_pix");
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"chop_crud");
    cpl_parameterlist_append(recipe->parameters,p);

    /* Flag to decide how we get the astrometric standard star information. 
       If "none", then use the local catalogues specified in the sof. If not 
       "none", then use one of the selection of catalogues available from CDS */

    p = cpl_parameter_new_enum("vimos.vimos_ima_science.cdssearch_astrom",
                               CPL_TYPE_STRING,
                               "CDS astrometric catalogue",
                               "vimos.vimos_ima_science",
                               "none",5,"none","2mass","usnob","ppmxl","wise");
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"cdssearch_astrom");
    cpl_parameterlist_append(recipe->parameters,p);

    /* The same for photometric standards */

    p = cpl_parameter_new_enum("vimos.vimos_ima_science.cdssearch_photom",
                               CPL_TYPE_STRING,
                               "CDS photometric catalogue",
                               "vimos.vimos_ima_science",
                               "none",2,"none","apass");
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"cdssearch_photom");
    cpl_parameterlist_append(recipe->parameters,p);

    /* Flag to ignore the provided fringe frame */

    p = cpl_parameter_new_value("vimos.vimos_ima_science.ignore_fringe",
                                CPL_TYPE_BOOL,"Ignore provided fringe frame?",
                                "vimos.vimos_ima_science",FALSE);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"ignore_fringe");
    cpl_parameterlist_append(recipe->parameters,p);

    /* ***Special parameters for object detection on stacks*** */

    /* Fill in the minimum object size */

    p = cpl_parameter_new_range("vimos.vimos_ima_science.stk_cat_ipix",
                                CPL_TYPE_INT,
                                "Minimum pixel area for each detected object",
                                "vimos.vimos_ima_science",10,1,100000);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"stk_cat_ipix");
    cpl_parameterlist_append(recipe->parameters,p);

    /* Fill in the detection threshold parameter */

    p = cpl_parameter_new_range("vimos.vimos_ima_science.stk_cat_thresh",
                                CPL_TYPE_DOUBLE,
                                "Detection threshold in sigma above sky",
                                "vimos.vimos_ima_science",1.5,1.0e-6,
                                1.0e10);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"stk_cat_thresh");
    cpl_parameterlist_append(recipe->parameters,p);

    /* Fill in flag to use deblending software or not */

    p = cpl_parameter_new_value("vimos.vimos_ima_science.stk_cat_icrowd",
                                CPL_TYPE_BOOL,"Use deblending?",
                                "vimos.vimos_ima_science",TRUE);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"stk_cat_icrowd");
    cpl_parameterlist_append(recipe->parameters,p);

    /* Fill in core radius */

    p = cpl_parameter_new_range("vimos.vimos_ima_science.stk_cat_rcore",
                                CPL_TYPE_DOUBLE,"Value of Rcore in pixels",
                                "vimos.vimos_ima_science",10.0,1.0e-6,
                                1024.0);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"stk_cat_rcore");
    cpl_parameterlist_append(recipe->parameters,p);

    /* Fill in background smoothing box size */

    p = cpl_parameter_new_range("vimos.vimos_ima_science.stk_cat_nbsize",
                                CPL_TYPE_INT,"Background smoothing box size",
                                "vimos.vimos_ima_science",128,1,2048);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"stk_cat_nbsize");
    cpl_parameterlist_append(recipe->parameters,p);

    /* ***Special parameters for stacking*** */

    /* The lower rejection threshold to be used during stacking */

    p = cpl_parameter_new_range("vimos.vimos_ima_science.stk_lthr",
                                CPL_TYPE_DOUBLE,"Low rejection threshold",
                                "vimos.vimos_ima_science",3.0,1.0e-6,
                                1.0e10);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"stk_lthr");
    cpl_parameterlist_append(recipe->parameters,p);

    /* The upper rejection threshold to be used during stacking */

    p = cpl_parameter_new_range("vimos.vimos_ima_science.stk_hthr",
                                CPL_TYPE_DOUBLE,"Upper rejection threshold",
                                "vimos.vimos_ima_science",3.0,1.0e-6,
                                1.0e10);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"stk_hthr");
    cpl_parameterlist_append(recipe->parameters,p);

    /* Stacking method */

    p = cpl_parameter_new_enum("vimos.vimos_ima_science.stk_method",
                               CPL_TYPE_STRING,"Stacking method",
                               "vimos.vimos_ima_science",
                               "linear",2,"nearest","linear");
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"stk_method");
    cpl_parameterlist_append(recipe->parameters,p);

    /* Flag to use seeing weighting */

    p = cpl_parameter_new_value("vimos.vimos_ima_science.stk_seeing",
                                CPL_TYPE_BOOL,"Weight by seeing?",
                                "vimos.vimos_ima_science",FALSE);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"stk_seeing");
    cpl_parameterlist_append(recipe->parameters,p);    

    /* Stacking speed method */

    p = cpl_parameter_new_enum("vimos.vimos_ima_science.stk_fast",
                               CPL_TYPE_STRING,"Stack using fast algorithm?",
                               "vimos.vimos_ima_science",
                               "auto",3,"fast","slow","auto");
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"stk_fast");
    cpl_parameterlist_append(recipe->parameters,p);

    /* Get the maximum number of frames we can stack using the fast
       algorithm */

    p = cpl_parameter_new_value("vimos.vimos_ima_science.stk_nfst",
                                CPL_TYPE_INT,
                                "Nframes to stack in fast mode",
                                "vimos.vimos_ima_science",16);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"stk_nfst");
    cpl_parameterlist_append(recipe->parameters,p);    

    /* Location for standard star cache */

    p = cpl_parameter_new_value("vimos.vimos_ima_science.cacheloc",
                                CPL_TYPE_STRING,
                                "Location for standard star cache",
                                "vimos.vimos_ima_science",".");
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"cacheloc");
    cpl_parameterlist_append(recipe->parameters,p);        

    /* Magnitude error cut */

    p = cpl_parameter_new_value("vimos.vimos_ima_science.magerrcut",
                                CPL_TYPE_DOUBLE,
                                "Magnitude error cut",
                                "vimos.vimos_science_magerrcut",100.0);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"magerrcut");
    cpl_parameterlist_append(recipe->parameters,p);        

    /* Get out of here */

    return(0);
}
        
/*---------------------------------------------------------------------------*/
/**
  @brief    Execute the plugin instance given by the interface
  @param    plugin  the plugin
  @return   0 if everything is ok
 */
/*---------------------------------------------------------------------------*/

static int vimos_ima_science_exec(cpl_plugin *plugin) {
    cpl_recipe  *recipe;

    /* Get the recipe out of the plugin */

    if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE)
        recipe = (cpl_recipe *)plugin;
    else 
        return(-1);

    return(vimos_ima_science(recipe->parameters,recipe->frames));
}
                                
/*---------------------------------------------------------------------------*/
/**
  @brief    Destroy what has been created by the 'create' function
  @param    plugin  the plugin
  @return   0 if everything is ok
 */
/*---------------------------------------------------------------------------*/

static int vimos_ima_science_destroy(cpl_plugin *plugin) {
    cpl_recipe *recipe ;

    /* Get the recipe out of the plugin */

    if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE)
        recipe = (cpl_recipe *)plugin;
    else 
        return(-1);

    cpl_parameterlist_delete(recipe->parameters);
    return(0);
}

/*---------------------------------------------------------------------------*/
/**
  @brief    The recipe data reduction part is implemented here 
  @param    parlist     the parameters list
  @param    framelist   the frames list
  @return   0 if everything is ok
 */
/*---------------------------------------------------------------------------*/

static int vimos_ima_science(cpl_parameterlist *parlist, 
                                 cpl_frameset *framelist) {
    const char *fctid="vimos_ima_science";
    configstruct cs;
    memstruct ps;
    cpl_parameter *p;
    cpl_size nlab;
    int nfail,status,i,j,isfirst,njsteps,n,ngood,fastslow,match,k;
    float gaincor_fac,gain,readnoise,fac;
    cpl_frame *catindex_a,*catindex_p,*frm,*template;
    cpl_frameset *chiplist;
    cpl_propertylist *pp,*epp;
    casu_fits *ff,*out,*outc,*outv,**good,**goodv,*ffv,*stack[VIMOS_NEXTN];
    casu_fits *stackv[VIMOS_NEXTN];
    casu_fits *stackconf[VIMOS_NEXTN];
    casu_tfits *tab,*stackcat[VIMOS_NEXTN],*mstds_a[VIMOS_NEXTN];
    casu_tfits *mstds_p[VIMOS_NEXTN];
    cpl_image *im,*imv;
    cpl_table *t,*stdscat;
    char filt[16],projid[16],*fname,bname[BUFSIZ],*junk1,*junk2,pcat[32];
    char *vimos_names[VIMOS_NEXTN],*assoc[2],photosys[8];

    /* Check validity of input frameset */

    if (framelist == NULL || cpl_frameset_get_size(framelist) <= 0) {
        cpl_msg_error(fctid,"Input framelist NULL or has no input data");
        return(-1);
    }

    /* Initialise some things */

    vimos_sci_init(&ps);

    /* Get the parameters */

    p = cpl_parameterlist_find(parlist,"vimos.vimos_ima_science.savecat");
    cs.savecat = cpl_parameter_get_bool(p);
    p = cpl_parameterlist_find(parlist,
                               "vimos.vimos_ima_science.savemstd");
    cs.savemstd = cpl_parameter_get_bool(p);
    p = cpl_parameterlist_find(parlist,
                               "vimos.vimos_ima_science.preview_only");
    cs.preview_only = cpl_parameter_get_bool(p);
    p = cpl_parameterlist_find(parlist,"vimos.vimos_ima_science.minphotom");
    cs.minphotom = cpl_parameter_get_int(p);
    p = cpl_parameterlist_find(parlist,
                               "vimos.vimos_ima_science.prettynames");
    cs.prettynames = (cpl_parameter_get_bool(p) ? 1 : 0);
    p = cpl_parameterlist_find(parlist,
                               "vimos.vimos_ima_science.chop_crud");
    if (! strcmp(cpl_parameter_get_string(p),"none")) 
        cs.chop_crud = 0;
    else if (! strcmp(cpl_parameter_get_string(p),"lowconf_block")) 
        cs.chop_crud = 1;
    else if (! strcmp(cpl_parameter_get_string(p),"hardcoded")) 
        cs.chop_crud = 2;
    else if (! strcmp(cpl_parameter_get_string(p),"lowconf_pix")) 
        cs.chop_crud = 3;
    else if (! strcmp(cpl_parameter_get_string(p),"hardconf_pix")) 
        cs.chop_crud = 4;
    p = cpl_parameterlist_find(parlist,
                               "vimos.vimos_ima_science.cdssearch_astrom");
    if (! strcmp(cpl_parameter_get_string(p),"none")) 
        cs.cdssearch_astrom = 0;
    else if (! strcmp(cpl_parameter_get_string(p),"2mass")) 
        cs.cdssearch_astrom = 1;
    else if (! strcmp(cpl_parameter_get_string(p),"usnob")) 
        cs.cdssearch_astrom = 2;
    else if (! strcmp(cpl_parameter_get_string(p),"ppmxl")) 
        cs.cdssearch_astrom = 3;
    else if (! strcmp(cpl_parameter_get_string(p),"wise")) 
        cs.cdssearch_astrom = 5;
    p = cpl_parameterlist_find(parlist,
                               "vimos.vimos_ima_science.cdssearch_photom");
    if (! strcmp(cpl_parameter_get_string(p),"none")) 
        cs.cdssearch_photom = 0;
    else if (! strcmp(cpl_parameter_get_string(p),"apass")) 
        cs.cdssearch_photom = 6;
    p = cpl_parameterlist_find(parlist,
                               "vimos.vimos_ima_science.ignore_fringe");
    cs.ignore_fringe = cpl_parameter_get_bool(p);

    /* Stacking special parameters */

    p = cpl_parameterlist_find(parlist,"vimos.vimos_ima_science.stk_lthr");
    cs.stk_lthr = cpl_parameter_get_double(p);
    p = cpl_parameterlist_find(parlist,"vimos.vimos_ima_science.stk_hthr");
    cs.stk_hthr = cpl_parameter_get_double(p);
    p = cpl_parameterlist_find(parlist,
                               "vimos.vimos_ima_science.stk_method");
    cs.stk_method = (! strcmp(cpl_parameter_get_string(p),"nearest") ? 0 : 1);
    p = cpl_parameterlist_find(parlist,
                               "vimos.vimos_ima_science.stk_seeing");
    cs.stk_seeing = cpl_parameter_get_bool(p);
    p = cpl_parameterlist_find(parlist,
                               "vimos.vimos_ima_science.stk_fast");
    if (! strcmp(cpl_parameter_get_string(p),"auto"))
        cs.stk_fast = -1;
    else if (! strcmp(cpl_parameter_get_string(p),"slow"))
        cs.stk_fast = 0;
    else if (! strcmp(cpl_parameter_get_string(p),"fast"))
        cs.stk_fast = 1;
    p = cpl_parameterlist_find(parlist,
                               "vimos.vimos_ima_science.stk_nfst");
    cs.stk_nfst = cpl_parameter_get_int(p);

    /* Object detection on stacks special parameters */

    p = cpl_parameterlist_find(parlist,
                               "vimos.vimos_ima_science.stk_cat_ipix");
    cs.stk_cat_ipix = cpl_parameter_get_int(p);
    p = cpl_parameterlist_find(parlist,
                               "vimos.vimos_ima_science.stk_cat_thresh");
    cs.stk_cat_thresh = (float)cpl_parameter_get_double(p);
    p = cpl_parameterlist_find(parlist,
                               "vimos.vimos_ima_science.stk_cat_icrowd");
    cs.stk_cat_icrowd = cpl_parameter_get_bool(p);
    p = cpl_parameterlist_find(parlist,
                               "vimos.vimos_ima_science.stk_cat_rcore");
    cs.stk_cat_rcore = (float)cpl_parameter_get_double(p);
    p = cpl_parameterlist_find(parlist,
                               "vimos.vimos_ima_science.stk_cat_nbsize");
    cs.stk_cat_nbsize = cpl_parameter_get_int(p);

    /* Cache location */

    p = cpl_parameterlist_find(parlist,
                               "vimos.vimos_ima_science.cacheloc");
    cs.cacheloc = (char *)cpl_parameter_get_string(p);

    /* Magnitude error cut */

    p = cpl_parameterlist_find(parlist,
                               "vimos.vimos_ima_science.magerrcut");
    cs.magerrcut = (float)cpl_parameter_get_double(p);

    /* Sort out raw from calib frames */

    if (vimos_dfs_set_groups(framelist) != CASU_OK) {
        cpl_msg_error(fctid,"Cannot identify RAW and CALIB frames");
        vimos_sci_tidy(&ps,0);
        return(-1);
    }

    /* Label the input frames */

    if ((ps.labels = cpl_frameset_labelise(framelist,casu_compare_tags,
                                           &nlab)) == NULL) {
        cpl_msg_error(fctid,"Cannot labelise the input frames");
        vimos_sci_tidy(&ps,0);
        return(-1);
    }

    /* Get the input science frames */

    nfail = 0;
    if ((ps.science_frames = 
         casu_frameset_subgroup(framelist,ps.labels,nlab,
                                VIMOS_SCI_OBJECT_RAW)) == NULL) {
        cpl_msg_error(fctid,"No science images to process!");
        vimos_sci_tidy(&ps,0);
        return(-1);
    }
    nfail += vimos_testfrms(ps.science_frames,1,1,1);

    /* See what chip set this is */

    if (vimos_load_names(cpl_frameset_get_position(ps.science_frames,0),vimos_names) != CASU_OK) {
        cpl_msg_error(fctid,"Cannot get the name set");
        vimos_sci_tidy(&ps,0);
        return(-1);
    }

    /* Now divide the science images up so that we have groups of 4
       images per exposure */

    vimos_get_groupsof4(ps.science_frames,vimos_names,&(ps.gr4),&(ps.ngr4));

    /* Check to see if there is a master bias frame */

    if ((ps.master_bias = 
         casu_frameset_subgroup_1(framelist,ps.labels,nlab,
                                  VIMOS_CAL_BIAS)) == NULL) {
        cpl_msg_error(fctid,"No master bias found");
        vimos_sci_tidy(&ps,0);
        return(-1);
    }
    nfail += vimos_testfrm_1(ps.master_dark,VIMOS_NEXTN,1,0);
        
    /* Check to see if there is a master dark frame */

    if ((ps.master_dark = 
         casu_frameset_subgroup_1(framelist,ps.labels,nlab,
                                  VIMOS_CAL_DARK)) == NULL) {
        cpl_msg_error(fctid,"No master dark found");
        vimos_sci_tidy(&ps,0);
        return(-1);
    }
    nfail += vimos_testfrm_1(ps.master_dark,VIMOS_NEXTN,1,0);
                
    /* Check to see if there is a master twilight flat frame */

    if ((ps.master_twilight_flat = 
         casu_frameset_subgroup_1(framelist,ps.labels,nlab,
                                  VIMOS_CAL_TWILIGHT_FLAT)) == NULL) {
        cpl_msg_error(fctid,"No master twilight flat found");
        vimos_sci_tidy(&ps,0);
        return(-1);
    }
    nfail += vimos_testfrm_1(ps.master_twilight_flat,VIMOS_NEXTN,1,0);

    /* Get the gain corrections */

    status = CASU_OK;
    if (casu_gaincor_calc(ps.master_twilight_flat,&i,&(ps.gaincors),
                          &status) != CASU_OK) {
        cpl_msg_error(fctid,"Error calculating gain corrections");
        vimos_sci_tidy(&ps,0);
        return(-1);
    }

    /* Check to see if there is a master confidence map. */

    if ((ps.master_conf = 
        casu_frameset_subgroup_1(framelist,ps.labels,nlab,
                                 VIMOS_CAL_CONF)) == NULL) {
        cpl_msg_error(fctid,"No master confidence map file found");
        vimos_sci_tidy(&ps,0);
        return(-1);
    }
    nfail += vimos_testfrm_1(ps.master_conf,VIMOS_NEXTN,1,0);
    ps.mask = casu_mask_define(framelist,ps.labels,nlab,VIMOS_CAL_CONF,
                               "NONE");

    /* Look for a master fringe map. It's not a problem if it doesn't 
       exist */

    ps.master_fringe = casu_frameset_subgroup_1(framelist,ps.labels,nlab,
                                                VIMOS_CAL_FRINGE);
    if (ps.master_fringe != NULL)
        nfail += vimos_testfrm_1(ps.master_fringe,VIMOS_NEXTN,1,0);

    /* Look for a master fringe variance map. It's not a problem if it doesn't 
       exist */

    ps.master_fringe_var = casu_frameset_subgroup_1(framelist,ps.labels,nlab,
                                                    VIMOS_CAL_FRINGE_VAR);
    if (ps.master_fringe_var != NULL)
        nfail += vimos_testfrm_1(ps.master_fringe_var,VIMOS_NEXTN,1,0);

    /* Check to see if there is a master photometric matched standards table */

    ps.master_mstd_phot = casu_frameset_subgroup_1(framelist,ps.labels,nlab,
                                                   VIMOS_CAL_MSTD_PHOT);
    if (ps.master_mstd_phot != NULL) 
        nfail += vimos_testfrm_1(ps.master_mstd_phot,VIMOS_NEXTN,0,0);
        
    /* Check to see if there is a master read/gain table */

    if ((ps.readgain = 
        casu_frameset_subgroup_1(framelist,ps.labels,nlab,
                                 VIMOS_CAL_READGAIN)) == NULL) {
        cpl_msg_error(fctid,"No master readgain table file found");
        vimos_sci_tidy(&ps,0);
        return(-1);
    }
    nfail += vimos_testrdgn(ps.master_conf,ps.readgain);

    /* Is the 2mass/ppmxl index file specified for astrom? */

    catindex_a = NULL;
    if (cs.cdssearch_astrom == 0) {
        if ((catindex_a = casu_frameset_subgroup_1(framelist,ps.labels,nlab,
                                                   VIMOS_CAL_2MASS_A)) == NULL) {
            if ((catindex_a = casu_frameset_subgroup_1(framelist,ps.labels,nlab,
                                                       VIMOS_CAL_PPMXL_A)) == NULL) {
                if ((catindex_a = casu_frameset_subgroup_1(framelist,ps.labels,nlab,
                                                           VIMOS_CAL_APASS_A)) == NULL) {
                    if ((catindex_a = casu_frameset_subgroup_1(framelist,ps.labels,
                                                               nlab,VIMOS_CAL_LOCCAT_A)) == NULL) {
        
                        cpl_msg_error(fctid,
                                      "No astrometric standard catalogue found -- cannot continue");
                        vimos_sci_tidy(&ps,0);
                        return(-1);
                    }
                }
            }
        }
        nfail += vimos_testfrm_1(catindex_a,1,0,0);
    } else {
        (void)casu_getstds_cdslist(cs.cdssearch_astrom,&junk1,&junk2,&status);
        freespace(junk1);
        freespace(junk2);
        if (status != CASU_OK) {
            status = CASU_OK;
            nfail++;
        }
    }

    /* See if there is a master standard star source catalogue. If there isn't
       then start looking for other ways to do the photometric solution */

    if ((ps.master_stdstar = casu_frameset_subgroup_1(framelist,ps.labels,nlab,
                                                      VIMOS_CAL_STDSTAR)) != NULL) {
        nfail += vimos_testfrm_1(ps.master_stdstar,VIMOS_NEXTN,0,0);
    } else {
        
        /* Is there an index file specified for photom? */

        catindex_p = NULL;
        pcat[0] = '\0';
        if (cs.cdssearch_photom == 0 && ps.master_mstd_phot == NULL) {
            if ((catindex_p = casu_frameset_subgroup_1(framelist,ps.labels,nlab,
                                                       VIMOS_CAL_APASS_P)) == NULL) {
                if ((catindex_p = casu_frameset_subgroup_1(framelist,ps.labels,
                                                           nlab,VIMOS_CAL_LOCCAT_P)) == NULL) {
                    cpl_msg_error(fctid,
                                  "No photometric standard catalogue found -- cannot continue");
                    vimos_sci_tidy(&ps,0);
                    return(-1);
                }
            }
            nfail += vimos_testfrm_1(catindex_p,1,0,0);
            pp = cpl_propertylist_load(cpl_frame_get_filename(catindex_p),1);
            if (cpl_propertylist_has(pp,"EXTNAME")) {
                strncpy(pcat,cpl_propertylist_get_string(pp,"EXTNAME"),32);
            } else {
                cpl_msg_error(fctid,"No extension name in %s",
                              cpl_frame_get_filename(catindex_p));
                nfail++;
            }
            cpl_propertylist_delete(pp);
        } else if (ps.master_mstd_phot != NULL) {
            pp = cpl_propertylist_load(cpl_frame_get_filename(ps.master_mstd_phot),0);
            strncpy(pcat,cpl_propertylist_get_string(pp,"PHOTCAT"),32);
            cpl_propertylist_delete(pp);
        } else {
            (void)casu_getstds_cdslist(cs.cdssearch_photom,&junk1,&junk2,&status);
            strncpy(pcat,junk1,32);
            freespace(junk1);
            freespace(junk2);
            if (status != CASU_OK) {
                status = CASU_OK;
                nfail++;
            }
        }
    }

    /* Check to see if there is a photometric table */

    if ((ps.phottab = casu_frameset_subgroup_1(framelist,ps.labels,nlab,
                                               VIMOS_CAL_PHOTTAB)) == NULL) {
        cpl_msg_error(fctid,"No photometric table found");
        vimos_sci_tidy(&ps,0);
        return(-1);
    }
    n = cpl_frame_get_nextensions(ps.phottab);
    nfail += vimos_testfrm_1(ps.phottab,n,0,0);

    /* Check the photometric table to see if it contains information about the
       photometric source we want to use */

    if (ps.master_stdstar == NULL) {
        match = 0;
        for (i = 1; i <= n; i++) {
            pp = cpl_propertylist_load(cpl_frame_get_filename(ps.phottab),i);
            if (cpl_propertylist_has(pp,"EXTNAME") &&
                strcmp(pcat,cpl_propertylist_get_string(pp,"EXTNAME")) == 0) {
                if (cpl_propertylist_has(pp,"PHOTOSYS")) {
                    strncpy(photosys,cpl_propertylist_get_string(pp,"PHOTOSYS"),8);
                } else {
                    cpl_msg_error(fctid,
                                  "Photcal table missing PHOTOSYS information in extension %d",i);
                    cpl_propertylist_delete(pp);
                    break;
                }
                ps.tphottab = cpl_table_load(cpl_frame_get_filename(ps.phottab),
                                             i,0);
                cpl_propertylist_delete(pp);
                match = 1;
                break;
            }
            cpl_propertylist_delete(pp);
        }
        if (! match) {
            cpl_msg_error(fctid,"Photcal table has no information on %s",pcat);
            nfail++;
        }

        /* Is the Northern Schlegel file specified? */

        if ((ps.schlf_n = casu_frameset_subgroup_1(framelist,ps.labels,nlab,
                                                   VIMOS_CAL_SCHL_N)) == NULL) {
            cpl_msg_error(fctid,"Schlegel North map not found -- cannot continue");
            vimos_sci_tidy(&ps,0);
            return(-1);
        }
        nfail += vimos_testfrm_1(ps.schlf_n,0,1,0);

        /* Is the Southern Schlegel file specified? */

        if ((ps.schlf_s = casu_frameset_subgroup_1(framelist,ps.labels,nlab,
                                                   VIMOS_CAL_SCHL_S)) == NULL) {
            cpl_msg_error(fctid,"Schlegel South map not found -- cannot continue");
            vimos_sci_tidy(&ps,0);
            return(-1);
        }
        nfail += vimos_testfrm_1(ps.schlf_s,0,1,0);
    } else {
        (void)strcpy(photosys,"VEGA");
    }

    /* Check the cache location is writable */

    if (access(cs.cacheloc,R_OK+W_OK+X_OK) != 0) {
        cpl_msg_error(fctid,"Cache location %s inacessible",cs.cacheloc);
        nfail++;
    }

    /* Ok if any of this failed, then get out of here. This makes it a bit 
       simpler later on since we now know that all of the specified files 
       have the correct number of extensions and will load properly */

    if (nfail > 0) {
        cpl_msg_error(fctid,
                      "There are %" CPL_SIZE_FORMAT " input file errors -- cannot continue",
                      (cpl_size)nfail);
        freeframe(catindex_a);
        freeframe(catindex_p);
        vimos_sci_tidy(&ps,0);
        return(-1);
    }

    /* Get catalogue parameters */

    if (cs.cdssearch_astrom == 0) {
        if (casu_catpars(catindex_a,&(ps.catpath_a),
                         &(ps.catname_a)) == CASU_FATAL) {
            vimos_sci_tidy(&ps,0);
            cpl_frame_delete(catindex_a);
            return(-1);
        }
        cpl_frame_delete(catindex_a);
    } else {
        ps.catpath_a = NULL;
        ps.catname_a = NULL;
    }
    if (cs.cdssearch_photom == 0 && ps.master_mstd_phot == NULL &&
        ps.master_stdstar == NULL) {
        if (casu_catpars(catindex_p,&(ps.catpath_p),
                         &(ps.catname_p)) == CASU_FATAL) {
            vimos_sci_tidy(&ps,0);
            cpl_frame_delete(catindex_p);
            return(-1);
        }
        cpl_frame_delete(catindex_p);
    } else {
        ps.catpath_p = NULL;
        ps.catname_p = NULL;
    }

    /* Get some info you might need */

    frm = cpl_frameset_get_position(ps.science_frames,0);
    ff = casu_fits_load(frm,CPL_TYPE_FLOAT,cpl_frame_get_nextensions(frm));
    pp = casu_fits_get_phu(ff);
    (void)vimos_pfits_get_filter(pp,filt);
    (void)vimos_pfits_get_projid(pp,projid);
    freefits(ff);

    /* Print out some stuff if you want a preview */

    if (cs.preview_only) {
        fprintf(stdout,"This is a paw group with %" CPL_SIZE_FORMAT " files\n",
                cpl_frameset_get_size(ps.science_frames));
        fprintf(stdout,"This is split into %" CPL_SIZE_FORMAT " individual observations\n",
                (cpl_size)ps.ngr4);
        njsteps = cpl_frameset_get_size(ps.science_frames);

        fprintf(stdout,"Filter: %8s, Project: %s\n",filt,projid);
        fprintf(stdout,"Bias: %s\n",cpl_frame_get_filename(ps.master_bias));
        fprintf(stdout,"Dark: %s\n",cpl_frame_get_filename(ps.master_dark));
        fprintf(stdout,"Flat: %s\n",
                cpl_frame_get_filename(ps.master_twilight_flat));
        fprintf(stdout,"CPM: %s\n",cpl_frame_get_filename(ps.master_conf));
        if (ps.master_fringe != NULL && ! cs.ignore_fringe) {
            fprintf(stdout,"Fringe: %s\n",
                    cpl_frame_get_filename(ps.master_fringe));
        } else {
            fprintf(stdout,"No defringing to be done\n");
        }
        fprintf(stdout,"A %" CPL_SIZE_FORMAT " point dither will be done\n",
                (cpl_size)njsteps);
        n = cpl_frameset_get_size(ps.science_frames);
        for (j = 0; j < n; j++) {
            frm = cpl_frameset_get_position(ps.science_frames,j);
            fprintf(stdout,"%s\n",cpl_frame_get_filename(frm));
        }
        if (cs.cdssearch_astrom != 0) {
            p = cpl_parameterlist_find(parlist,
                                       "vimos.vimos_ima_science.cdssearch_astrom");
            fprintf(stdout,"Astrometry will be done using CDS catalogue: %s\n",
                    cpl_parameter_get_string(p));
        } else {
            fprintf(stdout,"Astrometry will be done using local catalogue %s in %s\n",
                    ps.catname_a,ps.catpath_a);
        }
        if (ps.master_mstd_phot != NULL) {
            fprintf(stdout,"Photometry will be done using matched cat %s\n",
                    cpl_frame_get_filename(ps.master_mstd_phot));
        } else if (cs.cdssearch_photom != 0) {
            p = cpl_parameterlist_find(parlist,
                                       "vimos.vimos_ima_science.cdssearch_photom");
            fprintf(stdout,"Photometry will be done using CDS catalogue: %s\n",
                    cpl_parameter_get_string(p));
        } else {
            fprintf(stdout,"Photometry will be done using local catalogue %s in %s\n",
                    ps.catname_p,ps.catpath_p);
        }
        vimos_sci_tidy(&ps,0);
        return(0);
    }

    /* Get some workspace for the product frames */

    ps.product_frames_simple = cpl_malloc(ps.ngr4*sizeof(cpl_frame *));
    ps.product_frames_simple_var = cpl_malloc(ps.ngr4*sizeof(cpl_frame *));
    ps.product_frames_simple_cat = cpl_malloc(ps.ngr4*sizeof(cpl_frame *));
    ps.product_frames_simple_mstd_a = cpl_malloc(ps.ngr4*sizeof(cpl_frame *));
    ps.sci_fits = cpl_malloc(ps.ngr4*sizeof(casu_fits *));
    ps.sci_vars = cpl_malloc(ps.ngr4*sizeof(casu_fits *));
    ps.sci_cats = cpl_malloc(ps.ngr4*sizeof(casu_tfits *));
    ps.sci_mstd_a = cpl_malloc(ps.ngr4*sizeof(casu_tfits*));
    for (i = 0; i < ps.ngr4; i++) {
        ps.product_frames_simple[i] = NULL;
        ps.product_frames_simple_var[i] = NULL;
        ps.product_frames_simple_cat[i] = NULL;
        ps.product_frames_simple_mstd_a[i] = NULL;
        ps.sci_fits[i] = NULL;
        ps.sci_vars[i] = NULL;
        ps.sci_cats[i] = NULL;
        ps.sci_mstd_a[i] = NULL;
    }

    /* Now loop for all detectors... */

    template = NULL;
    assoc[0] = NULL;
    assoc[1] = NULL;
    ps.prod_im = cpl_frameset_new();
    ps.prod_tab = cpl_frameset_new();
    for (j = 1; j <= VIMOS_NEXTN; j++) {
        isfirst = (j == 1);

        /* Right, we want to bias, dark and flat correct, so we need to load 
           the mean calibration frames */

        ps.fbias = casu_fits_load(ps.master_bias,CPL_TYPE_FLOAT,j);
        ps.fdark = casu_fits_load(ps.master_dark,CPL_TYPE_FLOAT,j);
        ps.fflat = casu_fits_load(ps.master_twilight_flat,CPL_TYPE_FLOAT,j);
        ps.fconf = casu_fits_load(ps.master_conf,CPL_TYPE_INT,j);
        casu_mask_load(ps.mask,j,
                       (int)cpl_image_get_size_x(casu_fits_get_image(ps.fconf)),
                       (int)cpl_image_get_size_y(casu_fits_get_image(ps.fconf)));
        if (ps.master_fringe != NULL && ! cs.ignore_fringe) {
            ps.ffringe = casu_fits_load(ps.master_fringe,CPL_TYPE_FLOAT,j);
            if (ps.master_fringe_var != NULL) 
                ps.ffringe_var = casu_fits_load(ps.master_fringe_var,
                                                CPL_TYPE_FLOAT,j);
            status = CASU_OK;
            switch (cs.chop_crud) {
            case 1:
                vimos_chop_lowconfbands(ps.ffringe,ps.fconf,&status);
                vimos_chop_lowconfbands(ps.ffringe_var,ps.fconf,&status);
                break;
            case 2:
                vimos_chop_region(ps.ffringe,j,&status);
                vimos_chop_region(ps.ffringe_var,j,&status);
                break;
            case 4:
                vimos_chop_region(ps.ffringe,j,&status);
                vimos_chop_region(ps.ffringe_var,j,&status);
                break;
            }
        } else {
            ps.ffringe = NULL;
            ps.ffringe_var = NULL;
        }

        /* Gain corrections are in the same order as vimos_names so there's
           no trouble making sure we have the right one... */

        gaincor_fac = (ps.gaincors)[j-1];

        /* Get the readnoise and gain estimates */

        vimos_getrdgn(ps.readgain,casu_fits_get_extname(ps.fflat),
                      &readnoise,&gain);

        /* Create a duplicate flat field that is chopped like the
           science images */

        ps.fconf2 = casu_fits_duplicate(ps.fconf);
        status = CASU_OK;
        switch (cs.chop_crud) {
        case 1:
            vimos_chop_lowconfbands(ps.fconf2,ps.fconf,&status);
            break;
        case 2:
            vimos_chop_region(ps.fconf2,j,&status);
            break;
        case 3:
            vimos_chop_lowconfpix(ps.fconf2,&status);
            break;
        case 4:
            vimos_chop_region(ps.fconf2,j,&status);
            cpl_propertylist_erase(casu_fits_get_ehu(ps.fconf2),
                                   "ESO DRS CHOPCOR");
            vimos_chop_lowconfpix(ps.fconf2,&status);
            break;
        }

        /* Loop for each image and do the stage 1 corrections. */

        cpl_msg_info(fctid,"Stage 1 corrections for %s",vimos_names[j-1]);
        for (i = 0; i < ps.ngr4; i++) {
            status = CASU_OK;
            frm = ps.gr4[i].inf[j-1];
            if (frm == NULL) {
                im = casu_dummy_image(ps.fconf2);
                (ps.sci_fits)[i] = casu_fits_wrap(im,NULL,NULL,NULL);
                casu_fits_set_status((ps.sci_fits)[i],CASU_FATAL);
                im = casu_dummy_image(ps.fconf2);
                (ps.sci_vars)[i] = casu_fits_wrap(im,NULL,NULL,NULL);
                casu_fits_set_status((ps.sci_vars)[i],CASU_FATAL);
                cpl_msg_info(fctid,
                             "There is no image for detector %s in group %s",
                             vimos_names[j-1],ps.gr4[i].name);
                continue;
            }
            ff = casu_fits_load(frm,CPL_TYPE_FLOAT,
                                cpl_frame_get_nextensions(frm));
            (ps.sci_fits)[i] = ff;
            vimos_biascor(ff,ps.fbias,1,1,&status);
            float texp = 0.0f;
            vimos_pfits_get_exptime(casu_fits_get_ehu(ff),&texp);
            casu_darkcor(ff,ps.fdark,texp,&status);
            ffv = vimos_var_create(ff,ps.mask,readnoise,gain);
            casu_flatcor(ff,ps.fflat,&status);
            vimos_var_div_im(ffv,ps.fflat);
            casu_gaincor(ff,gaincor_fac,&status);
            vimos_var_divk(ffv,gaincor_fac);
            switch (cs.chop_crud) {
            case 1:
                vimos_chop_lowconfbands(ff,ps.fconf,&status);
                vimos_chop_lowconfbands(ffv,ps.fconf,&status);
                break;
            case 2:
                vimos_chop_region(ff,j,&status);
                vimos_chop_region(ffv,j,&status);
                break;
            case 4:
                vimos_chop_region(ff,j,&status);
                vimos_chop_region(ffv,j,&status);
                break;
            }
            casu_fits_set_error(ff,status);
            cpl_propertylist_update_float(casu_fits_get_ehu(ff),"READNOIS",
                                          readnoise);
            cpl_propertylist_set_comment(casu_fits_get_ehu(ff),"READNOIS",
                                         "[ADU] Read noise used in reduction");
            cpl_propertylist_update_float(casu_fits_get_ehu(ff),"GAIN",
                                          gain);
            cpl_propertylist_set_comment(casu_fits_get_ehu(ff),"GAIN",
                                         "[e-/ADU] gain used in reduction");
            (void)casu_imcore(ff,ps.fconf2,cs.stk_cat_ipix,cs.stk_cat_thresh,
                              cs.stk_cat_icrowd,cs.stk_cat_rcore,
                              cs.stk_cat_nbsize,6,2.0,&tab,gain,&status);
            if (status == CASU_OK) {
                vimos_wcsfit(&ff,&tab,1,ps.catname_a,ps.catpath_a,
                             cs.cdssearch_astrom,cs.cacheloc,cs.savemstd,
                             0,&(ps.sci_mstd_a[i]));
                vimos_copywcs(casu_fits_get_ehu(ff),casu_fits_get_ehu(ffv));
            } else {
                if (cs.savecat) {
                    t = casu_dummy_catalogue(2);
                    tab = casu_tfits_wrap(t,NULL,NULL,NULL);
                } else {
                    tab = NULL;
                }
                if (cs.savemstd) 
                    ps.sci_mstd_a[i] = casu_tfits_wrap(casu_dummy_catalogue(2),
                                                       NULL,NULL,NULL);
                else
                    ps.sci_mstd_a[i] = NULL;
            }

            /* Store away the variance map and the catalogue if we want 
               to keep it */

            (ps.sci_vars)[i] = ffv;
            if (cs.savecat) 
                (ps.sci_cats)[i] = tab;
            else
                freetfits(tab);
        }

        /* Defringe the data */

        if (ps.ffringe != NULL && ! cs.ignore_fringe) {
            casu_defringe(ps.sci_fits,ps.ngr4,&(ps.ffringe),1,ps.mask,
                          512,&status);
            if (ps.ffringe_var != NULL) {
                for (i = 0; i < ps.ngr4; i++) {
                    fac = cpl_propertylist_get_float(casu_fits_get_ehu(ps.sci_fits[i]),
                                                     "ESO DRS FRNGSC1");
                    vimos_var_scaledsubt(ps.sci_vars[i],ps.ffringe_var,fac);
                }
            }
        }

        /* Redo the loop for now. Because later we'll need to add 
           defringing. If the casu_fits status for an image is bad
           then a dummy image is saved. */

        status = CASU_OK;
        cpl_msg_info(fctid,"Saving simple images for %s",vimos_names[j-1]);
        for (i = 0; i < ps.ngr4; i++) {
            for (k = 0; k < 4; k++) {
                frm = ps.gr4[i].inf[k];
                if (frm != NULL)
                    break;
            }
            if (template == NULL)
                template = cpl_frame_duplicate(frm);

            /* Save the variance map */

            if (j == 1) {
                fname = cpl_strdup(cpl_frame_get_filename(frm));
                vimos_sci_product_name(fname,SIMPLE_VAR,cs.prettynames,i+1,
                                       bname);
                freespace(fname);
            } else {
                strcpy(bname,
                       cpl_frame_get_filename((ps.product_frames_simple_var)[i]));
            }
            if (vimos_sci_save_simple((ps.sci_vars)[i],framelist,parlist,
                                      frm,isfirst,VIMOS_PRO_VAR_SCI,bname,
                                      vimos_names[j-1],assoc[0],
                                      ps.product_frames_simple_var+i) != CASU_OK) {
                cpl_msg_error(fctid,"Unable to save product %s[%d]",bname,j);
                cpl_frameset_delete(chiplist);
                vimos_sci_tidy(&ps,0);
                return(-1);
            }

            /* Save the image */

            if (j == 1) {
                fname = cpl_strdup(cpl_frame_get_filename(frm));
                vimos_sci_product_name(fname,SIMPLE_FILE,cs.prettynames,i+1,
                                       bname);
                freespace(fname);
                assoc[0] = cpl_strdup(cpl_frame_get_filename(ps.product_frames_simple_var[i]));
            } else {
                strcpy(bname,
                       cpl_frame_get_filename((ps.product_frames_simple)[i]));
            }
            if (vimos_sci_save_simple((ps.sci_fits)[i],framelist,parlist,
                                      frm,isfirst,VIMOS_PRO_SIMPLE_SCI,bname,
                                      vimos_names[j-1],assoc[0],
                                      ps.product_frames_simple+i) != CASU_OK) {
                cpl_msg_error(fctid,"Unable to save product %s[%d]",bname,j);
                cpl_frameset_delete(chiplist);
                vimos_sci_tidy(&ps,0);
                return(-1);
            }
            if (j == 1) {
                freespace(assoc[0]);
                cpl_frameset_insert(ps.prod_im,cpl_frame_duplicate(ps.product_frames_simple[i]));
            }

            /* Now the catalogue if that's what we want */

            if (cs.savecat) {
            	cpl_propertylist_update_int(casu_tfits_get_phu((ps.sci_cats)[i]), "NSTACK", ps.ngr4 * 4);
                if (vimos_sci_save_cat((ps.sci_cats)[i],framelist,parlist,
                                       frm,cs.prettynames,SIMPLE_CAT,
                                       i+1,vimos_names[j-1],photosys,
                                       ps.product_frames_simple_cat+i) != CASU_OK) {
                    cpl_msg_error(fctid,"Unable to save product %s[%d]",
                                  cpl_frame_get_filename(ps.product_frames_simple_cat[i]),j);
                    cpl_frameset_delete(chiplist);
                    vimos_sci_tidy(&ps,0);
                    return(-1);
                }
                if (j == 1) 
                    cpl_frameset_insert(ps.prod_tab,cpl_frame_duplicate(ps.product_frames_simple_cat[i]));
            }

            /* And the astrometric matched standards catalogue if requested */

            if (cs.savemstd) {
            	cpl_propertylist_update_int(casu_tfits_get_phu((ps.sci_mstd_a)[i]), "NSTACK", ps.ngr4 * 4);
                if (vimos_sci_save_cat((ps.sci_mstd_a)[i],framelist,parlist,
                                       frm,cs.prettynames,SIMPLE_MSTD_A,
                                       i+1,vimos_names[j-1],photosys,
                                       ps.product_frames_simple_mstd_a+i) != CASU_OK) {

                    cpl_msg_error(fctid,"Unable to save product %s[%d]",
                                  cpl_frame_get_filename(ps.product_frames_simple_mstd_a[i]),j);
                    cpl_frameset_delete(chiplist);
                    vimos_sci_tidy(&ps,0);
                    return(-1);
                }
            }
        }

        /* Redo the list by removing any dummies. If this is the
           first extension, then get the mjd info you need */

        ngood = 0;
        good = cpl_malloc(ps.ngr4*sizeof(casu_fits *));
        goodv = cpl_malloc(ps.ngr4*sizeof(casu_fits *));
        double sum = 0.0;
        float sum_exp = 0.0f;
        double mjdmax = -1.0e15;
        float tmax = 0.0;
        for (i = 0; i < ps.ngr4; i++) {
            if (casu_fits_get_status((ps.sci_fits)[i]) == CASU_OK) {
                goodv[ngood] = (ps.sci_vars)[i];
                good[ngood++] = (ps.sci_fits)[i];
                double mjd = 0.0;
                vimos_pfits_get_mjd(casu_fits_get_phu((ps.sci_fits)[i]),&mjd);
                sum += mjd;
                float t_exp = 0.0;
                vimos_pfits_get_exptime(casu_fits_get_phu((ps.sci_fits)[i]),&t_exp);
                sum_exp += t_exp;
                if (mjd > mjdmax) {
                    mjdmax = mjd;
                    tmax = t_exp;
                }
            }
        }

        /* Do the stacking. If there are no good entries, the create a
           dummy image */

        if (ngood > 0) {
            status = CASU_OK;
            cpl_msg_info(fctid,"Stacking images for %s",vimos_names[j-1]);
            if (cs.stk_fast == -1)
                fastslow = (ngood <= cs.stk_nfst);
            else
                fastslow = cs.stk_fast;
            casu_imstack(good,&(ps.fconf2),goodv,NULL,ngood,1,cs.stk_lthr,
                         cs.stk_hthr,cs.stk_method,cs.stk_seeing,fastslow,0,
                         "ESO DET DIT",&out,&outc,&outv,&status);
            casu_fits_set_status(out,status);
            casu_fits_set_status(outc,status);
        } else {
            cpl_msg_warning(fctid,"No good images for stacking %s",
                            vimos_names[j-1]);
            im = casu_dummy_image((ps.sci_fits)[0]);
            out = casu_fits_wrap(im,(ps.sci_fits)[0],NULL,NULL);
            casu_dummy_property(casu_fits_get_ehu(out));
            casu_fits_set_status(out,CASU_FATAL);
            im = casu_dummy_image(ps.fconf2);
            outc = casu_fits_wrap(im,ps.fconf2,NULL,NULL);
            casu_dummy_property(casu_fits_get_ehu(outc));
            casu_fits_set_status(outc,CASU_FATAL);
            imv = casu_dummy_image(ps.sci_vars[0]);
            outv = casu_fits_wrap(imv,ps.fconf,NULL,NULL);
        }

        double eff_exp = 0.0;
		if (ngood > 0) {
			sum /= (double)ngood;
			eff_exp = sum_exp / (double)ngood;
		}

        cpl_propertylist_update_double(casu_fits_get_ehu(out),"MJD_MEAN",
									   sum);
		cpl_propertylist_set_comment(casu_fits_get_ehu(out),"MJD_MEAN",
									 "Mean MJD of the input images");

		cpl_propertylist_update_double(casu_fits_get_ehu(outc),"MJD_MEAN",
									   sum);
		cpl_propertylist_set_comment(casu_fits_get_ehu(outc),"MJD_MEAN",
									 "Mean MJD of the input images");

		cpl_propertylist_update_double(casu_fits_get_ehu(outv),"MJD_MEAN",
									   sum);
		cpl_propertylist_set_comment(casu_fits_get_ehu(outv),"MJD_MEAN",
									 "Mean MJD of the input images");

		cpl_propertylist_update_double(casu_fits_get_phu(out),"MJD-END",
									   mjdmax+tmax/86400.0);
		cpl_propertylist_set_comment(casu_fits_get_phu(out),"MJD-END",
									 "End of observations");

		cpl_propertylist_update_double(casu_fits_get_phu(outc),"MJD-END",
									   mjdmax+tmax/86400.0);
		cpl_propertylist_set_comment(casu_fits_get_phu(outc),"MJD-END",
									 "End of observations");

		cpl_propertylist_update_double(casu_fits_get_phu(outv),"MJD-END",
									   mjdmax+tmax/86400.0);
		cpl_propertylist_set_comment(casu_fits_get_phu(outv),"MJD-END",
									 "End of observations");

		cpl_propertylist_update_int(casu_fits_get_phu(out),"NCOMBINE",
			4*ps.ngr4);
		cpl_propertylist_set_comment(casu_fits_get_phu(out),"NCOMBINE",
									 "Number of input raw files");

		cpl_propertylist_update_int(casu_fits_get_phu(outv),"NCOMBINE",
									4*ps.ngr4);
		cpl_propertylist_set_comment(casu_fits_get_phu(outv),"NCOMBINE",
									 "Number of input raw files");

		cpl_propertylist_update_int(casu_fits_get_phu(outc),"NCOMBINE",
									4*ps.ngr4);
		cpl_propertylist_set_comment(casu_fits_get_phu(outc),"NCOMBINE",
									 "Number of input raw files");

		cpl_propertylist_update_int(casu_fits_get_phu(out),"NSTACK",
			ps.ngr4);
		cpl_propertylist_set_comment(casu_fits_get_phu(out),"NSTACK",
									 "Number of stacked images");

		cpl_propertylist_update_int(casu_fits_get_phu(outv),"NSTACK",
									ps.ngr4);
		cpl_propertylist_set_comment(casu_fits_get_phu(outv),"NSTACK",
									 "Number of stacked images");

		cpl_propertylist_update_int(casu_fits_get_phu(outc),"NSTACK",
									ps.ngr4);
		cpl_propertylist_set_comment(casu_fits_get_phu(outc),"NSTACK",
									 "Number of stacked images");

		cpl_propertylist_update_double(casu_fits_get_phu(out),  "EXPTIME",  sum_exp);
		cpl_propertylist_update_double(casu_fits_get_phu(outc), "EXPTIME",  sum_exp);
		cpl_propertylist_update_double(casu_fits_get_phu(outv), "EXPTIME",  sum_exp);
		cpl_propertylist_update_double(casu_fits_get_phu(out),  "TEXPTIME", sum_exp);
		cpl_propertylist_update_double(casu_fits_get_phu(outc), "TEXPTIME", sum_exp);
		cpl_propertylist_update_double(casu_fits_get_phu(outv), "TEXPTIME", sum_exp);
		cpl_propertylist_update_double(casu_fits_get_phu(out),  "EFF_EXPT", eff_exp);
		cpl_propertylist_update_double(casu_fits_get_phu(outc), "EFF_EXPT", eff_exp);
		cpl_propertylist_update_double(casu_fits_get_phu(outv), "EFF_EXPT", eff_exp);

		cpl_propertylist_update_double(casu_fits_get_ehu(out),  "EXPTIME",  sum_exp);
		cpl_propertylist_update_double(casu_fits_get_ehu(outc), "EXPTIME",  sum_exp);
		cpl_propertylist_update_double(casu_fits_get_ehu(outv), "EXPTIME",  sum_exp);
		cpl_propertylist_update_double(casu_fits_get_ehu(out),  "TEXPTIME", sum_exp);
		cpl_propertylist_update_double(casu_fits_get_ehu(outc), "TEXPTIME", sum_exp);
		cpl_propertylist_update_double(casu_fits_get_ehu(outv), "TEXPTIME", sum_exp);
		cpl_propertylist_update_double(casu_fits_get_ehu(out),  "EFF_EXPT", eff_exp);
		cpl_propertylist_update_double(casu_fits_get_ehu(outc), "EFF_EXPT", eff_exp);
		cpl_propertylist_update_double(casu_fits_get_ehu(outv), "EFF_EXPT", eff_exp);

        freespace(good);
        freespace(goodv);

        /* Do the catalogue */

        mstds_a[j-1] = NULL;
        if (casu_fits_get_status(out) == CASU_OK) {
            cpl_msg_info(fctid,"Source extraction for %s",vimos_names[j-1]);
            vimos_pfits_get_gain(casu_fits_get_ehu(out),&gain);
            (void)casu_imcore(out,outc,cs.stk_cat_ipix,cs.stk_cat_thresh,
                              cs.stk_cat_icrowd,cs.stk_cat_rcore,
                              cs.stk_cat_nbsize,6,2.0,&tab,gain,&status);
            if (status == CASU_OK) {
                vimos_wcsfit(&out,&tab,1,ps.catname_a,ps.catpath_a,
                             cs.cdssearch_astrom,cs.cacheloc,cs.savemstd,
                             0,&(mstds_a[j-1]));
                vimos_copywcs(casu_fits_get_ehu(out),casu_fits_get_ehu(outv));
                vimos_copywcs(casu_fits_get_ehu(out),casu_fits_get_ehu(outc));
                if (ps.master_mstd_phot == NULL && ps.master_stdstar == NULL) {
                    (void)casu_getstds(casu_fits_get_ehu(ff),1,ps.catpath_p,
                                       ps.catname_p,cs.cdssearch_photom,
                                       cs.cacheloc,&stdscat,&status);
                    if (status != CASU_OK) {
                        cpl_msg_warning(fctid,"No photometric standards found");
                        status = CASU_OK;
                    }
                } else 
                    stdscat = NULL;
            } else {
                stdscat = NULL;
            }
        } else {
            cpl_msg_info(fctid,"No stack, source extraction skipped %s",
                         vimos_names[j-1]);
            t = casu_dummy_catalogue(2);
            tab = casu_tfits_wrap(t,NULL,
                                  cpl_propertylist_duplicate(casu_fits_get_phu(out)),
                                  cpl_propertylist_duplicate(casu_fits_get_ehu(out)));
            casu_tfits_set_status(tab,CASU_FATAL);
            stdscat = NULL;
            if (cs.savemstd) 
                mstds_a[j-1] = casu_tfits_wrap(casu_dummy_catalogue(2),NULL,
                                               cpl_propertylist_duplicate(casu_fits_get_phu(out)),
                                               cpl_propertylist_duplicate(casu_fits_get_ehu(out)));
            else
                mstds_a[j-1] = NULL;
        }
        if (ps.master_stdstar != NULL) {
            mstds_p[j-1] = NULL;
        } else if (ps.master_mstd_phot != NULL) {
            t = cpl_table_load(cpl_frame_get_filename(ps.master_mstd_phot),j,1);
            pp = cpl_propertylist_load(cpl_frame_get_filename(ps.master_mstd_phot),0);
            epp = cpl_propertylist_load(cpl_frame_get_filename(ps.master_mstd_phot),j);
            mstds_p[j-1] = casu_tfits_wrap(t,NULL,pp,epp);
        } else if (stdscat != NULL) {
            (void)casu_matchstds(casu_tfits_get_table(tab),
                                 stdscat,300.0,&t,&status);
            if (status == CASU_OK) {
                mstds_p[j-1] = casu_tfits_wrap(t,tab,NULL,NULL);
            } else {
                status = CASU_OK;
                if (cs.savemstd) 
                    mstds_p[j-1] = casu_tfits_wrap(casu_dummy_catalogue(2),NULL,
                                                   cpl_propertylist_duplicate(casu_fits_get_phu(out)),
                                                   cpl_propertylist_duplicate(casu_fits_get_ehu(out)));
                else
                    mstds_p[j-1] = NULL;
            }
            cpl_propertylist_update_string(casu_tfits_get_phu(mstds_p[j-1]),
                                      
     "PHOTCAT",pcat);
            cpl_propertylist_set_comment(casu_tfits_get_phu(mstds_p[j-1]),
                                         "PHOTCAT",
                                         "Originating photometry source");
        } else {
            if (cs.savemstd) 
                mstds_p[j-1] = casu_tfits_wrap(casu_dummy_catalogue(2),NULL,
                                               cpl_propertylist_duplicate(casu_fits_get_phu(out)),
                                               cpl_propertylist_duplicate(casu_fits_get_ehu(out)));
            else
                mstds_p[j-1] = NULL;
        }
        freetable(stdscat);
        stack[j-1] = out;
        stackv[j-1] = outv;
        stackcat[j-1] = tab;
        stackconf[j-1] = outc;

        /* Clean up and go on to the next extension */

        vimos_sci_tidy(&ps,1);
    }

    /* Do the photometry on the stack. If we're using a master standard
       star as reference, then we can get away with scaling a few things.
       Otherwise we need to do the photometry properly. */

    if (ps.master_stdstar != NULL) {
        status = vimos_sci_phot_stdstar(ps.master_stdstar,stack,stackcat, 
                                        "ESO TEL AIRM START");
        if(status != CASU_OK)
        {
            cpl_msg_error(fctid,"Unable to compute photometry");
            vimos_sci_tidy(&ps,0);
            return(-1);
        }

    } else {
        status = CASU_OK;
        status = casu_photcal_extinct(stack,mstds_p,stackcat,VIMOS_NEXTN,filt,
                             ps.tphottab,cs.minphotom,ps.schlf_n,ps.schlf_s,
                             "EFF_EXPT","ESO TEL AIRM START",cs.magerrcut,
                             &status);
        if(status != CASU_OK)
        {
            cpl_msg_error(fctid,"Unable to compute extinction");
            vimos_sci_tidy(&ps,0);
            return(-1);
        }
    }

    /* Update the RA/DEC keywords in the primary of the stack products */

    vimos_sci_update_hdr_stack(stack,stackconf[0],stackv[0],stackcat[0]);

    /* Now save them all */

    for (j = 1; j <= VIMOS_NEXTN; j++) {
        out = stack[j-1];
        outv = stackv[j-1];
        tab = stackcat[j-1];
        outc = stackconf[j-1];

        /* Save the stack variance map */

        cpl_msg_info(fctid,"Saving stack variance for %s",vimos_names[j-1]);
        if (vimos_sci_save_stack(outv,framelist,parlist,template,cs.prettynames,
                                 vimos_names[j-1],STACK_VAR,assoc,photosys,
                                 &(ps.product_frame_stack_var)) != CASU_OK) {
            freefits(outv);
            cpl_msg_error(fctid,"Unable to save stack");
            vimos_sci_tidy(&ps,0);
            return(-1);
        }
        freefits(outv);

        /* And the confidence map of the stack */

        cpl_msg_info(fctid,"Saving stack confidence map for %s",
                     vimos_names[j-1]);
        if (vimos_sci_save_stack_conf(outc,framelist,parlist,template,
                                      cs.prettynames,vimos_names[j-1],
                                      &(ps.product_frame_stack_conf)) != CASU_OK) {
            freefits(outc);
            cpl_msg_error(fctid,"Unable to save stack conf map");
            vimos_sci_tidy(&ps,0);
            return(-1);
        }
        freefits(outc);

        /* Save the stack */

        if (j == 1) {
            assoc[0] = cpl_strdup(cpl_frame_get_filename(ps.product_frame_stack_conf));
            assoc[1] = cpl_strdup(cpl_frame_get_filename(ps.product_frame_stack_var));
        }
        cpl_msg_info(fctid,"Saving stack for %s",vimos_names[j-1]);
        if (vimos_sci_save_stack(out,framelist,parlist,template,cs.prettynames,
                                 vimos_names[j-1],STACK_FILE,assoc,photosys,
                                 &(ps.product_frame_stack)) != CASU_OK) {
            freefits(out);
            cpl_msg_error(fctid,"Unable to save stack");
            vimos_sci_tidy(&ps,0);
            return(-1);
        }
        freefits(out);
        if (j == 1) {
            freespace(assoc[0]);
            freespace(assoc[1]);
        }

        /* And the source catalogue of the stack */

        cpl_msg_info(fctid,"Saving stack source catalogue for %s",
                     vimos_names[j-1]);
        cpl_propertylist_erase_regexp(casu_tfits_get_phu(tab),"PROV[0-9]*",0);
        cpl_propertylist_update_string(casu_tfits_get_phu(tab),"PROV1",
                                       cpl_frame_get_filename(ps.product_frame_stack));
        if (vimos_sci_save_cat(tab,framelist,parlist,template,
                               cs.prettynames,STACK_CAT,1,
                               vimos_names[j-1],photosys,
                               &(ps.product_frame_stack_cat)) != CASU_OK) {
            freetfits(tab);
            cpl_msg_error(fctid,"Unable to save stack catalogue");
            vimos_sci_tidy(&ps,0);
            return(-1);
        }
        freetfits(tab);

        /* Now save the matched standards catalogues if that's been requested */

        if (cs.savemstd) {
            tab = mstds_a[j-1];
            if (vimos_sci_save_cat(tab,framelist,parlist,template,
                                   cs.prettynames,STACK_MSTD_A,1,
                                   vimos_names[j-1],photosys,
                                   &(ps.product_frame_stack_mstd_a)) != CASU_OK) {
                cpl_msg_error(fctid,"Unable to save stack astrom mstds catalogue");
                vimos_sci_tidy(&ps,0);
                return(-1);
            }
            if (ps.master_stdstar == NULL) {
                tab = mstds_p[j-1];
                if (vimos_sci_save_cat(tab,framelist,parlist,template,
                                       cs.prettynames,STACK_MSTD_P,1,
                                       vimos_names[j-1],photosys,
                                       &(ps.product_frame_stack_mstd_p)) != CASU_OK) {
                    cpl_msg_error(fctid,"Unable to save stack photom mstds catalogue");
                    vimos_sci_tidy(&ps,0);
                    return(-1);
                }
            }
        }
        freetfits(mstds_a[j-1]);
        freetfits(mstds_p[j-1]);
    }
    freeframe(template);

    /* Update the headers of the simple images */

    vimos_sci_update_hdr(ps.prod_im,ps.product_frame_stack,1);
    if (cs.savecat)
        vimos_sci_update_hdr(ps.prod_tab,ps.product_frame_stack,2);

    /* Get out of here */

    vimos_sci_tidy(&ps,0);
    return(0);
}

/*---------------------------------------------------------------------------*/
/**
    \par Name:
        vimos_sci_save_simple
    \par Purpose:
        Save the simple image products
    \par Description:
        An extension of a simple image product is saved here
    \par Language:
        C
    \param obj
        The casu_fits structure for the image to be saved
    \param framelist
        The input recipe frameset
    \param parlist
        The input recipe parameter list
    \param template
        A template to be used in the cpl saving routine
    \param isfirst
        Set of this is the first extension to be written to the file
    \param tag
        The product tag for the saved file
    \param fname
        The output file name
    \param cname
        The detector name
    \param assoc
        A character string with the name of the associated variance map
    \param product_frame
        The output product frame
    \retval CASU_OK if everything is ok
    \retval CASU_FATAL if something failed
    \author
        Jim Lewis, CASU
 */
/*---------------------------------------------------------------------------*/

static int vimos_sci_save_simple(casu_fits *obj, cpl_frameset *framelist,
                                 cpl_parameterlist *parlist, 
                                 cpl_frame *template, int isfirst,
                                 const char *tag, char *fname, char *cname,
                                 char *assoc, cpl_frame **product_frame) {
    cpl_propertylist *plist;
    int isdummy;
    float texp;
    char filt[32];
    const char *fctid = "vimos_sci_save_simple";

    /* Get some information about this guy */

    isdummy = (casu_fits_get_status(obj) != CASU_OK);

    /* If we need to make a PHU then do that now based on the first frame
       in the input frame list. Get rid of the file if it already exists */

    if (isfirst) {
        if (access(fname,F_OK)) 
            remove(fname);

        /* Create a new product frame object and define some tags */

        *product_frame = cpl_frame_new();
        cpl_frame_set_filename(*product_frame,fname);

        /* Tag the product */

        cpl_frame_set_tag(*product_frame,tag);
        cpl_frame_set_type(*product_frame,CPL_FRAME_TYPE_IMAGE);
        cpl_frame_set_group(*product_frame,CPL_FRAME_GROUP_PRODUCT);
        cpl_frame_set_level(*product_frame,CPL_FRAME_LEVEL_FINAL);

        /* Set some PhaseIII parameters */

        plist = casu_fits_get_phu(obj);
        cpl_propertylist_update_string(plist,"RADECSYS","ICRS");
        cpl_propertylist_update_string(plist,"ORIGIN","ESO-PARANAL");
        cpl_propertylist_set_comment(plist,"ORIGIN",
                                     "European Southern Observatory");
        cpl_propertylist_update_string(plist,"TELESCOP","ESO-VLT-U3");
        cpl_propertylist_set_comment(plist,"TELESCOP","ESO telescope name");
        cpl_propertylist_update_string(plist,"INSTRUME","VIMOS");
        cpl_propertylist_set_comment(plist,"INSTRUME","Instrument used");
        cpl_propertylist_update_string(plist,"OBSTECH","IMAGE");
        cpl_propertylist_set_comment(plist,"OBSTECH","Observation Technique");
        cpl_propertylist_update_string(plist,"IMATYPE","PAWPRINT");
        cpl_propertylist_update_bool(plist,"ISAMP",1);
        cpl_propertylist_set_comment(plist,"ISAMP",
                                     "TRUE if image represents partially sampled sky");
        cpl_propertylist_update_bool(plist,"SINGLEXP",1);
        cpl_propertylist_set_comment(plist,"SINGLEXP",
                                     "TRUE if resulting from a single exposure");
        cpl_propertylist_update_string(plist,"PROV1",
                                       cpl_propertylist_get_string(plist,"ARCFILE"));
        cpl_propertylist_set_comment(plist,"PROV1","Originating raw science file");
        cpl_propertylist_update_int(plist,"NCOMBINE",4);
        cpl_propertylist_set_comment(plist,"NCOMBINE","Number of raw files");
        cpl_propertylist_update_int(plist,"NSTACK",1);
        cpl_propertylist_set_comment(plist,"NSTACK","Number of stacked images");
        if (strcmp(tag,VIMOS_PRO_SIMPLE_SCI) == 0) {
            cpl_propertylist_update_string(plist,"PRODCATG","SCIENCE.MEFIMAGE");
            cpl_propertylist_set_comment(plist,"PRODCATG","Data product category");
            cpl_propertylist_update_string(plist,"ASSON1",assoc);
            cpl_propertylist_set_comment(plist,"ASSON1","Associated file");
            cpl_propertylist_update_string(plist,"ASSOC1","ANCILLARY.VARMAP");
            cpl_propertylist_set_comment(plist,"ASSOC1","Associated file category");
        } else {
            cpl_propertylist_erase(plist,"PRODCATG");
            cpl_propertylist_erase(plist,"ASSON1");
            cpl_propertylist_erase(plist,"ASSON2");
            cpl_propertylist_erase(plist,"ASSOC1");
            cpl_propertylist_erase(plist,"ASSOC2");
        }
        vimos_pfits_get_filter(plist,filt);
        cpl_propertylist_update_string(plist,"FILTER",filt);
        cpl_propertylist_set_comment(plist,"FILTER","Filter used in observation"); 
        vimos_pfits_get_exptime(plist,&texp);
        cpl_propertylist_update_double(plist,"EXPTIME",(double)texp);
        cpl_propertylist_update_double(plist,"TEXPTIME",(double)texp);
        cpl_propertylist_update_double(plist,"EFF_EXPT",(double)texp);
        cpl_propertylist_update_double(plist,"MJD-END",
                                       cpl_propertylist_get_double(plist,"MJD-OBS") +
                                       (double)texp/86400.0);
        cpl_propertylist_set_comment(plist,"MJD-END","End of observations");
        cpl_propertylist_update_string(plist,"PROG_ID",
                                       cpl_propertylist_get_string(plist,"ESO OBS PROG ID"));
        cpl_propertylist_set_comment(plist,"PROG_ID","ESO programme identification");
        cpl_propertylist_update_int(plist,"OBID1",
                                    cpl_propertylist_get_int(plist,"ESO OBS ID"));
        cpl_propertylist_set_comment(plist,"OBID1","Observation block ID");
        cpl_propertylist_update_bool(plist,"M_EPOCH",0);
        cpl_propertylist_set_comment(plist,"M_EPOCH",
                                     "TRUE if resulting from multiple epochs");
        cpl_propertylist_update_string(plist,"REFERENC","");
        cpl_propertylist_set_comment(plist,"REFERENC","Bibliographic Reference");
        cpl_propertylist_update_string(plist,"FLUXCAL","UNCALIBRATED");
        cpl_propertylist_set_comment(plist,"FLUXCAL","Certifies the validity of PHOTZP");
        cpl_propertylist_update_double(plist,"DIT",(double)texp);
        cpl_propertylist_set_comment(plist,"DIT","Detector integration time");

        /* Set up the PHU header */

        vimos_dfs_set_product_primary_header(plist,*product_frame,framelist,
                                             parlist,VIMOS_RECIPENAME,
                                             "PRO-1.15",template,1);

        /* 'Save' the PHU image */                         

        if (cpl_image_save(NULL,fname,CPL_TYPE_UCHAR,plist,CPL_IO_DEFAULT) != 
            CPL_ERROR_NONE) {
            cpl_msg_error(fctid,"Cannot save product PHU");
            cpl_frame_delete(*product_frame);
            return(CASU_FATAL);
        }
        cpl_frameset_insert(framelist,*product_frame);
    }

    /* Get the extension property list */

    plist = casu_fits_get_ehu(obj);
    if (isdummy)
        casu_dummy_property(plist);

    /* Add some PhaseIII stuff */

    cpl_propertylist_update_string(plist,"BUNIT","ADU");
    cpl_propertylist_set_comment(plist,"BUNIT","Physical unit of array values");

    /* Fiddle with the header now */

    cpl_propertylist_update_string(plist,"EXTNAME",cname);
    vimos_dfs_set_product_exten_header(plist,*product_frame,framelist,parlist,
                                       VIMOS_RECIPENAME,"PRO-1.15",template);
    if (cpl_image_save(casu_fits_get_image(obj),fname,CPL_TYPE_FLOAT,plist,
                       CPL_IO_EXTEND) != CPL_ERROR_NONE) {
        cpl_msg_error(fctid,"Cannot save product image extension -- %s",
                      cpl_error_get_message());
        return(CASU_FATAL);
    }

    /* Get out of here */

    return(CASU_OK);
}

/*---------------------------------------------------------------------------*/
/**
    \par Name:
        vimos_sci_save_stack
    \par Purpose:
        Save an extension of stacked image data products
    \par Description:
        An extension of a stacked image data product is saved here. 
    \par Language:
        C
    \param stack
        The casu_fits structure for the stack to be saved
    \param framelist
        The input recipe frameset
    \param parlist
        The input recipe parameter list
    \param template
        A template to be used in the cpl saving routine
    \param fnametype
        The type of file name. This could either be 0 -> ESO 'predictable'
        names, 1 -> pretty names based on input file names, 2 -> temporary
        file names.   
    \param cname
        The detector name
    \param prodtype
        Either STACK_FILE or STACK_VAR depending on whether this is the
        stacked image or its variance map.
    \param assoc
        A list of the associated confidence maps and variance frames for the
        stacked images.
    \param product_frame
        The output product frame
    \retval CASU_OK if everything is ok
    \retval CASU_FATAL if something failed
    \author
        Jim Lewis, CASU
 */
/*---------------------------------------------------------------------------*/

static int vimos_sci_save_stack(casu_fits *stack, cpl_frameset *framelist,
                                cpl_parameterlist *parlist, 
                                cpl_frame *template, int fnametype,
                                char *cname, int prodtype, char *assoc[],
                                char *photosys, cpl_frame **product_frame) {
    cpl_propertylist *plist,*prov;
    int isfirst,isdummy,n,np;
    char bname[BUFSIZ],filt[32],*base,*tname,prop[9];
    float texp;
    const char *fctid = "vimos_sci_save_stack";

    /* Get some information about this guy */

    isdummy = (casu_fits_get_status(stack) != CASU_OK);
    isfirst = (*product_frame == NULL);
    tname = cpl_strdup(cpl_frame_get_filename(template));
    base = basename(tname);
    vimos_sci_product_name(base,prodtype,fnametype,1,bname);
    freespace(tname);

    /* If we need to make a PHU then do that now based on the first frame
       in the input frame list. Get rid of the file if it already exists */

    if (isfirst) {
        if (access(bname,F_OK)) 
            remove(bname);

        /* Create a new product frame object and define some tags */

        *product_frame = cpl_frame_new();
        cpl_frame_set_filename(*product_frame,bname);

        /* Tag this product */

        const cpl_boolean is_sci = prodtype == STACK_FILE;
        if (is_sci)
            cpl_frame_set_tag(*product_frame,VIMOS_PRO_JITTERED_SCI);
        else
            cpl_frame_set_tag(*product_frame,VIMOS_PRO_JITTERED_VAR);
        cpl_frame_set_type(*product_frame,CPL_FRAME_TYPE_IMAGE);
        cpl_frame_set_group(*product_frame,CPL_FRAME_GROUP_PRODUCT);
        cpl_frame_set_level(*product_frame,CPL_FRAME_LEVEL_FINAL);

        /* Set some PhaseIII headers */

        plist = casu_fits_get_phu(stack);
        cpl_propertylist_update_string(plist,"RADECSYS","ICRS");
        cpl_propertylist_update_string(plist,"ORIGIN","ESO-PARANAL");
        cpl_propertylist_set_comment(plist,"ORIGIN",
                                     "European Southern Observatory");
        cpl_propertylist_update_string(plist,"TELESCOP","ESO-VLT-U3");
        cpl_propertylist_set_comment(plist,"TELESCOP","ESO telescope name");
        cpl_propertylist_update_string(plist,"INSTRUME","VIMOS");
        cpl_propertylist_set_comment(plist,"INSTRUME","Instrument used");
        cpl_propertylist_update_string(plist,"OBSTECH","IMAGE");
        cpl_propertylist_set_comment(plist,"OBSTECH","Observation Technique");
        cpl_propertylist_update_string(plist,"IMATYPE","PAWPRINT");
        cpl_propertylist_update_bool(plist,"ISAMP",1);
        cpl_propertylist_set_comment(plist,"ISAMP",
                                     "TRUE if image represents partially sampled sky");
        n = cpl_propertylist_get_int(plist,"NSTACK");
        if (n != 1) 
            cpl_propertylist_update_bool(plist,"SINGLEXP",0);
        else
            cpl_propertylist_update_bool(plist,"SINGLEXP",1);
        cpl_propertylist_set_comment(plist,"SINGLEXP",
                                     "TRUE if resulting from a single exposure");        

        const char * prod_catg = NULL;
        if (is_sci) {
        	prod_catg = "SCIENCE.MEFIMAGE";
            cpl_propertylist_update_string(plist,"ASSON1",assoc[0]);
            cpl_propertylist_set_comment(plist,"ASSON1","Associated file");
            cpl_propertylist_update_string(plist,"ASSON2",assoc[1]);
            cpl_propertylist_set_comment(plist,"ASSON2","Associated file");
        } else {
            prod_catg = "ANCILLARY.VARMAP";
            cpl_propertylist_erase(plist,"ASSON1");
            cpl_propertylist_erase(plist,"ASSON2");
        }
        cpl_propertylist_update_string(plist,"PRODCATG",prod_catg);
        cpl_propertylist_set_comment(plist,"PRODCATG","Data product category");
        cpl_propertylist_erase(plist,"ASSOC1");
        cpl_propertylist_erase(plist,"ASSOC2");

        vimos_pfits_get_filter(casu_fits_get_ehu(stack),filt);
        cpl_propertylist_update_string(plist,"FILTER",filt);
        cpl_propertylist_set_comment(plist,"FILTER","Filter used in observation");
        vimos_pfits_get_exptime(plist,&texp);
        cpl_propertylist_update_double(plist,"DIT",(double)texp / (double)n);
        cpl_propertylist_update_double(plist,"EFF_EXPT",(double)texp / (double)n);
        cpl_propertylist_set_comment(plist,"DIT","Detector integration time");
        cpl_propertylist_update_double(plist,"EXPTIME",(double)texp);
        cpl_propertylist_set_comment(plist,"EXPTIME","Total integration time.");
        cpl_propertylist_update_double(plist,"TEXPTIME",(double)texp);
        cpl_propertylist_update_string(plist,"PROG_ID",
                                       cpl_propertylist_get_string(plist,"ESO OBS PROG ID"));
        cpl_propertylist_set_comment(plist,"PROG_ID","ESO programme identification");
        cpl_propertylist_update_int(plist,"OBID1",
                                    cpl_propertylist_get_int(plist,"ESO OBS ID"));
        cpl_propertylist_set_comment(plist,"OBID1","Observation block ID");
        cpl_propertylist_update_bool(plist,"M_EPOCH",0);
        cpl_propertylist_set_comment(plist,"M_EPOCH",
                                     "TRUE if resulting from multiple epochs");
        cpl_propertylist_update_string(plist,"REFERENC","");
        cpl_propertylist_set_comment(plist,"REFERENC","Bibliographic Reference");
        cpl_propertylist_update_string(plist,"FLUXCAL","ABSOLUTE");
        cpl_propertylist_set_comment(plist,"FLUXCAL",
                                     "Certifies the validity of PHOTZP");

        /* Set up the PHU header */

        vimos_dfs_set_product_primary_header(plist,*product_frame,
                                             framelist,parlist,
                                             VIMOS_RECIPENAME,
                                             "PRO-1.15",template,1);

        /* Fix provenance */
        
        prov = cpl_propertylist_new();
        cpl_propertylist_copy_property_regexp(prov,plist,
                                              "ESO PRO REC1 RAW[0-9]* NAME",0);
        np = cpl_propertylist_get_size(prov);
        for (n = 0; n < np; n++) {
            (void)snprintf(prop,8,"PROV%d",n+1);
            cpl_property_set_name(cpl_propertylist_get(prov,n),prop);
        }
        casu_merge_propertylists(plist,prov);
        cpl_propertylist_delete(prov);

        /* 'Save' the PHU image */                         

        if (cpl_image_save(NULL,bname,CPL_TYPE_UCHAR,plist,
                           CPL_IO_DEFAULT) != CPL_ERROR_NONE) {
            cpl_msg_error(fctid,"Cannot save product PHU");
            cpl_frame_delete(*product_frame);
            return(CASU_FATAL);
        }
        cpl_frameset_insert(framelist,*product_frame);
    }

    /* Get the extension property list */

    plist = casu_fits_get_ehu(stack);
    if (isdummy)
        casu_dummy_property(plist);

    /* Add some PhaseIII stuff */

    cpl_propertylist_update_string(plist,"BUNIT","ADU");
    cpl_propertylist_set_comment(plist,"BUNIT","Physical unit of array values");
    if (cpl_propertylist_has(plist,"PHOTSYS") == 0) {
        cpl_propertylist_update_string(plist,"PHOTSYS",photosys);
        cpl_propertylist_set_comment(plist,"PHOTSYS","Photometric System");
    }

    if (cpl_propertylist_has(plist,"EXPTIME")) {
    	cpl_propertylist_set_comment(plist, "EXPTIME", "Total integration time.");
    }

    /* Fiddle with the header now */

    cpl_propertylist_update_string(plist,"EXTNAME",cname);
    vimos_dfs_set_product_exten_header(plist,*product_frame,framelist,
                                       parlist,VIMOS_RECIPENAME,
                                       "PRO-1.15",template);
    if (cpl_image_save(casu_fits_get_image(stack),bname,CPL_TYPE_FLOAT,
                       plist,CPL_IO_EXTEND) != CPL_ERROR_NONE) {
        cpl_msg_error(fctid,"Cannot save product image extension -- %s",
                      cpl_error_get_message());
        return(CASU_FATAL);
    }

    /* Get out of here */

    return(CASU_OK);
}

/*---------------------------------------------------------------------------*/
/**
    \par Name:
        merge_with_phu
    \par Purpose:
		Append one list at the end of the other avoiding duplications
    \par Description:
 		Appends other_list to out_list, if out_list already contains a property with a given name,
 		the property is removed from out_list before appending the corresponding property from other_list
    \par Language:
        C
    \param header
        The destination list
    \param phu
        The source list
    \author
        Matteo Salmistraro, ESO
 */
/*---------------------------------------------------------------------------*/
static void merge_with_phu(cpl_propertylist * header, const cpl_propertylist * phu){

	cpl_size i = 0;
	for(i = 0; i < cpl_propertylist_get_size(phu); ++i){
		const cpl_property * p = cpl_propertylist_get_const(phu, i);
		const char * name = cpl_property_get_name(p);
		if(cpl_propertylist_has(header, name)){
			cpl_propertylist_erase(header, name);
		}
		cpl_propertylist_append_property(header, p);
	}
}

/*---------------------------------------------------------------------------*/
/**
    \par Name:
        vimos_sci_save_stack_conf
    \par Purpose:
        Save an extension of stacked image confidence map data products
    \par Description:
        An extension of a stacked image confidence map  data product is saved 
        here. 
    \par Language:
        C
    \param stack
        The casu_fits structure for the stack confidence map to be saved
    \param framelist
        The input recipe frameset
    \param parlist
        The input recipe parameter list
    \param template
        A template to be used in the cpl saving routine
    \param fnametype
        The type of file name. This could either be 0 -> ESO 'predictable'
        names, 1 -> pretty names based on input file names, 2 -> temporary
        file names.   
    \param cname
        The detector name
    \param product_frame
        The output product frame
    \retval CASU_OK if everything is ok
    \retval CASU_FATAL if something failed
    \author
        Jim Lewis, CASU
 */
/*---------------------------------------------------------------------------*/

static int vimos_sci_save_stack_conf(casu_fits *stack, cpl_frameset *framelist,
                                     cpl_parameterlist *parlist, 
                                     cpl_frame *template, int fnametype,
                                     char *cname, cpl_frame **product_frame) {
    cpl_propertylist *plist,*prov;
    int isfirst,isdummy,np,n;
    char bname[BUFSIZ],*base,*tname,filt[32],prop[9];
    const char *fctid = "vimos_sci_save_stack_conf";

    /* Get some information about this guy */

    isdummy = (casu_fits_get_status(stack) != CASU_OK);
    isfirst = (*product_frame == NULL);
    tname = cpl_strdup(cpl_frame_get_filename(template));
    base = basename(tname);
    vimos_sci_product_name(base,STACK_CONF,fnametype,1,bname);
    freespace(tname);

    /* If we need to make a PHU then do that now based on the first frame
       in the input frame list. Get rid of the file if it already exists */

    if (isfirst) {
        if (access(bname,F_OK)) 
            remove(bname);

        /* Create a new product frame object and define some tags */

        *product_frame = cpl_frame_new();
        cpl_frame_set_filename(*product_frame,bname);

        /* Tag the output frame */

        cpl_frame_set_tag(*product_frame,VIMOS_PRO_CONF_SCI);
        cpl_frame_set_type(*product_frame,CPL_FRAME_TYPE_IMAGE);
        cpl_frame_set_group(*product_frame,CPL_FRAME_GROUP_PRODUCT);
        cpl_frame_set_level(*product_frame,CPL_FRAME_LEVEL_FINAL);

        /* Set up the PHU header */
        const char * fname = cpl_frame_get_filename(template);
        plist = cpl_propertylist_load(fname,
                                      cpl_frame_get_nextensions(template));

        merge_with_phu(plist, casu_fits_get_phu(stack));

        vimos_dfs_set_product_primary_header(plist,*product_frame,framelist,
                                             parlist,VIMOS_RECIPENAME,
                                             "PRO-1.15",template,1);
        vimos_pfits_get_filter(casu_fits_get_ehu(stack),filt);
        cpl_propertylist_update_string(plist,"FILTER",filt);
        cpl_propertylist_set_comment(plist,"FILTER","Filter used in observation");
        cpl_propertylist_erase(plist,"PRODCATG");
        cpl_propertylist_erase(plist,"ASSON1");
        cpl_propertylist_erase(plist,"ASSON2");
        cpl_propertylist_erase(plist,"ASSOC1");
        cpl_propertylist_erase(plist,"ASSOC2");
        cpl_propertylist_erase(plist,"NCOMBINE");
        cpl_propertylist_erase(plist,"SINGLEXP");
        cpl_propertylist_update_string(plist, "PRODCATG", "ANCILLARY.WEIGHTMAP");
        cpl_propertylist_set_comment(plist,"PRODCATG","Data product category");

        /* Fix provenance */
        
        prov = cpl_propertylist_new();
        cpl_propertylist_copy_property_regexp(plist,prov,
                                              "ESO PRO REC1 RAW[0-9]* NAME",0);
        np = cpl_propertylist_get_size(prov);
        for (n = 0; n < np; n++) {
            (void)snprintf(prop,8,"PROV%d",n);
            cpl_property_set_name(cpl_propertylist_get(prov,n),prop);
        }
        casu_merge_propertylists(plist,prov);
        cpl_propertylist_delete(prov);

        /* 'Save' the PHU image */                         

        if (cpl_image_save(NULL,bname,CPL_TYPE_UCHAR,plist,
                           CPL_IO_DEFAULT) != CPL_ERROR_NONE) {
            cpl_msg_error(fctid,"Cannot save product PHU");
            cpl_propertylist_delete(plist);
            cpl_frame_delete(*product_frame);
            return(CASU_FATAL);
        }
        cpl_frameset_insert(framelist,*product_frame);
        cpl_propertylist_delete(plist);
    }

    /* Get the extension property list */

    plist = casu_fits_get_ehu(stack);
    if (isdummy)
        casu_dummy_property(plist);

    /* Fiddle with the header now */

    cpl_propertylist_update_string(plist,"EXTNAME",cname);
    vimos_dfs_set_product_exten_header(plist,*product_frame,framelist,
                                       parlist,VIMOS_RECIPENAME,
                                       "PRO-1.15",template);
    if (cpl_image_save(casu_fits_get_image(stack),bname,CPL_TYPE_FLOAT,
                       plist,CPL_IO_EXTEND) != CPL_ERROR_NONE) {
        cpl_msg_error(fctid,"Cannot save product image extension -- %s",
                      cpl_error_get_message());
        return(CASU_FATAL);
    }

    /* Get out of here */

    return(CASU_OK);
}

/*---------------------------------------------------------------------------*/
/**
    \par Name:
        vimos_sci_save_cat
    \par Purpose:
        Save an extension of image source catalogue data products
    \par Description:
        An extension of an image source catalogue data product is 
        saved here. 
    \par Language:
        C
    \param scat
        The casu_tfits structure for the source catalogue to be saved
    \param framelist
        The input recipe frameset
    \param parlist
        The input recipe parameter list
    \param template
        A template to be used in the cpl saving routine
    \param fnametype
        The type of file name. This could either be 0 -> ESO 'predictable'
        names, 1 -> pretty names based on input file names, 2 -> temporary
        file names.   
    \param ptype
        The product type. Either SIMPLE_CAT, STACK_CAT, SIMPLE_MSTD_A,
        STACK_MSTD_A or STACK_MSTD_P
    \param fnumber
        The file number if the ESO product name option is chosen
    \param cname
        The detector name
    \param photosys
        The photometric system
    \param product_frame
        The output product frame
    \retval CASU_OK if everything is ok
    \retval CASU_FATAL if something failed
    \author
        Jim Lewis, CASU
 */
/*---------------------------------------------------------------------------*/

static int vimos_sci_save_cat(casu_tfits *scat, cpl_frameset *framelist,
                              cpl_parameterlist *parlist, 
                              cpl_frame *template, int fnametype, 
                              int ptype, int fnumber, char *cname, 
                              char *photosys, cpl_frame **product_frame) {
    cpl_propertylist *plist;
    int isfirst,isdummy,n;
    char bname[BUFSIZ],*base,*tname;
    float texp;
    const char *fctid = "vimos_sci_save_cat";

    /* Get some information about this guy */

    isdummy = (casu_tfits_get_status(scat) != CASU_OK);
    isfirst = (*product_frame == NULL);
    tname = cpl_strdup(cpl_frame_get_filename(template));
    base = basename(tname);
    vimos_sci_product_name(base,ptype,fnametype,fnumber,bname);
    freespace(tname);

    /* If we need to make a PHU then do that now based on the first frame
       in the input frame list. Get rid of the file if it already exists */

    if (isfirst) {
        if (access(bname,F_OK)) 
            remove(bname);

        /* Create a new product frame object and define some tags */

        *product_frame = cpl_frame_new();
        cpl_frame_set_filename(*product_frame,bname);

        /* Tag the output frame */

        switch (ptype) {
        case SIMPLE_CAT:
            cpl_frame_set_tag(*product_frame,VIMOS_PRO_OBJCAT_SCI);
            break;
        case STACK_CAT:
            cpl_frame_set_tag(*product_frame,VIMOS_PRO_OBJCAT_SCI_JITTERED);
            break;
        case SIMPLE_MSTD_A:
            cpl_frame_set_tag(*product_frame,VIMOS_PRO_MATCHSTD_ASTROM);
            break;
        case STACK_MSTD_A:
            cpl_frame_set_tag(*product_frame,VIMOS_PRO_MATCHSTD_ASTROM);
            break;
        case STACK_MSTD_P:
            cpl_frame_set_tag(*product_frame,VIMOS_PRO_MATCHSTD_PHOTOM);
            break;
        }        
        cpl_frame_set_type(*product_frame,CPL_FRAME_TYPE_TABLE);
        cpl_frame_set_group(*product_frame,CPL_FRAME_GROUP_PRODUCT);
        cpl_frame_set_level(*product_frame,CPL_FRAME_LEVEL_FINAL);

        /* Set some PhaseIII parameters */

        plist = casu_tfits_get_phu(scat);
        cpl_propertylist_update_string(plist,"RADECSYS","ICRS");
        cpl_propertylist_update_string(plist,"ORIGIN","ESO-PARANAL");
        cpl_propertylist_set_comment(plist,"ORIGIN",
                                     "European Southern Observatory");
        cpl_propertylist_update_string(plist,"TELESCOP","ESO-VLT-U3");
        cpl_propertylist_set_comment(plist,"TELESCOP","ESO telescope name");
        cpl_propertylist_update_string(plist,"INSTRUME","VIMOS");
        cpl_propertylist_set_comment(plist,"INSTRUME","Instrument used");
        cpl_propertylist_update_string(plist,"PRODCATG","SCIENCE.SRCTBL");
        cpl_propertylist_set_comment(plist,"PRODCATG","Data product category");
        cpl_propertylist_update_string(plist,"OBSTECH","IMAGE");
        cpl_propertylist_set_comment(plist,"OBSTECH","Observation Technique");
        cpl_propertylist_update_string(plist,"PROG_ID",
                                       cpl_propertylist_get_string(plist,"ESO OBS PROG ID"));
        cpl_propertylist_set_comment(plist,"PROG_ID","ESO programme identification");
        cpl_propertylist_update_int(plist,"OBID1",
                                    cpl_propertylist_get_int(plist,"ESO OBS ID"));
        cpl_propertylist_set_comment(plist,"OBID1","Observation block ID");
        cpl_propertylist_update_bool(plist,"M_EPOCH",0);
        cpl_propertylist_set_comment(plist,"M_EPOCH",
                                     "TRUE if resulting from multiple epochs");
        cpl_propertylist_update_string(plist,"REFERENC","");
        cpl_propertylist_set_comment(plist,"REFERENC","Bibliographic Reference");
        n = cpl_propertylist_get_int(plist,"NSTACK");
        if (n != 1) 
            cpl_propertylist_update_bool(plist,"SINGLEXP",0);
        else
            cpl_propertylist_update_bool(plist,"SINGLEXP",1);
        cpl_propertylist_set_comment(plist,"SINGLEXP",
                                     "TRUE if resulting from a single exposure");        
        if (ptype == STACK_CAT) {
            vimos_pfits_get_exptime(plist,&texp);
            cpl_propertylist_update_double(plist,"EFF_EXPT",(double)texp / (double)n);
            cpl_propertylist_update_double(plist,"EXPTIME",(double)texp);
            cpl_propertylist_set_comment(plist,"EXPTIME", "Total integration time.");
            cpl_propertylist_update_double(plist,"TEXPTIME",(double)texp);
        }
        cpl_propertylist_erase(plist,"ASSON1");
        cpl_propertylist_erase(plist,"ASSON2");
        cpl_propertylist_erase(plist,"ASSOC1");
        cpl_propertylist_erase(plist,"ASSOC2");       

        /* Set up the PHU header */

        vimos_dfs_set_product_primary_header(plist,*product_frame,framelist,
                                             parlist,VIMOS_RECIPENAME,
                                             "PRO-1.15",NULL,0);

        /* 'Save' the PHU image */                         

        if (cpl_image_save(NULL,bname,CPL_TYPE_UCHAR,plist,
                           CPL_IO_DEFAULT) != CPL_ERROR_NONE) {
            cpl_msg_error(fctid,"Cannot save product PHU");
            cpl_frame_delete(*product_frame);
            return(CASU_FATAL);
        }
        cpl_frameset_insert(framelist,*product_frame);
    }

    /* Get the extension property list */

    plist = casu_tfits_get_ehu(scat);
    if (isdummy)
        casu_dummy_property(plist);

    /* Fiddle with the header now */

    if (cpl_propertylist_has(plist,"PHOTSYS") == 0) {
        cpl_propertylist_update_string(plist,"PHOTSYS",photosys);
        cpl_propertylist_set_comment(plist,"PHOTSYS","Photometric System");
    }
    cpl_propertylist_erase(plist,"BUNIT");
    cpl_propertylist_update_string(plist,"EXTNAME",cname);
    vimos_dfs_set_product_exten_header(plist,*product_frame,framelist,
                                       parlist,VIMOS_RECIPENAME,
                                       "PRO-1.15",NULL);
    if (cpl_table_save(casu_tfits_get_table(scat),NULL,
                       plist,bname,CPL_IO_EXTEND) != CPL_ERROR_NONE) {
        cpl_msg_error(fctid,"Cannot save product table extension -- %s",
                      cpl_error_get_message());
        return(CASU_FATAL);
    }

    /* Get out of here */

    return(CASU_OK);
}

/*---------------------------------------------------------------------------*/
/**
    \par Name:
        vimos_sci_product_name
    \par Purpose:
        Set up an output product file name
    \par Description:
        An output file name is defined based on a template name,  a
        'predictable' name or a temporary file name.
    \par Language:
        C
    \param template
        If the pretty file name option is used or if this is a temporary
        file name, then the this will be used as a reference name.
    \param producttype
        The type of product you are writing out. These are enumerated above.
    \param nametype
        This is: 0 -> predictable names, 1 -> pretty names based on input
        file names or 2 -> temporary file names.
    \param fnumber
        If the predictable file names are going to be used, then this is 
        a number that is appended to a file name root.
    \param outfname
        The output file name
    \return
        Nothing
    \author
        Jim Lewis, CASU
 */
/*---------------------------------------------------------------------------*/

static void vimos_sci_product_name(char *template, int producttype, 
                                   int nametype, int fnumber, char *outfname) {
    const char *esonames[] = {"exp_","exp_var_","exp_cat_","stack_",
                              "stack_conf_","stack_var_","stack_cat_",
                              "exp_mstd_a_","stack_mstd_a_",
                              "stack_mstd_p_"};
    const char *suffix[] = {"_ex","_ex_var","_ex_cat","_st","_st_conf",
                            "_st_var","_st_cat","_ex_mstd_a",
                            "_st_mstd_a","_st_mstd_p"};
    char *fname,*bname,*dot;

    /* If the name type is the ESO predictable sort, then it's easy. Just
       use the esonames defined above and append the correct number */

    switch (nametype) {
    case 0:
        (void)sprintf(outfname,"%s%d.fits",esonames[producttype],fnumber);
        break;

    /* If this is a temporary file, then just append tmp_ to the template
       name and return that */

    case 2:
        fname = cpl_strdup(template);
        bname = basename(fname);
        (void)sprintf(outfname,"tmp_%s",bname);
        freespace(fname);
        break;

    /* Ok, we want a pretty name... */

    case 1:
        fname = cpl_strdup(template);
        bname = basename(fname);
        (void)sprintf(outfname,"%s",bname);
        dot = strrchr(outfname,'.');
        (void)sprintf(dot,"%s.fits",suffix[producttype]);
        freespace(fname);
        break;

    /* something else ?? */

    default:
        (void)strcpy(outfname,"");
        break;
    }
    return;
}   

/*---------------------------------------------------------------------------*/
/**
    \par Name:
        vimos_sci_update_hdr
    \par Purpose:
        Update header of simple products to give them photometric info
    \par Description:
        Simple products are not photometrically calibrated. This routine
        gives the simple products the photometric information that is
        calculated for the stacked images, which will be roughly right.
        This is provided as a stop-gap until CPL can provide a way to update
        the header of a file without having to re-write the whole file.
    \par Language:
        C
    \param simple
        A frameset of simple products
    \param stack
        The stacked image
    \param typef
        The type of file: 1 => MEF of float images, 2 => MEF of binary
        tables
    \return
        Nothing
    \author
        Jim Lewis, CASU
 */
/*---------------------------------------------------------------------------*/

static void vimos_sci_update_hdr(cpl_frameset *simple, cpl_frame *stack, 
                                 int typef) {
    int i,j,k,np_prim=1,np_extn=4,nsimple;
    const char *wanted_primary[] = {"FLUXCAL"};
    const char *wanted_extn[] = {"PHOTSYS","ABMAGSAT","PHOTZP","ABMAGLIM"};
    char tmpfname[BUFSIZ];
    cpl_propertylist *p,*new_primary,*new_extn[VIMOS_NEXTN];
    cpl_frame *fr;
    cpl_image *im;
    cpl_table *tab;

    /* First of all extract the stuff you really want from the primary header
       of the stack */

    new_primary = cpl_propertylist_new();
    p = cpl_propertylist_load(cpl_frame_get_filename(stack),0);
    for (i = 0; i < np_prim; i++) 
        cpl_propertylist_copy_property(new_primary,p,wanted_primary[i]); 
    cpl_propertylist_delete(p);

    /* Now from the extensions */

    for (j = 1; j <= VIMOS_NEXTN; j++) {
        p = cpl_propertylist_load(cpl_frame_get_filename(stack),(cpl_size)j);
        new_extn[j-1] = cpl_propertylist_new();
        for (i = 0; i < np_extn; i++) 
            cpl_propertylist_copy_property(new_extn[j-1],p,wanted_extn[i]);
        cpl_propertylist_delete(p);
    }
        
    /* Loop for each file. Define a temporary file name and load the 
       primary header */

    nsimple = cpl_frameset_get_size(simple);
    for (k = 0; k < nsimple; k++) {
        fr = cpl_frameset_get_position(simple,k);
        (void)sprintf(tmpfname,"scitmp_%s",cpl_frame_get_filename(fr));
        if (access(tmpfname,F_OK)) 
            remove(tmpfname);
        p = cpl_propertylist_load(cpl_frame_get_filename(fr),0);

        /* Update the header and write it out */

        casu_merge_propertylists(p,new_primary);
        cpl_image_save(NULL,tmpfname,CPL_TYPE_UCHAR,p,CPL_IO_DEFAULT);
        cpl_propertylist_delete(p);

        /* Loop for each extension now and do exactly the same thing. The
           only difference now is that we have to load either the image
           or table so that we can write it back out again */

        for (i = 1; i <= VIMOS_NEXTN; i++) {
            p = cpl_propertylist_load(cpl_frame_get_filename(fr),i);
            casu_merge_propertylists(p,new_extn[i-1]);
            if (typef == 1) {
                im = cpl_image_load(cpl_frame_get_filename(fr),
                                    CPL_TYPE_FLOAT,0,i);
                cpl_image_save(im,tmpfname,CPL_TYPE_FLOAT,p,CPL_IO_EXTEND);
                cpl_image_delete(im);
            } else {
                tab = cpl_table_load(cpl_frame_get_filename(fr),i,0);
                cpl_table_save(tab,NULL,p,tmpfname,CPL_IO_EXTEND);
                cpl_table_delete(tab);
            }
            cpl_propertylist_delete(p);
        }

        /* Now get rid of the old file and rename the new one */

        remove(cpl_frame_get_filename(fr));
        rename(tmpfname,cpl_frame_get_filename(fr));
    }

    /* Delete some stuff and get out of here */

    cpl_propertylist_delete(new_primary);
    for (i = 0; i < VIMOS_NEXTN; i++)
        cpl_propertylist_delete(new_extn[i]);
}

/*---------------------------------------------------------------------------*/
/**
    \par Name:
        vimos_sci_update_hdr_stack
    \par Purpose:
        Update header of stack products to give them correct central 
        coordinates
    \par Description:
        The WCS of the stack images are used to define a true central
        point for the coverage of the current stack. This is written 
        to the RA and Dec keywords in the primary of the products
    \par Language:
        C
    \param stack
        The pawprint structure for the current stack
    \return
        Nothing
    \author
        Jim Lewis, CASU
 */
/*---------------------------------------------------------------------------*/

static void vimos_sci_update_hdr_stack(casu_fits *stack[], casu_fits *conf,
                                       casu_fits *var, casu_tfits *cat) {
    double ra1,dec1,ra2,dec2,ramin,ramax,decmin,decmax;
    int i,status;
    cpl_propertylist *p[4];

    /* For each extension of the input stack, get the extension header 
       and work out the central RA and Dec. The NAXISn keywords haven't
       been updated yet as the stack hasn't yet been saved. */

    ramin = 370.0;
    ramax = -10.0;
    decmin = 95.0;
    decmax = -95.0;
    for (i = 1; i <= VIMOS_NEXTN; i++) {
        status = CASU_OK;
        cpl_propertylist_update_int(casu_fits_get_ehu(stack[i-1]),"NAXIS1",
                                    cpl_image_get_size_x(casu_fits_get_image(stack[i-1])));
        cpl_propertylist_update_int(casu_fits_get_ehu(stack[i-1]),"NAXIS2",
                                    cpl_image_get_size_y(casu_fits_get_image(stack[i-1])));
        cpl_propertylist_update_int(casu_fits_get_ehu(stack[i-1]),"ZNAXIS1",
                                    cpl_image_get_size_x(casu_fits_get_image(stack[i-1])));
        cpl_propertylist_update_int(casu_fits_get_ehu(stack[i-1]),"ZNAXIS2",
                                    cpl_image_get_size_y(casu_fits_get_image(stack[i-1])));
        (void)casu_coverage(casu_fits_get_ehu(stack[i-1]),0,&ra1,&ra2,
                                              &dec1,&dec2,&status);
        //The following formulas ensure that there are no issues due to wrapping.
        //Still, at the end of the process, ramin could wrap over ramax.
        ramin = ramin - fabs(min(fmod(ra1 - ramin + 180 + 360 , 360) - 180, 0));
        ramax = ramax + max(fmod(ra2 - ramax + 180 + 360, 360) - 180, 0);
        decmin = min(decmin,dec1);
        decmax = max(decmax,dec2);
    }
    //Compute the mean of the two angles (https://en.wikipedia.org/wiki/Mean_of_circular_quantities)
    ra1 = atan2(0.5*(sin(ramin*CPL_MATH_RAD_DEG) + sin(ramax*CPL_MATH_RAD_DEG)), 
                0.5*(cos(ramin*CPL_MATH_RAD_DEG) + cos(ramax*CPL_MATH_RAD_DEG)))*CPL_MATH_DEG_RAD;
    if (ra1 < 0.0)
        ra1 += 360.0;
    dec1 = 0.5*(decmin + decmax);

    /* Load up the propertylists to be updated */

    p[0] = casu_fits_get_phu(stack[0]);
    p[1] = casu_fits_get_phu(conf);
    p[2] = casu_fits_get_phu(var);
    p[3] = casu_tfits_get_phu(cat);

    /* Loop for each one now and update the RA and Dec value */
        
    for (i = 0; i < 4; i++) {
        cpl_propertylist_update_double(p[i],"RA",ra1);
        cpl_propertylist_update_double(p[i],"DEC",dec1);
    }

//Refactor ALL this code in a separate function
    double psf[VIMOS_NEXTN], psf_mean=0, psf_stddev=0;
    int ngood=0;
    for (i = 0; i < VIMOS_NEXTN; i++)
    {
        psf[i]=cpl_propertylist_get_float(casu_fits_get_ehu(stack[i]),"PSF_FWHM");
        if(psf[i]>0)
            ngood++;
    }
    if(ngood>0 && ngood<VIMOS_NEXTN)
    {
        for (i = 0; i < VIMOS_NEXTN; i++) 
        {
            if(psf[i]>0)
            {
                psf_mean+=psf[i];
                psf_stddev+=psf[i]*psf[i];
            }
        }
        psf_mean = psf_mean / ngood;
        psf_stddev= ngood > 1? sqrt(psf_stddev / ngood - psf_mean * psf_mean) : 0.0;
        for (i = 0; i < VIMOS_NEXTN; i++)
        {
            if(psf[i]<=0)
            {
                cpl_propertylist_update_float(casu_fits_get_ehu(stack[i]),"PSF_FWHM", psf_mean);
                cpl_propertylist_set_comment(casu_fits_get_ehu(stack[i]),
                     "PSF_FWHM", "spatial res computed from ext avg [arcsec]");
                if(psf_stddev != 0.0) 
                {
                    cpl_propertylist_update_float(casu_fits_get_ehu(stack[i]),"PSF_FERR", psf_stddev);
                    cpl_propertylist_set_comment(casu_fits_get_ehu(stack[i]),
                         "PSF_FERR", "error in avg-computed PSF [arcsec]");
                }
                else
                {
                    cpl_propertylist_update_float(casu_fits_get_ehu(stack[i]),"PSF_FERR", 0.3);
                    cpl_propertylist_set_comment(casu_fits_get_ehu(stack[i]),
                         "PSF_FERR", "error in avg-computed PSF [arcsec]");
                }
            }
        }
    }
    else
    {
        for (i = 0; i < VIMOS_NEXTN; i++)
        {
            double fwhm_dimm_start = cpl_propertylist_get_float(casu_fits_get_phu(stack[i]),
                "ESO TEL AMBI FWHM START");
            double fwhm_dimm_end = cpl_propertylist_get_float(casu_fits_get_phu(stack[i]),
                "ESO TEL AMBI FWHM END");
            if(fwhm_dimm_start == -1 || fwhm_dimm_end == -1)
            {
                cpl_propertylist_update_float(casu_fits_get_ehu(stack[i]),"PSF_FWHM", -1);
                cpl_propertylist_set_comment(casu_fits_get_ehu(stack[i]),
                     "PSF_FWHM", "no spatial res could be computed");
                cpl_propertylist_update_float(casu_fits_get_ehu(stack[i]),"PSF_FERR", -1);
                cpl_propertylist_set_comment(casu_fits_get_ehu(stack[i]),
                     "PSF_FERR", "no spatial res error could be computed");
                
            }
            else
            {
                double airmass_start = cpl_propertylist_get_float(casu_fits_get_phu(stack[i]),
                    "ESO TEL AIRM START");
                double airmass_end = cpl_propertylist_get_float(casu_fits_get_phu(stack[i]),
                    "ESO TEL AIRM END");
                double airmass = (airmass_start + airmass_end) /2;
                double fwhm_dimm = (fwhm_dimm_start + fwhm_dimm_end) /2;
                char filter[16];
                vimos_pfits_get_filter(casu_fits_get_phu(stack[i]),filter);
                double wavelength = 0; //in nm
                double wavelength_dimm = 0.5; //in micrometer
                // http://www.eso.org/sci/facilities/paranal/decommissioned/vimos/doc/VLT-MAN-ESO-14610-3509_v98.pdf (page 121)
                if(strncmp(filter, "U", 1)==0)
                    wavelength=370.0;
                else if(strncmp(filter, "B", 1)==0)
                    wavelength=430.0;
                else if(strncmp(filter, "V", 1)==0)
                    wavelength=546.0;
                else if(strncmp(filter, "R", 1)==0)
                    wavelength=648.5;
                else if(strncmp(filter, "I", 1)==0)
                    wavelength=830.0;
                else if(strncmp(filter, "z", 1)==0)
                    wavelength=950.0;
                double k = sqrt(1.-78.08 * pow(wavelength * 1e-9, 0.4) * 
                                pow(airmass, -0.6) * pow(fwhm_dimm, -1./3.));
                double fwhm_ins = (fwhm_dimm * pow(wavelength_dimm, 0.2) * pow(airmass, 0.6) * k) / pow(wavelength*1.e-3, 0.2);
                cpl_propertylist_update_float(casu_fits_get_ehu(stack[i]),"PSF_FWHM", fwhm_ins);
                cpl_propertylist_set_comment(casu_fits_get_ehu(stack[i]),
                     "PSF_FWHM", "spatial res computed from DIMM seeing [arcsec]");
                cpl_propertylist_update_float(casu_fits_get_ehu(stack[i]),"PSF_FERR", -1);
                cpl_propertylist_set_comment(casu_fits_get_ehu(stack[i]),
                     "PSF_FERR", "error in DIMM-computed PSF [arcsec]");
            }
        }
    }
}


/*---------------------------------------------------------------------------*/
/**
    \par Name:
        vimos_sci_phot_stdstar
    \par Purpose:
        Update header of stack products to give them a photometric calibration
        based on the values in the headers of a master standard star source
        catalogue.
    \par Description:
        Update header of stack products to give them a photometric calibration
        based on the values in the headers of a master standard star source
        catalogue.
    \par Language:
        C
    \param stdstar
        The frame for the standard star catalogue
    \param stack
        The list of images for the stack
    \param stackcat
        The list of tables for the stack source catalogues
    \param airmasskey
        The keyword for the airmass
    \return
        Nothing
    \author
        Jim Lewis, CASU
 */
/*---------------------------------------------------------------------------*/

static int vimos_sci_phot_stdstar(cpl_frame *stdstar, casu_fits **stack,
                                   casu_tfits **stackcat,
                                   const char *airmasskey) {
    int i;
    cpl_propertylist *p;
    float am_std,ext_std,am_obj,ext_obj,dext,magzpt_obj,magzpt_std,magzrr_std;
    float ap3_obj,dap3,dx,skybright;
    float extcoef,vals[VIMOS_NEXTN],vals2[VIMOS_NEXTN];
    double *cd;
    char *photsys;
    cpl_wcs *wcs;

    /* Go through the extension headers of the primary and get the value of
       the extension coefficient. Loop through until you find a non-zero
       one. */

    for (i = 1; i <= VIMOS_NEXTN; i++) {
        p = cpl_propertylist_load(cpl_frame_get_filename(stdstar),i);
        extcoef = (float)cpl_propertylist_get_double(p,"ESO DRS EXTCOEF");
        if (extcoef != 0.0) 
            break;
    }
    
    /* Read the necessary info from the primary headers */

    p = cpl_propertylist_load(cpl_frame_get_filename(stdstar),0);
    am_std = (float)cpl_propertylist_get_double(p,airmasskey);
    ext_std = extcoef*(am_std - 1.0);
    am_obj = (float)cpl_propertylist_get_double(casu_fits_get_phu(stack[0]),
                                                airmasskey);
    ext_obj = extcoef*(am_obj - 1.0);
    cpl_propertylist_delete(p);
    dext = ext_obj - ext_std;

    /* Add some things to the primaries */

    cpl_propertylist_update_string(casu_fits_get_phu(stack[0]),"FLUXCAL",
                                   "ABSOLUTE");
    cpl_propertylist_set_comment(casu_fits_get_phu(stack[0]),"FLUXCAL",
                                 "Certifies the validity of PHOTZP");
    cpl_propertylist_update_string(casu_tfits_get_phu(stackcat[0]),"FLUXCAL",
                                   "ABSOLUTE");
    cpl_propertylist_set_comment(casu_tfits_get_phu(stackcat[0]),"FLUXCAL",
                                 "Certifies the validity of PHOTZP");

    /* Now loop through the image extensions and get the magnitude zeropoints.
       Work out a median value to cover all extensions */

    int nvalid = 0;
    for (i = 1; i <= VIMOS_NEXTN; i++) {
        const char * fname = cpl_frame_get_filename(stdstar);
        p = cpl_propertylist_load(fname,i);
        if (i == 1) 
            photsys = cpl_strdup(cpl_propertylist_get_string(p,"PHOTSYS"));
        magzpt_std = (float)cpl_propertylist_get_double(p,"ESO QC MAGZPT");
        magzrr_std = (float)cpl_propertylist_get_double(p,"ESO QC MAGZERR");
        const float ap3_std = (float)cpl_propertylist_get_double(p,"APCOR3");
        ap3_obj = (float)cpl_propertylist_get_double(casu_tfits_get_ehu(stackcat[i-1]),
            "APCOR3");
        dap3 = ap3_obj - ap3_std;
        magzpt_obj = magzpt_std + dap3 + dext;
        if(magzpt_std != 0)
        {
            vals[nvalid] = magzpt_obj;
            vals2[nvalid] = magzrr_std;
            nvalid++;
        }
        cpl_propertylist_delete(p);
    }

    /* Get a median value */
    if(nvalid > 0)
    {
        magzpt_obj = casu_med(vals,NULL,nvalid);
        magzrr_std = casu_med(vals2,NULL,nvalid);
    }
    else
    {
        cpl_msg_error(cpl_func, "No valid zeropoints in "VIMOS_CAL_STDSTAR);
        freespace(photsys);
        return -1;
    }


    magzpt_obj = casu_med(vals,NULL,VIMOS_NEXTN);
    magzrr_std = casu_med(vals2,NULL,VIMOS_NEXTN);

    /* Loop through again with the median zeropoint value. Start by updating
       the header values of the zeropoint */

    for (i = 1; i <= VIMOS_NEXTN; i++) {
        cpl_propertylist_update_float(casu_fits_get_ehu(stack[i-1]),
            "ESO QC MAGZPT",magzpt_obj);
        cpl_propertylist_set_comment(casu_fits_get_ehu(stack[i-1]),
            "ESO QC MAGZPT","[mag] photometric zeropoint");
        cpl_propertylist_update_float(casu_tfits_get_ehu(stackcat[i-1]),
            "ESO QC MAGZPT",magzpt_obj);
        cpl_propertylist_set_comment(casu_tfits_get_ehu(stackcat[i-1]),
            "ESO QC MAGZPT","[mag] photometric zeropoint");
        cpl_propertylist_update_float(casu_fits_get_ehu(stack[i-1]),
            "ESO QC MAGZERR",magzrr_std);
        cpl_propertylist_set_comment(casu_fits_get_ehu(stack[i-1]),
            "ESO QC MAGZERR","[mag] photometric zeropoint error");
        cpl_propertylist_update_float(casu_tfits_get_ehu(stackcat[i-1]),
            "ESO QC MAGZERR",magzrr_std);
        cpl_propertylist_set_comment(casu_tfits_get_ehu(stackcat[i-1]),
            "ESO QC MAGZERR","[mag] photometric zeropoint error");
        const float eff_expt = (float)cpl_propertylist_get_double(casu_fits_get_phu(stack[0]),"EFF_EXPT");
        const float photzp = magzpt_obj + 2.5 * log10f(eff_expt);
        cpl_propertylist_update_float(casu_fits_get_ehu(stack[i-1]),
            "PHOTZP",photzp);
        cpl_propertylist_set_comment(casu_fits_get_ehu(stack[i-1]),
            "PHOTZP","[mag] photometric zeropoint");        
        cpl_propertylist_update_float(casu_tfits_get_ehu(stackcat[i-1]),
            "PHOTZP",photzp);
        cpl_propertylist_set_comment(casu_tfits_get_ehu(stackcat[i-1]),
            "PHOTZP","[mag] photometric zeropoint");        
        cpl_propertylist_update_float(casu_fits_get_ehu(stack[i-1]),
            "PHOTZPER",magzrr_std);
        cpl_propertylist_set_comment(casu_fits_get_ehu(stack[i-1]),
            "PHOTZPER","[mag] photometric zeropoint error");
        cpl_propertylist_update_float(casu_tfits_get_ehu(stackcat[i-1]),
            "PHOTZPER",magzrr_std);
        cpl_propertylist_set_comment(casu_tfits_get_ehu(stackcat[i-1]),
            "PHOTZPER","[mag] photometric zeropoint error");
        cpl_propertylist_update_string(casu_fits_get_ehu(stack[i-1]),
            "PHOTSYS",photsys);
        cpl_propertylist_set_comment(casu_fits_get_ehu(stack[i-1]),
            "PHOTSYS","Photometric System");        
        cpl_propertylist_update_string(casu_tfits_get_ehu(stackcat[i-1]),
            "PHOTSYS",photsys);
        cpl_propertylist_set_comment(casu_tfits_get_ehu(stackcat[i-1]),
            "PHOTSYS","Photometric System");        
        cpl_propertylist_update_bool(casu_fits_get_ehu(stack[0]),"ZPFUDGED",0);
        cpl_propertylist_update_bool(casu_tfits_get_ehu(stackcat[0]),"ZPFUDGED",
                                     0);
        /* Now the limiting magnitude */
        const float skynoise = (float)cpl_propertylist_get_double(casu_tfits_get_ehu(stackcat[i-1]),
            "ESO DRS SKYNOISE");
        const float rcore = (float)cpl_propertylist_get_double(casu_tfits_get_ehu(stackcat[i-1]),
            "ESO DRS RCORE");
        const float apcor3 = (float)cpl_propertylist_get_double(casu_tfits_get_ehu(stackcat[i-1]),
			"APCOR3");
        const float extinction = 0.0; //Assumed 0 in this case that the external catalog is good enough
        const float lim = casu_calculate_abmag_lim(magzpt_obj, skynoise, rcore, eff_expt, apcor3, extinction);
        cpl_propertylist_update_float(casu_fits_get_ehu(stack[i-1]),
            "ESO QC LIMITING_MAG",lim);
        cpl_propertylist_set_comment(casu_fits_get_ehu(stack[i-1]),
            "ESO QC LIMITING_MAG","[mag] 5 sigma limiting mag.");
        cpl_propertylist_update_float(casu_tfits_get_ehu(stackcat[i-1]),
            "ESO QC LIMITING_MAG",lim);
        cpl_propertylist_set_comment(casu_tfits_get_ehu(stackcat[i-1]),
            "ESO QC LIMITING_MAG","[mag] 5 sigma limiting mag.");
        cpl_propertylist_update_float(casu_fits_get_ehu(stack[i-1]),
            "ABMAGLIM",lim);
        cpl_propertylist_set_comment(casu_fits_get_ehu(stack[i-1]),
            "ABMAGLIM","[mag] 5 sigma limiting mag");
        cpl_propertylist_update_float(casu_tfits_get_ehu(stackcat[i-1]),
            "ABMAGLIM",lim);
        cpl_propertylist_set_comment(casu_tfits_get_ehu(stackcat[i-1]),
            "ABMAGLIM","[mag] 5 sigma limiting mag");

        /* And the sky brightness */

        const float mean_sky = (float)cpl_propertylist_get_double(casu_tfits_get_ehu(stackcat[i-1]),
            "ESO QC MEAN_SKY");
        wcs = cpl_wcs_new_from_propertylist(casu_fits_get_ehu(stack[i-1]));
        cd = cpl_matrix_get_data((cpl_matrix *)cpl_wcs_get_cd(wcs));
        dx = 3600.0*sqrt(fabs(cd[0]*cd[3] - cd[1]*cd[2]));
        cpl_wcs_delete(wcs);
        skybright = magzpt_obj - 2.5*log10(mean_sky/(eff_expt*dx*dx));
        cpl_propertylist_update_float(casu_fits_get_ehu(stack[i-1]),
            "ESO QC SKYBRIGHT",skybright);
        cpl_propertylist_set_comment(casu_fits_get_ehu(stack[i-1]),
            "ESO QC SKYBRIGHT","[mag/arcsec**2] sky brightness");
        cpl_propertylist_update_float(casu_tfits_get_ehu(stackcat[i-1]),
            "ESO QC SKYBRIGHT",skybright);
        cpl_propertylist_set_comment(casu_tfits_get_ehu(stackcat[i-1]),
            "ESO QC SKYBRIGHT","[mag/arcsec**2] sky brightness");

        /* And the saturation level */
        const float psf_fwhm = (float)cpl_propertylist_get_double(casu_tfits_get_ehu(stackcat[i-1]),
                "PSF_FWHM");
        const float pixel_scale = (float)cpl_propertylist_get_double(casu_tfits_get_ehu(stackcat[i-1]),
                "ESO QC WCS_SCALE");
        const float satlev = 65e3f;

        const float satmag = casu_calculate_abmag_sat(magzpt_obj, satlev, mean_sky, psf_fwhm, pixel_scale, eff_expt);
        cpl_propertylist_update_float(casu_fits_get_ehu(stack[i-1]),
            "ABMAGSAT",satmag);
        cpl_propertylist_set_comment(casu_fits_get_ehu(stack[i-1]),
            "ABMAGSAT","[mag] Saturation limit for point sources");
        cpl_propertylist_update_float(casu_tfits_get_ehu(stackcat[i-1]),
            "ABMAGSAT",satmag);
        cpl_propertylist_set_comment(casu_tfits_get_ehu(stackcat[i-1]),
            "ABMAGSAT","[mag] Saturation limit for point sources");
    }
    freespace(photsys);
    return CASU_OK;
}       
    
/**@}*/

/*---------------------------------------------------------------------------*/
/**
    \par Name:
        vimos_sci_init
    \par Purpose:
        Initialise pointers in the global memory structure.
    \par Description:
        Initialise pointers in the global memory structure.
    \par Language:
        C
    \param ps
        The memory structure
    \return
        Nothing
    \author
        Jim Lewis, CASU
 */
/*---------------------------------------------------------------------------*/

static void vimos_sci_init(memstruct *ps) {

    /* Level 0 stuff */

    ps->labels = NULL;
    ps->gr4 = NULL;
    ps->ngr4 = 0;
    ps->master_bias = NULL;
    ps->master_dark = NULL;
    ps->master_twilight_flat = NULL;
    ps->master_conf = NULL;
    ps->master_fringe = NULL;
    ps->master_fringe_var = NULL;
    ps->master_mstd_phot = NULL;
    ps->readgain = NULL;
    ps->mask = NULL;
    ps->phottab = NULL;
    ps->tphottab = NULL;
    ps->science_frames = NULL;
    ps->product_frames_simple = NULL;
    ps->product_frames_simple_var = NULL;
    ps->product_frames_simple_cat = NULL;
    ps->product_frames_simple_mstd_a = NULL;
    ps->product_frame_stack = NULL;
    ps->product_frame_stack_conf = NULL;
    ps->product_frame_stack_var = NULL;
    ps->product_frame_stack_cat = NULL;
    ps->product_frame_stack_mstd_a = NULL;
    ps->product_frame_stack_mstd_p = NULL;
    ps->gaincors = NULL;
    ps->catpath_a = NULL;
    ps->catname_a = NULL;
    ps->catpath_p = NULL;
    ps->catname_p = NULL;
    ps->schlf_n = NULL;
    ps->schlf_s = NULL;
    ps->prod_im = NULL;
    ps->prod_tab = NULL;

    ps->fbias = NULL;
    ps->fdark = NULL;
    ps->fflat = NULL;
    ps->fconf = NULL;
    ps->fconf2 = NULL;
    ps->ffringe = NULL;
    ps->ffringe_var = NULL;
    ps->sci_fits = NULL;
    ps->sci_vars = NULL;
    ps->sci_cats = NULL;
    ps->sci_mstd_a = NULL;
}

/*---------------------------------------------------------------------------*/
/**
    \par Name:
        vimos_sci_tidy
    \par Purpose:
        Free allocated workspace in the global memory structure
    \par Description:
        Free allocated workspace in the global memory structure. The tidy works
        on two levels. Level 1 is for items that are usually cleared up after
        each extension is processed. Level 2 is for cleaning up the whole
        recipe
    \par Language:
        C
    \param ps
        The memory structure
    \param level
        The level of the tidy to be done. 1: Tidy up after finishing an 
        extension, 0: Tidy up after finishing the recipe.
    \return
        Nothing
    \author
        Jim Lewis, CASU
 */
/*---------------------------------------------------------------------------*/

static void vimos_sci_tidy(memstruct *ps, int level) {
    int i;

    /* Level 1 stuff */

    freefits(ps->fbias); 
    freefits(ps->fdark); 
    freefits(ps->fflat); 
    freefits(ps->fconf);
    freefits(ps->fconf2);
    freefits(ps->ffringe);
    freefits(ps->ffringe_var);
    if (ps->sci_fits != NULL) {
        for (i = 0; i < ps->ngr4; i++) {
            freefits(ps->sci_fits[i]);
            freefits(ps->sci_vars[i]);
            freetfits(ps->sci_cats[i]);
            freetfits(ps->sci_mstd_a[i]);
        }
    }

    if (level == 1)
        return;
    
    /* Level 0 stuff */

    freespace(ps->labels);
    vimos_free_groupsof4(&(ps->ngr4),&(ps->gr4));
    freeframe(ps->master_bias);
    freeframe(ps->master_dark);
    freeframe(ps->master_twilight_flat);
    freeframe(ps->master_conf);
    freeframe(ps->master_fringe);
    freeframe(ps->master_fringe_var);
    freeframe(ps->master_mstd_phot);
    freeframe(ps->readgain);
    freemask(ps->mask);
    freeframe(ps->phottab);
    freetable(ps->tphottab);
    freeframeset(ps->science_frames);
    freespace(ps->product_frames_simple);        /* NB: We only have to delete */
    freespace(ps->product_frames_simple_var);    /* the arrays and not the */
    freespace(ps->product_frames_simple_mstd_a); /* frames as these get passed */
    freespace(ps->product_frames_simple_cat);    /* back to esorex */
    freespace(ps->gaincors);
    freespace(ps->catpath_a);
    freespace(ps->catname_a);
    freespace(ps->catpath_p);
    freespace(ps->catname_p);
    freeframe(ps->schlf_n);
    freeframe(ps->schlf_s);
    freespace(ps->sci_fits);
    freespace(ps->sci_vars);
    freespace(ps->sci_cats);
    freespace(ps->sci_mstd_a);
    freeframeset(ps->prod_im);
    freeframeset(ps->prod_tab);
}


/*

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

Revision 1.44  2015/11/19 12:44:44  jim
Fixed to remove offending keywords to keep ASG happy

Revision 1.43  2015/11/12 18:25:09  jim
Fixed provenance issues

Revision 1.42  2015/11/05 12:03:44  jim
Modified to fix issue with MJD-END in headers

Revision 1.41  2015/10/22 12:13:49  jim
A few touchups

Revision 1.40  2015/10/22 11:00:44  jim
fixed typo

Revision 1.39  2015/10/20 11:55:02  jim
Added stk_nfst parameter

Revision 1.38  2015/10/15 11:28:55  jim
modified definition of NCOMBINE and added NSTACK

Revision 1.37  2015/09/18 10:43:00  jim
fixed bug for mstd_p file with no rows

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

Revision 1.35  2015/09/15 10:39:26  jim
Fixed so that variance arrays get same wcs as their science images

Revision 1.34  2015/09/11 09:48:16  jim
changed declaration of some parameters to ranges where appropriate

Revision 1.33  2015/08/12 11:25:01  jim
Fixed bug affecting the loading of master photometric matched standards
catalogues

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

Revision 1.31  2015/08/03 12:49:02  jim
jittered catalogues have a different product label now from single
exposure catalogues. Fixed bug in update_hdr

Revision 1.30  2015/07/29 09:20:29  jim
Fixed some inconsistencies in the file naming in ugly mode and also
plugged a few memory leaks

Revision 1.29  2015/07/02 12:30:15  jim
Fixed so that we get the correct value of FLUXCAL now...

Revision 1.28  2015/06/29 10:36:19  jim
photosys now comes from photcal table header

Revision 1.27  2015/06/25 11:54:46  jim
Fixed call casu_imstack to fix expkey parameter

Revision 1.26  2015/06/08 09:46:31  jim
Added phase3 headers

Revision 1.25  2015/05/21 07:48:27  jim
Fixed bug where error was thrown if the photcal tab had more than one
extension

Revision 1.24  2015/05/13 11:59:51  jim
Fixed so that master fringe images are chopped properly

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

Revision 1.22  2015/04/30 12:24:55  jim
Fixed bug where 'predictable' names for catalogues are not unique. Added
defringing

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

Revision 1.20  2015/03/10 11:58:49  jim
Fixed another uninitialised value issue

Revision 1.19  2015/03/10 11:34:06  jim
Fixed missing initialisation

Revision 1.18  2015/03/10 08:42:46  jim
Plugged a few more memory leaks

Revision 1.17  2015/03/09 19:51:40  jim
Fixed memory leak

Revision 1.16  2015/03/09 15:39:26  jim
Fixed memory leak

Revision 1.15  2015/02/17 11:25:01  jim
Added ability to include a photometric matched standards catalogue as
an input calibration file

Revision 1.14  2015/02/14 12:39:11  jim
Now allow the user to save matched standards catalogues

Revision 1.13  2015/01/29 12:09:20  jim
Modified comments, moved some support routines to vimos_utils.c. Also
modified the grouping algorithm to help deal with missing files

Revision 1.12  2014/12/17 13:01:30  jim
Commented out testfrms on science images for the moment

Revision 1.11  2014/12/12 21:40:16  jim
changed so there is no inverse variance anymore

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

Revision 1.9  2014/05/07 05:23:46  jim
Bug fix so that trimmed master confidence map is used to create individual
confidence map

Revision 1.8  2014/05/06 15:55:28  jim
Fixed problem with cdssearch

Revision 1.7  2014/05/06 07:35:19  jim
modified to use new casu_getstds, to use readgain table and to create
modified confidence maps

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

Revision 1.5  2013/11/05 09:45:59  jim
Added different chop_crud options

Revision 1.4  2013/11/04 05:58:25  jim
Added chop_lowconf option

Revision 1.3  2013/10/25 09:31:22  jim
Fixed docs


*/