vircam_linearity_analyse.c

/* $Id: vircam_linearity_analyse.c,v 1.51 2008/01/22 19:47:56 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/01/22 19:47:56 $
 * $Revision: 1.51 $
 * $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_pfits.h"
#include "vircam_dfs.h"
#include "vircam_mods.h"
#include "vircam_channel.h"
#include "vircam_stats.h"
#include "vircam_paf.h"

/* Function prototypes */

static int vircam_linearity_analyse_create(cpl_plugin *) ;
static int vircam_linearity_analyse_exec(cpl_plugin *) ;
static int vircam_linearity_analyse_destroy(cpl_plugin *) ;
static int vircam_linearity_analyse(cpl_parameterlist *, cpl_frameset *) ;
static int vircam_linearity_analyse_lastbit(int jext, cpl_frameset *framelist,
                                            cpl_parameterlist *parlist);
static int vircam_linearity_analyse_save(cpl_frameset *framelist,
                                         cpl_parameterlist *parlist);
static int vircam_linearity_analyse_domedark_groups(void);
static double *vircam_linearity_analyse_genstat(vir_fits *fframe, int *bpm,
                                                parquet *p, int np);
static double *vircam_linearity_tweakfac(double **fdata, double *mjd, int nim,
                                         int nchan, double *facrng, 
                                         double *maxdiff);
static void vircam_mjdsort(double **fdata, double *mjd, int n);
static cpl_table *vircam_linearity_analyse_diagtab_init(int np, int nrows);
static void vircam_linearity_analyse_init(void);
static void vircam_linearity_analyse_tidy(int level);

/* Static global variables */

static struct {

    /* Input */

    int     norder;
    float   lthr;
    float   hthr;
    int     maxbpmfr;
    int     adjust;
    int     diagnostic;
    int     extenum;

    /* Output */

    float   linearity;
    float   linerror;
    float   bad_pixel_stat;
    int     bad_pixel_num;
    float   facrng;
    float   maxdiff;

} vircam_linearity_analyse_config;

typedef struct {
    cpl_frameset   *darks;
    cpl_frameset   *domes;
    int            ndarks;
    int            ndomes;
    float          exptime;
    unsigned char  flag;
    vir_fits       **proc;
} ddgrp;

#define OK_FLAG       0
#define SATURATE_FLAG 1

#define SUBSET 128
#define SUBSET2 (SUBSET/2)

static struct {
    int              *labels;
    cpl_frameset     *domelist;
    int              ndomes;
    cpl_frameset     *darklist;
    int              ndarks;
    cpl_frameset     *domecheck;
    int              ndomecheck;
    cpl_frameset     *darkcheck;
    int              ndarkcheck;
    cpl_frameset     *sorted;
    cpl_frame        *chanfrm;
    vir_tfits        *chantab;
    cpl_table        *lchantab;
    cpl_array        *bpm_array;
    ddgrp            *ddg;
    int              nddg;
    vir_fits         **flatlist;
    int              nflatlist;
    cpl_propertylist *plist;
    cpl_propertylist *elist;
    int              nx;
    int              ny;
    cpl_propertylist *phupaf;
    cpl_table        *diag1;
    cpl_table        *diag2;
} ps;

static int isfirst;
static int dummy;
static cpl_frame *product_frame_chantab = NULL;
static cpl_frame *product_frame_bpm = NULL;
static cpl_frame *product_frame_diag1 = NULL;
static cpl_frame *product_frame_diag2 = NULL;

static char vircam_linearity_analyse_description[] =
"vircam_linearity_analyse -- VIRCAM linearity mapping recipe.\n\n"
"Form master dark images from the input raw frames and use these to\n"
"dark correct a series of dome flat exposures Using the dark\n"
"corrected dome flat series, work out linearity coefficients for\n"
"each data channel. The program expects the following files in the SOF\n"
"    Tag                   Description\n"
"    -----------------------------------------------------------------------\n"
"    %-21s A list of raw dome flat images\n"
"    %-21s A list of raw dark images\n"
"    %-21s The channel table\n"
"    %-21s A list of raw dome flat images at the monitor exposure time"
"    %-21s A list of raw dark images at the monitor exposure time"
"The first three of these are required. The last two are only required if"
"the light source monitoring algorithm is to be used"
"\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_linearity_analyse_description,
                   VIRCAM_LIN_DOME_RAW,VIRCAM_LIN_DARK_RAW,
                   VIRCAM_CAL_CHANTAB_INIT,VIRCAM_LIN_DOME_CHECK,
                   VIRCAM_LIN_DARK_CHECK);

    cpl_plugin_init(plugin,
                    CPL_PLUGIN_API,
                    VIRCAM_BINARY_VERSION,
                    CPL_PLUGIN_TYPE_RECIPE,
                    "vircam_linearity_analyse",
                    "VIRCAM linearity analysis recipe",
                    alldesc,
                    "Jim Lewis",
                    "jrl@ast.cam.ac.uk",
                    vircam_get_license(),
                    vircam_linearity_analyse_create,
                    vircam_linearity_analyse_exec,
                    vircam_linearity_analyse_destroy);

    cpl_pluginlist_append(list,plugin);

    return(0);
}


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

static int vircam_linearity_analyse_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 parameters. First the polynomial order */

    p = cpl_parameter_new_range("vircam.vircam_linearity_analyse.norder",
                                CPL_TYPE_INT,
                                "Order of polynomial fit",
                                "vircam.vircam_linearity_analyse",
                                2,1,6);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"nord");
    cpl_parameterlist_append(recipe->parameters,p);

    /* The lower threshold */

    p = cpl_parameter_new_value("vircam.vircam_linearity_analyse.lthr",
                                CPL_TYPE_DOUBLE,
                                "Lower bad pixel threshold",
                                "vircam.vircam_linearity_analyse",8.0);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"lthr");
    cpl_parameterlist_append(recipe->parameters,p);

    /* The upper threshold */

    p = cpl_parameter_new_value("vircam.vircam_linearity_analyse.hthr",
                                CPL_TYPE_DOUBLE,
                                "Upper bad pixel threshold",
                                "vircam.vircam_linearity_analyse",8.0);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"hthr");
    cpl_parameterlist_append(recipe->parameters,p);

    /* The maximum number of frames to be used in forming the bad pixel mask */

    p = cpl_parameter_new_value("vircam.vircam_linearity_analyse.maxbpmfr",
                                CPL_TYPE_INT,
                                "Maximum # frames used in bpm analysis",
                                "vircam.vircam_linearity_analyse",10);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"maxbpmfr");
    cpl_parameterlist_append(recipe->parameters,p);

    /* The flag to allow statistics of the dome flat frames to be adjusted
       using a set of monitor exposures */

    p = cpl_parameter_new_value("vircam.vircam_linearity_analyse.adjust",
                                CPL_TYPE_BOOL,
                                "Adjust stats with monitor set",
                                "vircam.vircam_linearity_analyse",1);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"adjust");
    cpl_parameterlist_append(recipe->parameters,p);

    /* The flag to allow diagnostic curves to be written out */

    p = cpl_parameter_new_value("vircam.vircam_linearity_analyse.diagnostic",
                                CPL_TYPE_BOOL,
                                "Write out diagnostic tables",
                                "vircam.vircam_linearity_analyse",0);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"diagnostic");
    cpl_parameterlist_append(recipe->parameters,p);

    /* Extension number of input frames to use */

    p = cpl_parameter_new_range("vircam.vircam_linearity_analyse.extenum",
                                CPL_TYPE_INT,
                                "Extension number to be done, 0 == all",
                                "vircam.vircam_linearity_analyse",
                                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_linearity_analyse_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_linearity_analyse(recipe->parameters,recipe->frames));
}

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

