vircam_genlincur.c

/* $Id: vircam_genlincur.c,v 1.24 2008/01/22 19:45:24 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:45:24 $
 * $Revision: 1.24 $
 * $Name:  $
 */

/* Includes */

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

#include <cpl.h>

#include <math.h>

#include "vircam_mods.h"
#include "vircam_utils.h"
#include "vircam_stats.h"
#include "vircam_pfits.h"
#include "vircam_channel.h"


#define SZCOLNAME 16

static double nom_val = 10000;
static double nom_kfac = 0.1;


static double linval(double inval, double k, 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, vir_tfits *chantab, 
                            int norder, cpl_table **lchantab, 
                            double **lindata, int *status) {

    const char *fctid = "vircam_genlincur";
    int retval,i,j,nbad,oldnorder,nullval,k;
    long np;
    double *meds,sigfit,**aco,c0,val,lin_nom,*temp,*polyco,pt,*work,kfac;
    double a0,a1,sum,sumsq;
    parquet *p,*pp;
    cpl_table *ctab,*lc;
    cpl_array *exparray,*medsarray,*polyfitco;
    char colname[SZCOLNAME];

    /* Inherited status */

    *lchantab = NULL;
    if (*status != VIR_OK)
        return(*status);

    /* Check that you have enough images in the list */

    if (nimages < norder+1) {
        cpl_msg_error(fctid,"Not enought images (%d) for fit order (%d)",
                      nimages,norder);
        FATAL_ERROR
    }

    /* Open the parquet structure for the channel table */

    ctab = vircam_tfits_get_table(chantab);
    retval = vircam_chan_fill(ctab,&pp,&np);
    if (retval != VIR_OK) {
        *status = retval;
        return(retval);
    }

    /* 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++) {
            snprintf(colname,SZCOLNAME,"coeff_%d",i);
            cpl_table_erase_column(lc,colname);
        }
    } else if (oldnorder < norder) {
        for (i = oldnorder+1; i <= norder; i++) {
            snprintf(colname,SZCOLNAME,"coeff_%d",i);
            if (cpl_table_has_column(lc,colname)) 
                continue;
            cpl_table_new_column(lc,colname,CPL_TYPE_DOUBLE);
        }
    }

    /* Get some memory for the fitting arrays */

    exparray = cpl_array_wrap_double(exps,nimages);
    medsarray = cpl_array_new(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));

    /* 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++) {
        p = pp + 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,1,3.0,3.0,&polyfitco,
                           &sigfit) != VIR_OK) {
            nbad++;
            cpl_table_set_int(lc,"norder",i,norder);
            cpl_table_set_double(lc,"coeff_1",i,1.0);
            for (k = 1; k < norder; k++) {
                snprintf(colname,SZCOLNAME,"coeff_%d",k+1);
                cpl_table_set_double(lc,colname,i,0.0);
            }
            continue;
        }
        polyco = cpl_array_get_data_double(polyfitco);

        /* 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) != VIR_OK) {
            nbad++;
            cpl_table_set_int(lc,"norder",i,norder);
            cpl_table_set_double(lc,"coeff_1",i,1.0);
            for (k = 1; k < norder; k++) {
                snprintf(colname,SZCOLNAME,"coeff_%d",k+1);
                cpl_table_set_double(lc,colname,i,0.0);
            }
            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));
            snprintf(colname,SZCOLNAME,"coeff_%d",j+1);
            cpl_table_set_double(lc,colname,i,polyco[j]);
        }
        cpl_table_set_int(lc,"norder",i,norder);

        /* Work out the percentage of non-linearity for nominal kfac and
           intensity values */

        val = linval(nom_val,nom_kfac,0.5,10,polyco,norder);
        lin_nom = 100*fabs(val - nom_val)/nom_val;

        /* 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. */

        work = cpl_malloc(nimages*sizeof(double));
        for (j = 0; j < nimages; j++) {
            kfac = mindit/exps[j];
            work[j] = linval(meds[j],kfac,0.5,10,polyco,norder);
            (*lindata)[j*np+i] = work[j];
        }
        vircam_linfit(nimages,exps,work,&a0,&a1,&sigfit);
        sigfit *= 100.0/nom_val;
        freearray(polyfitco)
                       
        /* Put the nominal linearity and the fit quality into the header */

        cpl_table_set_double(lc,"lin_10000_err",i,sigfit);
        cpl_table_set_double(lc,"lin_10000",i,lin_nom);
        freespace(work);
    }

    /* Tidy and get out of here */

    *lchantab = cpl_table_duplicate(lc);
    cpl_array_unwrap(exparray);
    freearray(medsarray);
    freespace2(aco,norder);
    freespace(temp);
    freetable(lc);
    vircam_chan_free(np,&pp);
    if (nbad != 0) {
        cpl_msg_warning(fctid,"%d channels have a failed solution",nbad);
        WARN_RETURN
    }
    GOOD_STATUS
}


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

