vircam_getstds.c

/* $Id: vircam_getstds.c,v 1.21 2008/05/06 08:40:10 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/05/06 08:40:10 $
 * $Revision: 1.21 $
 * $Name:  $
 */

/* Includes */

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

#include <sys/stat.h>
#include <sys/types.h>
#include <unistd.h>
#include <string.h>
#include <cpl.h>

#include <math.h>

#include "vircam_mods.h"
#include "vircam_utils.h"
#include "vircam_wcsutils.h"

#define CACHEDIR "catcache"
#define CACHEIND "catcache/index"
#define SZBUF 1024

static cpl_table *vircam_2mass_extract(char *path, float ramin, float ramax,
                                       float decmin, float decmax);
static cpl_table *check_cache(char *catname, float ra1_im, float ra2_im, 
                              float dec1_im, float dec2_im);
static void addto_cache(cpl_table *stds, char *catname, float ramin, 
                        float ramax, float decmin, float decmax);

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

extern int vircam_getstds(cpl_propertylist *plist, int cache, char *path,
                          char *catname, cpl_table **stds, int *status) {
    const char *fctid = "vircam_getstds";
    double x1,x2,y1,y2,r,d,xx,yy;
    float ramin,ramax,decmin,decmax,*ra,*dec,*x,*y;
    cpl_wcs *wcs;
    int n,i;
    cpl_propertylist *p;

    /* Inherited status */

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

    /* Get the coverage of the input WCS */

    (void)vircam_coverage(plist,0,&x1,&x2,&y1,&y2,status);
    ramin = (float)x1;
    ramax = (float)x2;
    decmin = (float)y1;
    decmax = (float)y2;

    /* If using the cache, then check if for a previously used table */

    *stds = NULL;
    if (cache)
        *stds = check_cache(catname,ramin,ramax,decmin,decmax);

    /* If there was nothing in the cache (or you didn't use it) then search
       the standard 2mass catalogues */

    if (*stds == NULL) {        

        /* Read the standards from the catalogue */

        *stds = vircam_2mass_extract(path,ramin,ramax,decmin,decmax);
        if (*stds == NULL) {
            cpl_msg_error(fctid,"Unable to extract data in %s\n",path);
            FATAL_ERROR
        }
        
        /* Add this table to the cache if you want to */

        if (cache)
            addto_cache(*stds,catname,ramin,ramax,decmin,decmax);
    }

    /* If there are no rows in the table, this may be a cause for concern.
       So add the columns into the table and then set a warning return 
       status */

    n = cpl_table_get_nrow(*stds);
    if (n == 0) {
        cpl_table_new_column(*stds,"xpredict",CPL_TYPE_FLOAT);
        cpl_table_new_column(*stds,"ypredict",CPL_TYPE_FLOAT);
        WARN_RETURN
    }

    /* Now fill the coordinates in */

    wcs = cpl_wcs_new_from_propertylist((const cpl_propertylist *)plist);
    ra = cpl_table_get_data_float(*stds,"RA");
    dec = cpl_table_get_data_float(*stds,"Dec");
    x = cpl_malloc(n*sizeof(*x));
    y = cpl_malloc(n*sizeof(*y));
    for (i = 0; i < n; i++) {
        r = (double)ra[i];
        d = (double)dec[i];
        vircam_radectoxy(wcs,r,d,&xx,&yy);
        x[i] = (float)xx;
        y[i] = (float)yy;
    }
    cpl_wcs_delete(wcs);
    
    /* Add the predicted x,y coordinates columns to the table */

    cpl_table_wrap_float(*stds,x,"xpredict");
    cpl_table_set_column_unit(*stds,"xpredict","pixels");
    cpl_table_wrap_float(*stds,y,"ypredict");
    cpl_table_set_column_unit(*stds,"ypredict","pixels");

    /* Finally sort this by ypredict */

    p = cpl_propertylist_new();
    cpl_propertylist_append_bool(p,"ypredict",0);
    cpl_table_sort(*stds,p);
    cpl_propertylist_delete(p);

    /* Get out of here */

    GOOD_STATUS
}


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

