fors_data.c

00001 /* $Id: fors_data.c,v 1.41 2009/05/14 11:40:51 cizzo Exp $
00002  *
00003  * This file is part of the FORS Library
00004  * Copyright (C) 2002-2006 European Southern Observatory
00005  *
00006  * This program is free software; you can redistribute it and/or modify
00007  * it under the terms of the GNU General Public License as published by
00008  * the Free Software Foundation; either version 2 of the License, or
00009  * (at your option) any later version.
00010  *
00011  * This program is distributed in the hope that it will be useful,
00012  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014  * GNU General Public License for more details.
00015  *
00016  * You should have received a copy of the GNU General Public License
00017  * along with this program; if not, write to the Free Software
00018  * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301 USA
00019  */
00020 
00021 /*
00022  * $Author: cizzo $
00023  * $Date: 2009/05/14 11:40:51 $
00024  * $Revision: 1.41 $
00025  * $Name:  $
00026  */
00027 
00028 #ifdef HAVE_CONFIG_H
00029 #include <config.h>
00030 #endif
00031 
00032 #include <fors_data.h>
00033 
00034 #include <cpl_wcs.h>
00035 #include <fors_star.h>
00036 #include <fors_instrument.h>
00037 #include <fors_std_star.h>
00038 #include <fors_utils.h>
00039 
00040 #include <math.h>
00041 #include <stdbool.h>
00042 #include <string.h>
00043 
00044 
00045 const char *const FORS_DATA_PHOT_FILTER    = "FILTER";
00046 const char *const FORS_DATA_PHOT_EXTCOEFF  = "EXT";
00047 const char *const FORS_DATA_PHOT_DEXTCOEFF  = "DEXT";
00048 const char *const FORS_DATA_PHOT_ZEROPOINT = "ZPOINT";
00049 const char *const FORS_DATA_PHOT_DZEROPOINT = "DZPOINT";
00050 const char *const FORS_DATA_PHOT_COLORTERM = "COL";
00051 const char *const FORS_DATA_PHOT_DCOLORTERM = "DCOL";
00052 
00053 
00054 #undef cleanup
00055 #define cleanup \
00056 do { \
00057     cpl_wcs_delete(wcs); \
00058     cpl_matrix_delete(from); \
00059     cpl_matrix_delete(to); \
00060     cpl_array_delete(status); \
00061 } while(0)
00062 
00067 void
00068 fors_std_star_list_apply_wcs(               fors_std_star_list      *stars,
00069                                             const cpl_propertylist  *header)
00070 {
00071     cpl_wcs *wcs = NULL;
00072 
00073     cpl_matrix *from = NULL;
00074     cpl_matrix *to = NULL;
00075     cpl_array *status = NULL;
00076     
00077     cassure_automsg(                        stars != NULL,
00078                                             CPL_ERROR_NULL_INPUT,
00079                                             return);
00080     cassure_automsg(                        header != NULL,
00081                                             CPL_ERROR_NULL_INPUT,
00082                                             return);
00083 
00084     if (fors_std_star_list_size(stars) == 0) {
00085         cleanup;
00086         return;
00087     }
00088     
00089     wcs = cpl_wcs_new_from_propertylist(header);
00090     assure( !cpl_error_get_code(), return,
00091             "Failed to get WCS from header");
00092     
00093     {
00094         fors_std_star *star;
00095         int i;
00096 
00097         from = cpl_matrix_new(fors_std_star_list_size(stars), 2);
00098         
00099         for (star = fors_std_star_list_first(stars), i = 0;
00100              star != NULL;
00101              star = fors_std_star_list_next(stars), i++) {
00102             
00103             cpl_matrix_set(from, i, 0, star->ra);
00104             cpl_matrix_set(from, i, 1, star->dec);
00105         }
00106         
00107         cpl_wcs_convert(wcs, from,
00108                              &to, &status,
00109                              CPL_WCS_WORLD2PHYS);
00110 
00111         /* The WCSLIB call sometimes fails with the error message
00112            "WCSLIB One or more input coordinates invalid",
00113            while all status flags are 0.
00114         */
00115         if (cpl_error_get_code() == CPL_ERROR_UNSPECIFIED) {
00116             cpl_msg_warning(cpl_func,
00117                             "Ignoring WCSLIB unspecified error");
00118             cpl_error_reset();
00119         }
00120 
00121         assure( !cpl_error_get_code(), return,
00122                 "Failed to convert from world to physical coordinates");
00123 
00124         assure( cpl_matrix_get_ncol(to) == 2,
00125                 return, "%d columns, 2 expected",
00126                 cpl_matrix_get_ncol(to));
00127 
00128         assure( cpl_matrix_get_nrow(to) == fors_std_star_list_size(stars),
00129                 return, "%d rows, %d expected",
00130                 cpl_matrix_get_nrow(to), fors_std_star_list_size(stars));
00131         
00132         assure( status != NULL, return, NULL );
00133 
00134         assure( cpl_array_get_size(status) == fors_std_star_list_size(stars),
00135                 return, "Status array size is %d, %d expected",
00136                 cpl_array_get_size(status), fors_std_star_list_size(stars));
00137         
00138         for (star = fors_std_star_list_first(stars), i = 0;
00139              star != NULL;
00140              star = fors_std_star_list_next(stars), i++) {
00141             
00142             if (cpl_array_get_int(status, i, NULL) != 0) {
00143                 cpl_msg_warning(cpl_func, "Catalogue star %d: "
00144                                 "non-zero status = %d from WCSLIB",
00145                                 i, cpl_array_get_int(status, i, NULL));
00146             }
00147             star->pixel->x = cpl_matrix_get(to, i, 0);
00148             star->pixel->y = cpl_matrix_get(to, i, 1);
00149         }
00150     }
00151 
00152 #if 0 /* pre WCSLIB code */
00153     double deg_pr_pixel = 0.1/3600;
00154     fors_std_star *star;
00155 
00156     cpl_msg_warning(cpl_func,
00157                     "WCSLIB not available, "
00158                     "applying fake transformation");
00159     
00160     for (star = fors_std_star_list_first(stars);
00161          star != NULL;
00162          star = fors_std_star_list_next(stars)) {
00163         
00164         star->pixel->x = (star->ra  / deg_pr_pixel - 293000)/10;
00165         star->pixel->y = (star->dec / deg_pr_pixel / 10 + 169800)/10;
00166     }
00167 #endif
00168     
00169     cleanup;
00170     return;
00171 }
00172 
00173 
00174 #undef cleanup
00175 #define cleanup \
00176 do { \
00177     cpl_table_delete(t); \
00178 } while (0)
00179 
00192 void fors_phot_table_load(const cpl_frame *phot_table_frame,
00193                           const fors_setting *setting,
00194                           double *color_term,
00195                           double *dcolor_term,
00196                           double *ext_coeff,
00197                           double *dext_coeff,
00198                           double *expected_zeropoint,
00199                           double *dexpected_zeropoint)
00200 {
00201     cpl_table *t = NULL;
00202     
00203     assure( setting != NULL, return, NULL );
00204 /* %%%
00205     assure( setting->filter_name != NULL, return, "No filter name provided");
00206 */
00207     assure( phot_table_frame != NULL, return, NULL );
00208     /* color_term may be NULL */
00209     assure( ext_coeff  != NULL, return, NULL );
00210     /* expected_zeropoint may be NULL */
00211 
00212     assure( (color_term == NULL) == (dcolor_term == NULL), return, NULL );
00213     assure( (ext_coeff == NULL) == (dext_coeff == NULL), return, NULL );
00214     assure( (expected_zeropoint == NULL) == (dexpected_zeropoint == NULL), 
00215             return, NULL );
00216 
00217     assure( cpl_frame_get_filename(phot_table_frame) != NULL, return, NULL );
00218 
00219     if (setting->filter_name == NULL) {
00220         cpl_msg_warning(cpl_func, "Zeropoint computation is not supported "
00221                         "for non-standard filters");
00222         *color_term = 0.0;
00223         *dcolor_term = 0.0;
00224         *ext_coeff = 0.0;
00225         *dext_coeff = 0.0;
00226         *expected_zeropoint = 0.0;
00227         *dexpected_zeropoint = 0.0;
00228         return;
00229     }
00230 
00231     /* Verify instrument setting */
00232     //fixme: for user-friendliness we should verify the setting, except the filter
00233     //        (because this table contains data for all filters)
00234     if (0) {
00235     fors_setting_verify(setting, phot_table_frame, NULL);
00236     assure( !cpl_error_get_code(), return, 
00237             "Could not verify %s setting", 
00238             cpl_frame_get_filename(phot_table_frame));
00239     }
00240 
00241     t = cpl_table_load(cpl_frame_get_filename(phot_table_frame), 1, 1);
00242     
00243     assure( !cpl_error_get_code(), return, "Could not load %s",
00244             cpl_frame_get_filename(phot_table_frame));
00245 
00246     assure( cpl_table_get_nrow(t) > 0, return,
00247             "Empty table %s", cpl_frame_get_filename(phot_table_frame));
00248 
00249     {
00250         const char *const colname[] = {
00251             FORS_DATA_PHOT_FILTER,
00252             FORS_DATA_PHOT_EXTCOEFF,
00253             FORS_DATA_PHOT_DEXTCOEFF,
00254             FORS_DATA_PHOT_ZEROPOINT,
00255             FORS_DATA_PHOT_DZEROPOINT,
00256             FORS_DATA_PHOT_COLORTERM,
00257             FORS_DATA_PHOT_DCOLORTERM};
00258         
00259         const cpl_type coltype[] = {CPL_TYPE_STRING,
00260                                     CPL_TYPE_DOUBLE,
00261                                     CPL_TYPE_DOUBLE,
00262                                     CPL_TYPE_DOUBLE,
00263                                     CPL_TYPE_DOUBLE,
00264                                     CPL_TYPE_DOUBLE,
00265                                     CPL_TYPE_DOUBLE};
00266 
00267         unsigned i;
00268         for (i = 0; i < sizeof(colname) / sizeof(*colname); i++) {
00269             
00270             assure( cpl_table_has_column(t, colname[i]), return,
00271                     "%s: Missing column %s",
00272                     cpl_frame_get_filename(phot_table_frame), colname[i]);
00273             
00274             assure( cpl_table_get_column_type(t, colname[i]) == coltype[i], return,
00275                     "%s column %s type is %s, %s expected",
00276                     cpl_frame_get_filename(phot_table_frame),
00277                     colname[i],
00278                     fors_type_get_string(cpl_table_get_column_type(t, colname[i])),
00279                     fors_type_get_string(coltype[i]));
00280             
00281             assure( cpl_table_count_invalid(t, colname[i]) == 0, return,
00282                     "%s column %s has invalid values",
00283                     cpl_frame_get_filename(phot_table_frame),
00284                     colname[i]);
00285         }
00286     }
00287     /* Input validation done */
00288     
00289     cpl_msg_debug(cpl_func, "Searching for filter: %s", setting->filter_name);
00290             
00291     bool found = false;
00292     int i;
00293     for (i = 0; i < cpl_table_get_nrow(t) && !found; i++) {
00294         const char *phot_filter = cpl_table_get_string(t, FORS_DATA_PHOT_FILTER, i);
00295 
00296         assure( phot_filter != NULL, return, "%s, row %d: Null %s",
00297                 cpl_frame_get_filename(phot_table_frame), i+1, FORS_DATA_PHOT_FILTER);
00298         
00299         if (strcmp(setting->filter_name, phot_filter) == 0) {
00300 
00301             found = true;
00302             
00303             if (color_term != NULL) {
00304                 *color_term  = cpl_table_get_double(t, FORS_DATA_PHOT_COLORTERM, i, NULL);
00305                 *dcolor_term = cpl_table_get_double(t, FORS_DATA_PHOT_DCOLORTERM, i, NULL);
00306             }
00307             
00308             *ext_coeff  = cpl_table_get_double(t, FORS_DATA_PHOT_EXTCOEFF, i, NULL);
00309             *dext_coeff = cpl_table_get_double(t, FORS_DATA_PHOT_DEXTCOEFF, i, NULL);
00310             if (expected_zeropoint != NULL) {
00311                 *expected_zeropoint =
00312                     cpl_table_get_double(t, FORS_DATA_PHOT_ZEROPOINT, i, NULL);
00313                 *dexpected_zeropoint =
00314                     cpl_table_get_double(t, FORS_DATA_PHOT_DZEROPOINT, i, NULL);
00315             }
00316         }
00317     }
00318 
00319     if (found == false) {
00320         cpl_msg_warning(cpl_func, "Entry for filter %s missing in input "
00321                         "photometric table (%s): assuming all photometric "
00322                         "coefficients Z, E, and color term, equal to zero.",
00323                         setting->filter_name, 
00324                         cpl_frame_get_filename(phot_table_frame));
00325         *color_term          = 0.0;
00326         *dcolor_term         = 0.0;
00327         *ext_coeff           = 0.0;
00328         *dext_coeff          = 0.0;
00329         *expected_zeropoint  = 0.0;
00330         *dexpected_zeropoint = 0.0;
00331     }
00332     
00333     /*
00334     assure( found, return, "%s: Missing entry for filter '%s'",
00335             cpl_frame_get_filename(phot_table_frame), setting->filter_name);
00336     */
00337     
00338     cleanup;
00339     return;
00340 }
00341 
00342 
00343 #undef cleanup
00344 #define cleanup \
00345 do { \
00346     cpl_table_delete(table); \
00347 } while(0)
00348 
00365 cpl_table *fors_phot_coeff_create(const fors_setting *setting,
00366                                   double color_term,
00367                                   double dcolor_term,
00368                                   double ext_coeff,
00369                                   double dext_coeff,
00370                                   double zeropoint,
00371                                   double dzeropoint)
00372 {
00373     cpl_table *table = cpl_table_new(1);
00374 
00375     if (table == NULL) {
00376         return NULL;
00377     }
00378 
00379     if (dcolor_term > 0.0 || dext_coeff > 0.0 || dzeropoint > 0.0) {
00380         cpl_table_new_column(table, FORS_DATA_PHOT_FILTER, CPL_TYPE_STRING);
00381         cpl_table_set_string(table, FORS_DATA_PHOT_FILTER, 
00382                              0, setting->filter_name);
00383     }
00384     else {
00385         cleanup;
00386         return NULL;
00387     }
00388 
00389     /*
00390      * Create only the necessary columns for the new table
00391      */
00392 
00393     if (dext_coeff > 0.0) {
00394         cpl_table_new_column(table, FORS_DATA_PHOT_EXTCOEFF, CPL_TYPE_DOUBLE);
00395         cpl_table_new_column(table, FORS_DATA_PHOT_DEXTCOEFF, CPL_TYPE_DOUBLE);
00396         cpl_table_set_double(table, FORS_DATA_PHOT_EXTCOEFF, 0, ext_coeff);
00397         cpl_table_set_double(table, FORS_DATA_PHOT_DEXTCOEFF, 0, dext_coeff);
00398     }
00399 
00400     if (dzeropoint > 0.0) {
00401         cpl_table_new_column(table, FORS_DATA_PHOT_ZEROPOINT, CPL_TYPE_DOUBLE);
00402         cpl_table_new_column(table, FORS_DATA_PHOT_DZEROPOINT, CPL_TYPE_DOUBLE);
00403         cpl_table_set_double(table, FORS_DATA_PHOT_ZEROPOINT, 0, zeropoint);
00404         cpl_table_set_double(table, FORS_DATA_PHOT_DZEROPOINT, 0, dzeropoint);
00405     }
00406 
00407     if (dcolor_term > 0.0) {
00408         cpl_table_new_column(table, FORS_DATA_PHOT_COLORTERM, CPL_TYPE_DOUBLE);
00409         cpl_table_new_column(table, FORS_DATA_PHOT_DCOLORTERM, CPL_TYPE_DOUBLE);
00410         cpl_table_set_double(table, FORS_DATA_PHOT_COLORTERM, 0, color_term);
00411         cpl_table_set_double(table, FORS_DATA_PHOT_DCOLORTERM, 0, dcolor_term);
00412     }
00413 
00414     return table;
00415 }

Generated on Wed Jun 24 14:11:17 2009 for FORS Pipeline Reference Manual by  doxygen 1.4.7