extern int vircam_lincor(vir_fits *infile, vir_tfits *lchantab, int kconst, 
                         int ndit, int *status) {
    int retval,i,norder;
    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;
    const char *fctid = "vircam_lincor";
    parquet *pp;
    cpl_propertylist *plist;
    cpl_table *lctab;
    parquet *p;

    /* Inherited status */

    if (*status != VIR_OK)
        return(*status);

    /* Open the parquet structure for the channel table */

    lctab = vircam_tfits_get_table(lchantab);
    retval = vircam_chan_fill(lctab,&p,&np);
    if (retval != VIR_OK)
        return(retval);

    /* Get the data array for the image */

    data = cpl_image_get_data(vircam_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(vircam_fits_get_image(infile));
    naxis[1] = (long)cpl_image_get_size_y(vircam_fits_get_image(infile));

    /* Get the required parameters from the header */

    plist = vircam_fits_get_ehu(infile);
    if (vircam_pfits_get_exptime(plist,&texp) != VIR_OK) {
        vircam_chan_free(np,&p);
        cpl_msg_error(fctid,"No exposure time in %s",
                      vircam_fits_get_fullname(infile));
        FATAL_ERROR
    }
    if (vircam_pfits_get_mindit(plist,&reset_time) != VIR_OK) {
        vircam_chan_free(np,&p);
        cpl_msg_error(fctid,"No mindit time in %s",
                      vircam_fits_get_fullname(infile));
        FATAL_ERROR
    }
    read_time = reset_time;
    if (vircam_pfits_get_ditdelay(plist,&delay_time) != VIR_OK) {
        vircam_chan_free(np,&p);
        cpl_msg_error(fctid,"No dit delay time in %s",
                      vircam_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;

        /* Now for each pixel */

        for (j = 0; j < ncpts; j++) {

            /* Get the k-factor */

            rind = vircam_chan_d2r(pp,j);
            aind = vircam_chan_r2a(pp,naxis,rind);
            if (kconst == 1)
                lkfac = kfac_nom;
            else
                lkfac = getkfac(j,ncpts,reset_time,read_time,delay_time,texp);

            /* Calculate the linearised value now */

            inval = ((double)data[aind])/dnd;
            outval = linval(inval,lkfac,0.5,10,lbb,norder);
            data[aind] = (float)(dnd*outval);
        }
    }

    /* Add the linearity table to the DRS header */

    cpl_propertylist_update_string(vircam_fits_get_ehu(infile),
                                   "ESO DRS LINCOR",
                                   vircam_tfits_get_filename(lchantab));

    /* Right, get out of here */

    vircam_chan_free(np,&p);
    GOOD_STATUS
}

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

static double linval(double inval, double k, double tolerance, int niter,
                     double *b, int norder) {
    int jj,iter;
    double sc1,val_old,val,tol;

    val = inval;
    iter = 0;
    tol = tolerance + 1.0;
    while (iter < niter && tol > tolerance) {
        val_old = val;
        iter++;
        sc1 = inval;
        for (jj = 1; jj < norder; jj++)
            sc1 -= b[jj]*pow(val,(double)jj+1)*
                ((pow((1.0+k),(double)jj+1) - pow(k,(double)(jj+1))));
        val = sc1;
        tol = fabs(val - val_old);
    }
    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: vircam_genlincur.c,v $
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


*/

---