vircam_genlincur.c

/* $Id: vircam_genlincur.c,v 1.34 2012-01-15 17:40:09 jim Exp $
 *
 * This file is part of the VIRCAM 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: 2012-01-15 17:40:09 $
 * $Revision: 1.34 $
 * $Name: not supported by cvs2svn $
 */
 
/* Includes */
 
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
 
#include <cpl.h>
 
#include <math.h>
 
#include <casu_utils.h>
#include <casu_stats.h>
#include "vircam_utils.h"
#include "vircam_mods.h"
#include "vircam_pfits.h"
#include "vircam_channel.h"
 
 
#define SZCOLNAME 16
 
static double nom_val = 10000;
 
 
static double linval(double inval, double *kfacs, double tolerance,
                     int niter, double *b, int norder);
static double getkfac(long index, long ncpts, float reset_time, 
                      float read_time, float delay_time, float exptime);
static void getco(double *a, int nord, int m);

 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
extern int vircam_genlincur(double **fdata, int nimages, double *exps, 
                            double mindit, casu_tfits *chantab, 
                            int norder, cpl_table **lchantab, 
                            double **lindata, int *status) {
 
    const char *fctid = "vircam_genlincur";
    int retval,i,j,nbad,oldnorder,k,ii,nullval;
    long np;
    double *meds,sigfit,**aco,c0,lin_nom,*temp,*polyco,pt,*work,kfac;
    double sum,t10000,*kfacs;
    parquet *pp;
    cpl_table *ctab,*lc;
    cpl_array *exparray,*medsarray,*polyfitco,*workarray;
 
    /* Inherited status */
 
    *lchantab = NULL;
    if (*status != CASU_OK)
        return(*status);
 
    /* Check that you have enough images in the list */
 
    if (nimages < norder+1) {
        cpl_msg_error(fctid,
                      "Not enought images (%" CPL_SIZE_FORMAT ") for fit order (%" CPL_SIZE_FORMAT ")",
                      (cpl_size)nimages,(cpl_size)norder);
        FATAL_ERROR
    }
 
    /* Open the parquet structure for the channel table; for now we
       only need the number of channels in each structure, rather than
       the structure itself */
 
    ctab = casu_tfits_get_table(chantab);
    retval = vircam_chan_fill(ctab,&pp,&np);
    if (retval != CASU_OK) {
        *status = retval;
        return(retval);
    }
    vircam_chan_free(np,&pp);
 
    /* Create an output channel table. Copy the input channel table and then
       massage the linearity part */
 
    lc = cpl_table_duplicate(ctab);
    oldnorder = cpl_table_get_int(lc,"norder",0,&nullval);
    if (oldnorder > norder) {
        for (i = norder+1; i <= oldnorder; i++) {
            char * colname = cpl_sprintf("coeff_%d", i);
            //snprintf(colname,SZCOLNAME,"coeff_%d", i);
            cpl_table_erase_column(lc,colname);
            cpl_free(colname);
        }
    } else if (oldnorder < norder) {
        for (i = oldnorder+1; i <= norder; i++) {
            char * colname = cpl_sprintf("coeff_%d", i);
            //snprintf(colname,SZCOLNAME,"coeff_%d", i);
            if (cpl_table_has_column(lc,colname)) {
                cpl_free(colname);
                continue;
            }
            cpl_table_new_column(lc,colname,CPL_TYPE_DOUBLE);
            cpl_free(colname);
        }
    }
 
    /* Get some memory for the fitting arrays */
 
    exparray = cpl_array_wrap_double(exps,nimages);
    medsarray = cpl_array_new((cpl_size)nimages,CPL_TYPE_DOUBLE);
    meds = cpl_array_get_data_double(medsarray);
    aco = cpl_malloc(norder*sizeof(double *));
    for (i = 0; i < norder; i++) 
        aco[i] = cpl_malloc(norder*sizeof(double));
    temp = cpl_malloc(norder*sizeof(double));
    kfacs = cpl_malloc(norder*sizeof(double));
 
    /* Get memory for output array of linearised stats */
 
    *lindata = cpl_malloc(nimages*np*sizeof(double));
 
    /* Loop for each channel */
 
    nbad = 0;
    for (i = 0; i < np; i++) {
 
        /* Load the data up for this channel */
 
        for (j = 0; j < nimages; j++) 
            meds[j] = fdata[j][i];
 
        /* Do the initial fit */
 
        if (vircam_polyfit(exparray,medsarray,norder,1,2,2.0,2.0,&polyfitco,
                           &sigfit) != CASU_OK) {
            nbad++;
            cpl_table_set_int(lc,"norder",(cpl_size)i,norder);
            cpl_table_set_double(lc,"coeff_1",(cpl_size)i,1.0);
            for (k = 1; k < norder; k++) {
                char * colname = cpl_sprintf("coeff_%d",k+1);
                //snprintf(colname,SZCOLNAME,"coeff_%d",k+1);
                cpl_table_set_double(lc,colname,(cpl_size)i,0.0);
                cpl_free(colname);
            }
            freearray(polyfitco);
            continue;
        }
        polyco = cpl_array_get_data_double(polyfitco);
 
        /* Work out linearity */
 
        for (j = 0; j < nimages; j++) 
            if (meds[j] > nom_val)
                break;
 
        if (j == 0) { /* The first frame already exeeds nom_val */
            nbad++;
            cpl_table_set_int(lc,"norder",(cpl_size)i,norder);
            cpl_table_set_double(lc,"coeff_1",(cpl_size)i,1.0);
            for (k = 1; k < norder; k++) {
                char * colname = cpl_sprintf("coeff_%d",k+1);
                //snprintf(colname,SZCOLNAME,"coeff_%d",k+1);
                cpl_table_set_double(lc,colname,(cpl_size)i,0.0);
                cpl_free(colname);
            }
            freearray(polyfitco);
            continue;
        }
 
        t10000 = exps[j-1] + (nom_val - meds[j-1])/(meds[j] - meds[j-1]);
        sum = 0.0;
        for (j = 0; j < norder; j++)
            sum += (double)(j+1)*polyco[j]*pow(t10000,(double)j);
        lin_nom = 100.0*fabs(sum - polyco[0])/polyco[0];
 
        /* Get intermediate coefficients for matrix */
 
        for (j = 0; j < norder; j++) {
            getco(temp,norder,j+1);
            for (k = 0; k < norder; k++) {
                pt = pow(mindit,(double)(j-k));
                aco[j][k] = pt*temp[k];
            }
        }
 
        /* Solve matrix equation to do the back substitution */
 
        if (vircam_solve_gauss(aco,polyco,norder) != CASU_OK) {
            nbad++;
            cpl_table_set_int(lc,"norder",(cpl_size)i,norder);
            cpl_table_set_double(lc,"coeff_1",(cpl_size)i,1.0);
            for (k = 1; k < norder; k++) {
                char * colname = cpl_sprintf("coeff_%d",k+1);
                //snprintf(colname,SZCOLNAME,"coeff_%d",k+1);
                cpl_table_set_double(lc,colname,(cpl_size)i,0.0);
                cpl_free(colname);
            }
            freearray(polyfitco);
            continue;
        }
 
        /* Now normalise to unit slope and write the result to the table*/
 
        c0 = polyco[0];
        for (j = 0; j < norder; j++) {
            polyco[j] /= pow(c0,(double)(j+1));
            char * colname = cpl_sprintf("coeff_%d",j+1);
            //snprintf(colname,SZCOLNAME,"coeff_%d",j+1);
            cpl_table_set_double(lc,colname,(cpl_size)i,polyco[j]);
            cpl_free(colname);
        }
        cpl_table_set_int(lc,"norder",(cpl_size)i,norder);
 
        /* Work out how well the solution creates a linear system. Loop 
           for each input image and work out a 'linearised' median. Then
           do a linear fit to the linearsed median vs exposure time. */
 
        workarray = cpl_array_new((cpl_size)nimages,CPL_TYPE_DOUBLE);
        work = cpl_array_get_data_double(workarray);
        for (j = 0; j < nimages; j++) {
            kfac = mindit/exps[j];
            kfacs[0] = 1.0;
            for (ii = 1; ii < norder; ii++) 
                kfacs[ii] = pow((kfac+1.0),(double)(ii+1)) - 
                    pow(kfac,(double)(ii+1));
            work[j] = linval(meds[j],kfacs,0.5,10,polyco,norder);
            (*lindata)[j*np+i] = work[j];
        }       
        freearray(polyfitco);
        (void)vircam_polyfit(exparray,workarray,2,0,2,2.0,2.0,&polyfitco,
                             &sigfit);
        polyco = cpl_array_get_data_double(polyfitco);
        sigfit *= 100.0/nom_val;
        freearray(workarray);
        freearray(polyfitco);
                       
        /* Put the nominal linearity and the fit quality into the header */
 
        cpl_table_set_double(lc,"lin_10000_err",(cpl_size)i,sigfit);
        cpl_table_set_double(lc,"lin_10000",(cpl_size)i,lin_nom);
    }
 
    /* Tidy and get out of here */
 
    *lchantab = cpl_table_duplicate(lc);
    cpl_array_unwrap(exparray);
    freearray(medsarray);
    freespace2(aco,norder);
    freespace(temp);
    freetable(lc);
    freespace(kfacs);
    if (nbad != 0) {
        cpl_msg_warning(fctid,
                        "%" CPL_SIZE_FORMAT " channels have a failed solution",
                        (cpl_size)nbad);
        WARN_RETURN
    }
    GOOD_STATUS
}
 
 
/*---------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/
 
extern int vircam_lincor(casu_fits *infile, casu_tfits *lchantab, int kconst, 
                         int ndit, int *status) {
    int retval,i,norder,ii;
    long naxis[2],j,rind,aind,ncpts,np;
    float *data,texp,reset_time,read_time,delay_time;
    double kfac_nom,lkfac,inval,outval,*lbb,dnd,*kfacs;
    const char *fctid = "vircam_lincor";
    parquet *pp;
    cpl_propertylist *plist;
    cpl_table *lctab;
    parquet *p;
 
    /* Inherited status */
 
    if (*status != CASU_OK)
        return(*status);
 
    /* Do we even need to be here? */
 
    if (cpl_propertylist_has(casu_fits_get_ehu(infile),"ESO DRS LINCOR"))
        return(*status);
 
    /* Open the parquet structure for the channel table */
 
    lctab = casu_tfits_get_table(lchantab);
    retval = vircam_chan_fill(lctab,&p,&np);
    if (retval != CASU_OK)
        return(retval);
 
    /* Get the data array for the image */
 
    data = cpl_image_get_data(casu_fits_get_image(infile));
    if (data == NULL) {
        vircam_chan_free(np,&p);
        cpl_msg_error(fctid,"Error mapping data in input image");
        FATAL_ERROR
    }
    naxis[0] = (long)cpl_image_get_size_x(casu_fits_get_image(infile));
    naxis[1] = (long)cpl_image_get_size_y(casu_fits_get_image(infile));
 
    /* Get the required parameters from the header */
 
    plist = casu_fits_get_ehu(infile);
    if (vircam_pfits_get_exptime(plist,&texp) != CASU_OK) {
        vircam_chan_free(np,&p);
        cpl_msg_error(fctid,"No exposure time in %s",
                      casu_fits_get_fullname(infile));
        FATAL_ERROR
    }
    if (vircam_pfits_get_mindit(plist,&reset_time) != CASU_OK) {
        vircam_chan_free(np,&p);
        cpl_msg_error(fctid,"No mindit time in %s",
                      casu_fits_get_fullname(infile));
        FATAL_ERROR
    }
    read_time = reset_time;
    if (vircam_pfits_get_ditdelay(plist,&delay_time) != CASU_OK) {
        vircam_chan_free(np,&p);
        cpl_msg_error(fctid,"No dit delay time in %s",
                      casu_fits_get_fullname(infile));
        FATAL_ERROR
    }
      
    /* If there is a constant k factor, then calculate it now */
 
    kfac_nom = (double)(read_time/texp);
 
    /* Factor to take the number of DITs into account */
 
    dnd = (double)ndit;
        
    /* Loop for each channel now */
 
    for (i = 0; i < np; i++) {
        pp = p + i;
        ncpts = (pp->delta_i)*(pp->delta_j);
 
        /* Load up the fit coefficients. If there is only one coefficient
           this is by definition 1 and therefore we can skip this channel */
 
        norder = pp->norder;
        if (norder == 1)
            continue;
        lbb = pp->bb;
 
        /* Get workspace for K array and fill it in for situation of constant
           k factor */
 
        kfacs = cpl_malloc(norder*sizeof(double));
        if (kconst) {
            kfacs[0] = 1.0;
            for (ii = 1; ii < norder; ii++) 
                kfacs[ii] = pow((kfac_nom+1.0),(double)(ii+1)) - 
                    pow(kfac_nom,(double)(ii+1));
        }
 
        /* Now for each pixel */
 
        for (j = 0; j < ncpts; j++) {
 
            /* Get the 'non-constant' k-factor and fill in the K array. */
 
            rind = vircam_chan_d2r(pp,j);
            aind = vircam_chan_r2a(pp,naxis,rind);
            if (! kconst) {
                lkfac = getkfac(j,ncpts,reset_time,read_time,delay_time,texp);
                kfacs[0] = 1.0;
                for (ii = 1; ii < norder; ii++) 
                     kfacs[ii] = pow((lkfac+1.0),(double)(ii+1)) - 
                         pow(lkfac,(double)(ii+1));
            }
 
            /* Calculate the linearised value now */
 
            inval = ((double)data[aind])/dnd;
            outval = linval(inval,kfacs,0.5,10,lbb,norder);
            data[aind] = (float)(dnd*outval);
        }
        freespace(kfacs);
    }
 
    /* Add the linearity table to the DRS header */
 
    cpl_propertylist_update_string(casu_fits_get_ehu(infile),
                                   "ESO DRS LINCOR",
                                   casu_tfits_get_filename(lchantab));
 
    /* Right, get out of here */
 
    vircam_chan_free(np,&p);
    GOOD_STATUS
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
static double linval(double inval, double *kfacs, double tolerance,
                     int niter, double *b, int norder) {
    int jj,iter;
    double val_old,val,tol,sum;
 
    val = inval;
    iter = 0;
    tol = tolerance + 1.0;
    while (iter < niter && tol > tolerance) {
        val_old = val;
        iter++;
        sum = 0.0;
        for (jj = norder - 1; jj >= 1; jj--) 
            sum = (sum + b[jj]*kfacs[jj])*val;
        sum *= val;
        val = inval - sum;
        tol = fabs(val - val_old);
        if (val > 65535.0) {
            val = 65535.0;
            break;
        } else if (val < -1000.0) {
            val = -1000.0;
            break;
        }
    }
    return(val);
}
 
    
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
static double getkfac(long index, long npts, float reset_time, 
                      float read_time, float delay_time, float exptime) {
    double tkfac,dt1,dt2,dt3,dt4,df;
 
    df = ((double)index/(double)npts);
    dt1 = (double)exptime;
    dt2 = (double)read_time;
    dt3 = (double)reset_time;
    dt4 = (double)delay_time;
    tkfac = (dt3 + dt4 + (dt2 - dt3)*df)/dt1;
    return(tkfac);
}
 
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/
 
