vircam_standard_process.c

/* $Id: vircam_standard_process.c,v 1.22 2008/05/06 12:15:21 jim Exp $
 *
 * This file is part of the VIRCAM Pipeline
 * Copyright (C) 2005 Cambridge Astronomy Survey Unit
 *
 * 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: 2008/05/06 12:15:21 $
 * $Revision: 1.22 $
 * $Name:  $
 */

/* Includes */

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

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

#include "vircam_utils.h"
#include "vircam_mask.h"
#include "vircam_pfits.h"
#include "vircam_dfs.h"
#include "vircam_mods.h"
#include "vircam_fits.h"
#include "vircam_tfits.h"
#include "vircam_jmp_utils.h"

/* Function prototypes */

static int vircam_standard_process_create(cpl_plugin *);
static int vircam_standard_process_exec(cpl_plugin *);
static int vircam_standard_process_destroy(cpl_plugin *);
static int vircam_standard_process(cpl_parameterlist *,cpl_frameset *);
static cpl_propertylist *vircam_standard_process_dummyqc(int type);


static char vircam_standard_process_description[] =
"vircam_standard_process -- VIRCAM standard field processing recipe.\n\n"
"Process a complete pawprint for standard fields in VIRCAM data.\n"
"Remove instrumental signature, interleave microstep sequences (if done),\n"
"combine jitters, photometrically and astrometrically calibrate the pawprint\n"
"image. Photometric calibration can be done with 2mass and 1 other source\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 dark frame\n"
"    %-21s A master twilight flat frame\n"
"    %-21s A channel table\n"
"    %-21s A photometric calibration table\n"
"    %-21s A readnoise/gain file\n"
"    %-21s A master confidence map or\n"
"    %-21s A master bad pixel mask\n"
"    %-21s A master 2mass catalogue index\n"
"    %-21s A second reference catalogue index\n"
"All of the above are required\n"
"\n";

/* Function code */

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

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,vircam_standard_process_description,
                   VIRCAM_STD_OBJECT_RAW,VIRCAM_CAL_DARK,
                   VIRCAM_CAL_TWILIGHT_FLAT,VIRCAM_CAL_CHANTAB,
                   VIRCAM_CAL_PHOTTAB,VIRCAM_CAL_READGAINFILE,VIRCAM_CAL_CONF,
                   VIRCAM_CAL_BPM,VIRCAM_CAL_2MASS,VIRCAM_CAL_REFCAT);

    cpl_plugin_init(plugin,
                    CPL_PLUGIN_API,
                    VIRCAM_BINARY_VERSION,
                    CPL_PLUGIN_TYPE_RECIPE,
                    "vircam_standard_process",
                    "VIRCAM standard field processing recipe",
                    alldesc,
                    "Jim Lewis",
                    "jrl@ast.cam.ac.uk",
                    vircam_get_license(),
                    vircam_standard_process_create,
                    vircam_standard_process_exec,
                    vircam_standard_process_destroy);

    cpl_pluginlist_append(list,plugin);

    return(0);
}

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

