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_tools.h>
00033
00034 #include <fors_pfits.h>
00035 #include <fors_utils.h>
00036
00037 #include <cpl.h>
00038 #include <math.h>
00039 #include <stdbool.h>
00040
00041
00045
00046
00049 #undef cleanup
00050 #define cleanup \
00051 do { \
00052 cpl_propertylist_delete(header); \
00053 } while(0)
00054
00063 double
00064 fors_star_ext_corr(fors_star_list *stars,
00065 const fors_setting *setting,
00066 double ext_coeff,
00067 double dext_coeff,
00068 const cpl_frame *raw_frame)
00069 {
00070 cpl_propertylist *header = NULL;
00071
00072
00073 cpl_msg_info(cpl_func, "Extinction correction");
00074
00075 assure( cpl_frame_get_filename(raw_frame) != NULL, return -1, NULL );
00076
00077 header = cpl_propertylist_load(cpl_frame_get_filename(raw_frame), 0);
00078 assure( !cpl_error_get_code(), return -1,
00079 "Failed to load %s primary header",
00080 cpl_frame_get_filename(raw_frame));
00081
00082
00083 double avg_airmass = fors_get_airmass(header);
00084 assure( !cpl_error_get_code(), return -1,
00085 "%s: Could not read airmass",
00086 cpl_frame_get_filename(raw_frame));
00087
00088 cpl_msg_indent_more();
00089 cpl_msg_info(cpl_func, "Exposure time = %f s", setting->exposure_time);
00090 cpl_msg_info(cpl_func, "Gain = %f ADU/e-", setting->average_gain);
00091 cpl_msg_info(cpl_func, "Ext. coeff. = %f +- %f mag/airmass",
00092 ext_coeff, dext_coeff);
00093 cpl_msg_info(cpl_func, "Avg. airmass = %f airmass", avg_airmass);
00094
00095
00096 cpl_msg_indent_less();
00097
00098 {
00099 fors_star *star;
00100
00101 for (star = fors_star_list_first(stars);
00102 star != NULL;
00103 star = fors_star_list_next(stars)) {
00104 star->magnitude_corr = star->magnitude
00105 + 2.5*log(setting->average_gain)/log(10)
00106 + 2.5*log(setting->exposure_time)/log(10)
00107 - ext_coeff * avg_airmass;
00108
00109
00110
00111
00112 star->dmagnitude_corr = sqrt(star->dmagnitude * star->dmagnitude
00113 + dext_coeff*dext_coeff * avg_airmass*avg_airmass);
00114 }
00115 }
00116
00117 cleanup;
00118 return avg_airmass;
00119 }
00120
00128 cpl_table *
00129 fors_create_sources_table(fors_star_list *sources)
00130 {
00131 cpl_table *t = NULL;
00132
00133 t = cpl_table_new(fors_star_list_size(sources));
00134 cpl_table_new_column(t, "X", CPL_TYPE_DOUBLE);
00135 cpl_table_new_column(t, "Y", CPL_TYPE_DOUBLE);
00136 cpl_table_new_column(t, "FWHM", CPL_TYPE_DOUBLE);
00137 cpl_table_new_column(t, "A", CPL_TYPE_DOUBLE);
00138 cpl_table_new_column(t, "B", CPL_TYPE_DOUBLE);
00139 cpl_table_new_column(t, "THETA", CPL_TYPE_DOUBLE);
00140 cpl_table_new_column(t, "ELL", CPL_TYPE_DOUBLE);
00141 cpl_table_new_column(t, "INSTR_MAG", CPL_TYPE_DOUBLE);
00142 cpl_table_new_column(t, "DINSTR_MAG", CPL_TYPE_DOUBLE);
00143 cpl_table_new_column(t, "INSTR_CMAG", CPL_TYPE_DOUBLE);
00144 cpl_table_new_column(t, "DINSTR_CMAG", CPL_TYPE_DOUBLE);
00145 cpl_table_new_column(t, "CLASS_STAR", CPL_TYPE_DOUBLE);
00146
00147 cpl_table_new_column(t, "OBJECT", CPL_TYPE_STRING);
00148 cpl_table_new_column(t, "RA", CPL_TYPE_DOUBLE);
00149 cpl_table_new_column(t, "DEC", CPL_TYPE_DOUBLE);
00150 cpl_table_new_column(t, "MAG", CPL_TYPE_DOUBLE);
00151 cpl_table_new_column(t, "DMAG", CPL_TYPE_DOUBLE);
00152 cpl_table_new_column(t, "CAT_MAG", CPL_TYPE_DOUBLE);
00153 cpl_table_new_column(t, "DCAT_MAG", CPL_TYPE_DOUBLE);
00154 cpl_table_new_column(t, "COLOR", CPL_TYPE_DOUBLE);
00155 cpl_table_new_column(t, "USE_CAT", CPL_TYPE_INT);
00156
00157
00158 cpl_table_new_column(t, "SHIFT_X", CPL_TYPE_DOUBLE);
00159 cpl_table_new_column(t, "SHIFT_Y", CPL_TYPE_DOUBLE);
00160 cpl_table_new_column(t, "ZEROPOINT", CPL_TYPE_DOUBLE);
00161 cpl_table_new_column(t, "DZEROPOINT", CPL_TYPE_DOUBLE);
00162 cpl_table_new_column(t, "WEIGHT", CPL_TYPE_DOUBLE);
00163
00164 {
00165 fors_star *s;
00166 int i;
00167 for (s = fors_star_list_first(sources), i = 0;
00168 s != NULL;
00169 s = fors_star_list_next(sources), i++) {
00170
00171 const fors_std_star *id = s->id;
00172
00173 cpl_table_set_double(t, "X", i, s->pixel->x);
00174 cpl_table_set_double(t, "Y", i, s->pixel->y);
00175 cpl_table_set_double(t, "FWHM", i, s->fwhm);
00176 cpl_table_set_double(t, "A", i, s->semi_major);
00177 cpl_table_set_double(t, "B", i, s->semi_minor);
00178 cpl_table_set_double(t, "THETA", i, s->orientation);
00179 cpl_table_set_double(t, "ELL", i, fors_star_ellipticity(s, NULL));
00180 cpl_table_set_double(t, "INSTR_MAG", i, s->magnitude);
00181 cpl_table_set_double(t, "DINSTR_MAG", i, s->dmagnitude);
00182 cpl_table_set_double(t, "INSTR_CMAG", i, s->magnitude_corr);
00183 cpl_table_set_double(t, "DINSTR_CMAG", i, s->dmagnitude_corr);
00184 cpl_table_set_double(t, "CLASS_STAR", i, s->stellarity_index);
00185 cpl_table_set_double(t, "WEIGHT", i, s->weight);
00186
00187 if (id != NULL) {
00188 cpl_table_set_string(t, "OBJECT" , i, id->name);
00189 cpl_table_set_double(t, "RA" , i, id->ra);
00190 cpl_table_set_double(t, "DEC" , i, id->dec);
00191 cpl_table_set_double(t, "MAG" , i, id->magnitude);
00192 cpl_table_set_double(t, "DMAG" , i, id->dmagnitude);
00193 cpl_table_set_double(t, "CAT_MAG", i, id->cat_magnitude);
00194 cpl_table_set_double(t, "DCAT_MAG", i, id->dcat_magnitude);
00195 cpl_table_set_double(t, "COLOR" , i, id->color);
00196 cpl_table_set_double(t, "SHIFT_X", i, s->pixel->x - id->pixel->x);
00197 cpl_table_set_double(t, "SHIFT_Y", i, s->pixel->y - id->pixel->y);
00198 cpl_table_set_double(t, "ZEROPOINT", i,
00199 fors_star_get_zeropoint(s, NULL));
00200 cpl_table_set_double(t, "DZEROPOINT", i,
00201 fors_star_get_zeropoint_err(s, NULL));
00202 cpl_table_set_int (t, "USE_CAT", i, 1);
00203
00204
00205 }
00206 else {
00207 cpl_table_set_invalid(t, "RA" , i);
00208 cpl_table_set_invalid(t, "DEC", i);
00209 cpl_table_set_invalid(t, "MAG", i);
00210 cpl_table_set_invalid(t, "DMAG", i);
00211 cpl_table_set_invalid(t, "SHIFT_X", i);
00212 cpl_table_set_invalid(t, "SHIFT_Y", i);
00213 cpl_table_set_invalid(t, "ZEROPOINT", i);
00214 cpl_table_set_invalid(t, "DZEROPOINT", i);
00215 }
00216 }
00217 }
00218
00219 return t;
00220 }
00221
00222 #undef cleanup
00223 #define cleanup \
00224 do { \
00225 fors_image_delete(&image); \
00226 fors_image_delete(&image2); \
00227 } while(0)
00228
00235 double
00236 fors_fixed_pattern_noise(const fors_image *master,
00237 double convert_ADU,
00238 double master_noise)
00239 {
00240 double master_fixed_pattern_noise;
00241 fors_image *image = NULL;
00242 fors_image *image2 = NULL;
00243
00244 assure( master != NULL, return -1, NULL );
00245
00246
00247
00248
00249 if (fors_image_get_size_x(master) >= 121 &&
00250 fors_image_get_size_y(master) >= 121) {
00251
00252 int mid_x = (fors_image_get_size_x(master) + 1) / 2;
00253 int mid_y = (fors_image_get_size_y(master) + 1) / 2;
00254
00255 image = fors_image_duplicate(master);
00256 fors_image_crop(image,
00257 mid_x - 50, mid_y - 50,
00258 mid_x + 50, mid_y + 50);
00259
00260 image2 = fors_image_duplicate(master);
00261 fors_image_crop(image2,
00262 mid_x + 10 - 50, mid_y + 10 - 50,
00263 mid_x + 10 + 50, mid_y + 10 + 50);
00264
00265 fors_image_subtract(image, image2);
00266
00267 master_fixed_pattern_noise =
00268 fors_image_get_stdev(image, NULL) / sqrt(2);
00269
00270
00271 master_fixed_pattern_noise *= convert_ADU;
00272
00273
00274 if (master_fixed_pattern_noise >= master_noise) {
00275
00276 master_fixed_pattern_noise = sqrt(master_fixed_pattern_noise*
00277 master_fixed_pattern_noise
00278 -
00279 master_noise*
00280 master_noise);
00281 }
00282 else {
00283 cpl_msg_warning(cpl_func,
00284 "Zero-shift noise (%f ADU) is greater than "
00285 "accumulated zero-shift and fixed pattern noise (%f ADU), "
00286 "setting fixed pattern noise to zero",
00287 master_noise,
00288 master_fixed_pattern_noise);
00289 master_fixed_pattern_noise = 0;
00290 }
00291 }
00292 else {
00293 cpl_msg_warning(cpl_func,
00294 "Master flat too small (%dx%d), "
00295 "need size 121x121 to compute master flat "
00296 "fixed pattern noise",
00297 fors_image_get_size_x(master),
00298 fors_image_get_size_y(master));
00299 master_fixed_pattern_noise = -1;
00300 }
00301
00302 cleanup;
00303 return master_fixed_pattern_noise;
00304 }
00305
00306
00307 #undef cleanup
00308 #define cleanup \
00309 do { \
00310 fors_image_delete(&image); \
00311 fors_image_delete(&image2); \
00312 } while(0)
00313
00320 double
00321 fors_fixed_pattern_noise_bias(const fors_image *first_raw,
00322 const fors_image *second_raw,
00323 double ron)
00324 {
00325 double bias_fixed_pattern_noise;
00326 fors_image *image = NULL;
00327 fors_image *image2 = NULL;
00328 int nx, ny;
00329
00330 assure( first_raw != NULL, return -1, NULL );
00331 assure( second_raw != NULL, return -1, NULL );
00332
00333
00334
00335
00336
00337 nx = fors_image_get_size_x(first_raw);
00338 ny = fors_image_get_size_y(first_raw);
00339
00340 image = fors_image_duplicate(first_raw);
00341 fors_image_crop(image,
00342 1, 1,
00343 nx - 10, ny - 10);
00344
00345 image2 = fors_image_duplicate(second_raw);
00346 fors_image_crop(image2,
00347 11, 11,
00348 nx, ny);
00349
00350 fors_image_subtract(image, image2);
00351
00352 bias_fixed_pattern_noise = fors_image_get_stdev_robust(image, 50, NULL)
00353 / sqrt(2);
00354
00355
00356
00357
00358
00359 if (bias_fixed_pattern_noise > ron) {
00360
00361 bias_fixed_pattern_noise = sqrt(bias_fixed_pattern_noise *
00362 bias_fixed_pattern_noise
00363 -
00364 ron * ron);
00365 }
00366 else {
00367 cpl_msg_warning(cpl_func,
00368 "Zero-shift noise (%f ADU) is greater than "
00369 "accumulated zero-shift and fixed pattern "
00370 "noise (%f ADU), "
00371 "setting fixed pattern noise to zero",
00372 ron,
00373 bias_fixed_pattern_noise);
00374 bias_fixed_pattern_noise = 0;
00375 }
00376
00377 cleanup;
00378 return bias_fixed_pattern_noise;
00379 }
00380
00381
00382 #undef cleanup
00383 #define cleanup
00384
00389 double
00390 fors_get_airmass(const cpl_propertylist *header)
00391 {
00392 double airmass_start, airmass_end;
00393 airmass_start = cpl_propertylist_get_double(header, FORS_PFITS_AIRMASS_START);
00394 assure( !cpl_error_get_code(), return -1,
00395 "Could not read %s from header",
00396 FORS_PFITS_AIRMASS_START);
00397
00398 airmass_end = cpl_propertylist_get_double(header, FORS_PFITS_AIRMASS_END);
00399 assure( !cpl_error_get_code(), return -1,
00400 "Could not read %s from header",
00401 FORS_PFITS_AIRMASS_END);
00402
00403 return 0.5 * (airmass_start + airmass_end);
00404 }
00405