static void getco(double *a, int nord, int m) {
    int i,j,start;
 
    for (i = 0; i < nord; i++)
        a[i] = 0.0;
    start = m-1;
    a[start] = 1.0;
    j = 1;
    for (i = start-1; i >= 0; i--) {
        j++;
        a[i] = a[i+1]*(double)(m - j + 2)/(double)(j-1);
    }
}
 

 
 
/*
 
$Log: not supported by cvs2svn $
Revision 1.33  2010/06/03 12:15:31  jim
A few mods to get rid of compiler warnings
 
Revision 1.32  2009/12/11 06:53:35  jim
Minor changes to documentation
 
Revision 1.31  2009/09/09 09:45:36  jim
Modified to speed things up
 
Revision 1.30  2009/06/08 08:08:15  jim
Fixed memory leak and changed clipping parameters to vircam_polyfit
 
Revision 1.29  2009/05/20 12:18:42  jim
Modified so that if the operation is already done, then it just returns
 
Revision 1.28  2009/02/20 10:50:55  jim
Removed superfluous declarations
 
Revision 1.27  2008/10/21 08:38:48  jim
Final estiamte on linearity error done with vircam_polyfit now
 
Revision 1.26  2008/09/29 11:26:18  jim
Modified linval to provide a lower limit in case the fit goes crazy
 
Revision 1.25  2008/08/28 09:05:37  jim
Fixed bug where QC was being duplicated from master BPM on rare occasions.
Fixed bug where ARCFILE wasn't being written to the paf file for illum_cor
tables. Sky combine is done with medians
 
Revision 1.24  2008/01/22 19:45:24  jim
New version of genlincur to take into account the equality of readout and
reset time
 
Revision 1.23  2007/11/26 09:57:14  jim
Linearity correction routines now take account of ndit
 
Revision 1.22  2007/11/22 12:36:15  jim
Modified to return linearised values in an array
 
Revision 1.21  2007/11/20 09:38:57  jim
changed definition of fit quality to percentage error at 10000 counts
 
Revision 1.20  2007/11/14 14:47:32  jim
Modified the qualfit definition to be back in line with DRLD
 
Revision 1.19  2007/11/14 12:34:21  jim
Fixed header comments
 
Revision 1.18  2007/11/14 10:48:29  jim
Major rewrite to incorporate simpler and more robust algorithm
 
Revision 1.17  2007/03/29 12:19:39  jim
Little changes to improve documentation
 
Revision 1.16  2007/03/01 12:42:41  jim
Modified slightly after code checking
 
Revision 1.15  2006/11/27 12:08:18  jim
Modified lincor to get a better answer. Also modified the way the fit quality
is calculated
 
Revision 1.14  2006/09/29 11:19:31  jim
changed aliases on parameter names
 
Revision 1.13  2006/07/04 09:19:05  jim
replaced all sprintf statements with snprintf
 
Revision 1.12  2006/06/09 11:26:26  jim
Small changes to keep lint happy
 
Revision 1.11  2006/04/20 11:23:15  jim
Now medians the data before accumulation in the event that k is constant.
 
Revision 1.10  2006/03/23 21:18:47  jim
Minor changes mainly to comment headers
 
Revision 1.9  2006/03/22 13:36:50  jim
cosmetic changes to keep lint happy
 
Revision 1.8  2006/03/15 10:43:41  jim
Fixed a few things
 
Revision 1.7  2006/03/08 14:32:21  jim
Lots of little modifications
 
Revision 1.6  2006/03/03 14:29:46  jim
Modified definition of vir_fits and channel table
 
Revision 1.5  2006/03/01 10:31:28  jim
Now uses new vir_fits objects
 
Revision 1.4  2006/02/18 11:45:59  jim
Fixed a couple of memory bugs
 
Revision 1.3  2006/01/23 22:58:14  jim
Added vircam_lincor module
 
Revision 1.2  2006/01/23 16:06:03  jim
Added in header comments and section to evaluate the goodness of fit
 
Revision 1.1  2006/01/23 10:31:56  jim
New file
 
 
*/