static int vircam_standard_process_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 the minimum object size */

    p = cpl_parameter_new_value("vircam.vircam_standard_process.ipix",
                                CPL_TYPE_INT,
                                "Minimum pixel area for each detected object",
                                "vircam.vircam_standard_process",5);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"ipix");
    cpl_parameterlist_append(recipe->parameters,p);

    /* Fill in the detection threshold parameter */

    p = cpl_parameter_new_value("vircam.vircam_standard_process.thresh",
                                CPL_TYPE_DOUBLE,
                                "Detection threshold in sigma above sky",
                                "vircam.vircam_standard_process",2.0);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"thr");
    cpl_parameterlist_append(recipe->parameters,p);

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

    p = cpl_parameter_new_value("vircam.vircam_standard_process.icrowd",
                                CPL_TYPE_BOOL,"Use deblending?",
                                "vircam.vircam_standard_process",0);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"icrowd");
    cpl_parameterlist_append(recipe->parameters,p);

    /* Fill in core radius */

    p = cpl_parameter_new_value("vircam.vircam_standard_process.rcore",
                                CPL_TYPE_DOUBLE,"Value of Rcore in pixels",
                                "vircam.vircam_standard_process",4.0);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"rcore");
    cpl_parameterlist_append(recipe->parameters,p);

    /* Fill in background smoothing box size */

    p = cpl_parameter_new_value("vircam.vircam_standard_process.nbsize",
                                CPL_TYPE_INT,"Background smoothing box size",
                                "vircam.vircam_standard_process",64);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"nb");
    cpl_parameterlist_append(recipe->parameters,p);

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

    p = cpl_parameter_new_value("vircam.vircam_standard_process.savecat",
                                CPL_TYPE_BOOL,"Save catalogue?",
                                "vircam.vircam_standard_process",1);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"savecat");
    cpl_parameterlist_append(recipe->parameters,p);


    /* Fill in flag to destripe the images */

    p = cpl_parameter_new_value("vircam.vircam_standard_process.destripe",
                                CPL_TYPE_BOOL,"Destripe images?",
                                "vircam.vircam_standard_process",1);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"destripe");
    cpl_parameterlist_append(recipe->parameters,p);

    /* Fill in flag to correct sky background */

    p = cpl_parameter_new_value("vircam.vircam_standard_process.skycor",
                                CPL_TYPE_BOOL,"Sky correct images?",
                                "vircam.vircam_standard_process",1);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"skycor");
    cpl_parameterlist_append(recipe->parameters,p);

    /* Extension number of input frames to use */

    p = cpl_parameter_new_range("vircam.vircam_standard_process.extenum",
                                CPL_TYPE_INT,
                                "Extension number to be done, 0 == all",
                                "vircam.vircam_standard_process",
                                1,0,16);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"ext");
    cpl_parameterlist_append(recipe->parameters,p);
        
    /* Get out of here */

    return(0);
}
    
    
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

static int vircam_standard_process_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(vircam_standard_process(recipe->parameters,recipe->frames));
}
                                
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

static int vircam_standard_process_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);
}

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

