vircam_twilight_flat_combine.c

/* $Id: vircam_twilight_flat_combine.c,v 1.40 2007/11/26 09:59:06 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: 2007/11/26 09:59:06 $
 * $Revision: 1.40 $
 * $Name:  $
 */

/* Includes */

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

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

#include "vircam_utils.h"
#include "vircam_mask.h"
#include "vircam_dfs.h"
#include "vircam_mods.h"
#include "vircam_stats.h"
#include "vircam_fits.h"
#include "vircam_pfits.h"
#include "vircam_channel.h"
#include "vircam_paf.h"
#include "vircam_wcsutils.h"

/* Define values for bit mask that flags dummy results */

#define MEANTWI     1
#define CONFMAP     2
#define RATIMG      4
#define STATS_TAB   8

/* Function prototypes */

static int vircam_twilight_flat_combine_create(cpl_plugin *) ;
static int vircam_twilight_flat_combine_exec(cpl_plugin *) ;
static int vircam_twilight_flat_combine_destroy(cpl_plugin *) ;
static int vircam_twilight_flat_combine(cpl_parameterlist *, cpl_frameset *) ;
static int vircam_twilight_flat_combine_save(cpl_frameset *framelist, 
                                             cpl_parameterlist *parlist);
static void vircam_twilight_flat_combine_dummy_products(void);
static void vircam_twilight_flat_combine_normal(int jext);
static int vircam_twilight_flat_combine_lastbit(int jext, 
                                                cpl_frameset *framelist,
                                                cpl_parameterlist *parlist);
static void vircam_twilight_flat_combine_init(void);
static void vircam_twilight_flat_combine_tidy(int level);

/* Static global variables */

static struct {

    /* Input */

    float       lthr;
    float       hthr;
    int         combtype;
    int         scaletype;
    int         xrej;
    float       thresh;
    int         ncells;
    int         extenum;

    /* Output */

    float       flatrms;
    float       flatratio_med;
    float       flatratio_rms;

} vircam_twilight_flat_combine_config;


static struct {
    vir_fits         **good;
    int              ngood;
    int              *labels;
    cpl_frameset     *twilightlist;
    cpl_frame        *master_dark;
    cpl_frame        *master_twilight_flat;
    vir_mask         *master_mask;
    cpl_frame        *chantab;

    cpl_image        *outimage;
    cpl_image        *outconf;
    vir_fits         **twilights;
    int              ntwilights;
    cpl_propertylist *drs;
    cpl_propertylist *drs2;
    unsigned char    *rejmask;
    unsigned char    *rejplus;
    vir_fits         *mfimage;
    cpl_image        *ratioimg;
    cpl_table        *ratioimstats;
    vir_tfits        *ctable;
    vir_fits         *mdark;
    cpl_propertylist *phupaf;
} ps;

static int isfirst;
static cpl_frame *product_frame_mean_twi = NULL;
static cpl_frame *product_frame_conf = NULL;
static cpl_frame *product_frame_ratioimg = NULL;
static cpl_frame *product_frame_ratioimg_stats = NULL;
static int we_expect;
static int we_get;

static char vircam_twilight_flat_combine_description[] =
"vircam_twilight_flat_combine -- VIRCAM twilight flat combine recipe.\n\n"
"Combine a list of twilight flat frames into a mean frame. Optionally\n"
"compare the output frame to a master twilight flat frame\n\n"
"The program accepts the following files in the SOF:\n\n"
"    Tag                   Description\n"
"    -----------------------------------------------------------------------\n"
"    %-21s A list of raw twilight flat images\n"
"    %-21s A master dark frame\n"
"    %-21s Optional reference twilight flat frame\n"
"    %-21s Optional channel table or\n"
"    %-21s Optional initial channel table\n"
"    %-21s Optional master bad pixel map or\n"
"    %-21s Optional master confidence map\n"
"If no master twilight flat is made available, then no comparison will be\n"
"done. This means there will be no output ratio image. If a master twilight\n"
"is available, but no channel table is, then a ratio image will be formed\n"
"but no stats will be written."
"\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_twilight_flat_combine_description,
                   VIRCAM_TWI_RAW,VIRCAM_CAL_DARK,VIRCAM_REF_TWILIGHT_FLAT,
                   VIRCAM_CAL_CHANTAB,VIRCAM_CAL_CHANTAB_INIT,VIRCAM_CAL_BPM,
                   VIRCAM_CAL_CONF);

    cpl_plugin_init(plugin,
                    CPL_PLUGIN_API,
                    VIRCAM_BINARY_VERSION,
                    CPL_PLUGIN_TYPE_RECIPE,
                    "vircam_twilight_flat_combine",
                    "VIRCAM twilight combination recipe",
                    alldesc,
                    "Jim Lewis",
                    "jrl@ast.cam.ac.uk",
                    vircam_get_license(),
                    vircam_twilight_flat_combine_create,
                    vircam_twilight_flat_combine_exec,
                    vircam_twilight_flat_combine_destroy);

    cpl_pluginlist_append(list,plugin);

    return(0);
}

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