static cpl_table *vircam_2mass_extract(char *path, float ramin1, float ramax1, 
                                       float decmin, float decmax) {
    cpl_table *t,*s,*o;
    int i,nrows,start,finish,first_index,last_index,irow,init,j;
    int first_index_ra,last_index_ra,wrap,iwrap;
    float dectest,ratest,ramin,ramax;
    char fullname[SZBUF];
    cpl_array *a;
    char *deccol[] = {"Dec"};
    cpl_propertylist *p;

    /* Create an output table */

    o = cpl_table_new(0);
    init = 1;

    /* Create a cpl array */

    a = cpl_array_wrap_string(deccol,1);

    /* Is there a wrap around problem? */

    wrap = (ramin1 < 0.0 && ramax1 > 0.0) ? 2 : 1;

    /* Loop for each query. If there is a wrap around problem then we need 2
       queries. If not, then we only need 1 */

    for (iwrap = 0; iwrap < wrap; iwrap++) {
        if (wrap == 2) {
            if (iwrap == 0) {
                ramin = ramin1 + 360.0;
                ramax = 360.0;
            } else {
                ramin = 0.000001;
                ramax = ramax1;
            }
        } else {
            ramin = ramin1;
            ramax = ramax1;
        }

        /* Find out where in the index to look */

        first_index_ra = (int)ramin;
        last_index_ra = min((int)ramax,359);
        
        /* Look at the min and max RA and decide which files need to be 
           opened. */

        for (i = first_index_ra; i <= last_index_ra; i++) {

            /* Ok, we've found one that needs opening. Read the file with 
               the relevant CPL call */

            (void)snprintf(fullname,SZBUF,"%s/npsc%03d.fits",path,i);

            /* Read the propertylist so that you know how many rows there
               are in the table */

            p = cpl_propertylist_load(fullname,1);
            if (p == NULL) {
                freetable(o);
                cpl_array_unwrap(a);
                return(NULL);
            }
            nrows = cpl_propertylist_get_int(p,"NAXIS2");
            cpl_propertylist_delete(p);

            /* Load various rows until you find the Dec range that you 
               have specified. First the minimum Dec */

            start = 0;
            finish = nrows;
            first_index = nrows/2;
            while (finish - start >= 2) {
                t = cpl_table_load_window(fullname,1,0,a,first_index,1);
                dectest = cpl_table_get_float(t,"Dec",0,NULL);
                cpl_table_delete(t);
                if (dectest < decmin) {
                    start = first_index;
                    first_index = (first_index + finish)/2;
                } else {
                    finish = first_index;
                    first_index = (first_index + start)/2;
                }
            }

            /* Load various rows until you find the Dec range that you 
               have specified. Now the maximum Dec */

            start = first_index;
            finish = nrows;
            last_index = start + (finish - start)/2;
            while (finish - start >= 2) {
                t = cpl_table_load_window(fullname,1,0,a,last_index,1);
                dectest = cpl_table_get_float(t,"Dec",0,NULL);
                cpl_table_delete(t);
                if (dectest < decmax) {
                    start = last_index;
                    last_index = (last_index + finish)/2;
                } else {
                    finish = last_index;
                    last_index = (last_index + start)/2;
                }
            }    
            if (last_index < first_index)
                last_index = first_index;

            /* Ok now now load all the rows in the relevant dec limits */

            nrows = last_index - first_index + 1;
            if ((t = cpl_table_load_window(fullname,1,0,NULL,first_index,
                                           nrows)) == NULL) {
                freetable(o);
                cpl_array_unwrap(a);
                return(NULL);
            }
            cpl_table_unselect_all(t);

            /* Right, we now know what range of rows to search. Go through 
               these and pick the ones that are in the correct range of RA.
               If a row qualifies, then 'select' it. */

            for (j = 0; j < nrows; j++) {
                ratest = cpl_table_get_float(t,"RA",j,NULL);
                if (cpl_error_get_code() != CPL_ERROR_NONE) {
                    cpl_table_delete(t);
                    cpl_array_unwrap(a);
                    freetable(o);
                    return(NULL);
                }
                if (ratest >= ramin && ratest <= ramax)
                    cpl_table_select_row(t,j);
            }

            /* Extract the rows that have been selected now and append them
               onto the output table */

            s = cpl_table_extract_selected(t);
            if (init == 1) {
                cpl_table_copy_structure(o,t);
                init = 0;
            }
            irow = cpl_table_get_nrow(o) + 1;
            cpl_table_insert(o,s,irow);

            /* Tidy up */

            cpl_table_delete(t);
            cpl_table_delete(s);
        }
    }

    /* Ok, now just return the table and get out of here */

    cpl_array_unwrap(a);
    return(o);
}
        
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