static int vircam_standard_process(cpl_parameterlist *parlist, 
                                   cpl_frameset *framelist) {
    const char *fctid="vircam_standard_process";
    cpl_parameter *p;
    int nlab,jst,jfn,status,j,i,retval,nusteps,isconf,live,ndit;
    float readnoise,gain,gaincor_fac;
    vir_fits *ff;
    cpl_propertylist *ehu,*pp;
    cpl_frame *catindex;

    /* 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\n");
        return(-1);
    }

    /* Initialise some things */

    vircam_jmp_init();
    (void)strncpy(vircam_recipename,fctid,VIRCAM_PATHSZ);
    (void)snprintf(vircam_recipepaf,VIRCAM_PATHSZ,"VIRCAM/%s",fctid);
    recflag = RECSTD;

    /* Get the parameters */

    p = cpl_parameterlist_find(parlist,
                               "vircam.vircam_standard_process.ipix");
    vircam_jmp_config.ipix = cpl_parameter_get_int(p);
    p = cpl_parameterlist_find(parlist,
                               "vircam.vircam_standard_process.thresh");
    vircam_jmp_config.threshold = (float)cpl_parameter_get_double(p);
    p = cpl_parameterlist_find(parlist,
                               "vircam.vircam_standard_process.icrowd");
    vircam_jmp_config.icrowd = cpl_parameter_get_bool(p);
    p = cpl_parameterlist_find(parlist,
                               "vircam.vircam_standard_process.rcore");
    vircam_jmp_config.rcore = (float)cpl_parameter_get_double(p);
    p = cpl_parameterlist_find(parlist,
                               "vircam.vircam_standard_process.nbsize");
    vircam_jmp_config.nbsize = cpl_parameter_get_int(p);
    p = cpl_parameterlist_find(parlist,
                               "vircam.vircam_standard_process.savecat");
    vircam_jmp_config.savecat = cpl_parameter_get_bool(p);
    p = cpl_parameterlist_find(parlist,
                               "vircam.vircam_standard_process.destripe");
    vircam_jmp_config.destripe = cpl_parameter_get_bool(p);
    p = cpl_parameterlist_find(parlist,
                               "vircam.vircam_standard_process.skycor");
    vircam_jmp_config.skycor = cpl_parameter_get_bool(p);
    p = cpl_parameterlist_find(parlist,
                               "vircam.vircam_standard_process.extenum");
    vircam_jmp_config.extenum = cpl_parameter_get_int(p);

    /* Sort out raw from calib frames */

    if (vircam_dfs_set_groups(framelist) != VIR_OK) {
        cpl_msg_error(fctid,"Cannot identify RAW and CALIB frames");
        vircam_jmp_tidy(0);
        return(-1);
    }

    /* Label the input frames */

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

    /* Get the input science frames */

    if ((ps.science_frames = 
         vircam_frameset_subgroup(framelist,ps.labels,nlab,
                                  VIRCAM_STD_OBJECT_RAW)) == NULL) {
        cpl_msg_error(fctid,"No science images to process!");
        vircam_jmp_tidy(0);
        return(-1);
    }

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

    if ((ps.master_dark = 
         vircam_frameset_subgroup_1(framelist,ps.labels,nlab,
                                    VIRCAM_CAL_DARK)) == NULL) {
        cpl_msg_error(fctid,"No master dark found");
        vircam_jmp_tidy(0);
        return(-1);
    }
        
    /* Check to see if there is a master twilight flat frame */

    if ((ps.master_twilight_flat = 
         vircam_frameset_subgroup_1(framelist,ps.labels,nlab,
                                    VIRCAM_CAL_TWILIGHT_FLAT)) == NULL) {
        cpl_msg_error(fctid,"No master twilight flat found");
        vircam_jmp_tidy(0);
        return(-1);
    }
        
    /* Get the gain corrections */

    status = VIR_OK;
    if (vircam_gaincor_calc(ps.master_twilight_flat,&i,&(ps.gaincors),
                            &status) != VIR_OK) {
        cpl_msg_error(fctid,"Error calculating gain corrections");
        vircam_jmp_tidy(0);
        return(-1);
    }
        
    /* Check to see if there is a readgain file */

    if ((ps.readgain_file = 
         vircam_frameset_subgroup_1(framelist,ps.labels,nlab,
                                    VIRCAM_CAL_READGAINFILE)) == NULL) {
        cpl_msg_error(fctid,"No master readnoise/gain file found");
        vircam_jmp_tidy(0);
        return(-1);
    }
        
    /* Check to see if there is a master confidence map. If there isn't
       then look for a bad pixel mask that can be converted into a 
       confidence map (in an emergency) */

    isconf = 1;
    if ((ps.master_conf = 
         vircam_frameset_subgroup_1(framelist,ps.labels,nlab,
                                    VIRCAM_CAL_CONF)) == NULL) {
        isconf = 0;
        if ((ps.master_conf = 
             vircam_frameset_subgroup_1(framelist,ps.labels,nlab,
                                        VIRCAM_CAL_BPM)) == NULL) {
            cpl_msg_error(fctid,"No master confidence map found");
            vircam_jmp_tidy(0);
            return(-1);
        }
    }
    ps.mask = vircam_mask_define(framelist,ps.labels,nlab);
        
    /* Check to see if there is a channel table */

    if ((ps.chantab = vircam_frameset_subgroup_1(framelist,ps.labels,nlab,
                                                 VIRCAM_CAL_CHANTAB)) == NULL) {
        cpl_msg_error(fctid,"No channel table found");
        vircam_jmp_tidy(0);
        return(-1);
    }

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

    if ((ps.phottab = vircam_frameset_subgroup_1(framelist,ps.labels,nlab,
                                                 VIRCAM_CAL_PHOTTAB)) == NULL) {
        cpl_msg_error(fctid,"No photometric table found");
        vircam_jmp_tidy(0);
        return(-1);
    }
    if ((ps.tphottab = cpl_table_load(cpl_frame_get_filename(ps.phottab),1,0)) == NULL) {
        cpl_msg_error(fctid,"Unable to load photometric table");
        vircam_jmp_tidy(0);
        return(-1);
    }

    /* Is the 2mass index file specified? */

    if ((catindex = vircam_frameset_subgroup_1(framelist,ps.labels,nlab,
                                               VIRCAM_CAL_2MASS)) == NULL) {
        cpl_msg_info(fctid,"No 2MASS index found -- cannot continue");
        vircam_jmp_tidy(0);
        return(-1);
    }
    
    /* Get catalogue parameters */

    if (vircam_catpars(catindex,&(ps.catpath),&(ps.catname)) == VIR_FATAL) {
        vircam_jmp_tidy(0);
        cpl_frame_delete(catindex);
        return(-1);
    }
    cpl_frame_delete(catindex);

    /* Is there a second photometric catalogue index file specified? */

    if ((catindex = vircam_frameset_subgroup_1(framelist,ps.labels,nlab,
                                               VIRCAM_CAL_2MASS)) != NULL) {
    
        /* Get catalogue parameters */

        if (vircam_catpars(catindex,&(ps.catpath2),&(ps.catname2)) == VIR_FATAL) {
            vircam_jmp_tidy(0);
            cpl_frame_delete(catindex);
            return(-1);
        }
        cpl_frame_delete(catindex);
    }

    /* Get the number of DITs */

    pp = cpl_propertylist_load(cpl_frame_get_filename(cpl_frameset_get_frame(ps.science_frames,0)),0);
    if (vircam_pfits_get_ndit(pp,&ndit) != VIR_OK) {
        cpl_msg_error(fctid,"No value for NDIT available");
        freepropertylist(pp);
        vircam_jmp_tidy(0);
        return(-1);
    }
    cpl_propertylist_delete(pp);

    /* Now, how many image extensions do we want to do? If the extension
       number is zero, then we loop for all possible extensions. If it
       isn't then we just do the extension specified */

    vircam_exten_range(vircam_jmp_config.extenum,
                       (const cpl_frame *)cpl_frameset_get_frame(ps.science_frames,0),
                       &jst,&jfn);
    if (jst == -1 || jfn == -1) {
        cpl_msg_error(fctid,"Unable to continue");
        vircam_jmp_tidy(0);
        return(-1);
    }

    /* Now loop for all the extensions... */

    status = VIR_OK;
    for (j = jst; j <= jfn; j++) {
        isfirst = (j == jst);
        gaincor_fac = (ps.gaincors)[j-1];

        /* Load up the calibration frames into vir_fits/vir_tfits structures 
           It is a fatal error if any one of them can't load properly */
        
        ps.fdark = vircam_fits_load(ps.master_dark,CPL_TYPE_FLOAT,j);
        if (ps.fdark == NULL) {
            cpl_msg_error(fctid,"Error loading master dark %s[%d]\n%s",
                          cpl_frame_get_filename(ps.master_dark),j,
                          cpl_error_get_message());
            vircam_jmp_tidy(0);
            return(-1);
        }
        ps.fflat = vircam_fits_load(ps.master_twilight_flat,CPL_TYPE_FLOAT,j);
        if (ps.fflat == NULL) {
            cpl_msg_error(fctid,"Error loading master flat %s[%d]\n%s",
                          cpl_frame_get_filename(ps.master_twilight_flat),j,
                          cpl_error_get_message());
            vircam_jmp_tidy(0);
            return(-1);
        }
        ps.fconf = vircam_fits_load(ps.master_conf,CPL_TYPE_INT,j);
        if (ps.fconf == NULL) {
            cpl_msg_error(fctid,"Error loading master conf %s[%d]\n%s",
                          cpl_frame_get_filename(ps.master_conf),j,
                          cpl_error_get_message());
            vircam_jmp_tidy(0);
            return(-1);
        }
        if (! isconf) 
            vircam_jmp_bpm2conf();
        if (vircam_mask_load(ps.mask,j,
                             cpl_image_get_size_x(vircam_fits_get_image(ps.fconf)),
                             cpl_image_get_size_y(vircam_fits_get_image(ps.fconf))) != VIR_OK) {
            cpl_msg_error(fctid,"Error loading mask from master conf %s[%d]\n%s",
                          cpl_frame_get_filename(ps.master_conf),j,
                          cpl_error_get_message());
            vircam_jmp_tidy(0);
            return(-1);
        }
        ps.fchantab = vircam_tfits_load(ps.chantab,j);
        if (ps.fchantab == NULL) {
            cpl_msg_error(fctid,"Error loading channel table %s[%d]\n%s",
                          cpl_frame_get_filename(ps.chantab),j,
                          cpl_error_get_message());
            vircam_jmp_tidy(0);
            return(-1);
        }

        /* Load up the vir_fits structures for the science images */

        ps.nscience = cpl_frameset_get_size(ps.science_frames);
        ps.sci_fits = vircam_fits_load_list(ps.science_frames,CPL_TYPE_FLOAT,j);
        if (ps.sci_fits == NULL) {
            cpl_msg_error(fctid,"Error loading science frames extension %d: %s",
                          j,cpl_error_get_message());
            vircam_jmp_tidy(0);
            return(-1);
        }

        /* Loop through and mark the frames where the header says the detector
           wasn't live */

        for (i = 0; i < ps.nscience; i++) {
            ff = ps.sci_fits[i];
            vircam_pfits_get_detlive(vircam_fits_get_ehu(ff),&live);
            if (! live) 
                vircam_fits_set_error(ff,VIR_FATAL);
        }

        /* Get the readnoise and gain estimate for this extension */

        vircam_jmp_get_readnoise_gain(j,&readnoise,&gain);

        /* Loop for all the science frames and do the 2d corrections */

        cpl_msg_info(fctid,"Doing stage1 corrections on %s",
                     vircam_fits_get_extname(ps.sci_fits[0]));
        for (i = 0; i < ps.nscience; i++) {
            ff = ps.sci_fits[i];
            cpl_propertylist_update_float(vircam_fits_get_ehu(ff),"READNOIS",
                                          readnoise);
            cpl_propertylist_set_comment(vircam_fits_get_ehu(ff),"READNOIS",
                                         "[e-] Readnoise used in processing");
            cpl_propertylist_update_float(vircam_fits_get_ehu(ff),"GAIN",gain);
            cpl_propertylist_set_comment(vircam_fits_get_ehu(ff),"GAIN",
                                         "[e-/adu] Gain used in processing");
            if (vircam_fits_get_status(ff) == VIR_FATAL) {
                cpl_msg_info(fctid,"Detector is flagged dead in %s",
                             vircam_fits_get_fullname(ff));
                continue;
            }
            status = VIR_OK;
            (void)vircam_darkcor(ff,ps.fdark,1.0,&status);
            (void)vircam_lincor(ff,ps.fchantab,1,ndit,&status);
            (void)vircam_flatcor(ff,ps.fflat,&status);
            (void)vircam_gaincor(ff,gaincor_fac,&status);
            if (vircam_jmp_config.destripe) 
                (void)vircam_destripe(ff,ps.mask,&status);
            vircam_fits_set_error(ff,status);
        }

        /* Do a simple sky correction if requested */

        if (vircam_jmp_config.skycor) {
            cpl_msg_info(fctid,"Doing sky correction");
            vircam_jmp_skycor();
        }

        /* Calculate the illumination correction */

        cpl_msg_info(fctid,"Doing illumination correction");
        (void)strcpy(current_cat,ps.catname);
        (void)strcpy(current_catpath,ps.catpath);
        vircam_jmp_illum();     
        
        /* Save the simple images */ 

        cpl_msg_info(fctid,"Saving simple images");
        vircam_jmp_save_simple(framelist,parlist);
        vircam_mask_clear(ps.mask);

        /* Save the illumination correction table */

        cpl_msg_info(fctid,"Saving illumination correction table");
        dummyqc = vircam_standard_process_dummyqc(3);
        vircam_jmp_save_illum(framelist,parlist);
        freepropertylist(dummyqc);

        /* Look at the first frame in the list and see if the number of
           microsteps is greater than 1. If so, then we'll have to go
           through interleaving. */

        retval = vircam_pfits_get_nusteps(vircam_fits_get_phu(ps.sci_fits[0]),
                                          &nusteps);
        if (retval != VIR_OK) {
            cpl_msg_warning(fctid,"Unable to get nusteps from header.\nAssuming no microstepping");
            nusteps = 1;
        }
        
        /* If the number of microsteps is 1 then copy over the good frames
           into a fits list for dithering. */

        ps.ndith = 0;
        if (nusteps == 1) {
            cpl_msg_info(fctid,"No interleaving will be done");
            ps.dith_input = cpl_malloc(ps.nscience*sizeof(vir_fits *));
            ps.dithc_input = cpl_malloc(sizeof(vir_fits *));
            for (i = 0; i < ps.nscience; i++) {
                if (vircam_fits_get_status(ps.sci_fits[i]) == VIR_OK) 
                    ps.dith_input[ps.ndith++] = ps.sci_fits[i];
            }
            ps.dithc_input[0] = ps.fconf;
            ps.ndithc = 1;
            interlv = 0;

        /* If the number of microsteps is more than 1, then we need
           to do interleaving. The interleaving routine define
           ps.dith_input. */

        } else {
            cpl_msg_info(fctid,"Interleaving");
            vircam_jmp_interleave();
            cpl_msg_info(fctid,"Saving superframe images");
            vircam_jmp_save_super(framelist,parlist);
            interlv = 1;
        }

        /* Work out the jitter offsets and the stack the jitter frame */

        cpl_msg_info(fctid,"Working out jitter offsets");
        vircam_jmp_dither_offsets();
        cpl_msg_info(fctid,"Stacking jittered frame");
        vircam_jmp_dither_images();

        /* Do a catalogue generation */

        cpl_msg_info(fctid,"Doing object extraction");
        vircam_jmp_catalogue();

        /* Create a matched standards table */

        cpl_msg_info(fctid,"Matching objects with 2mass standards");
        vircam_jmp_matched_stds();

        /* Do a WCS fit for the dithered image */

        cpl_msg_info(fctid,"Fitting a WCS");
        vircam_jmp_wcsfit();

        /* Finally do the photometric zeropoint fit */

        cpl_msg_info(fctid,"Doing 2mass photometric zeropoint calculation");
        vircam_jmp_photcal();
        if (vircam_fits_get_status(ps.stack_frame) == VIR_OK) {
            ehu = vircam_fits_get_ehu(ps.stack_frame);
            cpl_propertylist_update_float(ehu,"ESO QC ZPT_2MASS",
                                          cpl_propertylist_get_float(ehu,"ESO QC MAGZPT"));
            cpl_propertylist_set_comment(ehu,"ESO QC ZPT_2MASS",
                                         cpl_propertylist_get_comment(ehu,"ESO QC MAGZPT"));
        }

        /* Redo the photometric zeropoint with another catalogue if we
           have one */

        if (ps.catpath2 != NULL) {
            freetable(ps.matchstds);
            (void)strcpy(current_cat,ps.catname2);
            (void)strcpy(current_catpath,ps.catpath2);
            cpl_msg_info(fctid,"Matching objects with %s standards",
                         ps.catname2);
            vircam_jmp_matched_stds();
            cpl_msg_info(fctid,"Doing %s photometric zeropoint calculation",
                         ps.catname2);
            vircam_jmp_photcal();
            if (vircam_fits_get_status(ps.stack_frame) == VIR_OK) {
                ehu = vircam_fits_get_ehu(ps.stack_frame);
                cpl_propertylist_update_float(ehu,"ESO QC ZPT_STDS",
                                              cpl_propertylist_get_float(ehu,"ESO QC MAGZPT"));
                cpl_propertylist_set_comment(ehu,"ESO QC ZPT_STDS",
                                             cpl_propertylist_get_comment(ehu,"ESO QC MAGZPT"));
                cpl_propertylist_update_string(ehu,"ESO QC ZPT_STDS_CAT",
                                               ps.catname2);
                cpl_propertylist_set_comment(ehu,"ESO QC ZPT_STDS_CAT",
                                             "Catalogue used in zeropoint calc");
            }
        }

        /* Save the dithered images */

        cpl_msg_info(fctid,"Saving stacked images");    
        dummyqc = vircam_standard_process_dummyqc(1);
        vircam_jmp_save_stack(framelist,parlist);
        freepropertylist(dummyqc);
        if (vircam_jmp_config.savecat) {
            cpl_msg_info(fctid,"Saving stacked image catalogues");
            dummyqc = vircam_standard_process_dummyqc(2);
            vircam_jmp_save_catalogue(framelist,parlist);
            freepropertylist(dummyqc);
        }

        /* Clean up on aisle 12! */

        vircam_jmp_tidy(1);
    }

    /* Final cleanup */

    vircam_jmp_tidy(0);
    return(0);
}