static int vircam_linearity_analyse_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_linearity_analyse(cpl_parameterlist *parlist,
                                    cpl_frameset *framelist) {
    const char *fctid="vircam_linearity_analyse";
    char colname[16];
    int i,jst,jfn,j,status,nlab,k,nbad,ntimes,init,ngood,nuse,krem,irem,n;
    int ndarks,ndomes,kk,live,nbmax,ngood_flats,*bpm,nalloc,nfdata,adjust,ndit;
    long np,npi;
    float *idata,low,high,med,sig,mindit,sat,*exps,badfrac,exp;
    unsigned char *rejmask,*rejplus;
    double *dexps,**fdata,*d,*mjds,*mjdcheck,mjd,**cdata,*cf,fac;
    double facrng,maxdiff,*lindata;
    vir_fits **darks,**domes,*outfits,*test,*outdark,*fframe;
    parquet *pp;
    cpl_parameter *p;
    cpl_propertylist *drs,*plist;
    cpl_image *master_img,*im,*outimage;
    cpl_frame *frame;
    
    /* 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 a few things */

    vircam_linearity_analyse_init();

    /* Get the parameters */

    p = cpl_parameterlist_find(parlist,
                               "vircam.vircam_linearity_analyse.norder");
    vircam_linearity_analyse_config.norder = cpl_parameter_get_int(p);
    p = cpl_parameterlist_find(parlist,
                               "vircam.vircam_linearity_analyse.lthr");
    vircam_linearity_analyse_config.lthr = (float)cpl_parameter_get_double(p);
    p = cpl_parameterlist_find(parlist,
                               "vircam.vircam_linearity_analyse.hthr");
    vircam_linearity_analyse_config.hthr = (float)cpl_parameter_get_double(p);
    p = cpl_parameterlist_find(parlist,
                               "vircam.vircam_linearity_analyse.maxbpmfr");
    vircam_linearity_analyse_config.maxbpmfr = cpl_parameter_get_int(p);
    p = cpl_parameterlist_find(parlist,
                               "vircam.vircam_linearity_analyse.adjust");
    vircam_linearity_analyse_config.adjust = cpl_parameter_get_bool(p);
    p = cpl_parameterlist_find(parlist,
                               "vircam.vircam_linearity_analyse.diagnostic");
    vircam_linearity_analyse_config.diagnostic = cpl_parameter_get_bool(p);
    p = cpl_parameterlist_find(parlist,
                               "vircam.vircam_linearity_analyse.extenum");
    vircam_linearity_analyse_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_linearity_analyse_tidy(2);
        return(-1);
    }

    /* Get framelist labels */

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

    /* Get the dome flat frames */

    if ((ps.domelist = vircam_frameset_subgroup(framelist,ps.labels,nlab,
                                                VIRCAM_LIN_DOME_RAW)) == NULL) {
        cpl_msg_error(fctid,"Cannot find dome flat frames in input frameset");
        vircam_linearity_analyse_tidy(2);
        return(-1);
    }
    ps.ndomes = cpl_frameset_get_size(ps.domelist);

    /* Check to make sure that NDIT == 1 */

    plist = cpl_propertylist_load(cpl_frame_get_filename(cpl_frameset_get_frame(ps.domelist,0)),0);
    (void)vircam_pfits_get_ndit(plist,&ndit);
    freepropertylist(plist);
    if (ndit != 1) {
        cpl_msg_error(fctid,"NDIT=%d. Recipe requires that ndit == 1",ndit);
        vircam_linearity_analyse_tidy(2);
        return(-1);
    }

    /* Get the dark frames */

    if ((ps.darklist = vircam_frameset_subgroup(framelist,ps.labels,nlab,
                                                VIRCAM_LIN_DARK_RAW)) == NULL) {
        cpl_msg_error(fctid,"Cannot find dark frames in input frameset");
        vircam_linearity_analyse_tidy(2);
        return(-1);
    }
    ps.ndarks = cpl_frameset_get_size(ps.darklist);

    /* If you are planning to adjust the stats for the dome frames by using
       a monitor exposure set, then read that frameset in now */

    if (vircam_linearity_analyse_config.adjust) {
        if ((ps.domecheck = vircam_frameset_subgroup(framelist,ps.labels,nlab,
                                                     VIRCAM_LIN_DOME_CHECK)) == NULL) {
            cpl_msg_info(fctid,"No monitor frames found in sof. No adjustments made to stats");
            vircam_linearity_analyse_config.adjust = 0;
            ps.ndomecheck = 0;
        } else {
            ps.ndomecheck = cpl_frameset_get_size(ps.domecheck);
            if ((ps.darkcheck = vircam_frameset_subgroup(framelist,ps.labels,
                                                         nlab,VIRCAM_LIN_DARK_CHECK)) == NULL) {
                cpl_msg_info(fctid,"No darks for monitor frames found in sof. No adjustments made to stats");
                vircam_linearity_analyse_config.adjust = 0;
                ps.ndomecheck = 0;
                freeframeset(ps.domecheck);
                ps.ndarkcheck = 0;
            } else {
                ps.ndarkcheck = cpl_frameset_get_size(ps.darkcheck);
            }
        }       
    }
    
    /* Check to see if there is a channel table. If so, then read it */

    if ((ps.chanfrm = vircam_frameset_subgroup_1(framelist,ps.labels,nlab,
                                                 VIRCAM_CAL_CHANTAB_INIT)) == NULL) {
        cpl_msg_error(fctid,"No initial channel table found");
        vircam_linearity_analyse_tidy(2);
        return(-1);
    }

    /* Create a sorted frameset (needed for creating output header with 
       properties copied from a dome rather than from a dark) */

    ps.sorted = cpl_frameset_new();
    for (j = 0; j < ps.ndomes; j++) 
        cpl_frameset_insert(ps.sorted,cpl_frame_duplicate(cpl_frameset_get_frame(ps.domelist,j)));
    for (j = 0; j < ps.ndarks; j++) 
        cpl_frameset_insert(ps.sorted,cpl_frame_duplicate(cpl_frameset_get_frame(ps.darklist,j)));
    cpl_frameset_insert(ps.sorted,cpl_frame_duplicate(ps.chanfrm));

    /* Group the domes and darks by exposure time */

    if (vircam_linearity_analyse_domedark_groups() != 0) {
        vircam_linearity_analyse_tidy(2);
        return(-1);
    }

    /* See if the number of exposure times is too small for the order of the
       polynomial you want to fit. If it is, the adjust the order */

    if (ps.nddg < vircam_linearity_analyse_config.norder+1) {
        cpl_msg_warning(fctid,
                        "Number of exposure times is too small: %d, order: %d\nTaking fit down to order %d",
                        ps.nddg,vircam_linearity_analyse_config.norder,
                        ps.nddg-1);
        vircam_linearity_analyse_config.norder = ps.nddg - 1;
    }

    /* 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_linearity_analyse_config.extenum,
                       (const cpl_frame *)cpl_frameset_get_frame(ps.ddg[0].darks,0),
                       &jst,&jfn);
    if (jst == -1 || jfn == -1) {
        cpl_msg_error(fctid,"Unable to continue");
        vircam_linearity_analyse_tidy(2);
        return(-1);
    }

    /* Now loop for all the extensions and do the BPM analysis... */

    for (j = jst; j <= jfn; j++) {
        cpl_msg_info(fctid,"Beginning BPM work on extension %d",j);
        isfirst = (j == jst);
        dummy = 0;
        vircam_linearity_analyse_config.bad_pixel_stat = 0.0;
        vircam_linearity_analyse_config.bad_pixel_num = 0;
        vircam_linearity_analyse_config.linerror = 0.0;
        vircam_linearity_analyse_config.linearity = 0.0;
        vircam_linearity_analyse_config.facrng = 0.0;
        vircam_linearity_analyse_config.maxdiff = 0.0;

        /* Get some standard info in case we need it for dummy products */

        test = vircam_fits_load(cpl_frameset_get_first(ps.ddg[0].domes),
                                CPL_TYPE_FLOAT,j);
        ps.plist = cpl_propertylist_duplicate(vircam_fits_get_phu(test));
        ps.elist = cpl_propertylist_duplicate(vircam_fits_get_ehu(test));
        ps.nx = cpl_image_get_size_x(vircam_fits_get_image(test));
        ps.ny = cpl_image_get_size_y(vircam_fits_get_image(test));
        vircam_fits_delete(test);

        /* Load up the channel table for this detector and verify it. Get 
           saturation level for this detector from FITS header */

        ps.chantab = vircam_tfits_load(ps.chanfrm,j);
        if (ps.chantab == NULL) {
            cpl_msg_error(fctid,"Channel table extension %d failed to load",j);
            dummy = 1;
            if (vircam_linearity_analyse_lastbit(j,framelist,parlist) != 0)
                return(-1);
            continue;
        } else if (vircam_chantab_verify(vircam_tfits_get_table(ps.chantab)) 
                   != VIR_OK) {
            cpl_msg_error(fctid,"Channel table extension %d has errors",j);
            dummy = 1;
            if (vircam_linearity_analyse_lastbit(j,framelist,parlist) != 0)
                return(-1);
            continue;
        }
        if (vircam_pfits_get_saturation(vircam_tfits_get_ehu(ps.chantab),
                                        &sat) != VIR_OK) {
            cpl_msg_error(fctid,"Channel table extension header %d missing saturation info",j);
            dummy = 1;
            if (vircam_linearity_analyse_lastbit(j,framelist,parlist) != 0)
                return(-1);
            vircam_linearity_analyse_tidy(1);
            continue;
        }

        /* Get the channel structure */

        vircam_chan_fill(vircam_tfits_get_table(ps.chantab),&pp,&np);

        /* If doing diagnostics, then create the tables now */

        if (vircam_linearity_analyse_config.diagnostic) {
            ps.diag1 = vircam_linearity_analyse_diagtab_init(np,ps.ndomes);
            if (vircam_linearity_analyse_config.adjust)
                ps.diag2 = vircam_linearity_analyse_diagtab_init(np,ps.ndomecheck);
        }

        /* Check DETLIVE for this extension */

        if (vircam_pfits_get_detlive((const cpl_propertylist *)ps.elist,&live) 
            != VIR_OK) {
            cpl_msg_error(fctid,"No DET LIVE keyword in this extension");
            dummy = 1;
            if (vircam_linearity_analyse_lastbit(j,framelist,parlist) != 0) 
                return(-1);
            vircam_linearity_analyse_tidy(1);
            continue;
        }
        if (! live) {
            cpl_msg_info(fctid,"Detector flagged dead");
            dummy = 1;
            if (vircam_linearity_analyse_lastbit(j,framelist,parlist) != 0) 
                return(-1);
            vircam_linearity_analyse_tidy(1);
            continue;
        }

        /* Get the value of MINDIT so that you can work out the rough level
           of the flats before the reset frame was subtracted off */

        if (vircam_pfits_get_mindit((const cpl_propertylist *)ps.elist,
                                    &mindit) != VIR_OK) {
            cpl_msg_error(fctid,"No value of MINDIT found in extension %d",j);
            dummy = 1;
            if (vircam_linearity_analyse_lastbit(j,framelist,parlist) != 0)
                return(-1);
            continue;
        }

        /* Get a representative from each exposure group and check if it
           is saturated. If it is, then reject the group from further 
           analysis */

        ngood = 0;
        exps = cpl_malloc(ps.nddg*sizeof(float));
        ngood_flats = 0;
        for (i = 0; i < ps.nddg; i++) {
            test = vircam_fits_load(cpl_frameset_get_first(ps.ddg[i].domes),
                                    CPL_TYPE_FLOAT,j);
            med = cpl_image_get_median((const cpl_image*)vircam_fits_get_image(test));
            med *= (1.0 + mindit/ps.ddg[i].exptime);
            if (med > sat) {
                ps.ddg[i].flag = SATURATE_FLAG;
            } else {
                ngood++;
                exps[ngood-1] = ps.ddg[i].exptime;
                ngood_flats += ps.ddg[i].ndomes;
            }
            vircam_fits_delete(test);
        }
        exps = cpl_realloc(exps,ngood*sizeof(float));
               
        /* Are there enough non-saturated exposures for linearity fit? */

        if (ngood <  vircam_linearity_analyse_config.norder+1) {
            cpl_msg_info(fctid,"Too few unsaturated flats for linearity fit for extension %d",j);
            dummy = 1;
            cpl_free(exps);
            if (vircam_linearity_analyse_lastbit(j,framelist,parlist) != 0) 
                return(-1);
            continue;
        }           

        /* Sort the exposure array */

        vircam_sort(&exps,ngood,1);

        /* Loop for each exposure time. When you have enough flats, then
           you can quit this loop */

        nuse = min(vircam_linearity_analyse_config.maxbpmfr,ngood_flats);
        ps.nflatlist = 0;
        ps.flatlist = cpl_malloc(nuse*sizeof(vir_fits *));
        for (i = ngood-1; i >= 0; i--) {
            for (k = 0; k <= ps.nddg; k++) {
                if (ps.ddg[k].exptime == exps[i]) {
                    krem = k;
                    break;
                }
            }

            /* Load the dark frames from this exposure time */

            darks = vircam_fits_load_list(ps.ddg[krem].darks,CPL_TYPE_FLOAT,j);
            ndarks = ps.ddg[krem].ndarks;
            if (darks == NULL) {
                cpl_msg_error(fctid,"Error loading darks extension %d, exptime %g",
                              j,ps.ddg[krem].exptime);
                continue;
            }

            /* Form a mean dark for this exposure time. If there is only one
               then don't bother, just load up that one frame. */

            if (ndarks == 1) {
                outdark = vircam_fits_duplicate(darks[0]);
            } else {
                status = VIR_OK;
                (void)vircam_imcombine(darks,ndarks,1,1,1,5.0,&outimage,
                                       &rejmask,&rejplus,&drs,&status);
                freespace(rejmask);
                freespace(rejplus);
                freepropertylist(drs);
                if (status != VIR_OK) {
                    cpl_msg_error(fctid,"Dark combine failure extension %d exposure %g",
                                  j,ps.ddg[krem].exptime);
                    freefitslist(darks,ndarks);
                    continue;
                }
                outdark = vircam_fits_wrap(outimage,darks[0],NULL,NULL);
            }
            freefitslist(darks,ndarks);

            /* Load the flats for this group */

            domes = vircam_fits_load_list(ps.ddg[krem].domes,CPL_TYPE_FLOAT,j);
            ndomes = ps.ddg[krem].ndomes;
            if (domes == NULL) {
                cpl_msg_error(fctid,"Error loading domes extension %d, exptime %g",
                              j,ps.ddg[i].exptime);
                freefits(outdark);
                continue;
            }
            
            /* Now loop for each flat in this group or until you have
               filled the flats buffer */

            for (kk = 0; kk < ndomes; kk++) {
                status = VIR_OK;
                vircam_darkcor(domes[kk],outdark,1.0,&status);
                ps.flatlist[ps.nflatlist] = vircam_fits_duplicate(domes[kk]);
                ps.ddg[krem].proc[kk] = ps.flatlist[ps.nflatlist];
                ps.nflatlist++;
                if (ps.nflatlist == nuse)
                    break;
            }

            /* Tidy up a bit */

            freefitslist(domes,ndomes);
            freefits(outdark);

            /* Do we have enough yet? */

            if (ps.nflatlist == nuse) 
                break;
        }
        freespace(exps);
        ps.flatlist = cpl_realloc(ps.flatlist,
                                  (ps.nflatlist)*sizeof(vir_fits *));
        
        /* Generate a bad pixel mask now */

        status = VIR_OK;
        (void)vircam_genbpm(ps.flatlist,ps.nflatlist,
                            vircam_linearity_analyse_config.lthr,
                            vircam_linearity_analyse_config.hthr,
                            &(ps.bpm_array),&nbad,&badfrac,&status);
        bpm = cpl_array_get_data_int(ps.bpm_array);

        /* Store away some useful info */

        vircam_linearity_analyse_config.bad_pixel_num = nbad;
        vircam_linearity_analyse_config.bad_pixel_stat = badfrac;

        /* Right. Free the pointer for flatlist, but don't delete the 
           vir_fits structures in the list because they are copied into
           the ddg structure */

        freespace(ps.flatlist);
        ps.nflatlist = 0;

        /* Get an initial allocation of space to hold the stats */

        nalloc = 16;
        fdata = cpl_malloc(nalloc*sizeof(double *));
        dexps = cpl_malloc(nalloc*sizeof(double));
        mjds = cpl_malloc(nalloc*sizeof(double));

        /* Loop through the ddg structure, missing out any that have
           overexposed flats. Work out the stats for the images that 
           remain. */

        cpl_msg_info(fctid,"Beginning linearity work on extension %d",j);
        nfdata = 0;
        outdark = NULL;
        for (i = 0; i < ps.nddg; i++) {
            if (ps.ddg[i].flag == SATURATE_FLAG)
                continue;
            for (k = 0; k < ps.ddg[i].ndomes; k++) {

                /* If this particular frame wasn't processed, then you need
                   to form a dark for this group. */

                if (ps.ddg[i].proc[k] == NULL) {
                    if (outdark == NULL) {

                        /* Load the dark frames from this exposure time */

                        darks = vircam_fits_load_list(ps.ddg[i].darks,
                                                      CPL_TYPE_FLOAT,j);
                        ndarks = ps.ddg[i].ndarks;
                        if (darks == NULL) {
                            cpl_msg_error(fctid,"Error loading darks extension %d, exptime %g",
                                          j,ps.ddg[i].exptime);
                            continue;
                        }

                        /* Form a mean dark for this exposure time. If there 
                           is only one then don't bother, just load up that 
                           one frame. */

                        if (ps.ddg[i].ndarks == 1) {
                            outdark = vircam_fits_duplicate(darks[0]);
                        } else {
                            status = VIR_OK;
                            (void)vircam_imcombine(darks,ndarks,1,1,1,5.0,
                                                   &outimage,&rejmask,&rejplus,
                                                   &drs,&status);
                            freespace(rejmask);
                            freespace(rejplus);
                            freepropertylist(drs);
                            if (status != VIR_OK) {
                                cpl_msg_error(fctid,
                                              "Dark combine failure extension %d exposure %g",
                                              j,ps.ddg[i].exptime);
                                freefitslist(darks,ndarks);
                                continue;
                            }
                            outdark = vircam_fits_wrap(outimage,darks[0],
                                                       NULL,NULL);
                        }
                        freefitslist(darks,ndarks);
                    }

                    /* Load the flat and dark correct it */

                    frame = cpl_frameset_get_frame(ps.ddg[i].domes,k);
                    fframe = vircam_fits_load(frame,CPL_TYPE_FLOAT,j);
                    vircam_darkcor(fframe,outdark,1.0,&status);

                /* If this frame has already been corrected, then use it */

                } else {
                    fframe = ps.ddg[i].proc[k];
                }

                /* Generate the stats for this frame and store it */

                d = vircam_linearity_analyse_genstat(fframe,bpm,pp,np);
                if (nfdata >= nalloc) {
                    nalloc += 16;
                    fdata = cpl_realloc(fdata,nalloc*sizeof(double *));
                    dexps = cpl_realloc(dexps,nalloc*sizeof(double));
                    mjds = cpl_realloc(mjds,nalloc*sizeof(double));
                }
                (void)vircam_pfits_get_mjd(vircam_fits_get_phu(fframe),&mjd);
                mjds[nfdata] = mjd;
                dexps[nfdata] = (double)(ps.ddg[i].exptime);            
                fdata[nfdata] = d;              

                /* If doing diagnostic curves, then add the relevant
                   information to the first table now */

                if (ps.diag1 != NULL) {
                    cpl_table_set_string(ps.diag1,"filename",nfdata,
                                         vircam_fits_get_filename(fframe));
                    cpl_table_set_double(ps.diag1,"exptime",nfdata,
                                         dexps[nfdata]);
                    cpl_table_set_double(ps.diag1,"mjd",nfdata,mjd);
                    for (n = 1; n <= np; n++) {
                        snprintf(colname,16,"rawflux_%02d",n);
                        cpl_table_set_double(ps.diag1,colname,nfdata,d[n-1]);
                    }
                }
                if (ps.ddg[i].proc[k] != NULL) { 
                    freefits(ps.ddg[i].proc[k]);
                } else {
                    freefits(fframe);
                }
                nfdata++;
            }
            freefits(outdark);
        }
        freefits(outdark);
        if (ps.diag1 != NULL) 
            cpl_table_set_size(ps.diag1,nfdata);

        /* Now, if we are going to tweak the stats using the monitor exposure
           then we should do that now */

        if (vircam_linearity_analyse_config.adjust) {

            /* Get some workspace for the data array */

            cdata = cpl_malloc(ps.ndomecheck*sizeof(double *));

            /* Get the exposure time for the monitor set and make sure
               that the set isn't saturated */

            test = vircam_fits_load(cpl_frameset_get_first(ps.domecheck),
                                    CPL_TYPE_FLOAT,j);
            (void)vircam_pfits_get_exptime(vircam_fits_get_phu(test),&exp);
            med = cpl_image_get_median((const cpl_image*)vircam_fits_get_image(test));
            med *= (1.0 + mindit/exp);
            adjust = 1;
            if (med > sat) {
                cpl_msg_info(fctid,"Monitor exposures saturated. No drift adjustment made");
                adjust = 0;
            }
            vircam_fits_delete(test);

            /* Ok assuming all that's going well, then form a mean dark */

            if (adjust) {

                darks = vircam_fits_load_list(ps.darkcheck,CPL_TYPE_FLOAT,j);
                ndarks = ps.ndarkcheck;
                if (darks == NULL) {
                    cpl_msg_error(fctid,
                                  "Error loading check darks extension %d",j);
                    continue;
                }

                /* Form a mean dark for this exposure time. If there 
                   is only one then don't bother, just load up that 
                   one frame. */

                if (ndarks == 1) {
                    outdark = vircam_fits_duplicate(darks[0]);
                } else {
                    status = VIR_OK;
                    (void)vircam_imcombine(darks,ndarks,1,1,1,5.0,
                                           &outimage,&rejmask,&rejplus,
                                           &drs,&status);
                    freespace(rejmask);
                    freespace(rejplus);
                    freepropertylist(drs);
                    if (status != VIR_OK) {
                        cpl_msg_error(fctid,
                                      "Dark combine failure extension %d exposure %g",
                                      j,ps.ddg[irem].exptime);
                        freefitslist(darks,ndarks);
                        continue;
                    }
                    outdark = vircam_fits_wrap(outimage,darks[0],
                                               NULL,NULL);
                }
                freefitslist(darks,ndarks);             

                /* Now, loop through the monitor domes, dark correct and then
                   do the stats */
                
                mjdcheck = cpl_malloc(ps.ndomecheck*sizeof(double));
                for (i = 0; i < ps.ndomecheck; i++) {
                    frame = cpl_frameset_get_frame(ps.domecheck,i);
                    fframe = vircam_fits_load(frame,CPL_TYPE_FLOAT,j);
                    vircam_darkcor(fframe,outdark,1.0,&status);
                    d = vircam_linearity_analyse_genstat(fframe,bpm,pp,np);
                    cdata[i] = d;
                    vircam_pfits_get_mjd(vircam_fits_get_phu(fframe),&mjd);
                    vircam_pfits_get_exptime(vircam_fits_get_phu(fframe),&exp);
                    mjdcheck[i] = mjd;

                    /* If doing diagnostics then fill in the table now */

                    if (ps.diag2 != NULL) {
                        cpl_table_set_string(ps.diag2,"filename",i,
                                             vircam_fits_get_filename(fframe));
                        cpl_table_set_double(ps.diag2,"exptime",i,(double)exp);
                        cpl_table_set_double(ps.diag2,"mjd",i,mjd);
                        for (n = 1; n <= np; n++) {
                            snprintf(colname,16,"rawflux_%02d",n);
                            cpl_table_set_double(ps.diag2,colname,i,
                                                 d[n-1]);
                            snprintf(colname,16,"linflux_%02d",n);
                            cpl_table_set_double(ps.diag2,colname,i,
                                                 d[n-1]);
                        }
                    }
                    freefits(fframe);
                }
                freefits(outdark);

                /* Generate the correction factors now */

                cf = vircam_linearity_tweakfac(cdata,mjdcheck,ps.ndomecheck,
                                               np,&facrng,&maxdiff);
                vircam_linearity_analyse_config.facrng = 100.0*(float)facrng;
                vircam_linearity_analyse_config.maxdiff = 100.0*(float)maxdiff;
                if (ps.diag2 != NULL) {
                    for (i = 0; i < ps.ndomecheck; i++) 
                        cpl_table_set_double(ps.diag2,"adjust_fac",i,cf[i]);
                }

                /* Ok, now do the correction for each of the linearity
                   sequence frames */

                for (i = 0; i < nfdata; i++) {
                    mjd = mjds[i];
                    krem = -1;
                    for (k = 0; k < ps.ndomecheck; k++) {
                        if (mjd < mjdcheck[k]) {
                            krem = k;
                            break;
                        }
                    }
                    if (krem == -1) {
                        fac = cf[ps.ndomecheck-1];
                    } else if (krem == 0) {
                        fac = cf[0];
                    } else {
                        fac = 0.5*(cf[krem -1]  + cf[krem]);
                    }
                    for (k = 0; k < np; k++)
                        fdata[i][k] /= fac;
                    if (ps.diag1 != NULL) 
                        cpl_table_set_double(ps.diag1,"adjust_fac",i,fac);
                }

                /* Get rid of some stuff now */

                freespace2(cdata,ps.ndomecheck);
                freespace(cf);
                freespace(mjdcheck);
            }
        }

        /* Do some intermediate tidying */
        
        freespace(mjds);
        vircam_chan_free(np,&pp);

        /* Right, there should be no images left in memory now. Do the
           linearity analysis now. */

        (void)vircam_genlincur(fdata,nfdata,dexps,(double)mindit,ps.chantab,
                               vircam_linearity_analyse_config.norder,
                               &(ps.lchantab),&lindata,&status);
        if (ps.diag1 != NULL) {
            for (i = 0; i < nfdata; i++) {
                for (n = 0; n < np; n++) {
                    snprintf(colname,16,"linflux_%02d",n+1);
                    cpl_table_set_double(ps.diag1,colname,i,lindata[i*np+n]);
                }
            }
        }
        freespace2(fdata,nfdata);
        freespace(dexps);
        freespace(lindata);
        if (status != VIR_OK) {
            cpl_msg_error(fctid,"Linearity curve fit failed extension %d",j);
            dummy = 1;
            if (vircam_linearity_analyse_lastbit(j,framelist,parlist) != 0)
                return(-1);
            vircam_linearity_analyse_tidy(1);
            continue;
        }
        
        /* Get some QC1 parameters by finding the average fit error and
           the average percentage non-linearity */

        vircam_linearity_analyse_config.linearity = 
            (float)cpl_table_get_column_mean(ps.lchantab,"lin_10000");
        vircam_linearity_analyse_config.linerror = 
            (float)cpl_table_get_column_mean(ps.lchantab,"lin_10000_err");

        /* Save new linearity info */

        if (vircam_linearity_analyse_lastbit(j,framelist,parlist) != 0)
            return(-1);
    }

    /* Get out of here */

    vircam_linearity_analyse_tidy(2);
    return(0);
}


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

