fors_data.c

00001 /* $Id: fors_data.c,v 1.32 2008/04/01 10:45:29 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: 2008/04/01 10:45:29 $
00024  * $Revision: 1.32 $
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_std_star.h>
00037 #include <fors_utils.h>
00038 
00039 #include <math.h>
00040 #include <stdbool.h>
00041 #include <string.h>
00042 
00043 struct filtername FORS_DATA_FILTERLIST[] = {
00044     /* FORS1 */
00045     {"U_BESS", FILTER_U},
00046     {"u_HIGH", FILTER_U},
00047     {"B_BESS", FILTER_B},
00048     {"b_HIGH", FILTER_B},
00049     {"g_HIGH", FILTER_G},
00050     {"V_BESS", FILTER_V},
00051     {"v_HIGH", FILTER_V},
00052     {"R_BESS", FILTER_R},
00053     {"I_BESS", FILTER_I},
00054 
00055     /* FORS2 */
00056     {"U_SPECIAL", FILTER_U},
00057     {"B_BESS"   , FILTER_B},
00058     {"V_BESS"   , FILTER_V},
00059     {"R_SPECIAL", FILTER_R},
00060     {"I_BESS"   , FILTER_I}
00061 };
00062 
00063 const char *const FORS_DATA_STD_MAG[FORS_NUM_FILTER] =
00064 {"U",
00065  "B",
00066  //"G",
00067  "V",  /* G uses V */
00068  "V",
00069  "R",
00070  "I",
00071  "Z"};
00072 
00073 const char *const FORS_DATA_STD_DMAG[FORS_NUM_FILTER] =
00074 {"ERR_U",
00075  "ERR_B",
00076  //"ERR_G",
00077  "ERR_V", /* G uses V */
00078  "ERR_V",
00079  "ERR_R",
00080  "ERR_I",
00081  "ERR_Z"};
00082 
00083 const char *const FORS_DATA_STD_COL[FORS_NUM_FILTER] = 
00084 {"U_B",
00085  "B_V",
00086  "B_V",
00087  "B_V",
00088  "V_R",
00089  "V_R",
00090  "?Z?"};
00091 
00092 const char *const FORS_DATA_STD_RA   = "RA";
00093 const char *const FORS_DATA_STD_DEC  = "DEC";
00094 const char *const FORS_DATA_STD_NAME = "OBJECT";
00095 
00096 
00097 const char *const FORS_DATA_PHOT_FILTER    = "FILTER";
00098 const char *const FORS_DATA_PHOT_EXTCOEFF  = "EXT";
00099 const char *const FORS_DATA_PHOT_DEXTCOEFF  = "DEXT";
00100 const char *const FORS_DATA_PHOT_ZEROPOINT = "ZPOINT";
00101 const char *const FORS_DATA_PHOT_DZEROPOINT = "DZPOINT";
00102 const char *const FORS_DATA_PHOT_COLORTERM = "COL";
00103 const char *const FORS_DATA_PHOT_DCOLORTERM = "DCOL";
00104 
00105 
00106 #undef cleanup
00107 #define cleanup \
00108 do { \
00109     cpl_propertylist_delete(raw_header); \
00110     cpl_wcs_delete(wcs); \
00111     cpl_matrix_delete(from); \
00112     cpl_matrix_delete(to); \
00113     cpl_array_delete(status); \
00114 } while(0)
00115 
00120 static void
00121 apply_wcs(fors_std_star_list *stars, const cpl_frame *raw_frame)
00122 {
00123     cpl_propertylist *raw_header = NULL;
00124     cpl_wcs *wcs = NULL;
00125 
00126     cpl_matrix *from = NULL;
00127     cpl_matrix *to = NULL;
00128     cpl_array *status = NULL;
00129 
00130     assure( cpl_frame_get_filename(raw_frame) != NULL, return, NULL );
00131 
00132     if (fors_std_star_list_size(stars) == 0) {
00133         cleanup;
00134         return;
00135     }
00136 
00137     raw_header = cpl_propertylist_load(cpl_frame_get_filename(raw_frame), 0);
00138     
00139     assure( !cpl_error_get_code(), return,
00140             "Failed to load %s header",
00141             cpl_frame_get_filename(raw_frame));
00142     
00143     wcs = cpl_wcs_new_from_propertylist(raw_header);
00144     assure( !cpl_error_get_code(), return,
00145             "Failed to get WCS from %s header",
00146             cpl_frame_get_filename(raw_frame));
00147     
00148     {
00149         fors_std_star *star;
00150         int i;
00151 
00152         from = cpl_matrix_new(fors_std_star_list_size(stars), 2);
00153         
00154         for (star = fors_std_star_list_first(stars), i = 0;
00155              star != NULL;
00156              star = fors_std_star_list_next(stars), i++) {
00157             
00158             cpl_matrix_set(from, i, 0, star->ra);
00159             cpl_matrix_set(from, i, 1, star->dec);
00160         }
00161         
00162         cpl_wcs_convert(wcs, from,
00163                              &to, &status,
00164                              CPL_WCS_WORLD2PHYS);
00165 
00166         /* The WCSLIB call sometimes fails with the error message
00167            "WCSLIB One or more input coordinates invalid",
00168            while all status flags are 0.
00169         */
00170         if (cpl_error_get_code() == CPL_ERROR_UNSPECIFIED) {
00171             cpl_msg_warning(cpl_func,
00172                             "Ignoring WCSLIB unspecified error");
00173             cpl_error_reset();
00174         }
00175 
00176         assure( !cpl_error_get_code(), return,
00177                 "Failed to convert from world to physical coordinates");
00178 
00179         assure( cpl_matrix_get_ncol(to) == 2,
00180                 return, "%d columns, 2 expected",
00181                 cpl_matrix_get_ncol(to));
00182 
00183         assure( cpl_matrix_get_nrow(to) == fors_std_star_list_size(stars),
00184                 return, "%d rows, %d expected",
00185                 cpl_matrix_get_nrow(to), fors_std_star_list_size(stars));
00186         
00187         assure( status != NULL, return, NULL );
00188 
00189         assure( cpl_array_get_size(status) == fors_std_star_list_size(stars),
00190                 return, "Status array size is %d, %d expected",
00191                 cpl_array_get_size(status), fors_std_star_list_size(stars));
00192         
00193         for (star = fors_std_star_list_first(stars), i = 0;
00194              star != NULL;
00195              star = fors_std_star_list_next(stars), i++) {
00196             
00197             if (cpl_array_get_int(status, i, NULL) != 0) {
00198                 cpl_msg_warning(cpl_func, "Catalogue star %d: "
00199                                 "non-zero status = %d from WCSLIB",
00200                                 i, cpl_array_get_int(status, i, NULL));
00201             }
00202             star->pixel->x = cpl_matrix_get(to, i, 0);
00203             star->pixel->y = cpl_matrix_get(to, i, 1);
00204         }
00205     }
00206 
00207 #if 0 /* pre WCSLIB code */
00208     double deg_pr_pixel = 0.1/3600;
00209     fors_std_star *star;
00210 
00211     cpl_msg_warning(cpl_func,
00212                     "WCSLIB not available, "
00213                     "applying fake transformation");
00214     
00215     for (star = fors_std_star_list_first(stars);
00216          star != NULL;
00217          star = fors_std_star_list_next(stars)) {
00218         
00219         star->pixel->x = (star->ra  / deg_pr_pixel - 293000)/10;
00220         star->pixel->y = (star->dec / deg_pr_pixel / 10 + 169800)/10;
00221     }
00222 #endif
00223     
00224     cleanup;
00225     return;
00226 }
00227 
00228 
00229 #undef cleanup
00230 #define cleanup \
00231 do { \
00232     cpl_table_delete(t); \
00233 } while (0)
00234 
00247 void fors_phot_table_load(const cpl_frame *phot_table_frame,
00248                           const fors_setting *setting,
00249               double *color_term,
00250               double *dcolor_term,
00251               double *ext_coeff,
00252               double *dext_coeff,
00253                           double *expected_zeropoint,
00254                           double *dexpected_zeropoint)
00255 {
00256     cpl_table *t = NULL;
00257     
00258     assure( setting != NULL, return, NULL );
00259 /* %%%
00260     assure( setting->filter_name != NULL, return, "No filter name provided");
00261 */
00262     assure( phot_table_frame != NULL, return, NULL );
00263     /* color_term may be NULL */
00264     assure( ext_coeff  != NULL, return, NULL );
00265     /* expected_zeropoint may be NULL */
00266 
00267     assure( (color_term == NULL) == (dcolor_term == NULL), return, NULL );
00268     assure( (ext_coeff == NULL) == (dext_coeff == NULL), return, NULL );
00269     assure( (expected_zeropoint == NULL) == (dexpected_zeropoint == NULL), 
00270         return, NULL );
00271 
00272     assure( cpl_frame_get_filename(phot_table_frame) != NULL, return, NULL );
00273 
00274     if (setting->filter_name == NULL) {
00275         cpl_msg_warning(cpl_func, "Zeropoint computation is not supported "
00276                         "for non-standard filters");
00277         *color_term = 0.0;
00278         *dcolor_term = 0.0;
00279         *ext_coeff = 0.0;
00280         *dext_coeff = 0.0;
00281         *expected_zeropoint = 0.0;
00282         *dexpected_zeropoint = 0.0;
00283         return;
00284     }
00285 
00286     /* Verify instrument setting */
00287     //fixme: for user-friendliness we should verify the setting, except the filter
00288     //        (because this table contains data for all filters)
00289     if (0) {
00290     fors_setting_verify(setting, phot_table_frame, NULL);
00291     assure( !cpl_error_get_code(), return, 
00292             "Could not verify %s setting", 
00293             cpl_frame_get_filename(phot_table_frame));
00294     }
00295 
00296     t = cpl_table_load(cpl_frame_get_filename(phot_table_frame), 1, 1);
00297     
00298     assure( !cpl_error_get_code(), return, "Could not load %s",
00299             cpl_frame_get_filename(phot_table_frame));
00300 
00301     assure( cpl_table_get_nrow(t) > 0, return,
00302             "Empty table %s", cpl_frame_get_filename(phot_table_frame));
00303 
00304     {
00305         const char *const colname[] = {
00306             FORS_DATA_PHOT_FILTER,
00307             FORS_DATA_PHOT_EXTCOEFF,
00308             FORS_DATA_PHOT_DEXTCOEFF,
00309             FORS_DATA_PHOT_ZEROPOINT,
00310             FORS_DATA_PHOT_DZEROPOINT,
00311             FORS_DATA_PHOT_COLORTERM,
00312             FORS_DATA_PHOT_DCOLORTERM};
00313         
00314         const cpl_type coltype[] = {CPL_TYPE_STRING,
00315                                     CPL_TYPE_DOUBLE,
00316                                     CPL_TYPE_DOUBLE,
00317                                     CPL_TYPE_DOUBLE,
00318                                     CPL_TYPE_DOUBLE,
00319                                     CPL_TYPE_DOUBLE,
00320                                     CPL_TYPE_DOUBLE};
00321 
00322         unsigned i;
00323         for (i = 0; i < sizeof(colname) / sizeof(*colname); i++) {
00324             
00325             assure( cpl_table_has_column(t, colname[i]), return,
00326                     "%s: Missing column %s",
00327                     cpl_frame_get_filename(phot_table_frame), colname[i]);
00328             
00329             assure( cpl_table_get_column_type(t, colname[i]) == coltype[i], return,
00330                     "%s column %s type is %s, %s expected",
00331                     cpl_frame_get_filename(phot_table_frame),
00332                     colname[i],
00333                     fors_type_get_string(cpl_table_get_column_type(t, colname[i])),
00334                     fors_type_get_string(coltype[i]));
00335             
00336             assure( cpl_table_count_invalid(t, colname[i]) == 0, return,
00337                     "%s column %s has invalid values",
00338                     cpl_frame_get_filename(phot_table_frame),
00339                     colname[i]);
00340         }
00341     }
00342     /* Input validation done */
00343     
00344     cpl_msg_debug(cpl_func, "Searching for filter: %s", setting->filter_name);
00345             
00346     bool found = false;
00347     int i;
00348     for (i = 0; i < cpl_table_get_nrow(t) && !found; i++) {
00349         const char *phot_filter = cpl_table_get_string(t, FORS_DATA_PHOT_FILTER, i);
00350 
00351         assure( phot_filter != NULL, return, "%s, row %d: Null %s",
00352                 cpl_frame_get_filename(phot_table_frame), i+1, FORS_DATA_PHOT_FILTER);
00353         
00354         if (strcmp(setting->filter_name, phot_filter) == 0) {
00355 
00356             found = true;
00357             
00358             if (color_term != NULL) {
00359                 *color_term  = cpl_table_get_double(t, FORS_DATA_PHOT_COLORTERM, i, NULL);
00360                 *dcolor_term = cpl_table_get_double(t, FORS_DATA_PHOT_DCOLORTERM, i, NULL);
00361             }
00362             
00363             *ext_coeff  = cpl_table_get_double(t, FORS_DATA_PHOT_EXTCOEFF, i, NULL);
00364             *dext_coeff = cpl_table_get_double(t, FORS_DATA_PHOT_DEXTCOEFF, i, NULL);
00365             if (expected_zeropoint != NULL) {
00366                 *expected_zeropoint =
00367                     cpl_table_get_double(t, FORS_DATA_PHOT_ZEROPOINT, i, NULL);
00368                 *dexpected_zeropoint =
00369                     cpl_table_get_double(t, FORS_DATA_PHOT_DZEROPOINT, i, NULL);
00370             }
00371         }
00372     }
00373     
00374     assure( found, return, "%s: Missing entry for filter '%s'",
00375             cpl_frame_get_filename(phot_table_frame), setting->filter_name);
00376     
00377     cleanup;
00378     return;
00379 }
00380 
00381 
00382 #undef cleanup
00383 #define cleanup \
00384 do { \
00385     cpl_table_delete(t); \
00386     cpl_free((void *)ERR_B1); \
00387     cpl_free((void *)ERR_C1); \
00388     cpl_free((void *)ERR_C2); \
00389 } while(0)
00390 
00403 fors_std_star_list *
00404 fors_std_cat_load(const cpl_frameset *cat_frames,
00405           const cpl_frame *raw_frame,
00406           const fors_setting *setting,
00407                   double color_term, double dcolor_term)
00408 {
00409     fors_std_star_list *c = NULL;
00410     cpl_table *t = NULL;
00411     const char *filename;
00412     const char *ERR_B1 = NULL;
00413     const char *ERR_C1 = NULL;
00414     const char *ERR_C2 = NULL;
00415 
00416     assure( cat_frames != NULL, return c, NULL );
00417     assure( setting != NULL, return c, NULL );
00418 
00419     assure( cpl_frame_get_filename(raw_frame) != NULL, return c, NULL );
00420 
00421     /* For each input table
00422            if it contains the required column, then load 
00423     */
00424     c = fors_std_star_list_new();
00425 
00426     const cpl_frame *cat_frame;
00427 
00428     for (cat_frame = cpl_frameset_get_first_const(cat_frames);
00429          cat_frame != NULL;
00430          cat_frame = cpl_frameset_get_next_const(cat_frames)) {
00431         
00432         
00433         filename = cpl_frame_get_filename(cat_frame);
00434         assure( filename != NULL, return c, NULL );
00435 
00436         cpl_table_delete(t);
00437         t = cpl_table_load(filename, 1, 1);
00438         assure( !cpl_error_get_code(), return c, "Could not load FITS catalogue %s",
00439                 filename);
00440         
00441         assure( cpl_table_has_column(t, FORS_DATA_STD_RA), return c, 
00442                 "%s: Missing column %s", filename, FORS_DATA_STD_RA);
00443         
00444         assure( cpl_table_get_column_type(t, FORS_DATA_STD_RA) == CPL_TYPE_DOUBLE, 
00445                 return c,
00446                 "%s: Column %s type is %s, double expected", filename, FORS_DATA_STD_RA,
00447                 fors_type_get_string(cpl_table_get_column_type(t, FORS_DATA_STD_RA)));
00448         
00449         
00450         assure( cpl_table_has_column(t, FORS_DATA_STD_DEC), return c, 
00451                 "%s: Missing column %s", filename, FORS_DATA_STD_DEC);
00452         
00453         assure( cpl_table_get_column_type(t, FORS_DATA_STD_DEC) == CPL_TYPE_DOUBLE, 
00454                 return c, 
00455                 "%s: Column %s type is %s, double expected", filename, FORS_DATA_STD_DEC,
00456                 fors_type_get_string(cpl_table_get_column_type(t, FORS_DATA_STD_DEC)));
00457 
00458 
00459         assure( cpl_table_has_column(t, FORS_DATA_STD_NAME), return c, 
00460                 "%s: Missing column %s", filename, FORS_DATA_STD_NAME);
00461         
00462         assure( cpl_table_get_column_type(t, FORS_DATA_STD_NAME) == CPL_TYPE_STRING, 
00463                 return c,
00464                 "%s: Column %s type is %s, string expected", filename, 
00465                 FORS_DATA_STD_NAME,
00466                 fors_type_get_string(cpl_table_get_column_type(t, FORS_DATA_STD_NAME)));
00467         
00468 
00469         const char *B1 = FORS_DATA_STD_MAG[setting->filter];
00470         assure( B1 != NULL, return c, NULL );
00471         
00472         if ( cpl_table_has_column(t, FORS_DATA_STD_MAG[setting->filter]) ) {
00473             
00474             /* 
00475                 The error propagation depends on which
00476                 error bars are available, which are assumed to be
00477                 uncorrelated
00478 
00479                 Pseudo code:
00480 
00481                 given band B1 and color term C1-C2
00482 
00483                 If have ERR_B1 and ERR_C1 and ERR_C2
00484                     // Stetson like
00485                     B1   + c(C1 - B1)      =  (1-c)B1  + c C1
00486                    (1-c)^2errB1^2 + c^2 errC1^2
00487                     +errc^2 * (C1-B1)^2
00488 
00489                     B1   + c(B1 - C2)      =  (1+c)B1  - c C2     
00490 
00491                    (1+c)^2errB1^2 + c^2 errC2^2
00492                     +errc^2 (B1-C2)^2
00493 
00494                     B1   + c(C1 - C2)
00495                     errB1^2 + c^2 errC1^2 + c^2 errC2^2
00496                     +errc^2 (C1-C2)^2
00497 
00498                     Special case:
00499             G = V + 0.56*(B-V) - 0.12  (Fukugita et al. 1996, AJ 111, p1748)
00500 
00501             magG   = G + c*(B-V) = V + (0.56+c)*(B-V) - 0.12
00502                            = (1 - 0.56 - c)*V  + (0.56+c)*B - 0.12
00503 
00504                     errG^2 = (1 - 0.56 - c)^2*errV^2  + (0.56+c)^2*errB^2
00505                              +errc^2 (B-V)^2
00506                 else
00507                     // Landolt like
00508                     magU  = V      +    (B-V)   + (U-B) + c(U-B)
00509                     err^2 = errV^2 + err(B-V)^2 + (1+c)^2 err(U-B)^2
00510                           + errc^2 (U-B)^2
00511 
00512                     magB  = V      +    (B-V) + c(B-V)
00513                     err^2 = errV^2 + (1+c)^2 err(B-V)^2
00514                           + errc^2 (B-V)^2
00515 
00516                     magG  = V + (0.56+c)*(B-V) - 0.12  (Fukugita et al. 1996, AJ 111, p1748)
00517                     err^2 = errV^2 + (0.56+c)^2*err(B-V)^2
00518                           + errc^2 (B-V)^2
00519 
00520                     magV  = V      + c(B-V)
00521                     err^2 = errV^2 + c^2 err(B-V)^2
00522                           + errc^2 (B-V)^2
00523 
00524                     magR  = V - (V-R) + c (V-R) = V + (c-1)(V-R)
00525                     err^2 = errV^2 + (c-1)^2 err(V-R)^2
00526                           + errc^2 (V-R)^2
00527 
00528                     magI  = V      -    (V-I)   + c      (V-R)
00529                     err^2 = errV^2 + err(V-I)^2 + c^2 err(V-R)^2
00530                           + errc^2 (V-R)^2
00531         */
00532 
00533             /* Find other bands, C1, C2 */
00534             const char *col =  FORS_DATA_STD_COL[setting->filter];
00535             double coeff = -color_term;  /* per convention */
00536         double dcoeff = dcolor_term;
00537 
00538             assure( strlen(col) == strlen("X_Y"), return c,
00539                     "Color term column must have format 'X_Y', is '%s'",
00540                     col );
00541             
00542             assure( col[1] == '_', return c,
00543                     "Color term column must have format 'X_Y', is '%s'",
00544                     col );
00545             
00546             char C1[2] = {'\0', '\0'};
00547             char C2[2] = {'\0', '\0'};
00548 
00549             *C1 = col[0];
00550             *C2 = col[2];
00551 
00552             ERR_B1 = cpl_sprintf("ERR_%s", B1);
00553             ERR_C1 = cpl_sprintf("ERR_%s", C1);
00554             ERR_C2 = cpl_sprintf("ERR_%s", C2);
00555 
00556         /* Catalog data */
00557         double cat_mag, dcat_mag, color;
00558             
00559         /* Color corrected magnitude */
00560             double mag, dmag;
00561 
00562             int i;
00563             for (i = 0; i < cpl_table_get_nrow(t); i++) {
00564                 if (cpl_table_has_column(t, ERR_B1) &&
00565                     cpl_table_has_column(t, ERR_C1) &&
00566                     cpl_table_has_column(t, ERR_C2)) {
00567                     
00568                     /* Stetson, four cases */
00569             if (setting->filter == FILTER_G) {
00570             if (cpl_table_is_valid(t, "V", i) &&
00571                             cpl_table_is_valid(t, "B", i) &&
00572                             cpl_table_is_valid(t, "ERR_B", i) &&
00573                             cpl_table_is_valid(t, "ERR_V", i)) {
00574                 double     v = cpl_table_get_float(t, "V", i, NULL);
00575                 double     b = cpl_table_get_float(t, "B", i, NULL);
00576                             double err_b = cpl_table_get_float(t, "ERR_B", i, NULL);
00577                             double err_v = cpl_table_get_float(t, "ERR_V", i, NULL);
00578 
00579                 color = b-v;
00580                             mag     = (1-0.56-coeff) * v + (0.56+coeff)*b - 0.12;
00581                 cat_mag = (1-0.56      ) * v + (0.56      )*b - 0.12;
00582 
00583                             dmag = 
00584                 sqrt(
00585                                     (1-0.56-coeff)*(1-0.56-coeff)*err_v*err_v +
00586                     (0.56  +coeff)*(0.56  +coeff)*err_b*err_b +
00587                     +
00588                     dcoeff*dcoeff*color*color);
00589 
00590                 dcat_mag = sqrt(
00591                 (1-0.56-0)*(1-0.56-0)*err_v*err_v +
00592                 (0.56  +0)*(0.56  +0)*err_b*err_b);
00593                         }
00594             }
00595             else if (*B1 == *C2) {
00596                         if (cpl_table_is_valid(t, B1, i) &&
00597                             cpl_table_is_valid(t, C1, i) &&
00598                             cpl_table_is_valid(t, ERR_B1, i) &&
00599                             cpl_table_is_valid(t, ERR_C1, i)) {
00600                 double     b1 = cpl_table_get_float(t, B1, i, NULL);
00601                 double     c1 = cpl_table_get_float(t, C1, i, NULL);
00602                             double err_b1 = cpl_table_get_float(t, ERR_B1, i, NULL);
00603                             double err_c1 = cpl_table_get_float(t, ERR_C1, i, NULL);
00604 
00605                 color = c1-b1;
00606 
00607                             cat_mag = b1;
00608                 dcat_mag = err_b1;
00609 
00610                             mag = 
00611                                 (1-coeff) * b1
00612                                 +  coeff  * c1;
00613                 dmag = 
00614                                 sqrt(
00615                                     (1-coeff)*(1-coeff)*err_b1*err_b1 +
00616                                            coeff*coeff *err_c1*err_c1
00617                     +
00618                     dcoeff*dcoeff*color*color);
00619             }
00620                         else continue;
00621                     }
00622                     else if (*B1 == *C1) {
00623                         if (cpl_table_is_valid(t, B1, i) &&
00624                             cpl_table_is_valid(t, C2, i) &&
00625                             cpl_table_is_valid(t, ERR_B1, i) &&
00626                             cpl_table_is_valid(t, ERR_C2, i)) {
00627                 double     b1 = cpl_table_get_float(t, B1, i, NULL);
00628                 double     c2 = cpl_table_get_float(t, C2, i, NULL);
00629                             double err_b1 = cpl_table_get_float(t, ERR_B1, i, NULL);
00630                             double err_c2 = cpl_table_get_float(t, ERR_C2, i, NULL);
00631 
00632                 color = b1-c2;
00633                             
00634                 cat_mag = b1;
00635                 dcat_mag = err_b1;
00636                 
00637                             mag = 
00638                                 (1+coeff) * b1
00639                                 -  coeff  * c2;
00640                             dmag = 
00641                                 sqrt(
00642                                     (1+coeff)*(1+coeff)*err_b1*err_b1 +
00643                     coeff*coeff *err_c2*err_c2
00644                     +
00645                     dcoeff*dcoeff*color*color);
00646 
00647                         }
00648                         else continue;
00649                     }
00650                     else {
00651                         /* All different */
00652                         if (cpl_table_is_valid(t, B1, i) &&
00653                             cpl_table_is_valid(t, C1, i) &&
00654                             cpl_table_is_valid(t, C2, i) &&
00655                             cpl_table_is_valid(t, ERR_B1, i) &&
00656                             cpl_table_is_valid(t, ERR_C1, i) &&
00657                             cpl_table_is_valid(t, ERR_C2, i)) {
00658                 double     b1 = cpl_table_get_float(t, B1, i, NULL);
00659                 double     c1 = cpl_table_get_float(t, C1, i, NULL);
00660                 double     c2 = cpl_table_get_float(t, C2, i, NULL);
00661                             double err_b1 = cpl_table_get_float(t, ERR_B1, i, NULL);
00662                             double err_c1 = cpl_table_get_float(t, ERR_C1, i, NULL);
00663                             double err_c2 = cpl_table_get_float(t, ERR_C2, i, NULL);
00664 
00665                 color = c1-c2;
00666 
00667                 cat_mag = b1;
00668                 dcat_mag = err_b1;
00669                 
00670                             mag = b1
00671                                 + coeff * c1
00672                                 - coeff * c2;
00673                             dmag = sqrt(
00674                                 err_b1*err_b1 +
00675                                 coeff*coeff * err_c1*err_c1 +
00676                                 coeff*coeff * err_c2*err_c2 +
00677                 dcoeff*dcoeff * color*color);
00678                         }
00679                         else continue;
00680                     }
00681                 } /* If stetson, else */
00682         else if (setting->filter == FILTER_G) {
00683             if (cpl_table_is_valid(t, "V", i) &&
00684             cpl_table_is_valid(t, "B_V", i) &&
00685             cpl_table_is_valid(t, "ERR_V", i) &&
00686             cpl_table_is_valid(t, "ERR_B_V", i)) {
00687             double     v   = cpl_table_get_float(t, "V", i, NULL);
00688             double    bv   = cpl_table_get_float(t, "B_V", i, NULL);
00689             double err_v   = cpl_table_get_float(t, "ERR_V", i, NULL);
00690             double err_b_v = cpl_table_get_float(t, "ERR_B_V", i, NULL);
00691 
00692             color = bv;
00693 
00694             cat_mag = v + (0.56)*(bv) - 0.12;
00695 
00696             mag  = v + (0.56+coeff)*(bv) - 0.12;
00697             dmag = 
00698                 sqrt(
00699                 err_v*err_v +
00700                 (0.56+coeff)*(0.56+coeff)*err_b_v*err_b_v +
00701                 +
00702                 dcoeff*dcoeff*bv*bv);
00703             dcat_mag = 
00704                 sqrt(
00705                 err_v*err_v +
00706                 (0.56)*(0.56)*err_b_v*err_b_v);
00707             }
00708         }
00709                 else switch (*B1) {
00710                     /* Landolt, every band is a special case */
00711                 case 'U':
00712                     if (cpl_table_is_valid(t, "U", i) &&
00713                         cpl_table_is_valid(t, "U_B", i) &&
00714                         cpl_table_is_valid(t, "ERR_V", i) &&
00715                         cpl_table_is_valid(t, "ERR_B_V", i) &&
00716                         cpl_table_is_valid(t, "ERR_U_B", i)) {
00717                         
00718                         double err_v  = cpl_table_get_float(t, "ERR_V", i, NULL);
00719                         double err_bv = cpl_table_get_float(t, "ERR_B_V", i, NULL);
00720                         double err_ub = cpl_table_get_float(t, "ERR_U_B", i, NULL);
00721                         double ub     = cpl_table_get_float(t, "U_B", i, NULL);
00722 
00723             color = ub;
00724 
00725             cat_mag = cpl_table_get_float(t, "U", i, NULL);
00726             
00727                         mag = cat_mag + coeff * ub;
00728                         dmag = sqrt(err_v*err_v + err_bv*err_bv+
00729                                     (1+coeff)*(1+coeff)*err_ub*err_ub +
00730                     dcoeff*dcoeff*ub*ub);
00731             
00732                         dcat_mag = sqrt(err_v*err_v + err_bv*err_bv+
00733                     err_ub*err_ub);
00734                     } 
00735                     else continue; /* to next for(...) iteration */
00736                     break; /* out of switch */
00737                 case 'B':
00738                     if (cpl_table_is_valid(t, "B", i) &&
00739                         cpl_table_is_valid(t, "B_V", i) &&
00740                         cpl_table_is_valid(t, "ERR_V", i) &&
00741                         cpl_table_is_valid(t, "ERR_B_V", i)) {
00742                         
00743                         double err_v  = cpl_table_get_float(t, "ERR_V", i, NULL);
00744                         double err_bv = cpl_table_get_float(t, "ERR_B_V", i, NULL);
00745             double     bv = cpl_table_get_float(t, "B_V", i, NULL);
00746 
00747             color = bv;
00748             cat_mag = cpl_table_get_float(t, "B", i, NULL);
00749 
00750                         mag = cat_mag +
00751                             coeff * bv;
00752                         dmag = sqrt(err_v*err_v +
00753                                     (1+coeff)*(1+coeff)*err_bv*err_bv +
00754                     dcoeff*dcoeff*bv*bv);
00755             dcat_mag = sqrt(err_v*err_v +
00756                     err_bv*err_bv);
00757                     }
00758                     else continue;
00759                     break;
00760                 case 'V':
00761                     if (cpl_table_is_valid(t, "V", i) &&
00762                         cpl_table_is_valid(t, "B_V", i) &&
00763                         cpl_table_is_valid(t, "ERR_V", i) &&
00764                         cpl_table_is_valid(t, "ERR_B_V", i)) {
00765                         
00766                         double err_v  = cpl_table_get_float(t, "ERR_V", i, NULL);
00767                         double err_bv = cpl_table_get_float(t, "ERR_B_V", i, NULL);
00768             double     bv = cpl_table_get_float(t, "B_V", i, NULL);
00769 
00770             color = bv;
00771             cat_mag = cpl_table_get_float(t, "V", i, NULL);
00772                         dcat_mag = err_v;
00773             
00774                         mag = cat_mag + coeff * bv;
00775                         dmag = sqrt(err_v*err_v +
00776                                     coeff*coeff*err_bv*err_bv +
00777                     dcoeff*dcoeff*bv*bv);
00778                     }
00779                     else continue;
00780                     break;
00781                 case 'R':
00782                     if (cpl_table_is_valid(t, "R", i) &&
00783                         cpl_table_is_valid(t, "V_R", i) &&
00784                         cpl_table_is_valid(t, "ERR_V", i) &&
00785                         cpl_table_is_valid(t, "ERR_V_R", i)) {
00786                         
00787                         double err_v  = cpl_table_get_float(t, "ERR_V", i, NULL);
00788                         double err_vr = cpl_table_get_float(t, "ERR_V_R", i, NULL);
00789             double     vr = cpl_table_get_float(t, "V_R", i, NULL);
00790 
00791             color = vr;
00792             cat_mag = cpl_table_get_float(t, "R", i, NULL);
00793                         mag = cat_mag + coeff * vr;
00794                         dmag = sqrt(err_v*err_v +
00795                                     (1-coeff)*(1-coeff)*err_vr*err_vr +
00796                     dcoeff*dcoeff*vr*vr);
00797                         dcat_mag = sqrt(err_v*err_v + err_vr*err_vr);
00798                     }
00799                     else continue;
00800                     break;
00801                 case 'I':
00802                     if (cpl_table_is_valid(t, "I", i) &&
00803                         cpl_table_is_valid(t, "V_R", i) &&
00804                         cpl_table_is_valid(t, "ERR_V", i) &&
00805                         cpl_table_is_valid(t, "ERR_V_I", i) &&
00806                         cpl_table_is_valid(t, "ERR_V_R", i)) {
00807                         
00808                         double err_v  = cpl_table_get_float(t, "ERR_V", i, NULL);
00809                         double err_vi = cpl_table_get_float(t, "ERR_V_I", i, NULL);
00810                         double err_vr = cpl_table_get_float(t, "ERR_V_R", i, NULL);
00811                         double     vr = cpl_table_get_float(t, "V_R", i, NULL);
00812 
00813             color = vr;
00814             cat_mag = cpl_table_get_float(t, "I", i, NULL);
00815                         mag = cat_mag + coeff * vr;
00816                         dmag = sqrt(err_v*err_v + err_vi*err_vi+
00817                                     coeff*coeff*err_vr*err_vr +
00818                     dcoeff*dcoeff*vr*vr);
00819                         dcat_mag = sqrt(err_v*err_v + err_vi*err_vi);
00820                     }
00821                     else continue;
00822                     break;
00823                 default:
00824                     assure( false, return c, 
00825                             "Unknown filter: %s", B1);
00826                     break;
00827                 }
00828                 
00829                 double ra = cpl_table_get_double(t, FORS_DATA_STD_RA,
00830                                                  i, NULL);
00831                 double dec = cpl_table_get_double(t, FORS_DATA_STD_DEC,
00832                                                   i, NULL);
00833 
00834                 const char *std_name = cpl_table_get_string(t, FORS_DATA_STD_NAME,
00835                                                             i);
00836                 
00837                 fors_std_star_list_insert(
00838                     c, fors_std_star_new(ra, dec, mag, dmag, 
00839                      cat_mag, dcat_mag,
00840                      color, std_name));
00841 
00842             } /* for each table row */
00843 
00844 #if 0  /* This is old code which removes doublets but is slow (O(n^2)),
00845         * Doublets will be removed during identification, after
00846         * selecting the (relatively few) stars that fall inside the CCD
00847         */
00848                     double nearest_dist;
00849                     if (fors_std_star_list_size(c) > 0) {
00850                         fors_std_star *nearest = fors_std_star_list_min_val(
00851                             c,
00852                             (fors_std_star_list_func_eval)
00853                             fors_std_star_dist_arcsec,
00854                             std);
00855                         
00856                         nearest_dist = fors_std_star_dist_arcsec(std, nearest);
00857                         cpl_msg_debug(cpl_func, "min dist = %f arcseconds",
00858                                       nearest_dist);
00859                     }
00860                     
00861                     if (fors_std_star_list_size(c) == 0 || nearest_dist > 5) {
00862                         fors_std_star_list_insert(c, std);
00863                     }
00864                     else {
00865                         fors_std_star_delete(&std);
00866                     }
00867 #endif
00868         } /* if has column B1 */
00869     else {
00870             cpl_msg_info(cpl_func, "Skipping catalog %s, no column %s",
00871                          filename, B1);
00872         }
00873     } /* for each frame */
00874     
00875     cpl_msg_info(cpl_func, "Found %d catalogue standards",
00876                  fors_std_star_list_size(c));
00877     
00878     apply_wcs(c, raw_frame);
00879 
00880     /* fors_std_cat_print(c); */
00881 
00882     cleanup;
00883     return c;
00884 }

Generated on Wed Sep 10 07:31:51 2008 for FORS Pipeline Reference Manual by  doxygen 1.4.6