00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
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
00112
00113
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
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
00206
00207 assure( phot_table_frame != NULL, return, NULL );
00208
00209 assure( ext_coeff != NULL, return, NULL );
00210
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
00232
00233
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
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
00335
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
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 }