vircam_photcal.c

/* $Id: vircam_photcal.c,v 1.29 2007/10/25 17:34:01 jim Exp $
 *
 * This file is part of the VIRCAM Pipeline
 * Copyright (C) 2005 Cambridge Astronomy Survey Unit
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 */

/*
 * $Author: jim $
 * $Date: 2007/10/25 17:34:01 $
 * $Revision: 1.29 $
 * $Name:  $
 */

/* Includes */

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

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

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

typedef struct {
    char *filt;
    char **columns;
    int ncolumns;
    float extinct;
    float *coloureq;
    float offset;
    int nmags;
} photstrct;

static photstrct p;

#define NOMPIXSIZE 0.34
#define INITALLOC 1024
#define SZBUF 1024

static double pixsize (cpl_propertylist *plist);
static int vircam_phot_open(cpl_table *phottab, char *filt);
static void vircam_phot_close(void);
static int extract_columns(cpl_table *tab);
static int extract_coleq(cpl_table *tab);

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

extern int vircam_photcal(vir_fits **images, cpl_table **mstds, 
                          cpl_propertylist **pl, int nimages, char *filt, 
                          cpl_table *phottab, int *status) {
    float **stdmagptr,*resall3,*resall5,cdfudge,saturate,apcor3,apcor5,exptime;
    float airmass,*catcore3,*catcore5,*resim3,*resim5,cf,fluxmag3,fluxmag5;
    float refmag,extinct,dm3,dm5,med3,mad,cut,lcut,hcut,med5,sig3,sig5;
    float rcore,lim3,lim5;
    int nresall,nalloc_resall,i,j,k,ncat,nresim;
    char junk[SZBUF];
    const char *fctid = "vircam_photcal";
    vir_fits *im;
    cpl_propertylist *ehu_im,*ehu_cat;
    cpl_table *stds,*cl;

    /* Inherited status */

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

    /* Check for nonsense errors */

    if (nimages <= 0) {
        cpl_msg_error(fctid,"No images included in photometric calibration");
        FATAL_ERROR
    }

    /* Set up the structure that will give us the colour equations for
       this filter later on */

    if (vircam_phot_open(phottab,filt) != VIR_OK)
        FATAL_ERROR

    /* Get a workspace to hold the pointers to the magnitude columns */

    stdmagptr = cpl_malloc(p.ncolumns*sizeof(float *));

    /* Get some workspace to hold all the zeropoints for all the images
       in the input list. This is an initial allocation and more can be
       made available later if needed */

    resall3 = cpl_malloc(INITALLOC*sizeof(float));
    resall5 = cpl_malloc(INITALLOC*sizeof(float));
    nresall = 0;
    nalloc_resall = INITALLOC;

    /* Loop for the input images and catalogues. Create some shorthand
       variables and work out what the zeropoint fudge factor will be 
       based on the sampling of the current frame and the nominal sampling
       of a VIRCAM image */

    for (i = 0; i < nimages; i++) {
        im = images[i];
        ehu_im = vircam_fits_get_ehu(im);
        stds = mstds[i];
        ehu_cat = pl[i];
        cdfudge = 2.5*log10((double)pixsize(ehu_im)/NOMPIXSIZE);

        /* Now for the input catalogues. Start by getting some useful info
           from the header */

        saturate = cpl_propertylist_get_float(ehu_cat,"ESO QC SATURATION");
        apcor3 = cpl_propertylist_get_float(ehu_cat,"APCOR3");
        apcor5 = cpl_propertylist_get_float(ehu_cat,"APCOR5");
        rcore = cpl_propertylist_get_float(ehu_cat,"ESO DRS RCORE");
        (void)vircam_pfits_get_exptime(ehu_cat,&exptime);
        (void)vircam_pfits_get_airmass(vircam_fits_get_phu(im),&airmass);
        if (cpl_error_get_code() != CPL_ERROR_NONE) {
            cpl_msg_error(fctid,"Unable to get header info for %s",
                          vircam_fits_get_fullname(im));
            cpl_error_reset();
            continue;
        }

        /* Get objects that are not saturated and aren't too elliptical */

        cpl_table_select_all(stds);  
        cpl_table_and_selected_float(stds,"Ellipticity",CPL_LESS_THAN,0.5);
        cpl_table_and_selected_float(stds,"Peak_height",CPL_LESS_THAN,
                                     saturate);
        if (cpl_error_get_code() != CPL_ERROR_NONE) {
            cpl_msg_error(fctid,"Unable select data from matched stds tab %s",
                          vircam_fits_get_fullname(im));
            cpl_error_reset();
            continue;
        }

        /* Get all objects where the magnitude errors are less than 0.1 */

        for (j = 0; j < p.ncolumns; j++) {
            (void)snprintf(junk,SZBUF,"%ssig",(p.columns)[j]);
            cpl_table_and_selected_float(stds,junk,CPL_LESS_THAN,0.1);
        }

        /* Now extract all those that have passed the test. If none are left
           then issue a warning and move on to the next one */

        cl = cpl_table_extract_selected(stds);
        ncat = cpl_table_get_nrow(cl);
        if (ncat == 0) {
            cpl_msg_error(fctid,"No good standards available for %s",
                         vircam_fits_get_fullname(im));
            cpl_table_delete(cl);
            cpl_error_reset();
            continue;
        }

        /* Dereference some of the columns */

        catcore3 = cpl_table_get_data_float(cl,"Aper_flux_3");
        catcore5 = cpl_table_get_data_float(cl,"Aper_flux_5");
        for (j = 0; j < p.ncolumns; j++) 
            stdmagptr[j] = cpl_table_get_data_float(cl,(p.columns)[j]);

        /* Get some workspace for the results arrays for this image */

        resim3 = cpl_malloc(ncat*sizeof(float));
        resim5 = cpl_malloc(ncat*sizeof(float));
        nresim = 0;

        /* Loop for all the standards */

        extinct = p.extinct*(airmass - 1.0);
        for (j = 0; j < ncat; j++) {
            
            /* Do core magnitude calculation */

            cf = catcore3[j]/exptime;
            if (cf < 1.0)
                cf = 1.0;
            fluxmag3 = 2.5*log10((double)cf) + apcor3;
            cf = catcore5[j]/exptime;
            if (cf < 1.0)
                cf = 1.0;
            fluxmag5 = 2.5*log10((double)cf) + apcor5;

            /* Work out a reference magnitude */

            refmag = p.offset;
            for (k = 0; k < p.nmags; k++) 
                refmag += ((p.coloureq)[k]*stdmagptr[k][j]);

            /* Work out zero points and store them away for later */

            dm3 = refmag + fluxmag3 + extinct;
            dm5 = refmag + fluxmag5 + extinct;
            resim3[nresim] = dm3 + cdfudge;
            resim5[nresim++] = dm5 + cdfudge;
            resall3[nresall] = dm3 + cdfudge;
            resall5[nresall++] = dm5 + cdfudge;
            if (nresall == nalloc_resall) {
                nalloc_resall += INITALLOC;
                resall3 = cpl_realloc(resall3,nalloc_resall*sizeof(float));
                resall5 = cpl_realloc(resall5,nalloc_resall*sizeof(float));
            }
        }

        /* Ok, what is the mean zeropoint for this image? */

        (void)vircam_medmad(resim3,NULL,(long)nresim,&med3,&mad);
        cut = max(3.0*1.48*mad,0.1);
        lcut = med3 - cut;
        hcut = med3 + cut;
        (void)vircam_meansigcut(resim3,NULL,nresim,lcut,hcut,&med3,&sig3);
        (void)vircam_medmad(resim5,NULL,(long)nresim,&med5,&mad);
        cut = max(3.0*1.48*mad,0.1);
        lcut = med5 - cut;
        hcut = med5 + cut;
        (void)vircam_meansigcut(resim5,NULL,nresim,lcut,hcut,&med5,&sig5);

        /* Delete some workspace */

        freespace(resim3);
        freespace(resim5);
        freetable(cl);

        /* Calculate the limiting magnitudes */

        lim3 = med3 - 2.5*log10((5.0*sig3*rcore*sqrt(M_PI))/exptime) -
            apcor3 - extinct;
        lim5 = med5 - 2.5*log10((5.0*sig5*rcore*sqrt(M_PI))/exptime) -
            apcor5 - extinct;
        
        /* Write these results to the header of the image and the catalogue.
           First the QC headers for the image. These are not allowed to be
           copied into the catalogue headers...*/

        cpl_propertylist_update_float(ehu_im,"ESO QC MAGZPT",med3);
        cpl_propertylist_set_comment(ehu_im,"ESO QC MAGZPT",
                                     "[mag] photometric zeropoint");
        cpl_propertylist_update_float(ehu_im,"ESO QC MAGZERR",sig3);
        cpl_propertylist_set_comment(ehu_im,"ESO QC MAGZERR",
                                     "[mag] photometric zeropoint error");
        cpl_propertylist_update_int(ehu_im,"ESO QC MAGNZPT",nresim);
        cpl_propertylist_set_comment(ehu_im,"ESO QC MAGNZPT",
                                     "number of stars in magzpt calc");
        cpl_propertylist_update_float(ehu_im,"ESO QC LIMITING_MAG",lim3);
        cpl_propertylist_set_comment(ehu_im,"ESO QC LIMITING_MAG",
                                     "[mag] 5 sigma limiting mag.");
        
        /* Now the DRS headers for the image */

        cpl_propertylist_update_int(ehu_im,"ESO DRS MAGNZPTIM",nresim);
        cpl_propertylist_set_comment(ehu_im,"ESO DRS MAGNZPTIM",
                                     "number of stars in image magzpt calc");

        cpl_propertylist_update_float(ehu_im,"ESO DRS ZPIM1",med3);
        cpl_propertylist_set_comment(ehu_im,"ESO DRS ZPIM1",
                                     "[mag] zeropoint 1*rcore this image only");
        cpl_propertylist_update_float(ehu_im,"ESO DRS ZPSIGIM1",sig3);
        cpl_propertylist_set_comment(ehu_im,"ESO DRS ZPSIGIM1",
                                     "[mag] zeropoint sigma 1*rcore this image only");
        cpl_propertylist_update_float(ehu_im,"ESO DRS LIMIT_MAG1",lim3);
        cpl_propertylist_set_comment(ehu_im,"ESO DRS LIMIT_MAG1",
                                     "[mag] 5 sigma limiting mag 1*rcore.");
        cpl_propertylist_update_float(ehu_im,"ESO DRS ZPIM2",med5);
        cpl_propertylist_set_comment(ehu_im,"ESO DRS ZPIM2",
                                     "[mag] zeropoint 2*rcore this image only");
        cpl_propertylist_update_float(ehu_im,"ESO DRS ZPSIGIM2",sig5);
        cpl_propertylist_set_comment(ehu_im,"ESO DRS ZPSIGIM2",
                                     "[mag] zeropoint sigma 2*rcore this image only");
        cpl_propertylist_update_float(ehu_im,"ESO DRS LIMIT_MAG2",lim5);
        cpl_propertylist_set_comment(ehu_im,"ESO DRS LIMIT_MAG2",
                                     "[mag] 5 sigma limiting mag core5.");
        cpl_propertylist_update_float(ehu_im,"ESO DRS EXTINCT",p.extinct);
        cpl_propertylist_set_comment(ehu_im,"ESO DRS EXTINCT",
                                     "[mag] Assumed extinction.");

        /* And for the catalogue */

        cpl_propertylist_update_int(ehu_cat,"ESO DRS MAGNZPTIM",nresim);
        cpl_propertylist_set_comment(ehu_cat,"ESO DRS MAGNZPTIM",
                                     "number of stars in image magzpt calc");

        cpl_propertylist_update_float(ehu_cat,"ESO DRS ZPIM1",med3);
        cpl_propertylist_set_comment(ehu_cat,"ESO DRS ZPIM1",
                                     "[mag] zeropoint 1*rcore this image only");
        cpl_propertylist_update_float(ehu_cat,"ESO DRS ZPSIGIM1",sig3);
        cpl_propertylist_set_comment(ehu_cat,"ESO DRS ZPSIGIM1",
                                     "[mag] zeropoint sigma 1*rcore this image only");
        cpl_propertylist_update_float(ehu_cat,"ESO DRS LIMIT_MAG1",lim3);
        cpl_propertylist_set_comment(ehu_cat,"ESO DRS LIMIT_MAG1",
                                     "[mag] 5 sigma limiting mag 1*rcore.");
        cpl_propertylist_update_float(ehu_cat,"ESO DRS ZPIM2",med5);
        cpl_propertylist_set_comment(ehu_cat,"ESO DRS ZPIM2",
                                     "[mag] zeropoint 2*rcore this image only");
        cpl_propertylist_update_float(ehu_cat,"ESO DRS ZPSIGIM2",sig5);
        cpl_propertylist_set_comment(ehu_cat,"ESO DRS ZPSIGIM2",
                                     "[mag] zeropoint sigma 2*rcore this image only");
        cpl_propertylist_update_float(ehu_cat,"ESO DRS LIMIT_MAG2",lim5);
        cpl_propertylist_set_comment(ehu_cat,"ESO DRS LIMIT_MAG2",
                                     "[mag] 5 sigma limiting mag 2*rcore.");
        cpl_propertylist_update_float(ehu_im,"ESO DRS EXTINCT",p.extinct);
        cpl_propertylist_set_comment(ehu_im,"ESO DRS EXTINCT",
                                     "[mag] Assumed extinction.");

    }

    /* Ok, what is the mean zeropoint for all images? */

    if (nresall > 0) {
        (void)vircam_medmad(resall3,NULL,(long)nresall,&med3,&mad);
        cut = max(3.0*1.48*mad,0.1);
        lcut = med3 - cut;
        hcut = med3 + cut;
        (void)vircam_meansigcut(resall3,NULL,nresall,lcut,hcut,&med3,&sig3);
        (void)vircam_medmad(resall5,NULL,(long)nresall,&med5,&mad);
        cut = max(3.0*1.48*mad,0.1);
        lcut = med5 - cut;
        hcut = med5 + cut;
        (void)vircam_meansigcut(resall5,NULL,nresall,lcut,hcut,&med5,&sig5);
    }

    /* Delete some workspace */

    freespace(resall3);
    freespace(resall5);
    freespace(stdmagptr);
    vircam_phot_close();

    /* Write these results to the header of the images and the catalogues. 
       First the images */

    if (nresall > 0) {
        for (i = 0; i < nimages; i++) {
            im = images[i];
            ehu_im = vircam_fits_get_ehu(im);
            stds = mstds[i];
            ehu_cat = pl[i];
            cpl_propertylist_update_int(ehu_im,"ESO DRS MAGNZPTALL",nresall);
            cpl_propertylist_set_comment(ehu_im,"ESO DRS MAGNZPTALL",
                                         "number of stars in all magzpt calc");
            cpl_propertylist_update_float(ehu_im,"ESO DRS ZPALL1",med3);
            cpl_propertylist_set_comment(ehu_im,"ESO DRS ZPALL1",
                                         "[mag] zeropoint 1*rcore all group images");
            cpl_propertylist_update_float(ehu_im,"ESO DRS ZPSIGALL1",sig3);
            cpl_propertylist_set_comment(ehu_im,"ESO DRS ZPSIGALL1",
                                         "[mag] zeropoint sigma 1*rcore all group images");
            cpl_propertylist_update_float(ehu_im,"ESO DRS ZPALL2",med5);
            cpl_propertylist_set_comment(ehu_im,"ESO DRS ZPALL2",
                                         "[mag] zeropoint 2*rcore all group images");
            cpl_propertylist_update_float(ehu_im,"ESO DRS ZPSIGALL2",sig5);
            cpl_propertylist_set_comment(ehu_im,"ESO DRS ZPSIGALL2",
                                         "[mag] zeropoint sigma 2*rcore all group images");

            /* Now the catalogues */

            cpl_propertylist_update_int(ehu_cat,"ESO DRS MAGNZPTALL",nresall);
            cpl_propertylist_set_comment(ehu_cat,"ESO DRS MAGNZPTALL",
                                         "number of stars in all magzpt calc");
            cpl_propertylist_update_float(ehu_cat,"ESO DRS ZPALL1",med3);
            cpl_propertylist_set_comment(ehu_cat,"ESO DRS ZPALL1",
                                         "[mag] zeropoint 1*rcore all group images");
            cpl_propertylist_update_float(ehu_cat,"ESO DRS ZPSIGALL1",sig3);
            cpl_propertylist_set_comment(ehu_cat,"ESO DRS ZPSIGALL1",
                                         "[mag] zeropoint sigma 1*rcore all group images");
            cpl_propertylist_update_float(ehu_cat,"ESO DRS ZPALL2",med5);
            cpl_propertylist_set_comment(ehu_cat,"ESO DRS ZPALL2",
                                         "[mag] zeropoint 2*rcore all group images");
            cpl_propertylist_update_float(ehu_cat,"ESO DRS ZPSIGALL2",sig5);
            cpl_propertylist_set_comment(ehu_cat,"ESO DRS ZPSIGALL2",
                                         "[mag] zeropoint sigma 2*rcore all group images");
        }
    }

    /* Get out of here */

    GOOD_STATUS

}

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