static cpl_propertylist *vircam_standard_process_dummyqc(int type) {
    cpl_propertylist *p;

    /* Get an empty property list */

    p = cpl_propertylist_new();

    /* Now switch for the various products */

    switch (type) {

    /* Stack images */

    case 1:
        cpl_propertylist_update_double(p,"ESO QC WCS_DCRVAL1",0.0);
        cpl_propertylist_set_comment(p,"ESO QC WCS_DCRVAL1",
                                     "[deg] change in crval1");
        cpl_propertylist_update_double(p,"ESO QC WCS_DCRVAL2",0.0);
        cpl_propertylist_set_comment(p,"ESO QC WCS_DCRVAL2",
                                     "[deg] change in crval2");
        cpl_propertylist_update_double(p,"ESO QC WCS_DTHETA",0.0);
        cpl_propertylist_set_comment(p,"ESO QC WCS_DTHETA",
                                     "[deg] change in rotation");
        cpl_propertylist_update_double(p,"ESO QC WCS_SCALE",0.0);
        cpl_propertylist_set_comment(p,"ESO QC WCS_SCALE",
                                     "[arcsec] mean plate scale");
        cpl_propertylist_update_double(p,"ESO QC WCS_SHEAR",0.0);
        cpl_propertylist_set_comment(p,"ESO QC WCS_SHEAR",
                                     "[deg] abs(xrot) - abs(yrot)");
        cpl_propertylist_update_double(p,"ESO QC WCS_RMS",0.0);
        cpl_propertylist_set_comment(p,"ESO QC WCS_RMS",
                                     "[arcsec] Average error in WCS fit");
        cpl_propertylist_update_float(p,"ESO QC MAGZPT",0.0);
        cpl_propertylist_set_comment(p,"ESO QC MAGZPT",
                                     "[mag] photometric zeropoint");
        cpl_propertylist_update_float(p,"ESO QC MAGZERR",0.0);
        cpl_propertylist_set_comment(p,"ESO QC MAGZERR",
                                     "[mag] photometric zeropoint error");
        cpl_propertylist_update_int(p,"ESO QC MAGNZPT",0);
        cpl_propertylist_set_comment(p,"ESO QC MAGNZPT",
                                     "number of stars in magzpt calc");
        cpl_propertylist_update_float(p,"ESO QC LIMITING_MAG",0.0);
        cpl_propertylist_set_comment(p,"ESO QC LIMITING_MAG",
                                     "[mag] 5 sigma limiting mag");
        cpl_propertylist_update_float(p,"ESO QC ZPT_2MASS",0.0);
        cpl_propertylist_set_comment(p,"ESO QC ZPT_2MASS",
                                     "[mag] photometric zeropoint");
        cpl_propertylist_update_float(p,"ESO QC ZPT_STDS",0.0);
        cpl_propertylist_set_comment(p,"ESO QC ZPT_STDS",
                                     "[mag] photometric zeropoint");
        cpl_propertylist_update_string(p,"ESO QC ZPT_STDS_CAT","");
        cpl_propertylist_set_comment(p,"ESO QC ZPT_STDS_CAT",
                                     "Catalogue used in zeropoint calc");
        break;
        
    /* Catalogues */

    case 2:
        cpl_propertylist_update_float(p,"ESO QC SATURATION",0.0);
        cpl_propertylist_set_comment(p,"ESO QC SATURATION",
                                     "[adu] Saturation level");
        cpl_propertylist_update_float(p,"ESO QC MEAN_SKY",0.0);
        cpl_propertylist_set_comment(p,"ESO QC MEAN_SKY",
                                     "[adu] Median sky brightness");
        cpl_propertylist_update_float(p,"ESO QC SKY_NOISE",0.0);
        cpl_propertylist_set_comment(p,"ESO QC SKY_NOISE",
                                     "[adu] Pixel noise at sky level");
        cpl_propertylist_update_float(p,"ESO QC IMAGE_SIZE",0.0);
        cpl_propertylist_set_comment(p,"ESO QC IMAGE_SIZE",
                                     "[pixels] Average FWHM of stellar objects");
        cpl_propertylist_update_float(p,"ESO QC ELLIPTICITY",0.0);
        cpl_propertylist_set_comment(p,"ESO QC ELLIPTICITY",
                                     "Average stellar ellipticity (1-b/a)");
        cpl_propertylist_update_float(p,"ESO QC APERTURE_CORR",0.0);
        cpl_propertylist_set_comment(p,"ESO QC APERTURE_CORR",
                                     "Stellar ap-corr 1x core flux");
        cpl_propertylist_update_int(p,"ESO QC NOISE_OBJ",0);
        cpl_propertylist_set_comment(p,"ESO QC NOISE_OBJ",
                                     "Number of noise objects");
        break;

    /* Illumination tables */

    case 3:
        cpl_propertylist_update_float(p,"ESO QC ILLUMCOR_RMS",0.0);
        cpl_propertylist_set_comment(p,"ESO QC ILLUMCOR_RMS",
                                     "RMS of illumination correction map");
        break;
    default:
        break;
    }

    /* Get out of here */

    return(p);
}

