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

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

/* Includes */

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

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

#include "vmutils.h"
#include "vimos_imaging_utils.h"
#include "vimos_dfs.h"
#include "casu_stats.h"
#include "casu_fits.h"
#include "casu_mask.h"
#include "casu_utils.h"
#include "casu_mods.h"

/* Structure definitions */

typedef struct {

    /* Input */

    float       thresh;

    /* Output */

    char        extname[11];
    float       readnoise;
    float       gain;
    float       covar;

} configstruct;

typedef struct {
    cpl_size          *labels;
    cpl_size          *labels2;
    cpl_size          *labels3;
    cpl_frameset      *biaslist;
    cpl_frameset      *flatlist;
    cpl_image         *biasim1;
    cpl_image         *biasim2;
    cpl_image         *flatim1;
    cpl_image         *flatim2;
    casu_mask         *master_mask;
    cpl_table         *tab;
} memstruct;

static char vimos_ima_det_noise_description[] =
"vimos_ima_det_noise -- VIMOS recipe to find readnoise and gain.\n\n"
"Compute the readnoise and gain for the vimos detectors using two\n"
"bias frames and two flat field frames.\n"
"The program accepts the following files in the SOF:\n\n"
"    Tag                   Description\n"
"    -----------------------------------------------------------------------\n"
"    %-21s A list of raw bias images\n"
"    %-21s A list of raw flat images\n"
"    %-21s Optional master bad pixel mask or\n"
"    %-21s Optional master confidence map\n"
"If a master bad pixel mask or confidence map is included then the"
"bad pixels will be flagged during the statistical analysis\n";

/* Function prototypes */

static int vimos_ima_det_noise_create(cpl_plugin *);
static int vimos_ima_det_noise_exec(cpl_plugin *);
static int vimos_ima_det_noise_destroy(cpl_plugin *);
static int vimos_ima_det_noise(cpl_parameterlist *, cpl_frameset *) ;
static int vimos_ima_det_noise_save(cpl_frameset *framelist, 
                                     cpl_parameterlist *parlist, 
                                     configstruct *cs, memstruct *ps,
                                     int nexten, cpl_propertylist *phu,
                                     cpl_frame **product_frame);
static void vimos_ima_det_noise_init(memstruct *ps);
static void vimos_ima_det_noise_tidy(memstruct *ps, int level);
static void vimos_ima_det_noise_autocorr(float *data, unsigned char *bpm, 
                                          int nx, int ny, float sig, 
                                          float *covar);

/**@{*/

/**

    \ingroup recipelist
    \defgroup vimos_ima_det_noise vimos_ima_det_noise
    \brief Compute the detector readnoise and gain for VIMOS

    \par Name: 
        vimos_ima_det_noise
    \par Purpose: 
        Compute the detector readnoise and gain for VIMOS 
    \par Description: 
        A pair of bias frames and a pair of flat field frames are used to
        work out the readnoise and gain of the VIMOS detectors. 
    \par Language:
        C
    \par Parameters:
        - \b thresh (float): The rejection threshold in numbers of sigmas.
    \par Input File Types:
        The following list gives the file types that can appear in the SOF
        file. The word in bold is the DO CATG value.
        - \b BIAS (required): A list of raw bias images. Only 2 will be used
        - \b FLAT_TWILIGHT (required): A list of raw flat images. Only 2 
             will be used
        - \b MASTER_BPM or \b MASTER_CONF (optional): If either is given, then
             it will be used to mask out bad pixels during and statistical 
             analysis that is done.
    \par Output Products:
        The following list gives the output data products that are generated
        by this recipe. The word in bold gives the PRO CATG keyword value for
        each product:
        - The output readgain table (\b MASTER_READGAIN)
    \par Output QC Parameters:
        None
    \par Notes
        None.
    \par Fatal Error Conditions:
        - NULL input frameset
        - Input frameset headers incorrect meaning that RAW and CALIB frame
          can't be distinguished
        - No bias or flat frames in the input frameset
        - Unable to save data products
    \par Non-Fatal Error Conditions:
        - No master bad pixel or confidence map. All pixels considered to be 
          good.
    \par Conditions Leading To Dummy Products:
        - Image extensions wouldn't load.
        - The detector for the current image extension is flagged dead
    \par Author:
        Jim Lewis, CASU
    \par Code Reference: 
        vimos_ima_det_noise.c

*/