static int vircam_twilight_flat_combine_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();

    /* Lower threshold for rejecting underexposed images */

    p = cpl_parameter_new_value("vircam.vircam_twilight_flat_combine.lthr",
                                CPL_TYPE_DOUBLE,
                                "Low rejection threshold for underexpsed images",
                                "vircam.vircam_twilight_flat_combine",
                                0.0);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"lthr");
    cpl_parameterlist_append(recipe->parameters,p);

    /* Upper threshold for rejecting overexposed images */

    p = cpl_parameter_new_value("vircam.vircam_twilight_flat_combine.hthr",
                                CPL_TYPE_DOUBLE,
                                "High rejection threshold for overexposed images",
                                "vircam.vircam_twilight_flat_combine",
                                65535.0);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"hthr");
    cpl_parameterlist_append(recipe->parameters,p);

    /* Fill in the parameters. First the combination type */

    p = cpl_parameter_new_range("vircam.vircam_twilight_flat_combine.combtype",
                                CPL_TYPE_INT,
                                "1 == Median,\n 2 == Mean",
                                "vircam.vircam_twilight_flat_combine",
                                2,1,2);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"combtype");
    cpl_parameterlist_append(recipe->parameters,p);

    /* The requested scaling */

    p = cpl_parameter_new_range("vircam.vircam_twilight_flat_combine.scaletype",
                                CPL_TYPE_INT,
                                "0 == none,\n 1 == additive offset,\n 2 == multiplicative offset,\n 3 == exposure time scaling + additive offset",
                                "vircam.vircam_twilight_flat_combine",
                                2,0,3);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"scaletype");
    cpl_parameterlist_append(recipe->parameters,p);
    
    /* Extra rejection cycle */

    p = cpl_parameter_new_value("vircam.vircam_twilight_flat_combine.xrej",
                                CPL_TYPE_BOOL,
                                "True if using extra rejection cycle",
                                "vircam.vircam_twilight_flat_combine",
                                TRUE);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"xrej");
    cpl_parameterlist_append(recipe->parameters,p);

    /* Rejection threshold */

    p = cpl_parameter_new_value("vircam.vircam_twilight_flat_combine.thresh",
                                CPL_TYPE_DOUBLE,
                                "Rejection threshold in sigma above background",
                                "vircam.vircam_twilight_flat_combine",5.0);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"thresh");
    cpl_parameterlist_append(recipe->parameters,p);

    /* How many cells to divide each data channel */

    p = cpl_parameter_new_enum("vircam.vircam_twilight_flat_combine.ncells",
                               CPL_TYPE_INT,
                               "Number of cells for data channel stats",
                               "vircam.vircam_twilight_flat_combine",8,7,1,2,4,
                               8,16,32,64);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"ncells");
    cpl_parameterlist_append(recipe->parameters,p);     

    /* Extension number of input frames to use */

    p = cpl_parameter_new_range("vircam.vircam_twilight_flat_combine.extenum",
                                CPL_TYPE_INT,
                                "Extension number to be done, 0 == all",
                                "vircam.vircam_twilight_flat_combine",
                                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_twilight_flat_combine_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_twilight_flat_combine(recipe->parameters,recipe->frames));
}
                                
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