static int vircam_linearity_analyse_save(cpl_frameset *framelist,
                                         cpl_parameterlist *parlist) {
    cpl_propertylist *plist,*elist,*pafprop;
    cpl_image *outimg;
    const char *outtab = "lchantab.fits";
    const char *outbpm = "bpm.fits";
    const char *outtabpaf = "lchantab";
    const char *outbpmpaf = "bpm";
    const char *outdiag1 = "ldiag1.fits";
    const char *outdiag2 = "ldiag2.fits";
    const char *fctid = "vircam_linearity_analyse_save";
    const char *recipeid = "vircam_linearity_analyse";
    int nx,ny,nord,*bpm,i;

    /* Do some stuff for the first extension to set up the frame */

    nord = vircam_linearity_analyse_config.norder;
    if (isfirst) {

        /* Set up the output frame */

        product_frame_chantab = cpl_frame_new();
        cpl_frame_set_filename(product_frame_chantab,outtab);
        cpl_frame_set_tag(product_frame_chantab,VIRCAM_PRO_CHANTAB);
        cpl_frame_set_type(product_frame_chantab,CPL_FRAME_TYPE_TABLE);
        cpl_frame_set_group(product_frame_chantab,CPL_FRAME_GROUP_PRODUCT);
        cpl_frame_set_level(product_frame_chantab,CPL_FRAME_LEVEL_FINAL);

        /* Set up the PRO keywords for primary */
        
        ps.phupaf = vircam_paf_phu_items(ps.plist);
        plist = cpl_propertylist_duplicate(ps.plist);
        vircam_dfs_set_product_primary_header(plist,product_frame_chantab,
                                              ps.sorted,parlist,
                                              (char *)recipeid,
                                              "PRO-1.15");

        /* Now define the table propertylist and give it an extension name and 
           PRO keywords */

        elist = cpl_propertylist_duplicate(ps.elist);
        vircam_dfs_set_product_exten_header(elist,product_frame_chantab,
                                            ps.sorted,parlist,
                                            (char *)recipeid,
                                            "PRO-1.15");
        
        /* Add QC1 info */

        cpl_propertylist_update_float(elist,"ESO QC LINEARITY",
                                      vircam_linearity_analyse_config.linearity);
        cpl_propertylist_set_comment(elist,"ESO QC LINEARITY",
                                     "% non-linearity at 10000 ADU");
        cpl_propertylist_update_float(elist,"ESO QC LINERROR",
                                      vircam_linearity_analyse_config.linerror);
        cpl_propertylist_set_comment(elist,"ESO QC LINERROR",
                                     "% error non-linearity at 10000 ADU");
        cpl_propertylist_update_float(elist,"ESO QC SCREEN_TOTAL",
                                      vircam_linearity_analyse_config.facrng);
        cpl_propertylist_set_comment(elist,"ESO QC SCREEN_TOTAL",
                                     "total % range in screen variation");
        cpl_propertylist_update_float(elist,"ESO QC SCREEN_STEP",
                                      vircam_linearity_analyse_config.maxdiff);
        cpl_propertylist_set_comment(elist,"ESO QC SCREEN_STEP",
                                     "maximum % step in screen variation");

        /* Set up a dummy table if necessary */

        if (dummy == 1) {
            vircam_dummy_property(elist);
            if (ps.lchantab == NULL) 
                ps.lchantab = vircam_chantab_new(nord,vircam_tfits_get_table(ps.chantab));
        }

        /* And finally save the table */
 
        if (cpl_table_save(ps.lchantab,plist,elist,outtab,CPL_IO_DEFAULT)
            != CPL_ERROR_NONE) {
            cpl_msg_error(fctid,"Cannot save product table extension");
            freepropertylist(plist);
            freepropertylist(elist);
            return(-1);
        }
        cpl_frameset_insert(framelist,product_frame_chantab);

        /* Write PAF */

        pafprop = vircam_paf_req_items(elist);
        vircam_merge_propertylists(pafprop,ps.phupaf);
        vircam_paf_append(pafprop,ps.plist,"ESO INS FILT1 NAME");
        if (vircam_paf_print((char *)outtabpaf,"VIRCAM/vircam_linearity_analyse",
                             "QC file",pafprop) != VIR_OK)
            cpl_msg_warning(fctid,"Unable to save PAF for linearity table");
        cpl_propertylist_delete(pafprop);

        /* Quick tidy */

        freepropertylist(plist);
        freepropertylist(elist);

        /* Set up the output bad pixel mask primary */

        product_frame_bpm = cpl_frame_new();
        cpl_frame_set_filename(product_frame_bpm,outbpm);
        cpl_frame_set_tag(product_frame_bpm,VIRCAM_PRO_BPM);
        cpl_frame_set_type(product_frame_bpm,CPL_FRAME_TYPE_IMAGE);
        cpl_frame_set_group(product_frame_bpm,CPL_FRAME_GROUP_PRODUCT);
        cpl_frame_set_level(product_frame_bpm,CPL_FRAME_LEVEL_FINAL);

        /* Set up the PRO keywords for primary header in the bpm */

        plist = ps.plist;
        vircam_dfs_set_product_primary_header(plist,product_frame_bpm,
                                              ps.sorted,parlist,
                                              (char *)recipeid,
                                              "PRO-1.15");

        /* Now save the PHU 'image' */

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

        /* Section for saving diagnostic tables. First the linearity
           sequence diagnostics */

        if (ps.diag1 != NULL) {

            /* Set up the output frame */

            product_frame_diag1 = cpl_frame_new();
            cpl_frame_set_filename(product_frame_diag1,outdiag1);
            cpl_frame_set_tag(product_frame_diag1,VIRCAM_PRO_LIN_DIAG1);
            cpl_frame_set_type(product_frame_diag1,CPL_FRAME_TYPE_TABLE);
            cpl_frame_set_group(product_frame_diag1,CPL_FRAME_GROUP_PRODUCT);
            cpl_frame_set_level(product_frame_diag1,CPL_FRAME_LEVEL_FINAL);

            /* Set up the PRO keywords for primary */

            plist = cpl_propertylist_duplicate(ps.plist);
            vircam_dfs_set_product_primary_header(plist,product_frame_diag1,
                                                  ps.sorted,parlist,
                                                  (char *)recipeid,
                                                  "PRO-1.15");

            /* Now define the table propertylist and give it an extension name and 
               PRO keywords */

            elist = cpl_propertylist_duplicate(ps.elist);
            vircam_dfs_set_product_exten_header(elist,product_frame_diag1,
                                                ps.sorted,parlist,
                                                (char *)recipeid,
                                                    "PRO-1.15");

            /* Set up a dummy property if necessary */

            if (dummy == 1) 
                vircam_dummy_property(elist);

            /* And finally save the table */

            if (cpl_table_save(ps.diag1,plist,elist,outdiag1,CPL_IO_DEFAULT)
                != CPL_ERROR_NONE) {
                cpl_msg_error(fctid,"Cannot save product table extension");
                freepropertylist(plist);
                freepropertylist(elist);
                return(-1);
            }
            cpl_frameset_insert(framelist,product_frame_diag1);
            freepropertylist(plist);
            freepropertylist(elist);
        }

        /* Now the monitor sequence diagnostics */

        if (ps.diag2 != NULL) {

            /* Set up the output frame */

            product_frame_diag2 = cpl_frame_new();
            cpl_frame_set_filename(product_frame_diag2,outdiag2);
            cpl_frame_set_tag(product_frame_diag2,VIRCAM_PRO_LIN_DIAG2);
            cpl_frame_set_type(product_frame_diag2,CPL_FRAME_TYPE_TABLE);
            cpl_frame_set_group(product_frame_diag2,CPL_FRAME_GROUP_PRODUCT);
            cpl_frame_set_level(product_frame_diag2,CPL_FRAME_LEVEL_FINAL);

            /* Set up the PRO keywords for primary */

            plist = cpl_propertylist_duplicate(ps.plist);
            vircam_dfs_set_product_primary_header(plist,product_frame_diag2,
                                                  ps.sorted,parlist,
                                                  (char *)recipeid,
                                                  "PRO-1.15");

            /* Now define the table propertylist and give it an extension name and 
               PRO keywords */

            elist = cpl_propertylist_duplicate(ps.elist);
            vircam_dfs_set_product_exten_header(elist,product_frame_diag2,
                                                ps.sorted,parlist,
                                                (char *)recipeid,
                                                    "PRO-1.15");

            /* Set up a dummy property if necessary */

            if (dummy == 1) 
                vircam_dummy_property(elist);

            /* And finally save the table */

            if (cpl_table_save(ps.diag2,plist,elist,outdiag2,CPL_IO_DEFAULT)
                != CPL_ERROR_NONE) {
                cpl_msg_error(fctid,"Cannot save product table extension");
                freepropertylist(plist);
                freepropertylist(elist);
                return(-1);
            }
            cpl_frameset_insert(framelist,product_frame_diag2);
            freepropertylist(plist);
            freepropertylist(elist);
        }


    /* Section for all other extensions */

    } else {

        /* Do the table extension PRO keywords */

        elist = cpl_propertylist_duplicate(ps.elist);
        vircam_dfs_set_product_exten_header(elist,product_frame_chantab,
                                            ps.sorted,parlist,
                                            (char *)recipeid,
                                            "PRO-1.15");

        /* Add QC1 info */

        cpl_propertylist_update_float(elist,"ESO QC LINEARITY",
                                      vircam_linearity_analyse_config.linearity);
        cpl_propertylist_set_comment(elist,"ESO QC LINEARITY",
                                     "% non-linearity at 10000 ADU");
        cpl_propertylist_update_float(elist,"ESO QC LINERROR",
                                      vircam_linearity_analyse_config.linerror);
        cpl_propertylist_set_comment(elist,"ESO QC LINERROR",
                                     "% error non-linearity at 10000 ADU");
        cpl_propertylist_update_float(elist,"ESO QC SCREEN_TOTAL",
                                      vircam_linearity_analyse_config.facrng);
        cpl_propertylist_set_comment(elist,"ESO QC SCREEN_TOTAL",
                                     "total % range in screen variation");
        cpl_propertylist_update_float(elist,"ESO QC SCREEN_STEP",
                                      vircam_linearity_analyse_config.maxdiff);
        cpl_propertylist_set_comment(elist,"ESO QC SCREEN_STEP",
                                     "maximum % step in screen variation");

        /* Set up a dummy table if necessary */

        if (dummy == 1) {
            vircam_dummy_property(elist);
            if (ps.lchantab == NULL) 
                ps.lchantab = vircam_chantab_new(nord,vircam_tfits_get_table(ps.chantab));
        }

        /* And finally save the table */
 
        if (cpl_table_save(ps.lchantab,NULL,elist,outtab,CPL_IO_EXTEND)
            != CPL_ERROR_NONE) {
            cpl_msg_error(fctid,"Cannot save product table extension");
            freepropertylist(elist);
            return(-1);
        }

        /* Write PAF */

        pafprop = vircam_paf_req_items(elist);
        vircam_merge_propertylists(pafprop,ps.phupaf);
        vircam_paf_append(pafprop,ps.plist,"ESO INS FILT1 NAME");
        if (vircam_paf_print((char *)outtabpaf,"VIRCAM/vircam_linearity_analyse",
                             "QC file",pafprop) != VIR_OK)
            cpl_msg_warning(fctid,"Unable to save PAF for BPM");
        cpl_propertylist_delete(pafprop);

        /* Quick tidy */

        freepropertylist(elist);

        /* Now the diagnostic tables */

        if (ps.diag1 != NULL) {
            elist = cpl_propertylist_duplicate(ps.elist);
            vircam_dfs_set_product_exten_header(elist,product_frame_diag1,
                                                ps.sorted,parlist,
                                                (char *)recipeid,
                                                "PRO-1.15");

            /* Set up a dummy property if necessary */

            if (dummy == 1) 
                vircam_dummy_property(elist);

            /* And finally save the table */

            if (cpl_table_save(ps.diag1,NULL,elist,outdiag1,CPL_IO_EXTEND)
                != CPL_ERROR_NONE) {
                cpl_msg_error(fctid,"Cannot save product table extension");
                freepropertylist(elist);
                return(-1);
            }
            freepropertylist(elist);
        }
        if (ps.diag2 != NULL) {
            elist = cpl_propertylist_duplicate(ps.elist);
            vircam_dfs_set_product_exten_header(elist,product_frame_diag2,
                                                ps.sorted,parlist,
                                                (char *)recipeid,
                                                "PRO-1.15");

            /* Set up a dummy property if necessary */

            if (dummy == 1) 
                vircam_dummy_property(elist);

            /* And finally save the table */

            if (cpl_table_save(ps.diag2,NULL,elist,outdiag2,CPL_IO_EXTEND)
                != CPL_ERROR_NONE) {
                cpl_msg_error(fctid,"Cannot save product table extension");
                freepropertylist(elist);
                return(-1);
            }
            freepropertylist(elist);
        }

    }

    /* Save the bpm extension now */

    plist = ps.elist;
    nx = ps.nx;
    ny = ps.ny;
    if (dummy && ps.bpm_array == NULL) {
        ps.bpm_array = cpl_array_new(nx*ny,CPL_TYPE_INT);
        bpm = cpl_array_get_data_int(ps.bpm_array);
        for (i = 0; i < nx*ny; i++)
            bpm[i] = 0;
    }
    bpm = cpl_array_get_data_int(ps.bpm_array);
    vircam_dfs_set_product_exten_header(plist,product_frame_bpm,
                                        ps.sorted,parlist,
                                        (char *)recipeid,
                                        "PRO-1.15");
    cpl_propertylist_update_float(plist,"ESO QC BAD_PIXEL_STAT",
                                  vircam_linearity_analyse_config.bad_pixel_stat);
    cpl_propertylist_set_comment(plist,"ESO QC BAD_PIXEL_STAT",
                                 "Fraction of pixels that are bad");
    cpl_propertylist_update_int(plist,"ESO QC BAD_PIXEL_NUM",
                                vircam_linearity_analyse_config.bad_pixel_num);
    cpl_propertylist_set_comment(plist,"ESO QC BAD_PIXEL_NUM",
                                 "Number of pixels that are bad");
    if (dummy)
        vircam_dummy_property(plist);
    outimg = cpl_image_wrap_int(nx,ny,bpm);
    if (cpl_image_save(outimg,outbpm,CPL_BPP_8_UNSIGNED,plist,
                       CPL_IO_EXTEND) != CPL_ERROR_NONE) {
        cpl_msg_error(fctid,"Cannot save product image extension");
        return(-1);
    }

    /* Write PAF */

    pafprop = vircam_paf_req_items(plist);
    vircam_merge_propertylists(pafprop,ps.phupaf);
    vircam_paf_append(pafprop,ps.plist,"ESO INS FILT1 NAME");
    if (vircam_paf_print((char *)outbpmpaf,"VIRCAM/vircam_linearity_analyse",
                         "QC file",pafprop) != VIR_OK)
        cpl_msg_warning(fctid,"Unable to save PAF for linearity table");
    cpl_propertylist_delete(pafprop);

    /* Quick tidy */

    cpl_image_unwrap(outimg);

    return(0);
}

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

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

    /* Save the new channel table and bad pixel map */

    cpl_msg_info(fctid,"Saving linearity table and bpm for extension %d",jext);
    retval = vircam_linearity_analyse_save(framelist,parlist);
    if (retval != 0) {
        vircam_linearity_analyse_tidy(2);
        return(-1);
    }

    /* Do some intermediate tidying */

    vircam_linearity_analyse_tidy(1);
    return(0);
}

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