/* Function code */

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

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

int cpl_plugin_get_info(cpl_pluginlist *list) {
    cpl_recipe  *recipe = cpl_calloc(1,sizeof(*recipe));
    cpl_plugin  *plugin = &recipe->interface;
    char alldesc[SZ_ALLDESC];
    (void)snprintf(alldesc,SZ_ALLDESC,vimos_ima_det_noise_description,
                   VIMOS_BIAS_RAW,VIMOS_TWI_RAW,VIMOS_CAL_BPM,
                   VIMOS_CAL_CONF);

    cpl_plugin_init(plugin,
                    CPL_PLUGIN_API,
                    VIMOS_BINARY_VERSION,
                    CPL_PLUGIN_TYPE_RECIPE,
                    "vimos_ima_det_noise",
                    "Detector noise for imaging",
                    alldesc,
                    "Jim Lewis",
                    "jrl@ast.cam.ac.uk",
                    vimos_get_license(),
                    vimos_ima_det_noise_create,
                    vimos_ima_det_noise_exec,
                    vimos_ima_det_noise_destroy);

    cpl_pluginlist_append(list,plugin);

    return(0);
}

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

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

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

    /* Rejection threshold */

    p = cpl_parameter_new_range("vimos.vimos_ima_det_noise.thresh",
                                CPL_TYPE_DOUBLE,
                                "Rejection threshold in sigma above background",
                                "vimos.vimos_ima_det_noise",5.0,1.0e-6,1.0e10);
    cpl_parameter_set_alias(p,CPL_PARAMETER_MODE_CLI,"thresh");
    cpl_parameterlist_append(recipe->parameters,p);

    /* Get out of here */

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

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

    /* Get the recipe out of the plugin */

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

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

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

    /* Get the recipe out of the plugin */

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

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

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