static int vircam_twilight_flat_combine_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_twilight_flat_combine(cpl_parameterlist *parlist, 
                                        cpl_frameset *framelist) {
    const char *fctid="vircam_twilight_flat_combine";
    int nlab,j,jst,jfn,retval,status,live,nx,ny,ndit;
    long i;
    cpl_parameter *p;
    vir_fits *ff;
    cpl_propertylist *pp;

    /* 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_twilight_flat_combine_init();
    we_expect = MEANTWI + CONFMAP;

    /* Get the parameters */

    p = cpl_parameterlist_find(parlist,"vircam.vircam_twilight_flat_combine.lthr");
    vircam_twilight_flat_combine_config.lthr = 
        (float)cpl_parameter_get_double(p);
    p = cpl_parameterlist_find(parlist,
                               "vircam.vircam_twilight_flat_combine.hthr");
    vircam_twilight_flat_combine_config.hthr = 
        (float)cpl_parameter_get_double(p);
    p = cpl_parameterlist_find(parlist,
                               "vircam.vircam_twilight_flat_combine.combtype");
    vircam_twilight_flat_combine_config.combtype = cpl_parameter_get_int(p);
    p = cpl_parameterlist_find(parlist,
                               "vircam.vircam_twilight_flat_combine.scaletype");
    vircam_twilight_flat_combine_config.scaletype = cpl_parameter_get_int(p);
    p = cpl_parameterlist_find(parlist,
                               "vircam.vircam_twilight_flat_combine.xrej");
    vircam_twilight_flat_combine_config.xrej = cpl_parameter_get_bool(p);
    p = cpl_parameterlist_find(parlist,
                               "vircam.vircam_twilight_flat_combine.thresh");
    vircam_twilight_flat_combine_config.thresh = 
        (float)cpl_parameter_get_double(p);
    p = cpl_parameterlist_find(parlist,
                               "vircam.vircam_twilight_flat_combine.ncells");
    vircam_twilight_flat_combine_config.ncells = cpl_parameter_get_int(p);
    p = cpl_parameterlist_find(parlist,
                               "vircam.vircam_twilight_flat_combine.extenum");
    vircam_twilight_flat_combine_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_twilight_flat_combine_tidy(2);
        return(-1);
    }

    /* Get the twilight frames */

    if ((ps.labels = cpl_frameset_labelise(framelist,vircam_compare_tags,
                                           &nlab)) == NULL) {
        cpl_msg_error(fctid,"Cannot labelise the input frames");
        vircam_twilight_flat_combine_tidy(2);
        return(-1);
    }
    if ((ps.twilightlist = vircam_frameset_subgroup(framelist,ps.labels,nlab,
                                                    VIRCAM_TWI_RAW)) == NULL) {
        cpl_msg_error(fctid,"Cannot find twilight frames in input frameset");
        vircam_twilight_flat_combine_tidy(2);
        return(-1);
    }
    ps.ntwilights = cpl_frameset_get_size(ps.twilightlist);

    /* 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_twilight_flat_combine_tidy(2);
        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_REF_TWILIGHT_FLAT)) == NULL)
        cpl_msg_info(fctid,"No master twilight flat found -- no ratio image will be formed");
    else
        we_expect |= RATIMG;
        
    /* Check to see if there is a master bad pixel map. If there isn't one 
       then look for a confidence map */

    ps.master_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) {
        if ((ps.chantab = vircam_frameset_subgroup_1(framelist,ps.labels,nlab,
                                                     VIRCAM_CAL_CHANTAB_INIT)) == NULL) {
            cpl_msg_info(fctid,"No channel table found -- no ratio image stats and no linearisation will be done");
        } else {
            cpl_msg_info(fctid,"Channel table is labelled INIT -- no linearisation will be done");
            if (we_expect & RATIMG)
                we_expect |= STATS_TAB;
        }
    } else if (we_expect & RATIMG) {
        we_expect |= STATS_TAB;
    }

    /* Get the number of DITs */

    pp = cpl_propertylist_load(cpl_frame_get_filename(cpl_frameset_get_frame(ps.twilightlist,0)),0);
    if (vircam_pfits_get_ndit(pp,&ndit) != VIR_OK) {
        cpl_msg_error(fctid,"No value for NDIT available");
        freepropertylist(pp);
        vircam_twilight_flat_combine_tidy(2);
        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_twilight_flat_combine_config.extenum,
                       (const cpl_frame *)cpl_frameset_get_frame(ps.twilightlist,0),
                       &jst,&jfn);
    if (jst == -1 || jfn == -1) {
        cpl_msg_error(fctid,"Unable to continue");
        vircam_twilight_flat_combine_tidy(2);
        return(-1);
    }

    /* Get some space for the good frames */

    ps.good = cpl_malloc(ps.ntwilights*sizeof(vir_fits *));

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

    for (j = jst; j <= jfn; j++) {
        status = VIR_OK;
        we_get = 0;
        isfirst = (j == jst);

        /* Load the images and the master dark. */

        ps.twilights = vircam_fits_load_list(ps.twilightlist,CPL_TYPE_FLOAT,j);
        if (ps.twilights == NULL) {
            cpl_msg_info(fctid,"Extension %d twilights wouldn't load",j);
            retval = vircam_twilight_flat_combine_lastbit(j,framelist,parlist);
            if (retval != 0)
                return(-1);
            continue;
        }

        /* Are any of these twilight flats any good? */
        
        ps.ngood = 0;
        for (i = 0; i < ps.ntwilights; i++) {
            ff = ps.twilights[i];
            vircam_pfits_get_detlive(vircam_fits_get_ehu(ff),&live);
            if (! live) {
                cpl_msg_info(fctid,"Detector flagged dead %s",
                             vircam_fits_get_fullname(ff));
                vircam_fits_set_error(ff,VIR_FATAL);
            } else {
                ps.good[ps.ngood] = ff;
                ps.ngood += 1;
            }
        }

        /* If there are no good images, then signal that we need to create 
           dummy products and move on */

        if (ps.ngood == 0) {
            cpl_msg_info(fctid,"All images flagged bad for this extension");
            retval = vircam_twilight_flat_combine_lastbit(j,framelist,parlist);
            if (retval != 0)
                return(-1);
            continue;
        }

        /* Sort out the images that are either over or under exposed */

        vircam_overexp(ps.good,&(ps.ngood),
                       vircam_twilight_flat_combine_config.lthr,
                       vircam_twilight_flat_combine_config.hthr,0);

        /* Check to see how many are left. If there aren't any, then
           signal a major error */

        if (ps.ngood == 0) {
            cpl_msg_info(fctid,"All images either under or overexposed");
            retval = vircam_twilight_flat_combine_lastbit(j,framelist,parlist);
            if (retval != 0)
                return(-1);
            continue;
        }

        /* Load up the mask */

        nx = cpl_image_get_size_x(vircam_fits_get_image(ps.good[0]));
        ny = cpl_image_get_size_y(vircam_fits_get_image(ps.good[0]));
        if (vircam_mask_load(ps.master_mask,j,nx,ny) == VIR_FATAL) {
            cpl_msg_info(fctid,"Unable to load mask image %s[%d]",
                         vircam_mask_get_filename(ps.master_mask),j);
            cpl_msg_info(fctid,"Forcing all pixels to be good from now on");
            vircam_mask_force(ps.master_mask,nx,ny);
        }

        /* Right, we want to dark correct, so we need to load the mean 
           dark and make sure it isn't a dummy */

        ps.mdark = vircam_fits_load(ps.master_dark,CPL_TYPE_FLOAT,j);
        if (ps.mdark == NULL) {
            cpl_msg_info(fctid,"Can't load master dark for extension %d",j);
            retval = vircam_twilight_flat_combine_lastbit(j,framelist,parlist);
            if (retval != 0)
                return(-1);
            continue;
        } else if (vircam_is_dummy(vircam_fits_get_ehu(ps.mdark))) {
            cpl_msg_info(fctid,"Can't master dark extension %d is a dummy",j);
            retval = vircam_twilight_flat_combine_lastbit(j,framelist,parlist);
            if (retval != 0)
                return(-1);
            continue;
        }

        /* Loop for each image and dark correct */

        cpl_msg_info(fctid,"Dark correcting extension %d\n",j);
        for (i = 0; i < ps.ngood; i++)
            vircam_darkcor((ps.good)[i],ps.mdark,1.0,&status);

        /* We need to load the channel table (if it exists) for linearisation
           and for the ratio image stats table */

        if (ps.chantab != NULL) {
            ps.ctable = vircam_tfits_load(ps.chantab,j);
            if (ps.ctable == NULL) {
                cpl_msg_info(fctid,"Channel table extension %d won't load",j);
            } else if (vircam_chantab_verify(vircam_tfits_get_table(ps.ctable)) != VIR_OK) {
                cpl_msg_info(fctid,"Channel table extension %d has errors",j);
                freetfits(ps.ctable);
            } else { 
                pp = cpl_propertylist_load(cpl_frame_get_filename(ps.chantab),
                                           j);
                if (vircam_is_dummy(pp)) {
                    cpl_msg_info(fctid,
                                 "Channel table extensions %d is a dummy",j);
                    freetfits(ps.ctable);
                }
                freepropertylist(pp);
            }
        } else 
            ps.ctable = NULL;
                
        /* Loop for each of the input images and linearise it if there
           is a channel table */

        if (ps.ctable != NULL) {
            cpl_msg_info(fctid,"Linearising extension %d\n",j);
            for (i = 0; i < ps.ngood; i++)
                (void)vircam_lincor((ps.good)[i],ps.ctable,1,ndit,&status);
        }

        /* Call the combine module */

        cpl_msg_info(fctid,"Doing combination for extension %d\n",j);
        (void)vircam_imcombine(ps.good,ps.ngood,
                               vircam_twilight_flat_combine_config.combtype,
                               vircam_twilight_flat_combine_config.scaletype,
                               vircam_twilight_flat_combine_config.xrej,
                               vircam_twilight_flat_combine_config.thresh,
                               &(ps.outimage),&(ps.rejmask),&(ps.rejplus),
                               &(ps.drs),&status);

        /* If these correction and combination routines failed at any stage
           then get out of here */

        if (status == VIR_OK) {
            we_get |= MEANTWI;
            vircam_twilight_flat_combine_normal(j);
        } else {
            cpl_msg_info(fctid,"A processing step failed");
        }

        /* Create any dummies and save the products */
        
        retval = vircam_twilight_flat_combine_lastbit(j,framelist,parlist);
        if (retval != 0)
            return(-1);

    }
    vircam_twilight_flat_combine_tidy(2);
    return(0);
}


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