/*

$Log: vircam_standard_process.c,v $
Revision 1.22  2008/05/06 12:15:21  jim
Changed to use new version of vircam_catpars

Revision 1.21  2007/11/26 09:59:06  jim
Recipe now takes ndit into account when doing linearity correction

Revision 1.20  2007/10/25 18:39:22  jim
Altered to remove some lint messages

Revision 1.19  2007/10/19 06:55:06  jim
Modifications made to use new method for directing the recipes to the
standard catalogues using the sof

Revision 1.18  2007/07/09 13:21:56  jim
Modified to use new version of vircam_exten_range

Revision 1.17  2007/06/13 08:11:27  jim
Modified docs to reflect changes in DFS tags

Revision 1.16  2007/05/08 21:31:16  jim
fixed typo

Revision 1.15  2007/05/08 10:42:44  jim
Added gain correction

Revision 1.14  2007/05/02 12:53:11  jim
typo fixes in docs

Revision 1.13  2007/04/26 13:09:40  jim
fixed typos

Revision 1.12  2007/04/04 16:05:59  jim
Modified to make paf information a bit more correct

Revision 1.11  2007/04/04 10:36:18  jim
Modified to use new dfs tags

Revision 1.10  2007/03/29 12:19:39  jim
Little changes to improve documentation

Revision 1.9  2007/03/14 14:49:13  jim
Fixed problem with missing paf files in jmp recipes if detlive = F. Also
fixed problem where extra dummy products were being created

Revision 1.8  2007/03/06 12:00:48  jim
Fixed stupid typo in header

Revision 1.7  2007/03/02 12:37:51  jim
Fixed small memory leak

Revision 1.6  2007/03/01 12:41:49  jim
Modified slightly after code checking

Revision 1.5  2007/02/07 10:12:40  jim
Removed calls to vircam_ndit_correct as this is now no longer necessary

Revision 1.4  2006/12/19 13:32:03  jim
Images that are flagged as dead detectors now generate an INFO rather than
an ERROR message

Revision 1.3  2006/11/29 12:28:45  jim
Modified so that the correct recipe names would appear in the headers of
data products

Revision 1.2  2006/11/28 20:57:43  jim
Added illumination correction section

Revision 1.1  2006/11/27 12:15:43  jim
Initial entry


*/

---