static int vimos_ima_det_noise(cpl_parameterlist *parlist, 
                                cpl_frameset *framelist) {
    const char *fctid="vimos_ima_det_noise";
    cpl_parameter *p;
    cpl_propertylist *phu;
    cpl_frame *prod;
    cpl_frameset *chiplist_b,*chiplist_f;
    memstruct ps;
    configstruct cs;
    int j,nx,ny,npts,retval,nbias,nflats,loadextn,i,status,nbsize=64,k,i2;
    cpl_size nlab,nlab2,nlab3,nx1,nx2,ny1,ny2;
    unsigned char *bpm,*bpm2;
    cpl_binary *bpmcpl;
    float *databias1,*databias2,*dataflat1,*dataflat2,meanbias1,sigbias;
    float lcut,hcut,meanbias2,meanflat1,sigflat,meanflat2,meanbiasdiff;
    float sigbiasdiff,meanflatdiff,sigflatdiff,counts,avbackbias,*skymap;
    float avbackflat,covar;
    char *vimos_names[VIMOS_NEXTN];

    /* Check validity of input frameset */

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

    /* Initialise some things */

    vimos_ima_det_noise_init(&ps);

    /* Get the parameters */

    p = cpl_parameterlist_find(parlist,"vimos.vimos_ima_det_noise.thresh");
    cs.thresh = (float)cpl_parameter_get_double(p);

    /* Sort out raw from calib frames */

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

    /* Get a list of the frame labels */ 

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

    /* Get the bias frames */

    if ((ps.biaslist = casu_frameset_subgroup(framelist,ps.labels,nlab,
                                              VIMOS_BIAS_RAW)) == NULL) {
        cpl_msg_error(fctid,"Cannot find bias frames in input frameset");
        vimos_ima_det_noise_tidy(&ps,0);
        return(-1);
    }
    nbias = cpl_frameset_get_size(ps.biaslist)/4;
    if (nbias < 2) {
        cpl_msg_error(fctid,"Only %d bias images. Need at least 2 in each chip",
                      nbias);
        vimos_ima_det_noise_tidy(&ps,0);
        return(-1);
    }

    /* Get the flat frames */

    if ((ps.flatlist = casu_frameset_subgroup(framelist,ps.labels,nlab,
                                              VIMOS_TWI_RAW)) == NULL) {
        cpl_msg_error(fctid,"Cannot find flat frames in input frameset");
        vimos_ima_det_noise_tidy(&ps,0);
        return(-1);
    }
    nflats = cpl_frameset_get_size(ps.flatlist)/4;
    if (nflats < 2) {
        cpl_msg_error(fctid,"Only %d flat images. Need at least 2 in each chip",
                      nflats);
        vimos_ima_det_noise_tidy(&ps,0);
        return(-1);
    }

    /* Now labelise the input framesets so that we separate them by the
       detector name */

    if ((ps.labels2 = cpl_frameset_labelise(ps.biaslist,vimos_compare_names,
                                            &nlab2)) == NULL) {
        cpl_msg_error(fctid,"Cannot labelise the bias frames by detector");
        vimos_ima_det_noise_tidy(&ps,0);
        return(-1);
    }
    if ((ps.labels3 = cpl_frameset_labelise(ps.flatlist,vimos_compare_names,
                                            &nlab3)) == NULL) {
        cpl_msg_error(fctid,"Cannot labelise the flat frames by detector");
        vimos_ima_det_noise_tidy(&ps,0);
        return(-1);
    }

    /* See what chip set this is */

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

    /* 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 = casu_mask_define(framelist,ps.labels,nlab,VIMOS_CAL_CONF,
                                      VIMOS_CAL_BPM);

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

    for (j = 1; j <= VIMOS_NEXTN; j++) {
        cs.readnoise = 0.0;
        cs.gain = 0.0;

         /* Get the subset of images in the input frameset that belong to
           each detector */

        if ((chiplist_b = vimos_get_det_frameset(ps.biaslist,ps.labels2,nlab2,
                                                 vimos_names[j-1])) == NULL) {
            cpl_msg_info(fctid,"No bias in input list for detector %s\n",
                         vimos_names[j-1]);
            phu = cpl_propertylist_new();
            vimos_ima_det_noise_save(framelist,parlist,&cs,&ps,j,phu,&prod);
            vimos_ima_det_noise_tidy(&ps,1);
            cpl_propertylist_delete(phu);
            continue;
        }
        if ((chiplist_f = vimos_get_det_frameset(ps.flatlist,ps.labels3,nlab3,
                                                 vimos_names[j-1])) == NULL) {
            cpl_msg_info(fctid,"No flats in input list for detector %s\n",
                         vimos_names[j-1]);
            phu = cpl_propertylist_new();
            cpl_frameset_delete(chiplist_b);
            vimos_ima_det_noise_save(framelist,parlist,&cs,&ps,j,phu,&prod);
            vimos_ima_det_noise_tidy(&ps,1);
            cpl_propertylist_delete(phu);
            continue;
        }

        /* Work out the where the middle 1024x1024 square is. This is to 
           avoid the bad stuff at the edges. */

        loadextn = cpl_frame_get_nextensions(cpl_frameset_get_position(chiplist_b,0));
        ps.biasim1 = cpl_image_load(cpl_frame_get_filename(cpl_frameset_get_position(chiplist_b,0)),
                                    CPL_TYPE_FLOAT,0,loadextn);
        nx = (int)cpl_image_get_size_x(ps.biasim1);
        ny = (int)cpl_image_get_size_y(ps.biasim1);
        nx1 = (cpl_size)(nx/2 - 512);
        nx2 = nx1 + 1023;
        ny1 = (cpl_size)(ny/2 - 512);
        ny2 = ny1 + 1023;
        npts = (nx2 - nx1 + 1)*(ny2 - ny1 + 1);
        freeimage(ps.biasim1);

        /* Load up the images. If they won't load then signal an error,
           and save dummy results. */

        ps.biasim1 = cpl_image_load_window(cpl_frame_get_filename(cpl_frameset_get_position(chiplist_b,0)),
                                           CPL_TYPE_FLOAT,0,loadextn,nx1,ny1,nx2,ny2);
        ps.biasim2 = cpl_image_load_window(cpl_frame_get_filename(cpl_frameset_get_position(chiplist_b,1)),
                                           CPL_TYPE_FLOAT,0,loadextn,nx1,ny1,nx2,ny2);
        ps.flatim1 = cpl_image_load_window(cpl_frame_get_filename(cpl_frameset_get_position(chiplist_f,0)),
                                           CPL_TYPE_FLOAT,0,loadextn,nx1,ny1,nx2,ny2);
        ps.flatim2 = cpl_image_load_window(cpl_frame_get_filename(cpl_frameset_get_position(chiplist_f,1)),
                                           CPL_TYPE_FLOAT,0,loadextn,nx1,ny1,nx2,ny2);
        snprintf(cs.extname,11,"%s",vimos_names[j-1]);
        if (j == 1) {
            if (ps.biasim1 != NULL) 
                phu = cpl_propertylist_load(cpl_frame_get_filename(cpl_frameset_get_position(chiplist_b,0)),loadextn);
            else if (ps.biasim2 != NULL) 
                phu = cpl_propertylist_load(cpl_frame_get_filename(cpl_frameset_get_position(chiplist_b,1)),loadextn);
            else if (ps.flatim1 != NULL) 
                phu = cpl_propertylist_load(cpl_frame_get_filename(cpl_frameset_get_position(chiplist_f,0)),loadextn);
            else if (ps.flatim2 != NULL) 
                phu = cpl_propertylist_load(cpl_frame_get_filename(cpl_frameset_get_position(chiplist_f,1)),loadextn);
            else 
                phu = cpl_propertylist_new();
        }
        cpl_frameset_delete(chiplist_b);
        cpl_frameset_delete(chiplist_f);
        if (ps.biasim1 == NULL || ps.biasim2 == NULL || ps.flatim1 == NULL ||
            ps.flatim2 == NULL) {
            cpl_error_reset();
            cpl_msg_error(fctid,
                          "Extension %" CPL_SIZE_FORMAT " image wouldn't load",
                          (cpl_size)j);
            vimos_ima_det_noise_save(framelist,parlist,&cs,&ps,j,phu,&prod);
            vimos_ima_det_noise_tidy(&ps,1);
            freepropertylist(phu);
            continue;
        }

        /* Load the master mask extension */

        retval = casu_mask_load(ps.master_mask,j,nx,ny);
        if (retval == CASU_FATAL) {
            cpl_msg_info(fctid,
                         "Unable to load mask image %s[%" CPL_SIZE_FORMAT "]",
                         casu_mask_get_filename(ps.master_mask),(cpl_size)j);
            cpl_msg_info(fctid,"Forcing all pixels to be good from now on");
            casu_mask_force(ps.master_mask,nx,ny);
        }
        bpm2 = casu_mask_get_data(ps.master_mask);
        bpm = cpl_calloc(npts,sizeof(unsigned char));
        k = 0;
        for (i2 = ny1; i2 <= ny2; i2++) 
            for (i = nx1; i <= nx2; i++) 
                bpm[k++] = bpm2[(i2-1)*nx+i-1];            
            
        /* Do the statistical computations now... */

        cpl_msg_info(fctid,"Doing computation for extension %" CPL_SIZE_FORMAT,
                     (cpl_size)j);
        databias1 = cpl_image_get_data_float(ps.biasim1);
        databias2 = cpl_image_get_data_float(ps.biasim2);
        dataflat1 = cpl_image_get_data_float(ps.flatim1);
        dataflat2 = cpl_image_get_data_float(ps.flatim2);
        casu_qmedsig(databias1,bpm,(long)npts,cs.thresh,3,-1000.0,65535.0,
                     &meanbias1,&sigbias);
        lcut = meanbias1 - cs.thresh*sigbias;
        hcut = meanbias1 + cs.thresh*sigbias;
        casu_qmedsig(databias1,bpm,(long)npts,cs.thresh,3,lcut,hcut,
                     &meanbias1,&sigbias);
        casu_qmedsig(databias2,bpm,(long)npts,cs.thresh,3,-1000.0,65535.0,
                     &meanbias2,&sigbias);
        lcut = meanbias2 - cs.thresh*sigbias;
        hcut = meanbias2 + cs.thresh*sigbias;
        casu_qmedsig(databias2,bpm,(long)npts,cs.thresh,3,lcut,hcut,
                     &meanbias2,&sigbias);
        casu_qmedsig(dataflat1,bpm,(long)npts,cs.thresh,10,-1000.0,65535.0,
                     &meanflat1,&sigflat);
        lcut = meanflat1 - cs.thresh*sigflat;
        hcut = meanflat1 + cs.thresh*sigflat;
        casu_qmedsig(dataflat1,bpm,(long)npts,cs.thresh,10,lcut,hcut,
                     &meanflat1,&sigflat);
        casu_qmedsig(dataflat2,bpm,(long)npts,cs.thresh,10,-1000.0,65535.0,
                     &meanflat2,&sigflat);
        lcut = meanflat2 - cs.thresh*sigflat;
        hcut = meanflat2 + cs.thresh*sigflat;
        casu_qmedsig(dataflat2,bpm,(long)npts,cs.thresh,10,lcut,hcut,
                     &meanflat2,&sigflat);
        
        /* Form the difference images and do some more stats */

        cpl_image_subtract(ps.biasim1,ps.biasim2);
        cpl_image_subtract(ps.flatim1,ps.flatim2);

        /* Remove any residual large scale structure in both */

        bpmcpl = cpl_malloc(npts*sizeof(cpl_binary));
        for (i = 0; i < npts; i++)
            bpmcpl[i] = (cpl_binary)bpm[i];
        status = CASU_OK;
        casu_backmap(databias1,bpmcpl,1024,1024,nbsize,&avbackbias,&skymap,
                     &status);
        for (i = 0; i < npts; i++) 
            databias1[i] -= skymap[i];
        freespace(skymap);
        casu_backmap(dataflat1,bpmcpl,1024,1024,nbsize,&avbackflat,&skymap,
                     &status);
        for (i = 0; i < npts; i++) 
            dataflat1[i] -= skymap[i];
        freespace(skymap);
        freespace(bpmcpl);

        /* Do some more stats */
        
        casu_qmedsig(databias1,bpm,(long)npts,cs.thresh,3,-1000.0,65535.0,
                     &meanbiasdiff,&sigbiasdiff);
        lcut = meanbiasdiff - cs.thresh*sigbiasdiff;
        hcut = meanbiasdiff + cs.thresh*sigbiasdiff;
        casu_qmedsig(databias1,bpm,(long)npts,cs.thresh,3,lcut,hcut,
                     &meanbiasdiff,&sigbiasdiff);
        casu_qmedsig(dataflat1,bpm,(long)npts,cs.thresh,3,-65535.0,65535.0,
                     &meanflatdiff,&sigflatdiff);
        lcut = meanflatdiff - cs.thresh*sigflatdiff;
        hcut = meanflatdiff + cs.thresh*sigflatdiff;
        casu_qmedsig(dataflat1,bpm,(long)npts,cs.thresh,3,lcut,hcut,
                     &meanflatdiff,&sigflatdiff);

        /* Now remove the DC level in the flat difference image */

        for (i = 0; i < npts; i++)
            dataflat1[i] -= meanflatdiff;

        /* Do a autocorrelation and find out what the covariance is */
        
        vimos_ima_det_noise_autocorr(dataflat1,bpm,1024,1024,sigflatdiff,&covar);

        /* Now do the computation */

        counts = meanflat1 + meanflat2 - meanbias1 - meanbias2;
        cs.gain = counts/(sigflatdiff*sigflatdiff - sigbiasdiff*sigbiasdiff);
        cs.covar = covar;
        cs.readnoise = sigbiasdiff/sqrt(2.0);

        /* Save the result */

        vimos_ima_det_noise_save(framelist,parlist,&cs,&ps,j,phu,&prod);
        vimos_ima_det_noise_tidy(&ps,1);
        freepropertylist(phu);
        freespace(bpm);
    }
    vimos_ima_det_noise_tidy(&ps,0);
    return(0);
}