static int vircam_twilight_flat_combine_save(cpl_frameset *framelist, 
                                             cpl_parameterlist *parlist) {
    cpl_propertylist *plist,*elist,*p,*pafprop;
    int status;
    const char *fctid = "vircam_twilight_flat_combine_save";
    const char *outfile = "twilightcomb.fits";
    const char *outdiff = "twilightratio.fits";
    const char *outdimst = "twilightratiotab.fits";
    const char *outconf = "twilightconf.fits";
    const char *outfilepaf = "twilightcomb";
    const char *outdiffpaf = "twilightratio";
    const char *recipeid = "vircam_twilight_flat_combine";

    /* If we need to make a PHU then do that now based on the first frame
       in the input frame list */

    if (isfirst) {

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

        product_frame_mean_twi = cpl_frame_new();
        cpl_frame_set_filename(product_frame_mean_twi,outfile);
        cpl_frame_set_tag(product_frame_mean_twi,VIRCAM_PRO_TWILIGHT_FLAT);
        cpl_frame_set_type(product_frame_mean_twi,CPL_FRAME_TYPE_IMAGE);
        cpl_frame_set_group(product_frame_mean_twi,CPL_FRAME_GROUP_PRODUCT);
        cpl_frame_set_level(product_frame_mean_twi,CPL_FRAME_LEVEL_FINAL);

        /* Set up the PHU header */

        plist = vircam_fits_get_phu(ps.twilights[0]);
        ps.phupaf = vircam_paf_phu_items(plist);
        vircam_dfs_set_product_primary_header(plist,product_frame_mean_twi,
                                              framelist,parlist,
                                              (char *)recipeid,
                                              "PRO-1.15");

        /* 'Save' the PHU image */                       

        if (cpl_image_save(NULL,outfile,CPL_BPP_8_UNSIGNED,plist,
                           CPL_IO_DEFAULT) != CPL_ERROR_NONE) {
            cpl_msg_error(fctid,"Cannot save product PHU");
            cpl_frame_delete(product_frame_mean_twi);
            return(-1);
        }
        cpl_frameset_insert(framelist,product_frame_mean_twi);

        /* Create a new product frame object and define some tags for the
           output confidence map */

        product_frame_conf = cpl_frame_new();
        cpl_frame_set_filename(product_frame_conf,outconf);
        cpl_frame_set_tag(product_frame_conf,VIRCAM_PRO_CONF);
        cpl_frame_set_type(product_frame_conf,CPL_FRAME_TYPE_IMAGE);
        cpl_frame_set_group(product_frame_conf,CPL_FRAME_GROUP_PRODUCT);
        cpl_frame_set_level(product_frame_conf,CPL_FRAME_LEVEL_FINAL);

        /* Set up the PHU header */

        plist = vircam_fits_get_phu(ps.twilights[0]);
        vircam_dfs_set_product_primary_header(plist,product_frame_conf,
                                              framelist,parlist,
                                              (char *)recipeid,
                                              "PRO-1.15");

        /* 'Save' the PHU image */                       

        if (cpl_image_save(NULL,outconf,CPL_BPP_8_UNSIGNED,plist,
                           CPL_IO_DEFAULT) != CPL_ERROR_NONE) {
            cpl_msg_error(fctid,"Cannot save product PHU");
            cpl_frame_delete(product_frame_conf);
            return(-1);
        }
        cpl_frameset_insert(framelist,product_frame_conf);

        /* Create a new product frame object for the difference image */

        if (we_expect & RATIMG) {
            product_frame_ratioimg = cpl_frame_new();
            cpl_frame_set_filename(product_frame_ratioimg,outdiff);
            cpl_frame_set_tag(product_frame_ratioimg,
                              VIRCAM_PRO_RATIOIMG_TWILIGHT_FLAT);
            cpl_frame_set_type(product_frame_ratioimg,CPL_FRAME_TYPE_IMAGE);
            cpl_frame_set_group(product_frame_ratioimg,
                                CPL_FRAME_GROUP_PRODUCT);
            cpl_frame_set_level(product_frame_ratioimg,CPL_FRAME_LEVEL_FINAL);

            /* Set up the PHU header */

            plist = vircam_fits_get_phu(ps.twilights[0]);
            vircam_dfs_set_product_primary_header(plist,product_frame_ratioimg,
                                                  framelist,parlist,
                                                  (char *)recipeid,
                                                  "PRO-1.15");

            /* 'Save' the PHU image */                   

            if (cpl_image_save(NULL,outdiff,CPL_BPP_8_UNSIGNED,plist,
                               CPL_IO_DEFAULT) != CPL_ERROR_NONE) {
                cpl_msg_error(fctid,"Cannot save product PHU");
                cpl_frame_delete(product_frame_ratioimg);
                return(-1);
            }
            cpl_frameset_insert(framelist,product_frame_ratioimg);
        }

        /* Create a new product frame object for the difference image stats 
           table */

        if (we_expect & STATS_TAB) {
            product_frame_ratioimg_stats = cpl_frame_new();
            cpl_frame_set_filename(product_frame_ratioimg_stats,outdimst);
            cpl_frame_set_tag(product_frame_ratioimg_stats,
                              VIRCAM_PRO_RATIOIMG_TWILIGHT_FLAT_STATS);
            cpl_frame_set_type(product_frame_ratioimg_stats,
                               CPL_FRAME_TYPE_TABLE);
            cpl_frame_set_group(product_frame_ratioimg_stats,
                                CPL_FRAME_GROUP_PRODUCT);
            cpl_frame_set_level(product_frame_ratioimg_stats,
                                CPL_FRAME_LEVEL_FINAL);

            /* Set up PHU header */

            plist = vircam_fits_get_phu(ps.twilights[0]);
            vircam_dfs_set_product_primary_header(plist,
                                                  product_frame_ratioimg_stats,
                                                  framelist,parlist,
                                                  (char *)recipeid,
                                                  "PRO-1.15");

            /* Fiddle with the extension header now */

            elist = vircam_fits_get_ehu(ps.twilights[0]);
            p = cpl_propertylist_duplicate(elist);
            vircam_merge_propertylists(p,ps.drs);
            if (! (we_get & STATS_TAB))
                vircam_dummy_property(p);
            vircam_dfs_set_product_exten_header(p,product_frame_ratioimg_stats,
                                                framelist,parlist,
                                                (char *)recipeid,
                                                "PRO-1.15");
            status = VIR_OK;
            vircam_removewcs(p,&status);

            /* And finally the difference image stats table */

            if (cpl_table_save(ps.ratioimstats,plist,p,outdimst,
                               CPL_IO_DEFAULT) != CPL_ERROR_NONE) {
                cpl_msg_error(fctid,"Cannot save product table extension");
                cpl_propertylist_delete(p);
                return(-1);
            }
            cpl_propertylist_delete(p);
            cpl_frameset_insert(framelist,product_frame_ratioimg_stats);
        }
    }

    /* Get the extension property list */

    plist = vircam_fits_get_ehu(ps.twilights[0]);
    vircam_merge_propertylists(plist,ps.drs);

    /* Fiddle with the header now */

    p = cpl_propertylist_duplicate(plist);
    if (! (we_get & MEANTWI)) 
        vircam_dummy_property(p);
    vircam_dfs_set_product_exten_header(p,product_frame_mean_twi,framelist,
                                        parlist,(char *)recipeid,
                                        "PRO-1.15");
                
    /* Now save the mean twilight flat image extension */

    cpl_propertylist_update_float(p,"ESO QC FLATRMS",
                                  vircam_twilight_flat_combine_config.flatrms);
    cpl_propertylist_set_comment(p,"ESO QC FLATRMS","RMS of output flat");
    if (cpl_image_save(ps.outimage,outfile,CPL_BPP_IEEE_FLOAT,p,
                       CPL_IO_EXTEND) != CPL_ERROR_NONE) {
        cpl_msg_error(fctid,"Cannot save product image extension");
        cpl_propertylist_delete(p);
        return(-1);
    }

    /* Write out PAF for mean image */

    pafprop = vircam_paf_req_items(p);
    vircam_merge_propertylists(pafprop,ps.phupaf);
    vircam_paf_append(pafprop,vircam_fits_get_phu(ps.twilights[0]),
                      "ESO INS FILT1 NAME");
    if (vircam_paf_print((char *)outfilepaf,
                         "VIRCAM/vircam_twilight_flat_combine",
                         "QC file",pafprop) != VIR_OK)
        cpl_msg_warning(fctid,"Unable to save PAF for mean twilight");
    cpl_propertylist_delete(pafprop);
    cpl_propertylist_delete(p);

    /* Now save the twilight ratio image extension */

    if (we_expect & RATIMG) {
        p = cpl_propertylist_duplicate(plist);
        if (! (we_get & RATIMG))
            vircam_dummy_property(p);
        cpl_propertylist_update_float(p,"ESO QC FLATRATIO_MED",
                                      vircam_twilight_flat_combine_config.flatratio_med);
        cpl_propertylist_set_comment(p,"ESO QC FLATRATIO_MED",
                                     "Median of ratio map");
        cpl_propertylist_update_float(p,"ESO QC FLATRATIO_RMS",
                                      vircam_twilight_flat_combine_config.flatratio_rms);
        cpl_propertylist_set_comment(p,"ESO QC FLATRATIO_RMS",
                                     "RMS of ratio map");
        vircam_dfs_set_product_exten_header(p,product_frame_ratioimg,
                                            framelist,parlist,
                                            (char *)recipeid,
                                            "PRO-1.15");
        if (cpl_image_save(ps.ratioimg,outdiff,CPL_BPP_IEEE_FLOAT,p,
                           CPL_IO_EXTEND) != CPL_ERROR_NONE) {
            cpl_propertylist_delete(p);
            cpl_msg_error(fctid,"Cannot save product image extension");
            return(-1);
        }

        /* Write out PAF for difference image */

        pafprop = vircam_paf_req_items(p);
        vircam_merge_propertylists(pafprop,ps.phupaf);
        vircam_paf_append(pafprop,vircam_fits_get_phu(ps.twilights[0]),
                          "ESO INS FILT1 NAME");
        if (vircam_paf_print((char *)outdiffpaf,
                             "VIRCAM/vircam_twilight_flat_combine",
                             "QC file",pafprop) != VIR_OK)
            cpl_msg_warning(fctid,"Unable to save PAF for twilight ratio image");
        cpl_propertylist_delete(pafprop);
        cpl_propertylist_delete(p);
    }

    /* Now any further ratio image stats tables */

    if (! isfirst && (we_expect & STATS_TAB)) {
        p = cpl_propertylist_duplicate(plist);
        if (! (we_get & STATS_TAB))
            vircam_dummy_property(p);
        vircam_dfs_set_product_exten_header(p,product_frame_ratioimg_stats,
                                            framelist,parlist,
                                            (char *)recipeid,
                                            "PRO-1.15");
        status = VIR_OK;
        vircam_removewcs(p,&status);
        if (cpl_table_save(ps.ratioimstats,NULL,p,outdimst,CPL_IO_EXTEND)
                           != CPL_ERROR_NONE) {
            cpl_msg_error(fctid,"Cannot save product table extension");
            cpl_propertylist_delete(p);
            return(-1);
        }       
        cpl_propertylist_delete(p);
    }

    /* Fiddle with the header now */

    vircam_merge_propertylists(plist,ps.drs2);
    p = cpl_propertylist_duplicate(plist);
    if (! (we_get & CONFMAP)) 
        vircam_dummy_property(p);

    /* Now save the confidence map */

    vircam_dfs_set_product_exten_header(p,product_frame_conf,framelist,
                                        parlist,(char *)recipeid,
                                        "PRO-1.15");
    if (cpl_image_save(ps.outconf,outconf,CPL_BPP_16_SIGNED,p,
                       CPL_IO_EXTEND) != CPL_ERROR_NONE) {
        cpl_msg_error(fctid,"Cannot save product image extension");
        cpl_propertylist_delete(p);
        return(-1);
    }
    cpl_propertylist_delete(p);

    /* Get out of here */

    return(0);
}

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