static int vircam_linearity_analyse_domedark_groups(void) {
    int i,j,found;
    float texp;
    cpl_frame *frame;
    cpl_propertylist *plist;
    const char *fctid = "vircam_linearity_analyse_domedark_groups";

    /* Start by getting the memory for the domedark groups */

    ps.ddg = cpl_calloc(ps.ndomes,sizeof(ddgrp));
    ps.nddg = 0;

    /* Loop for each of the dome frames and get its exposure time. If this
       doesn't exist, then signal an error and go on */

    for (i = 0; i < ps.ndomes; i++) {
        frame = cpl_frameset_get_frame(ps.domelist,i);
        plist = cpl_propertylist_load(cpl_frame_get_filename(frame),0);
        if (vircam_pfits_get_exptime(plist,&texp) != VIR_OK) {
            cpl_msg_warning(fctid,"No exposure time found in %s",
                            cpl_frame_get_filename(frame));
            cpl_propertylist_delete(plist);
            continue;
        }
        cpl_propertylist_delete(plist);

        /* Search the domedark groups to see if this exposure time has already 
           been used. If not, then create a new group. If it has then just add 
           this frame to the correct group */

        found = 0;
        for (j = 0; j < ps.nddg; j++) {
            if (ps.ddg[j].exptime == texp) {
                found = 1;
                break;
            }
        }
        if (found) {
            cpl_frameset_insert(ps.ddg[j].domes,cpl_frame_duplicate(frame));
            ps.ddg[j].ndomes += 1;
        } else {
            ps.ddg[ps.nddg].exptime = texp;
            ps.ddg[ps.nddg].darks = cpl_frameset_new();
            ps.ddg[ps.nddg].domes = cpl_frameset_new();
            ps.ddg[ps.nddg].ndarks = 0;
            ps.ddg[ps.nddg].ndomes = 1;
            ps.ddg[ps.nddg].flag = OK_FLAG;
            cpl_frameset_insert(ps.ddg[ps.nddg].domes,
                               cpl_frame_duplicate(frame));
            ps.nddg += 1;
        }
    }
            
    /* Right, now loop through all the darks and get their exposure times */

    for (i = 0; i < ps.ndarks; i++) {
        frame = cpl_frameset_get_frame(ps.darklist,i);
        plist = cpl_propertylist_load(cpl_frame_get_filename(frame),0);
        if (vircam_pfits_get_exptime(plist,&texp) != VIR_OK) {
            cpl_msg_warning(fctid,"No exposure time found in %s",
                            cpl_frame_get_filename(frame));
            cpl_propertylist_delete(plist);
            continue;
        }
        cpl_propertylist_delete(plist);

        /* Search the domedark groups to see if this dark fits into one of
           the defined groups. If not, then ignore it. If it does, then
           add it into the dark frameset */

        found = 0;
        for (j = 0; j < ps.nddg; j++) {
            if (ps.ddg[j].exptime == texp) {
                found = 1;
                break;
            }
        }
        if (found) {
            cpl_frameset_insert(ps.ddg[j].darks,cpl_frame_duplicate(frame));
            ps.ddg[j].ndarks += 1;
        }
    }

    /* Now go through the domedark groups and ditch any that don't have any
       dark frames */

    i = 0;
    while (i < ps.nddg) {
        if (ps.ddg[i].ndarks == 0) {
            cpl_msg_error(fctid,"No dark frames exist for exposure %g\n",
                          ps.ddg[i].exptime);
            freeframeset(ps.ddg[i].darks);
            freeframeset(ps.ddg[i].domes);
            for (j = i+1; j < ps.nddg; j++)
                ps.ddg[j-1] = ps.ddg[j];
            ps.nddg -= 1;
        } else
            i++;
    }

    /* Allocate some space for vir_fits arrays for processed domes */

    for (i = 0; i < ps.nddg; i++) {
        ps.ddg[i].proc = cpl_malloc(ps.ddg[i].ndomes*sizeof(vir_fits *));
        for (j = 0; j < ps.ddg[i].ndomes; j++) 
            ps.ddg[i].proc[j] = NULL;
    }                

    /* Resize the output array and return so long as there is anything
       left. If there isn't then signal a major error */

    if (ps.nddg > 0) {
        ps.ddg = cpl_realloc(ps.ddg,ps.nddg*sizeof(ddgrp));
        return(0);
    } else {
        cpl_msg_error(fctid,"There are no darks defined for input domes");
        return(-1);
    }
}

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