/*---------------------------------------------------------------------------*/
/**
  @brief    Do a 3x3 autocorrelation of an image
  @param    data         The input map
  @param    bpm          The input bad pixel mask
  @param    nx           The input x length of the map
  @param    ny           The input y length of the map
  @param    sig          The background sigma used for clipping
  @param    covar        The output covariance
  @return   nothing
 */
/*---------------------------------------------------------------------------*/

static void vimos_ima_det_noise_autocorr(float *data, unsigned char *bpm, 
                                          int nx, int ny, float sig, 
                                          float *covar) {
    int i,j,ioff,joff,n,ist,jst,ifn,jfn,ind,i1,j1,indoff;
    float clip,sum,av,acor[3][3],renorm;

    /* Set a clipping threshold */

    clip = 3.0*sig;

    /* Start by getting the average in the middle half of the array */

    ist = nx/4;
    ifn = 3*nx/4;
    jst = ny/4;
    jfn = 3*ny/4;
    n = 0;
    sum = 0.0;
    for (j = jst; j < jfn; j++) {
        for (i = ist; i < ifn; i++) {
            ind = j*ny + i;
            if (fabs(data[ind]) < clip && ! bpm[ind]) {
                n++;
                sum += data[ind];
            }
        }
    }
    av = sum/(float)n;
    
    /* Clear the auto-correlation array */

    for (i = 0; i < 3; i++)
        for (j = 0; j < 3; j++)
            acor[i][j] = 0.0;

    /* Now do the auto-correlation */

    for (j1 = 0; j1 < 3; j1++) {
        joff = -1 + j1;
        for (i1 = 0; i1 < 3; i1++) {
            ioff = -1 + i1;
            n = 0;
            sum = 0.0;
            for (j = jst; j < jfn; j++) {
                for (i = ist; i < ifn; i++) {
                    ind = j*nx + i;
                    indoff = (j + joff)*nx + (i + ioff);
                    if (fabs(data[ind]) < clip && ! bpm[ind] &&
                        fabs(data[indoff]) < clip && ! bpm[indoff]) {
                        n++;
                        sum += (data[ind] - av)*(data[indoff] - av);
                    }
                }
            }
            acor[i1][j1] = sum/(float)n;
        }
    }
    
    /* Now renormalise by the central value */

    renorm = acor[1][1];
    sum = 0.0;
    for (j = 0; j < 3; j++)
        for (i = 0; i < 3; i++)
            sum += acor[i][j]/renorm;

    /* Get out of here */

    *covar = sum;
}