static void vircam_twilight_flat_combine_dummy_products(void) {

    /* See if you even need to be here */

    if (we_get == we_expect)
        return;

    /* First an output combined twilight frame */

    if (! (we_get & MEANTWI)) {
        ps.outimage = vircam_dummy_image(ps.twilights[0]);
        vircam_twilight_flat_combine_config.flatrms = 0.0;
    }

    /* Now the confidence map */

    if (! (we_get & CONFMAP)) {
        ps.outconf = vircam_dummy_image(ps.twilights[0]);
        vircam_twilight_flat_combine_config.flatrms = 0.0;
    }

    /* Do a ratio image */

    if ((we_expect & RATIMG) && ! (we_get & RATIMG)) {
        vircam_twilight_flat_combine_config.flatratio_med = 0.0;
        vircam_twilight_flat_combine_config.flatratio_rms = 0.0;
        ps.ratioimg = vircam_dummy_image(ps.twilights[0]);
    }

    /* If a ratio image stats table is required, then do that now */
   
    if ((we_expect & STATS_TAB) && ! (we_get & STATS_TAB)) 
        ps.ratioimstats = vircam_create_diffimg_stats(0);

    return;
}

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

static void vircam_twilight_flat_combine_normal(int jext) {
    int nx,ny,ncells,status;
    long npi;
    unsigned char *bpm;
    float *idata,med,sig,gdiff,grms;
    const char *fctid="vircam_twilight_flat_combine_normal";

    /* Load up a bad pixel mask */

    nx = cpl_image_get_size_x(ps.outimage);
    ny = cpl_image_get_size_y(ps.outimage);
    npi = nx*ny;
    vircam_mask_load(ps.master_mask,jext,nx,ny);
    bpm = vircam_mask_get_data(ps.master_mask);

    /* Create the bad pixel mask */
    
    status = VIR_OK;
    (void)vircam_mkconf(ps.outimage,(char *)"None Available",ps.master_mask,
                        &(ps.outconf),&(ps.drs2),&status);
    if (status == VIR_OK) 
        we_get |= CONFMAP;
    else {
        cpl_msg_info(fctid,"Confidence map creation failed extension %d",jext);
        status = VIR_OK;
    }
        
    /* Work out the RMS of the mean twilight frame */

    idata = cpl_image_get_data(ps.outimage);
    vircam_medsig(idata,bpm,npi,&med,&sig);

    /* Divide through by the median */

    cpl_propertylist_update_float(ps.drs,"ESO DRS MEDFLAT",med);
    cpl_propertylist_set_comment(ps.drs,"ESO DRS MEDFLAT",
                                 "Median value before normalisation");
    cpl_image_divide_scalar(ps.outimage,med);
    vircam_medmad(idata,bpm,npi,&med,&sig);
    sig *= 1.48;
    vircam_twilight_flat_combine_config.flatrms = sig;

    /* Load up the master twilight flat */

    if (ps.master_twilight_flat != NULL) {
        ps.mfimage = vircam_fits_load(ps.master_twilight_flat,CPL_TYPE_FLOAT,jext);
        if (ps.mfimage == NULL) {
            cpl_msg_info(fctid,"Master twilight extension %d won't load",jext);
        } else if (vircam_is_dummy(vircam_fits_get_ehu(ps.mfimage))) {
            cpl_msg_info(fctid,"Master twilight extension %d is a dummy",jext);
            freefits(ps.mfimage);
        }
    } else
        ps.mfimage = NULL;


    /* Create a ratio image. NB: the difference image routine copes if the 
       input mean image or the channel tables are null.  Thus if either or
       both are null because of a failure to load, then the routine will do
       as  much as it can and return, allowing you to fill in the rest with
       dummy products */

    vircam_twilight_flat_combine_config.flatratio_med = 0.0;
    vircam_twilight_flat_combine_config.flatratio_rms = 0.0;
    ncells = vircam_twilight_flat_combine_config.ncells;
    vircam_difference_image(vircam_fits_get_image(ps.mfimage),ps.outimage,bpm,
                            vircam_tfits_get_table(ps.ctable),ncells,2,
                            &gdiff,&grms,&(ps.ratioimg),
                            &(ps.ratioimstats));
    vircam_mask_clear(ps.master_mask);
    vircam_twilight_flat_combine_config.flatratio_med = gdiff;
    vircam_twilight_flat_combine_config.flatratio_rms = grms;
    if (ps.ratioimg != NULL)
        we_get |= RATIMG;
    if (ps.ratioimstats != NULL)
        we_get |= STATS_TAB;
    return;
}

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