static cpl_table *check_cache(char *catname, float ra1_im, float ra2_im, 
                              float dec1_im, float dec2_im) {
    int wrap1,wrap2;
    FILE *fd;
    char fname[BUFSIZ],catname2[SZBUF],cat_cache[SZBUF];
    float best,ra1_cat,ra2_cat,dec1_cat,dec2_cat,d1,d2,fra,fdec,ftot;
    cpl_table *out_cat;

    /* Open the index file.  NB the path and file name are hardcoded */

    fd = fopen(CACHEIND,"r");
    if (fd == NULL)
        return(NULL);

    /* Check to see if there is wrap around in the coordinates */

    wrap1 = (ra1_im < 0.0 ? 1 : 0);

    /* Now see if you have any matching entries */

    best = 0.0;
    while (fscanf(fd,"%s %s %g %g %g %g",fname,catname2,&ra1_cat,&ra2_cat,
                  &dec1_cat,&dec2_cat) != EOF) {
        wrap2 = (ra1_cat < 0.0 ? 1 : 0);
        if (wrap1 != wrap2)
            continue;
        if (strcmp(catname,catname2))
            continue;
        
        /* Check to see if there is at least some overlap */

        if (!(((ra1_im >= ra1_cat && ra1_im <= ra2_cat) ||
             (ra2_im >= ra1_cat && ra2_im <= ra2_cat)) &&
            ((dec1_im >= dec1_cat && dec1_im <= dec2_cat) ||
             (dec2_im >= dec1_cat && dec2_im <= dec2_cat))))
            continue;

        /* Work out exactly how much there is in each coordinate */

        d1 = max(0.0,ra1_cat-ra1_im);
        d2 = max(0.0,ra2_im-ra2_cat);
        fra = 1.0 - (d1 + d2)/(ra2_im - ra1_im);
        d1 = max(0.0,dec1_cat-dec1_im);
        d2 = max(0.0,dec2_im-dec2_cat);
        fdec = 1.0 - (d1 + d2)/(dec2_im - dec1_im);
        ftot = fra*fdec;

        /* Keep track of which is the best one */

        if (ftot > best) {
            snprintf(cat_cache,SZBUF,"%s/%s",CACHEDIR,fname);
            best = ftot;
        }
    }
    fclose(fd);

    /* Return a bad status if there isn't sufficient overlap */

    if (best < 0.9)
        return(NULL);

    /* If there is enough overlap, then try and read the FITS table. If it
       reads successfully, then return the table pointer */

    out_cat = cpl_table_load(cat_cache,1,0);
    return(out_cat);
}
    
/*---------------------------------------------------------------------------*/
/*---------------------------------------------------------------------------*/

static void addto_cache(cpl_table *stds, char *catname, float ramin, 
                        float ramax, float decmin, float decmax) {
    FILE *fd;
    char newname[SZBUF];
    int i;

    /* Check to see if the cache directory exists.  If it doesn't, then create
       it. */

    if (access(CACHEDIR,0) != 0)
        mkdir(CACHEDIR,0755);

    /* Open the index file with 'append' access */

    fd = fopen(CACHEIND,"a");

    /* Check the files in the directory to get a number that isn't already
       being used */

    i = 0;
    while (i >= 0) {
        i++;
        snprintf(newname,SZBUF,"%s/cch_%08d",CACHEDIR,i);
        if (access(newname,F_OK) != 0) 
            break;
    }

    /* Now write the current entry and make a copy of the file into the 
       directory */

    snprintf(newname,SZBUF,"%s/cch_%08d",CACHEDIR,i);
    cpl_table_save(stds,NULL,NULL,newname,CPL_IO_DEFAULT);
    snprintf(newname,SZBUF,"cch_%08d",i);
    if (cpl_error_get_code() == CPL_ERROR_NONE)
        fprintf(fd, "%s %s %g %g %g %g\n",newname,catname,ramin-0.0005,
                ramax+0.0005,decmin-0.0005,decmax+0.0005);
    fclose(fd);
}

/* 

$Log: vircam_getstds.c,v $
Revision 1.21  2008/05/06 08:40:10  jim
Modified to use cpl_wcs interface

Revision 1.20  2007/10/25 17:34:00  jim
Modified to remove lint warnings

Revision 1.19  2007/10/19 09:25:09  jim
Fixed problems with missing includes

Revision 1.18  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.17  2007/10/15 12:50:53  jim
Modified to read new version of 2mass catalogue files

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

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

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

Revision 1.13  2007/01/17 23:54:00  jim
Plugged some memory leaks

Revision 1.12  2006/12/18 17:14:17  jim
Another tidy of some error messages

Revision 1.11  2006/12/18 16:41:47  jim
Added bit to check for the existence of the index file rather than asking cpl
to do that check with loading the table

Revision 1.10  2006/12/18 12:51:20  jim
Tightened up some of the error reporting

Revision 1.9  2006/12/15 12:03:27  jim
Fixed bug in 2mass_extract where index was offset by -1

Revision 1.8  2006/11/27 12:07:22  jim
Argument list now contains catname. Modified caching routines to take this
into account.

Revision 1.7  2006/08/11 12:45:41  jim
Modified to use wcslib

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

Revision 1.5  2006/06/08 14:50:09  jim
Initialised output table to NULL and fixed a problem in check_cache that
resets the cpl error code

Revision 1.4  2006/03/23 21:18:47  jim
Minor changes mainly to comment headers

Revision 1.3  2006/03/22 13:58:31  jim
Cosmetic fixes to keep lint happy

Revision 1.2  2006/01/23 16:05:33  jim
tidied up error handling

Revision 1.1  2006/01/23 10:31:56  jim
New file


*/

---