extern int vircam_illum(vir_fits **images, cpl_table **mstds, 
                        cpl_propertylist **pl, int nimages, char *filt,
                        cpl_table *phottab, int nbsize, cpl_table **illcor, 
                        float *illcor_rms,int *status) {
    const char *fctid = "vircam_illum";
    char junk[SZBUF];
    float **stdmagptr,fracx,fracy,**results,cdfudge,saturate,apcor3,exptime;
    float airmass,*catcore3,*xx,*yy,extinct,cf,fluxmag3,refmag,dm3,med,mad;
    float xmin,xmax,ymin,ymax,*resall,medall,lcut,hcut,sig,cut,*good;
    int nx,ny,ifracx,ifracy,nbsizx,nbsizy,nbx,nby,*nres,*nall,i,j,ncat,k;
    int ix,iy,ind,nresall,nallocall,ngood;
    vir_fits *im;
    cpl_propertylist *ehu_im,*ehu_cat;
    cpl_table *stds,*cl;

    /* Inherited status */

    *illcor = NULL;
    *illcor_rms = 0.0;
    if (*status != VIR_OK)
        return(*status);

    /* Check for nonsense errors */

    if (nimages <= 0) {
        cpl_msg_error(fctid,"No images included in photometric calibration");
        FATAL_ERROR
    }

    /* Set up the structure that will give us the colour equations for
       this filter later on */

    if (vircam_phot_open(phottab,filt) != VIR_OK)
        FATAL_ERROR

    /* Get a workspace to hold the pointers to the magnitude columns */

    stdmagptr = cpl_malloc(p.ncolumns*sizeof(float *));

    /* Check to see if nbsize is close to exact divisor */

    nx = cpl_image_get_size_x(vircam_fits_get_image(images[0]));
    ny = cpl_image_get_size_y(vircam_fits_get_image(images[0]));
    fracx = ((float)nx)/((float)nbsize);
    fracy = ((float)ny)/((float)nbsize);
    ifracx = (int)(fracx + 0.1);
    ifracy = (int)(fracy + 0.1);
    nbsizx = nx/ifracx;
    nbsizy = ny/ifracy;
    nbsize = max(vircam_nint(0.9*nbsize),min(nbsize,min(nbsizx,nbsizy)));
    nbsize = min(nx,min(ny,nbsize)); /* trap for small maps */

    /* Divide the map into partitions */

    nbx = nx/nbsize;
    nby = ny/nbsize;

    /* Get some workspace to hold all the zeropoints for all the images
       in the input list. This is an initial allocation and more can be
       made available later if needed */

    results = cpl_malloc(nbx*nby*sizeof(float *));
    good = cpl_malloc(nbx*nby*sizeof(float));
    nres = cpl_calloc(nbx*nby,sizeof(int));
    nall = cpl_malloc(nbx*nby*sizeof(int));
    for (i = 0; i < nbx*nby; i++) {
        results[i] = cpl_malloc(INITALLOC*sizeof(float));
        nall[i] = INITALLOC;
    }
    resall = cpl_malloc(INITALLOC*sizeof(float));
    nresall = 0;
    nallocall = INITALLOC;

    /* Set up the output illumination correction table */

    *illcor = vircam_illcor_newtab(nbx*nby);
    
    /* Loop for the input images and catalogues. Create some shorthand
       variables and work out what the zeropoint fudge factor will be 
       based on the sampling of the current frame and the nominal sampling
       of a VIRCAM image */

    for (i = 0; i < nimages; i++) {
        im = images[i];
        ehu_im = vircam_fits_get_ehu(im);
        stds = mstds[i];
        if (stds == NULL)
            continue;
        ehu_cat = pl[i];
        cdfudge = 2.5*log10((double)pixsize(ehu_im)/NOMPIXSIZE);

        /* Now for the input catalogues. Start by getting some useful info
           from the header */

        saturate = cpl_propertylist_get_float(ehu_cat,"ESO QC SATURATION");
        apcor3 = cpl_propertylist_get_float(ehu_cat,"APCOR3");
        (void)vircam_pfits_get_exptime(vircam_fits_get_phu(im),&exptime);
        (void)vircam_pfits_get_airmass(vircam_fits_get_phu(im),&airmass);
        if (cpl_error_get_code() != CPL_ERROR_NONE) {
            cpl_msg_error(fctid,"Unable to get header info for %s",
                          vircam_fits_get_fullname(im));
            cpl_error_reset();
            continue;
        }

        /* Get objects that are not saturated and aren't too elliptical */

        cpl_table_select_all(stds);  
        cpl_table_and_selected_float(stds,"Ellipticity",CPL_LESS_THAN,0.5);
        cpl_table_and_selected_float(stds,"Peak_height",CPL_LESS_THAN,
                                     saturate);
        if (cpl_error_get_code() != CPL_ERROR_NONE) {
            cpl_msg_error(fctid,"Unable select data from matched stds tab %s",
                          vircam_fits_get_fullname(im));
            cpl_error_reset();
            continue;
        }

        /* Get all objects where the magnitude errors are less than 0.1 */

        for (j = 0; j < p.ncolumns; j++) {
            (void)snprintf(junk,SZBUF,"%ssig",(p.columns)[j]);
            cpl_table_and_selected_float(stds,junk,CPL_LESS_THAN,0.1);
        }

        /* Now extract all those that have passed the test. If none are left
           then issue a warning and move on to the next one */

        cl = cpl_table_extract_selected(stds);
        ncat = cpl_table_get_nrow(cl);
        if (ncat == 0) {
            cpl_msg_error(fctid,"No good standards available for %s",
                         vircam_fits_get_fullname(im));
            cpl_table_delete(cl);
            cpl_error_reset();
            continue;
        }

        /* Dereference some of the columns */

        catcore3 = cpl_table_get_data_float(cl,"Aper_flux_3");
        xx = cpl_table_get_data_float(cl,"X_coordinate");
        yy = cpl_table_get_data_float(cl,"Y_coordinate");
        for (j = 0; j < p.ncolumns; j++) 
            stdmagptr[j] = cpl_table_get_data_float(cl,(p.columns)[j]);

        /* Loop for all the standards */

        extinct = p.extinct*(airmass - 1.0);
        for (j = 0; j < ncat; j++) {
            
            /* Do core magnitude calculation */

            cf = catcore3[j]/exptime;
            if (cf < 1.0)
                cf = 1.0;
            fluxmag3 = 2.5*log10((double)cf) + apcor3;

            /* Work out a reference magnitude */

            refmag = p.offset;
            for (k = 0; k < p.nmags; k++) 
                refmag += ((p.coloureq)[k]*stdmagptr[k][j]);

            /* Work out zero point */

            dm3 = refmag + fluxmag3 + extinct + cdfudge;

            /* What cell does this belong to? */

            ix = (int)(xx[j]/(float)nbsize);
            iy = (int)(yy[j]/(float)nbsize);
            ind = iy*nbx + ix;
            if (nres[ind] == nall[ind] - 1) {
                results[ind] = cpl_realloc(results[ind],
                                        (nall[ind]+INITALLOC)*sizeof(float));
                nall[ind] += INITALLOC;
            }
            results[ind][nres[ind]] = dm3;
            nres[ind] += 1;

            /* Add this into the global solution too */

            if (nresall == nallocall) {
                resall = cpl_realloc(resall,(nallocall+INITALLOC)*sizeof(float));
                nallocall += INITALLOC;
            }
            resall[nresall++] = dm3;
        }

        /* Get rid of the table subset */

        freetable(cl);
    }

    /* What is the global zero point ? */

    if (nresall != 0) {
        (void)vircam_medmad(resall,NULL,(long)nresall,&medall,&mad);
        cut = max(3.0*1.48*mad,0.1);
        lcut = medall - cut;
        hcut = medall + cut;
        (void)vircam_meansigcut(resall,NULL,(long)nresall,lcut,hcut,
                                &medall,&sig);
    }

    /* Ok, what is the mean zeropoint for each cell. Do all of this stuff
       even if nresall == 0 so that we at least have a table with all
       zeros at the end of this */

    ngood = 0;
    for (i = 0; i < nbx*nby; i++) {
        if (nres[i] > 0) {
            (void)vircam_medmad(results[i],NULL,(long)nres[i],&med,&mad);
            cut = max(3.0*1.48*mad,0.3);
            lcut = med - cut;
            hcut = med + cut;
            (void)vircam_meansigcut(results[i],NULL,(long)nres[i],lcut,hcut,
                                    &med,&sig);
            med -= medall;
            good[ngood++] = med;
        } else {
            med = -99.0;
        }

        /* Where are we in the map? */

        iy = i/nbx;
        ix = i - iy*nbx;
        xmin = ix*nbsize;
        xmax = xmin + nbsize - 1;
        if (ix == nbx)
            xmax = nx;
        ymin = iy*nbsize;
        ymax = ymin + nbsize - 1;
        if (iy == nby)
            ymax = ny;

        /* Write out the results into illumination table */

        cpl_table_set_float(*illcor,"xmin",i,xmin);
        cpl_table_set_float(*illcor,"xmax",i,xmax);
        cpl_table_set_float(*illcor,"ymin",i,ymin);
        cpl_table_set_float(*illcor,"ymax",i,ymax);
        cpl_table_set_float(*illcor,"illcor",i,med);
    }

    /* Get the RMS of the map */

    if (ngood > 0)
        vircam_medsig(good,NULL,(long)ngood,&med,illcor_rms);
    else
        *illcor_rms = 0.0;
    
    /* Delete some workspace */

    for (i = 0; i < nbx*nby; i++)
        freespace(results[i]);
    freespace(results);
    freespace(nres);
    freespace(nall);
    freespace(stdmagptr);
    freespace(resall);
    freespace(good);
    vircam_phot_close();

    /* Set out the appropriate return status */

    if (nresall != 0)
        GOOD_STATUS
    else 
        WARN_RETURN

}

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