/*---------------------------------------------------------------------------*/
/**
  @brief    The recipe data products are saved here
  @param    framelist    The input frame list
  @param    parlist      The input recipe parameter list
  @param    cs           The configuration structure defined in the main routine
  @param    ps           The memory structure defined in the main routine
  @param    nexten       The number of the extension
  @param    phu          The propertylist for the phu
  @param    product_frame  The frame for the output readgain product
  @return   0 if everything is ok
 */
/*---------------------------------------------------------------------------*/

static int vimos_ima_det_noise_save(cpl_frameset *framelist, 
                                     cpl_parameterlist *parlist, 
                                     configstruct *cs, memstruct *ps,
                                     int nexten, cpl_propertylist *phu,
                                     cpl_frame **product_frame) {
    const char *recipeid = "vimos_ima_det_noise";
    const char *fname = "readgain.fits";
    const char *fctid = "vimos_ima_det_noise_save";
    cpl_table *tab;
    cpl_propertylist *p;

    /* If we need to make a PHU then do that now. */

    if (nexten == 1) {

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

        *product_frame = cpl_frame_new();
        cpl_frame_set_filename(*product_frame,fname);
        cpl_frame_set_tag(*product_frame,VIMOS_PRO_READGAIN);
        cpl_frame_set_type(*product_frame,CPL_FRAME_TYPE_IMAGE);
        cpl_frame_set_group(*product_frame,CPL_FRAME_GROUP_PRODUCT);
        cpl_frame_set_level(*product_frame,CPL_FRAME_LEVEL_FINAL);

        /* Set the header for the primary */

        vimos_dfs_set_product_primary_header(phu,*product_frame,
                                             framelist,parlist,
                                             (char *)recipeid,"PRO-1.15",
                                             NULL,0);

        /* 'Save' the PHU image */                         

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

        /* Create the table */

        tab = cpl_table_new(4);
        cpl_table_new_column(tab,"EXTNAME",CPL_TYPE_STRING);
        cpl_table_new_column(tab,"READNOISE",CPL_TYPE_FLOAT);
        cpl_table_set_column_unit(tab,"READNOISE","ADUs");
        cpl_table_new_column(tab,"GAIN",CPL_TYPE_FLOAT);
        cpl_table_set_column_unit(tab,"GAIN","e-/ADU");
        cpl_table_new_column(tab,"COVAR",CPL_TYPE_FLOAT);
        cpl_table_set_column_unit(tab,"COVAR","");
        ps->tab = tab;
        
    }

    /* Fill in the most recent values */

    cpl_table_set_string(ps->tab,"EXTNAME",(cpl_size)(nexten-1),cs->extname);
    cpl_table_set_float(ps->tab,"READNOISE",(cpl_size)(nexten-1),cs->readnoise);
    cpl_table_set_float(ps->tab,"GAIN",(cpl_size)(nexten-1),cs->gain);
    cpl_table_set_float(ps->tab,"COVAR",(cpl_size)(nexten-1),cs->covar);

    /* If this is the last extension to be done, then save the table */

    if (nexten == VIMOS_NEXTN) {
        p = cpl_propertylist_new();
        vimos_dfs_set_product_exten_header(p,*product_frame,
                                           framelist,parlist,
                                           (char *)recipeid,
                                           "PRO-1.15",NULL);
        if (cpl_table_save(ps->tab,NULL,p,fname,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);
    }

    return(0);
}

/*---------------------------------------------------------------------------*/
/**
  @brief    Initialise the pointers in the memory structure
  @param    ps           The memory structure defined in the main routine
 */
/*---------------------------------------------------------------------------*/

static void vimos_ima_det_noise_init(memstruct *ps) {
    ps->labels = NULL;
    ps->labels2 = NULL;
    ps->labels3 = NULL;
    ps->biaslist = NULL;
    ps->flatlist = NULL;
    ps->biasim1 = NULL;
    ps->biasim2 = NULL;
    ps->flatim1 = NULL;
    ps->flatim2 = NULL;
    ps->master_mask = NULL;
    ps->tab = NULL;
}

/*---------------------------------------------------------------------------*/
/**
  @brief    Free any allocated workspace in the memory structure
  @param    ps           The memory structure defined in the main routine
  @param    level        The level of tidying. 
 */
/*---------------------------------------------------------------------------*/

static void vimos_ima_det_noise_tidy(memstruct *ps, int level) {
    
    freeimage(ps->biasim1);
    freeimage(ps->biasim2);
    freeimage(ps->flatim1);
    freeimage(ps->flatim2);

    if (level == 1)
        return;
    freespace(ps->labels);
    freespace(ps->labels2);
    freespace(ps->labels3);
    freeframeset(ps->biaslist);
    freeframeset(ps->flatlist);
    freemask(ps->master_mask);
    freetable(ps->tab);
}

/**@}*/
/*

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

Revision 1.8  2015/09/15 14:55:01  jim
Fixed call to casu_backmap

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

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

Revision 1.5  2015/06/08 09:38:31  jim
change clipping limits

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

Revision 1.3  2015/04/30 12:24:07  jim
Added autocorrelation

Revision 1.2  2015/02/14 12:38:09  jim
removed spurious print

Revision 1.1  2015/01/29 12:08:01  jim
New routine


*/