static double *vircam_linearity_analyse_genstat(vir_fits *fframe, int *bpm,
                                                parquet *p, int np) {
    int i,ist,ifn,jst,jfn,n,jind2,iind2,jj,nx,ii;
    parquet *pp;
    double *d;
    float *tmp,*data;

    /* Get the workspace for the output result */

    d = cpl_malloc(np*sizeof(double));

    /* Get the input data array */

    nx = cpl_image_get_size_x(vircam_fits_get_image(fframe));
    data = cpl_image_get_data_float(vircam_fits_get_image(fframe));

    /* Get some workspace for doing the median calculations */

    tmp = cpl_malloc(SUBSET*SUBSET*sizeof(float));

    /* Loop for each channel in the parquet structure */

    for (i = 0; i < np; i++) {
        pp = p + i;

        /* Take the central part of the channel */

        ist = ((pp->delta_i)/2 - SUBSET2);
        ifn = ist + SUBSET - 1;
        jst = ((pp->delta_j)/2 - SUBSET2);
        jfn = jst + SUBSET - 1;

        /* Put the data into the workspace and do a median */

        n = 0;
        for (jj = jst; jj <= jfn; jj++) {
            jind2 = (jj + pp->iymin - 1)*nx;
            for (ii = ist; ii <= ifn; ii++) {
                iind2 = jind2 + ii + pp->ixmin - 1;
                if (bpm[iind2] == 0) 
                    tmp[n++] = data[iind2];
            }
        }
        d[i] = (double)vircam_med(tmp,NULL,(long)n);
    }

    /* Tidy and get out of here */

    freespace(tmp);
    return(d);
}

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