static double pixsize (cpl_propertylist *plist) {
    double cd1_1,cd1_2,pix;

    if (vircam_pfits_get_cd11(plist,&cd1_1) != VIR_OK ||
        vircam_pfits_get_cd12(plist,&cd1_2) != VIR_OK) {
        pix = NOMPIXSIZE;
    } else {
        pix = sqrt(cd1_1*cd1_1 + cd1_2*cd1_2);
    }
    return(pix);
}
    
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

static int vircam_phot_open(cpl_table *phottab, char *filt) {
    const char *fctid = "vircam_phot_open";
    int ns,null,nerr;
    cpl_table *subset;
    const char *req_cols[5] = {"filter","extinction","offset","columns",
                               "coleq"};

    /* Initialise a few things */

    p.coloureq = NULL;
    p.columns = NULL;
    p.nmags = 0;
    p.ncolumns = 0;
    p.offset = 0.0;

    /* Check the table and make sure it has all the required columns */

    nerr = 0;
    for (ns = 0; ns < 5; ns++) {
        if (! cpl_table_has_column(phottab,req_cols[ns])) {
            cpl_msg_error(fctid,"Photometry table missing column %s",
                          req_cols[ns]);
            nerr++;
        }
    }
    if (nerr > 0)
        return(VIR_FATAL);

    /* Search the table and find the rows that matches the filter */

    ns = cpl_table_and_selected_string(phottab,"filter",CPL_EQUAL_TO,filt);
    if (ns <= 0) {
        cpl_msg_error(fctid,"Unable to match photometry table to filter %s",
                      filt);
        return(VIR_FATAL);
    } else if (ns > 1) {
        cpl_msg_error(fctid,"More than one row matches filter %s",filt);
    }

    /* Read the information from the selected row */

    subset = cpl_table_extract_selected(phottab);
    p.filt = (char *)cpl_table_get_string(subset,"filter",0);
    p.extinct = cpl_table_get_float(subset,"extinction",0,&null);
    p.offset = cpl_table_get_float(subset,"offset",0,&null);
    if (extract_columns(subset) != VIR_OK) {
        freetable(subset);
        return(VIR_FATAL);
    }
    if (extract_coleq(subset) != VIR_OK) {
        freetable(subset);
        return(VIR_FATAL);
    }
    freetable(subset);

    /* Get out of here */

    return(VIR_OK);
}
    

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