static int vircam_twilight_flat_combine_lastbit(int jext, 
                                                cpl_frameset *framelist,
                                                cpl_parameterlist *parlist) {
    int retval;
    const char *fctid="vircam_twilight_flat_combine_lastbit";

    /* Make whatever dummy products you need */

    vircam_twilight_flat_combine_dummy_products();

    /* Save everything */

    cpl_msg_info(fctid,"Saving products for extension %d",jext);
    retval = vircam_twilight_flat_combine_save(framelist,parlist);
    if (retval != 0) {
        vircam_twilight_flat_combine_tidy(2);
        return(-1);
    }

    /* Free some stuff up */

    vircam_twilight_flat_combine_tidy(1);
    return(0);
}

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

static void vircam_twilight_flat_combine_init(void) {
    ps.labels = NULL;
    ps.twilightlist = NULL;
    ps.twilights = NULL;
    ps.good = NULL;
    ps.master_dark = NULL;
    ps.master_twilight_flat = NULL;
    ps.master_mask = NULL;
    ps.chantab = NULL;
    ps.ctable = NULL;
    ps.outimage = NULL;
    ps.outconf = NULL;
    ps.drs = NULL;
    ps.drs2 = NULL;
    ps.rejmask = NULL;
    ps.rejplus = NULL;
    ps.mfimage = NULL;
    ps.ratioimg = NULL;
    ps.ratioimstats = NULL;
    ps.phupaf = NULL;
}

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