static double *vircam_linearity_tweakfac(double **fdata, double *mjd, int nim,
                                         int nchan, double *facrng, 
                                         double *maxdiff) {
    int i,ist,ifn,j;
    double *factors,sum,midval,minfac,maxfac;

    /* Get some memory for the output array */

    factors = cpl_malloc(nim*sizeof(double));

    /* First sort the data into order of mjd */

    vircam_mjdsort(fdata,mjd,nim);

    /* Which index is the midpoint? */

    if (nim % 2 == 0) {
        ist = nim/2 - 1;
        ifn = ist + 1;
    } else {
        ist = nim/2;
        ifn = ist;
    }

    /* Loop for each channel */

    for (i = 0; i < nchan; i++) {
        
        /* Get midpoint value */

        midval = 0.5*(fdata[ist][i] + fdata[ifn][i]);

        /* Now normalise all the ith channels by this value */

        for (j = 0; j < nim; j++) 
            fdata[j][i] /= midval;
    }

    /* Now loop for each image and average the values for all the channels in
       in image */

    *maxdiff = 0.0;
    for (j = 0; j < nim; j++) {
        sum = 0.0;
        for (i = 0; i < nchan; i++)
            sum += fdata[j][i];
        factors[j] = sum/(double)nchan;
        if (j == 0) {
            maxfac = factors[j];
            minfac = factors[j];
        } else {
            minfac = min(minfac,factors[j]);
            maxfac = max(maxfac,factors[j]);
            *maxdiff = max(*maxdiff,fabs(factors[j]-factors[j-1]));
        }
    }
    *facrng = maxfac - minfac;

    /* Get out of here */

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

static void vircam_mjdsort(double **fdata, double *mjd, int n) {
    int iii,ii,i,ifin,j,it;
    double tmpmjd,*tmpdata;


    iii = 2;
    while (iii < n)
        iii *= 2;
    iii = min(n,(3*iii)/4 - 1);

    while (iii > 1) {
        iii /= 2;
        ifin = n - iii;
        for (ii = 0; ii < ifin; ii++) {
            i = ii;
            j = i + iii;
            if (mjd[i] > mjd[j]) {
                tmpmjd = mjd[j];
                tmpdata = fdata[j];
                while (1) {
                    mjd[j] = mjd[i];
                    fdata[j] = fdata[i];
                    j = i;
                    i = i - iii;
                    if (i < 0 || mjd[0] <= tmpmjd) 
                        break;
                }
                mjd[j] = tmpmjd;
                fdata[j] = tmpdata;
            }
        }
    }
}

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

static cpl_table *vircam_linearity_analyse_diagtab_init(int np, int nrows) {
    int i;
    char colname[16];
    cpl_table *t;

    /* Create a new table */

    t = cpl_table_new(nrows);

    /* Add the first few columns */

    cpl_table_new_column(t,"filename",CPL_TYPE_STRING);
    cpl_table_new_column(t,"exptime",CPL_TYPE_DOUBLE);
    cpl_table_set_column_unit(t,"exptime","seconds");
    cpl_table_new_column(t,"mjd",CPL_TYPE_DOUBLE);
    cpl_table_set_column_unit(t,"mjd","days");

    /* Add columns for each of the channels' raw median flux and linearised 
       median flux */

    for (i = 1; i <= 16; i++) {
        (void)snprintf(colname,16,"rawflux_%02d",i);
        cpl_table_new_column(t,colname,CPL_TYPE_DOUBLE);
        cpl_table_set_column_unit(t,colname,"ADU");
        (void)snprintf(colname,16,"linflux_%02d",i);
        cpl_table_new_column(t,colname,CPL_TYPE_DOUBLE);
        cpl_table_set_column_unit(t,colname,"ADU");
    }

    /* Finally add the correction factors that were used */

    cpl_table_new_column(t,"adjust_fac",CPL_TYPE_DOUBLE);

    /* Right, get out of here */

    return(t);
}

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

static void vircam_linearity_analyse_init(void) {
    ps.labels = NULL;
    ps.domelist = NULL;
    ps.darklist = NULL;
    ps.domecheck = NULL;
    ps.darkcheck = NULL;
    ps.ndomes = 0;
    ps.ndarks = 0;
    ps.ndomecheck = 0;
    ps.ndarkcheck = 0;
    ps.chanfrm = NULL;
    ps.chantab = NULL;
    ps.lchantab = NULL;
    ps.flatlist = NULL;
    ps.bpm_array = NULL;
    ps.ddg = NULL;
    ps.plist = NULL;
    ps.elist = NULL;
    ps.phupaf = NULL;
    ps.diag1 = NULL;
    ps.diag2 = NULL;
}

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

static void vircam_linearity_analyse_tidy(int level) {
    int i;

    freetfits(ps.chantab);
    freearray(ps.bpm_array);
    freefitslist(ps.flatlist,ps.nflatlist);
    freetable(ps.lchantab);
    freepropertylist(ps.plist);
    freepropertylist(ps.elist);
    freetable(ps.diag1);
    freetable(ps.diag2);
    if (level == 1)
        return;
    
    freespace(ps.labels);
    freeframeset(ps.domelist);
    freeframeset(ps.darklist);
    freeframeset(ps.domecheck);
    freeframeset(ps.darkcheck);
    freeframeset(ps.sorted);
    freeframe(ps.chanfrm);
    if (ps.ddg != NULL) {
        for (i = 0; i < ps.nddg; i++) {
            freeframeset(ps.ddg[i].darks);
            freeframeset(ps.ddg[i].domes);
            freefitslist(ps.ddg[i].proc,ps.ddg[i].ndomes);
        }
        freespace(ps.ddg);
    }
    freepropertylist(ps.phupaf);
}

/* 

$Log: vircam_linearity_analyse.c,v $
Revision 1.51  2008/01/22 19:47:56  jim
New version to implement new algorithm

Revision 1.50  2007/11/26 09:58:49  jim
Now fails if given observation files done with NDIT != 1

Revision 1.49  2007/11/23 18:34:28  jim
fixed memory allocation bug

Revision 1.48  2007/11/22 12:36:55  jim
Modified to create diagnostic tables

Revision 1.47  2007/11/20 09:41:13  jim
Added ability to alter dome sequence stats by using the monitoring exposures

Revision 1.46  2007/11/14 10:42:25  jim
Substantial changes to incorporate new linearity analysis algorithm and to
restrict the amount of memory required to do the analysis (especially
the BPM work)

Revision 1.45  2007/09/07 13:32:12  jim
uses a sorted framelist to ensure that the correct information is given
to the output product header

Revision 1.44  2007/09/06 21:37:53  jim
fixed call to vircam_dfs_setup_product_ routines to use the full input
frameset

Revision 1.43  2007/08/29 09:20:33  jim
Primary header is now derived from the same header that forms the PAF rather
than starting off empty and allowing CPL to copy everything it thinks you
want...

Revision 1.42  2007/08/23 09:02:03  jim
Modified to check domes for DETLIVE before checking darks

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

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

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

Revision 1.38  2007/04/04 10:36:07  jim
Fixed typo preventing output of main PAF. Also modified to use dfs tags

Revision 1.37  2007/03/29 12:19:38  jim
Little changes to improve documentation

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

Revision 1.35  2007/02/19 21:13:04  jim
added bad pixel number QC parameter

Revision 1.34  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.33  2007/02/15 06:59:38  jim
Added ability to write QC paf files

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

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

Revision 1.30  2006/12/13 11:45:36  jim
Fixed scaling of sigma error

Revision 1.29  2006/12/11 22:47:12  jim
Fixed subtle bug in the way that stats were being done.

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

Revision 1.27  2006/11/10 09:23:46  jim
Fixed save routine so to use a new version of vircam_chantab_new

Revision 1.26  2006/10/31 10:27:27  jim
Fixed a few bugs and modified to make sure than an extension name appear
in each fits extension

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

Revision 1.24  2006/09/08 09:20:22  jim
major upgrade to main processing routine: to deal with bad input better; to
write out dummy results in the case of failure; to combine raw darks on the
fly for use in dark correction, rather than using master darks;

Revision 1.23  2006/08/03 13:26:44  jim
fixed another typo

Revision 1.22  2006/08/03 10:36:32  jim
Fixed typo

Revision 1.21  2006/06/20 19:06:38  jim
Added correction for ndit. Now adjusts the value of norder if the number
of frames given is too small for the order of polynomial requested

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

Revision 1.19  2006/06/06 13:03:42  jim
Fixed scaling that was causing funny stats

Revision 1.18  2006/05/27 21:40:06  jim
Bad pixels are now defined by a number of sigma above or below the mean

Revision 1.17  2006/05/09 09:30:47  jim
Fixed _save routine so that bad pixel mask is saved with unsigned char
data type

Revision 1.16  2006/05/08 12:32:12  jim
Changed default calling parameters for vircam_imcombine

Revision 1.15  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.14  2006/05/03 12:55:17  jim
Fixed some memory leaks

Revision 1.13  2006/05/02 13:26:32  jim
fixed bug where the wrong amount of memory was being allocated for the dark
exposure times

Revision 1.12  2006/05/02 11:36:29  jim
fixed illegal propertylist_delete calls

Revision 1.11  2006/04/27 09:46:01  jim
Modified DFS frame types to conform to new dictionary

Revision 1.10  2006/04/25 13:45:57  jim
Fixed to adhere to new calling sequence for vircam_dfs routines

Revision 1.9  2006/04/24 12:12:59  jim
Fixed --help documentation and sorted out filename extension problem
(.fit -> .fits)

Revision 1.8  2006/04/20 11:31:34  jim
Added bad pixel masking

Revision 1.7  2006/03/23 21:18:45  jim
Minor changes mainly to comment headers

Revision 1.6  2006/03/22 14:02:51  jim
cosmetic changes to keep lint happy

Revision 1.5  2006/03/22 12:13:51  jim
Modified to use new vircam_mask capability

Revision 1.4  2006/03/15 10:43:40  jim
Fixed a few things

Revision 1.3  2006/03/03 14:29:06  jim
Now calls routines with vir_fits.

Revision 1.2  2006/02/22 10:01:22  jim
Added full documentation

Revision 1.1  2006/02/18 11:49:58  jim
new file


*/

---