static void vircam_phot_close(void) {
    int j;

    for (j = 0; j < p.ncolumns; j++)
        freespace((p.columns)[j]);
    freespace(p.columns);
    freespace(p.coloureq);
}

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

static int extract_columns(cpl_table *tab) {
    int nv,i,j;
    char *v,*w;

    /* Get the relevant value from the table */

    v = cpl_strdup(cpl_table_get_string(tab,"columns",0));

    /* Count the number of commas in the string to see how many 
       columns there are */

    nv = 1;
    j = strlen(v);
    for (i = 0; i < j; i++)
        if (v[i] == ',')
            nv++;
    p.ncolumns = nv;

    /* Now parse the string into the column names */

    p.columns = cpl_malloc(nv*sizeof(char *));
    for (i = 0; i < nv; i++) {
        if (i == 0) 
            w = strtok(v,",");
        else
            w = strtok(NULL,",");
        (p.columns)[i] = cpl_strdup(w);
    }
    freespace(v);
    return(VIR_OK);
}

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

static int extract_coleq(cpl_table *tab) {
    int nv,i,j;
    char *v,*w;

    /* Get the relevant value from the table */

    v = cpl_strdup(cpl_table_get_string(tab,"coleq",0));

    /* Count the number of commas in the string to see how many 
       columns there are */

    nv = 1;
    j = strlen(v);
    for (i = 0; i < j; i++)
        if (v[i] == ',')
            nv++;
    p.nmags = nv;

    /* Now parse the string into the colour equation values */

    p.coloureq = cpl_malloc(nv*sizeof(float));
    for (i = 0; i < nv; i++) {
        if (i == 0) 
            w = strtok(v,",");
        else
            w = strtok(NULL,",");
        (p.coloureq)[i] = (float)atof(w);
    }
    freespace(v);
    return(VIR_OK);
}