static void vircam_twilight_flat_combine_tidy(int level) {
    freeimage(ps.outimage);
    freeimage(ps.outconf);
    freefitslist(ps.twilights,ps.ntwilights);
    freepropertylist(ps.drs);
    freepropertylist(ps.drs2);
    freespace(ps.rejmask);
    freespace(ps.rejplus);
    freefits(ps.mfimage);
    freeimage(ps.ratioimg);
    freetable(ps.ratioimstats);
    freetfits(ps.ctable);
    freefits(ps.mdark);
    if (level == 1)
        return;

    freespace(ps.good);
    freespace(ps.labels);
    freeframeset(ps.twilightlist);
    freeframe(ps.master_dark);
    freeframe(ps.master_twilight_flat);
    freemask(ps.master_mask);
    freeframe(ps.chantab);
    freepropertylist(ps.phupaf);
}

/*

$Log: vircam_twilight_flat_combine.c,v $
Revision 1.40  2007/11/26 09:59:06  jim
Recipe now takes ndit into account when doing linearity correction

Revision 1.39  2007/10/19 09:25:09  jim
Fixed problems with missing includes

Revision 1.38  2007/10/15 12:53:26  jim
Modified for compatibiliity with cpl_4.0

Revision 1.37  2007/07/18 15:35:42  jim
Added better error handling for missing or corrupt mask extensions

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

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

Revision 1.34  2007/05/08 10:42:33  jim
Added DRS MEDFLAT keyword

Revision 1.33  2007/04/30 09:40:17  jim
Added more stuff to paf files

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

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

Revision 1.30  2007/03/02 12:37:16  jim
Removed WCS stuff from table headers

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

Revision 1.28  2007/02/25 06:27:41  jim
plugged a few memory leaks

Revision 1.27  2007/02/15 11:54:09  jim
Modified to make a distinction between initial channel table and one that
has the proper linearity information

Revision 1.26  2007/02/15 06:59:38  jim
Added ability to write QC paf files

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

Revision 1.24  2007/02/06 13:11:12  jim
Fixed entry for PRO dictionary in cpl_dfs_set_product_header

Revision 1.23  2007/02/05 14:14:05  jim
Input master frame is now tagged as REFERENCE. QC removed from stats table
headers

Revision 1.22  2007/01/08 19:09:11  jim
Fixed memory leak

Revision 1.21  2006/12/13 13:26:09  jim
Fxied bad sigma estimate

Revision 1.20  2006/11/27 12:15:08  jim
changed calls to cpl_propertylist_append to cpl_propertylist_update

Revision 1.19  2006/11/10 09:24:49  jim
Fixed bug in save routine where the wrong propertylist was being saved

Revision 1.18  2006/09/29 11:19:31  jim
changed aliases on parameter names

Revision 1.17  2006/09/10 20:47:03  jim
Small documentation fix

Revision 1.16  2006/09/09 16:49:40  jim
Header comment update

Revision 1.15  2006/08/27 20:30:02  jim
Major mods to structure of the main processing routine to deal with missing
and dummy frames. Deals better with lower level failures too

Revision 1.14  2006/07/11 14:55:12  jim
Now checks for zeros in the output flat and replaces them

Revision 1.13  2006/06/20 19:07:01  jim
Corrects for ndit != 1

Revision 1.12  2006/06/15 09:58:58  jim
Minor changes to docs

Revision 1.11  2006/06/09 11:26:25  jim
Small changes to keep lint happy

Revision 1.10  2006/06/06 13:01:40  jim
Fixed so that the QC parameters go into the correct headers

Revision 1.9  2006/05/26 19:34:18  jim
Fixed recipe to normalise the flats

Revision 1.8  2006/05/17 14:43:58  jim
Fixed problem in save routine which messed up the PRO CATG keywords

Revision 1.7  2006/05/16 13:58:47  jim
Fixed memory leaks that occur from not closing images at the end of
the image extension loop

Revision 1.6  2006/05/09 08:54:40  jim
Fixed _save routine so that the confidence map is saved as a signed 16 bit
integer

Revision 1.5  2006/05/08 14:54:03  jim
Fixed little bug where the confidence map file name wasn't being declared to
the header correctly

Revision 1.4  2006/05/08 13:20:23  jim
Fixed bug where the fitslist for the flats got deleted twice

Revision 1.3  2006/05/04 11:53:15  jim
Fixed the way the _save routine works to be more consistent with the
standard CPL way of doing things

Revision 1.2  2006/04/27 14:21:26  jim
Renamed undefined tag

Revision 1.1  2006/04/27 09:44:47  jim
new file

*/

---