/*

$Log: vircam_photcal.c,v $
Revision 1.29  2007/10/25 17:34:01  jim
Modified to remove lint warnings

Revision 1.28  2007/10/19 06:55:06  jim
Modifications made to use new method for directing the recipes to the
standard catalogues using the sof

Revision 1.27  2007/05/15 08:52:01  jim
Fixed little bug in the clipping for the final magnitude zeropoint

Revision 1.26  2007/05/02 09:12:02  jim
Added extinction to DRS

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

Revision 1.24  2007/03/06 11:56:56  jim
Removed some QC parameters from the catalogues

Revision 1.23  2007/03/01 12:42:42  jim
Modified slightly after code checking

Revision 1.22  2007/02/25 06:34:20  jim
Plugged memory leak

Revision 1.21  2007/01/08 19:11:34  jim
Fixed incorrect FITS key comments

Revision 1.20  2006/12/19 13:30:34  jim
Initialise the output

Revision 1.19  2006/12/04 21:07:47  jim
Null value in illumination correction is changed to -99.0

Revision 1.18  2006/12/01 14:10:06  jim
trivial change to docs

Revision 1.17  2006/11/28 22:10:22  jim
Added illcor_rms to calling sequence of vircam_illum

Revision 1.16  2006/11/28 20:53:08  jim
Added vircam_illum

Revision 1.15  2006/11/27 12:03:12  jim
changed calls from cpl_propertylist_append to cpl_propertylist_update

Revision 1.14  2006/11/10 10:13:58  jim
Fixed error in docs

Revision 1.13  2006/10/02 13:47:33  jim
Added missing .h file to include list

Revision 1.12  2006/07/04 09:19:05  jim
replaced all sprintf statements with snprintf

Revision 1.11  2006/07/03 09:33:18  jim
Fixed a few things to keep the compiler happy

Revision 1.10  2006/06/13 14:09:09  jim
Better explanation of photcal tables might be wrong

Revision 1.9  2006/06/12 11:30:55  jim
Removed zeropoint column from photcal table

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

Revision 1.7  2006/06/06 13:07:14  jim
Modified some DRS keyword names

Revision 1.6  2006/05/31 14:23:32  jim
Fixed error message reset call

Revision 1.5  2006/05/31 11:40:09  jim
Added more QC and DRS parameters. Also added calculation of limiting
magnitudes. Tidied up some more docs

Revision 1.4  2006/05/30 14:27:54  jim
Added illcor parameter to call for photcal

Revision 1.3  2006/05/26 15:05:09  jim
Fixed lots of little bugs

Revision 1.2  2006/05/18 09:48:28  jim
Added docs

Revision 1.1  2006/05/17 12:05:30  jim
new file


*/

---