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_zeropoint_impl.h>
00033
00034 #include <fors_extract.h>
00035 #include <fors_identify.h>
00036 #include <fors_tools.h>
00037 #include <fors_star.h>
00038 #include <fors_std_star.h>
00039 #include <fors_image.h>
00040 #include <fors_data.h>
00041 #include <fors_setting.h>
00042 #include <fors_dfs.h>
00043 #include <fors_qc.h>
00044 #include <fors_pfits.h>
00045 #include <fors_utils.h>
00046
00047 #include <cpl.h>
00048
00049 #include <math.h>
00050
00057 const char *const fors_zeropoint_name = "fors_zeropoint";
00058 const char *const fors_zeropoint_description_short = "Compute zeropoint";
00059 const char *const fors_zeropoint_author = "Jonas M. Larsen";
00060 const char *const fors_zeropoint_email = PACKAGE_BUGREPORT;
00061 const char *const fors_zeropoint_description =
00062 "Input files:\n"
00063 " DO category: Type: Explanation: Number:\n"
00064 " STANDARD_IMG FITS image Phot. standard field 1\n"
00065 " MASTER_BIAS FITS image Master bias 1\n"
00066 " MASTER_SKY_FLAT_IMG FITS image Master sky flatfield 1\n"
00067 " FLX_STD_IMG FITS table Standard star catalog 1+\n"
00068 " PHOT_TABLE FITS table Filter ext. coeff, color 1\n"
00069 "\n"
00070 "Output files:\n"
00071 " DO category: Data type: Explanation:\n"
00072 " SOURCES_STD_IMG FITS image Unfiltered SExtractor output\n"
00073 " ALIGNED_PHOT FITS table\n"
00074 " PHOT_BACKGROUND_STD_IMG FITS image Reduced science image background\n"
00075 " STANDARD_REDUCED_IMG FITS image Reduced std image\n";
00076
00077 static double
00078 get_zeropoint(fors_star_list *stars,
00079 double cutoffE,
00080 double cutoffk,
00081 double dext_coeff,
00082 double dcolor_term,
00083 double avg_airmass,
00084 double *dzeropoint,
00085 int *n);
00086
00091 void fors_zeropoint_define_parameters(cpl_parameterlist *parameters)
00092 {
00093 const char *context = cpl_sprintf("fors.%s", fors_zeropoint_name);
00094
00095 fors_extract_define_parameters(parameters, context);
00096
00097 fors_identify_define_parameters(parameters, context);
00098
00099 const char *name, *full_name;
00100 cpl_parameter *p;
00101
00102 name = "magcutE";
00103 full_name = cpl_sprintf("%s.%s", context, name);
00104 p = cpl_parameter_new_value(full_name,
00105 CPL_TYPE_DOUBLE,
00106 "Zeropoint absolute cutoff (magnitude)",
00107 context,
00108 1.0);
00109 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
00110 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00111 cpl_parameterlist_append(parameters, p);
00112 cpl_free((void *)full_name); full_name = NULL;
00113
00114 name = "magcutk";
00115 full_name = cpl_sprintf("%s.%s", context, name);
00116 p = cpl_parameter_new_value(full_name,
00117 CPL_TYPE_DOUBLE,
00118 "Zeropoint kappa rejection parameter",
00119 context,
00120 5.0);
00121 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
00122 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00123 cpl_parameterlist_append(parameters, p);
00124 cpl_free((void *)full_name); full_name = NULL;
00125
00126
00127 cpl_free((void *)context);
00128
00129 return;
00130 }
00131
00132 #undef cleanup
00133 #define cleanup \
00134 do { \
00135 cpl_frameset_delete(std_frame); \
00136 cpl_frameset_delete(master_bias_frame); \
00137 cpl_frameset_delete(master_flat_frame); \
00138 cpl_frameset_delete(std_cat); \
00139 cpl_frameset_delete(phot_table); \
00140 fors_image_delete(&std); \
00141 fors_image_delete_const(&master_bias); \
00142 fors_image_delete(&master_flat); \
00143 cpl_table_delete(aligned_phot); \
00144 cpl_image_delete(background); \
00145 cpl_table_delete(sources); \
00146 fors_extract_method_delete(&em); \
00147 fors_identify_method_delete(&im); \
00148 fors_std_star_list_delete(&cat, fors_std_star_delete); \
00149 fors_star_list_delete(&stars, fors_star_delete); \
00150 cpl_free((void *)context); \
00151 fors_setting_delete(&setting); \
00152 cpl_propertylist_delete(qc); \
00153 cpl_propertylist_delete(product_header); \
00154 } while (0)
00155
00164 void fors_zeropoint(cpl_frameset *frames, const cpl_parameterlist *parameters)
00165 {
00166
00167 cpl_frameset *std_frame = NULL;
00168 fors_image *std = NULL;
00169
00170
00171 cpl_frameset *master_bias_frame = NULL;
00172 const fors_image *master_bias = NULL;
00173
00174 cpl_frameset *master_flat_frame = NULL;
00175 fors_image *master_flat = NULL;
00176
00177 cpl_frameset *std_cat = NULL;
00178 fors_std_star_list *cat = NULL;
00179
00180 cpl_frameset *phot_table = NULL;
00181 double color_term, dcolor_term;
00182 double ext_coeff, dext_coeff;
00183 double expected_zeropoint, dexpected_zeropoint;
00184
00185
00186 int nzeropoint = -1;
00187 cpl_table *aligned_phot = NULL;
00188 cpl_propertylist *qc = cpl_propertylist_new();
00189 double zeropoint, dzeropoint;
00190 fors_extract_sky_stats sky_stats;
00191 cpl_image *background = NULL;
00192 cpl_table *sources = NULL;
00193 cpl_propertylist *product_header = cpl_propertylist_new();
00194
00195
00196 extract_method *em = NULL;
00197 identify_method *im = NULL;
00198 double cutoffE, cutoffk;
00199
00200
00201 const char *context = cpl_sprintf("fors.%s", fors_zeropoint_name);
00202 fors_star_list *stars = NULL;
00203 fors_setting *setting = NULL;
00204 double avg_airmass;
00205
00206
00207 em = fors_extract_method_new(parameters, context);
00208 assure( !cpl_error_get_code(), return,
00209 "Could not get extraction parameters" );
00210
00211 im = fors_identify_method_new(parameters, context);
00212 assure( !cpl_error_get_code(), return,
00213 "Could not get identification parameters" );
00214
00215
00216 cpl_msg_indent_more();
00217 const char *name = cpl_sprintf("%s.%s", context, "magcutE");
00218 cutoffE = dfs_get_parameter_double_const(parameters,
00219 name);
00220 cpl_free((void *)name); name = NULL;
00221 cpl_msg_indent_less();
00222 assure( !cpl_error_get_code(), return, NULL );
00223
00224 cpl_msg_indent_more();
00225 name = cpl_sprintf("%s.%s", context, "magcutk");
00226 cutoffk = dfs_get_parameter_double_const(parameters,
00227 name);
00228 cpl_free((void *)name); name = NULL;
00229 cpl_msg_indent_less();
00230 assure( !cpl_error_get_code(), return, NULL );
00231
00232
00233 std_frame = fors_frameset_extract(frames, STANDARD_IMG);
00234 assure( cpl_frameset_get_size(std_frame) == 1, return,
00235 "Exactly 1 %s required. %d found",
00236 STANDARD_IMG, cpl_frameset_get_size(std_frame) );
00237
00238
00239 master_bias_frame = fors_frameset_extract(frames, MASTER_BIAS);
00240 assure( cpl_frameset_get_size(master_bias_frame) == 1, return,
00241 "One %s required. %d found",
00242 MASTER_BIAS, cpl_frameset_get_size(master_bias_frame) );
00243
00244 master_flat_frame = fors_frameset_extract(frames, MASTER_SKY_FLAT_IMG);
00245 assure( cpl_frameset_get_size(master_flat_frame) == 1, return,
00246 "One %s required. %d found",
00247 MASTER_SKY_FLAT_IMG, cpl_frameset_get_size(master_flat_frame) );
00248
00249 std_cat = fors_frameset_extract(frames, FLX_STD_IMG);
00250 assure( cpl_frameset_get_size(std_cat) >= 1, return,
00251 "One or more %s required. %d found",
00252 FLX_STD_IMG, cpl_frameset_get_size(std_cat));
00253
00254 phot_table = fors_frameset_extract(frames, PHOT_TABLE);
00255 assure( cpl_frameset_get_size(phot_table) == 1, return,
00256 "One %s required. %d found",
00257 PHOT_TABLE, cpl_frameset_get_size(phot_table));
00258
00259
00260
00261
00262 setting = fors_setting_new(cpl_frameset_get_first(std_frame));
00263 assure( !cpl_error_get_code(), return, "Could not get instrument setting" );
00264
00265 master_bias = fors_image_load(cpl_frameset_get_first(master_bias_frame),
00266 NULL, setting, NULL);
00267 assure( !cpl_error_get_code(), return,
00268 "Could not load master bias");
00269
00270
00271 std = fors_image_load(cpl_frameset_get_first(std_frame),
00272 master_bias, setting, NULL);
00273 assure( !cpl_error_get_code(), return, "Could not load standard image");
00274 fors_image_delete_const(&master_bias);
00275
00276
00277 master_flat = fors_image_load(cpl_frameset_get_first(master_flat_frame),
00278 NULL, setting, NULL);
00279 assure( !cpl_error_get_code(), return, "Could not load master flat");
00280
00281
00282 fors_image_divide_scalar(master_flat,
00283 fors_image_get_median(master_flat, NULL), -1.0);
00284 fors_image_divide(std, master_flat);
00285 assure( !cpl_error_get_code(), return, "Could not divide by master flat");
00286 fors_image_delete(&master_flat);
00287
00288
00289 stars = fors_extract(std, setting, em,
00290 &sky_stats, &background, &sources);
00291 assure( !cpl_error_get_code(), return, "Could not extract objects");
00292
00293 if (setting->filter_name != NULL) {
00294
00295 fors_phot_table_load(cpl_frameset_get_first(phot_table), setting,
00296 &color_term, &dcolor_term,
00297 &ext_coeff, &dext_coeff,
00298 &expected_zeropoint, &dexpected_zeropoint);
00299 assure( !cpl_error_get_code(), return,
00300 "Could not load photometry table" );
00301
00302
00303 cat = fors_std_cat_load(std_cat,
00304 cpl_frameset_get_first(std_frame),
00305 setting,
00306 color_term, dcolor_term);
00307
00308 assure( !cpl_error_get_code(), return,
00309 "Failed to create catalogue patterns" );
00310
00311
00312
00313 fors_identify(stars, cat, im);
00314 assure( !cpl_error_get_code(), return, "Failed to identify sources");
00315
00316
00317 avg_airmass = fors_star_ext_corr(stars, setting, ext_coeff, dext_coeff,
00318 cpl_frameset_get_first(std_frame));
00319 assure( !cpl_error_get_code(), return,
00320 "Extinction correction failed");
00321
00322
00323 zeropoint = get_zeropoint(stars, cutoffE, cutoffk,
00324 dext_coeff, dcolor_term,
00325 avg_airmass,
00326 &dzeropoint, &nzeropoint);
00327 assure( !cpl_error_get_code(), return,
00328 "Failed to get zeropoint" );
00329 }
00330 else {
00331 cpl_msg_warning(cpl_func, "Zeropoint computation is not supported "
00332 "for non-standard filters");
00333 color_term = 0.0;
00334 dcolor_term = 9999.0;
00335 ext_coeff = 0.0;
00336 dext_coeff = 9999.0;
00337 expected_zeropoint = 0.0;
00338 dexpected_zeropoint = 9999.0;
00339 zeropoint = 0.0;
00340 dzeropoint = 0.0;
00341 nzeropoint = 0;
00342 }
00343
00344
00345 fors_qc_start_group(qc, fors_qc_dic_version, setting->instrument);
00346
00347 fors_qc_write_group_heading(cpl_frameset_get_first(std_frame),
00348 ALIGNED_PHOT,
00349 setting->instrument);
00350 assure( !cpl_error_get_code(), return, "Could not write %s QC parameters",
00351 ALIGNED_PHOT);
00352
00353
00354 fors_qc_write_qc_double(qc,
00355 zeropoint,
00356 "QC.ZPOINT",
00357 "mag",
00358 "Frame zeropoint",
00359 setting->instrument);
00360 fors_qc_write_qc_double(qc,
00361 dzeropoint,
00362 "QC.ZPOINTRMS",
00363 "mag",
00364 "Uncertainty of frame zeropoint",
00365 setting->instrument);
00366 fors_qc_write_qc_int(qc,
00367 nzeropoint,
00368 "QC.ZPOINT.NSTARS",
00369 NULL,
00370 "Number of stars used for zeropoint computation",
00371 setting->instrument);
00372
00373 double derived_ext_coeff, derived_ext_coeff_err;
00374
00375 if (setting->filter_name != NULL) {
00376 cpl_msg_info(cpl_func,
00377 "Computing extinction "
00378 "(assuming zeropoint = %.3f +- %.3f mag)",
00379 expected_zeropoint,
00380 dexpected_zeropoint);
00381
00382 cpl_msg_indent_more();
00383
00384 if (nzeropoint > 0) {
00385 derived_ext_coeff = ext_coeff +
00386 (expected_zeropoint - zeropoint) / avg_airmass;
00387
00388
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401 derived_ext_coeff_err = sqrt(
00402 dexpected_zeropoint*dexpected_zeropoint +
00403 dzeropoint*dzeropoint) / avg_airmass - dext_coeff*dext_coeff;
00404
00405 cpl_msg_info(cpl_func, "Atmospheric extinction = "
00406 "%f +- %f mag/airmass", derived_ext_coeff,
00407 derived_ext_coeff_err);
00408 }
00409 else {
00410 cpl_msg_warning(cpl_func, "Too few stars available, "
00411 "setting extinction to zero");
00412 derived_ext_coeff = 0;
00413 derived_ext_coeff_err = 9999;
00414 }
00415 }
00416 else {
00417 derived_ext_coeff = 0;
00418 derived_ext_coeff_err = 9999;
00419 }
00420
00421 fors_qc_write_qc_double(qc,
00422 derived_ext_coeff,
00423 "QC.EXTCOEFF",
00424 "mag/airmass",
00425 "Atmospheric extinction (mag/airmass)",
00426 setting->instrument);
00427
00428 fors_qc_write_qc_double(qc,
00429 derived_ext_coeff_err,
00430 "QC.EXTCOEFFERR",
00431 "mag/airmass",
00432 "Error on atmospheric extinction (mag/airmass)",
00433 setting->instrument);
00434 cpl_msg_indent_less();
00435
00436 fors_qc_end_group();
00437
00438
00439
00440 aligned_phot = fors_create_sources_table(stars);
00441 assure( !cpl_error_get_code(), return,
00442 "Failed to create aligned photometry table");
00443
00444
00445 fors_dfs_save_table(frames, sources, SOURCES_STD,
00446 NULL, parameters, fors_zeropoint_name,
00447 cpl_frameset_get_first(std_frame));
00448 assure( !cpl_error_get_code(), return, "Saving %s failed",
00449 SOURCES_STD);
00450
00451
00452
00453
00454 cpl_propertylist_update_double(qc, FORS_PFITS_EXPOSURE_TIME,
00455 setting->exposure_time);
00456 cpl_propertylist_update_double(qc, "AIRMASS", avg_airmass);
00457
00458 fors_dfs_save_table(frames, aligned_phot, ALIGNED_PHOT,
00459 qc, parameters, fors_zeropoint_name,
00460 cpl_frameset_get_first(std_frame));
00461 assure( !cpl_error_get_code(), return, "Saving %s failed",
00462 ALIGNED_PHOT);
00463
00464 cpl_propertylist_delete(qc); qc = NULL;
00465
00466
00467 fors_dfs_add_wcs(product_header,
00468 cpl_frameset_get_first(std_frame), setting);
00469 fors_dfs_add_exptime(product_header,
00470 cpl_frameset_get_first(std_frame), 0.);
00471
00472 fors_dfs_save_image(frames, std, STANDARD_REDUCED_IMG,
00473 product_header, parameters, fors_zeropoint_name,
00474 cpl_frameset_get_first(std_frame));
00475 assure( !cpl_error_get_code(), return, "Saving %s failed",
00476 STANDARD_REDUCED_IMG);
00477
00478 dfs_save_image(frames, background, PHOT_BACKGROUND_STD_IMG,
00479 product_header, parameters, fors_zeropoint_name,
00480 setting->version);
00481 assure( !cpl_error_get_code(), return, "Saving %s failed",
00482 PHOT_BACKGROUND_STD_IMG);
00483
00484 cpl_image_delete(background); background = NULL;
00485
00486 if (setting->filter_name == NULL) {
00487
00488
00489
00490 cleanup;
00491 return;
00492 }
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506
00507
00508
00509 double color = fors_image_get_min(std);
00510 {
00511 fors_star *s;
00512 for (s = fors_star_list_first(stars);
00513 s != NULL;
00514 s = fors_star_list_next(stars)) {
00515 fors_image_draw(std, 0,
00516 s->pixel->x,
00517 s->pixel->y,
00518 10, color);
00519
00520 if (s->id != NULL) {
00521 fors_image_draw(std, 2,
00522 s->pixel->x,
00523 s->pixel->y,
00524 10, color);
00525 }
00526 }
00527 }
00528 {
00529 fors_std_star *s;
00530 for (s = fors_std_star_list_first(cat);
00531 s != NULL;
00532 s = fors_std_star_list_next(cat)) {
00533 fors_image_draw(std, 1,
00534 s->pixel->x,
00535 s->pixel->y,
00536 10, color);
00537 }
00538 fors_dfs_save_image(frames, std, "DEBUG",
00539 product_header, parameters, fors_zeropoint_name,
00540 cpl_frameset_get_first(std_frame));
00541 assure( !cpl_error_get_code(), return, "Saving %s failed",
00542 "DEBUG");
00543 }
00544
00545 cleanup;
00546 return;
00547 }
00548
00557 static bool
00558 zeropoint_inside(const fors_star *s,
00559 void *data)
00560 {
00561 struct {
00562 double hi, lo;
00563 double z, kappa;
00564 } *cuts = data;
00565
00566 double z = fors_star_get_zeropoint(s, NULL);
00567 double dz = fors_star_get_zeropoint_err(s, NULL);
00568
00569 return
00570 (cuts->lo <= z && z <= cuts->hi) ||
00571 (cuts->z - cuts->kappa*dz <= z && z <= cuts->z + cuts->kappa*dz);
00572 }
00573
00574 #undef cleanup
00575 #define cleanup \
00576 do { \
00577 fors_star_list_delete(&subset, fors_star_delete); \
00578 fors_star_list_delete(&identified, fors_star_delete); \
00579 } while(0)
00580
00592 static double
00593 get_zeropoint(fors_star_list *stars,
00594 double cutoffE,
00595 double cutoffk,
00596 double dext_coeff,
00597 double dcolor_term,
00598 double avg_airmass,
00599 double *dzeropoint,
00600 int *n)
00601 {
00602 fors_star_list *subset = NULL;
00603 fors_star_list *identified =
00604 fors_star_list_extract(stars, fors_star_duplicate,
00605 fors_star_is_identified, NULL);
00606
00607 assure( stars != NULL, return 0, NULL );
00608 assure( dzeropoint != NULL, return 0, NULL );
00609 assure( n != NULL, return 0, NULL );
00610
00611 if ( fors_star_list_size(identified) == 0 ) {
00612 cpl_msg_warning(cpl_func,
00613 "No identified stars for zeropoint computation");
00614 *n = 0;
00615 *dzeropoint = 0;
00616 cleanup;
00617 return 0;
00618 }
00619
00620 cpl_msg_info(cpl_func, "Computing zeropoint (assuming extinction)");
00621 cpl_msg_indent_more();
00622
00623 double zeropoint;
00624 double red_chisq = -1.0;
00625
00626
00627
00628
00629 zeropoint = fors_star_list_mean_optimal(identified,
00630 fors_star_get_zeropoint, NULL,
00631 fors_star_get_zeropoint_err, NULL,
00632 dzeropoint,
00633 fors_star_list_size(identified) >= 2 ? &red_chisq : NULL);
00634
00635 cpl_msg_info(cpl_func, "Optimal zeropoint (no rejection, %d stars) = %f mag",
00636 fors_star_list_size(identified), zeropoint);
00637
00638
00639
00640
00641
00642
00643 struct {
00644 double hi, lo;
00645 double z, kappa;
00646 } cuts;
00647 cuts.hi = zeropoint + 5*cutoffE;
00648 cuts.lo = zeropoint - 5*cutoffE;
00649 cuts.z = zeropoint;
00650 cuts.kappa = cutoffk;
00651
00652 subset = fors_star_list_extract(identified, fors_star_duplicate,
00653 zeropoint_inside, &cuts);
00654
00655
00656 if ( fors_star_list_size(subset) == 0 ) {
00657 cpl_msg_warning(cpl_func,
00658 "All stars rejected (%f mag). Cannot "
00659 "compute zeropoint", 5*cutoffE);
00660 *n = 0;
00661 *dzeropoint = 0;
00662 cleanup;
00663 return 0;
00664 }
00665
00666 zeropoint = fors_star_list_mean_optimal(subset,
00667 fors_star_get_zeropoint, NULL,
00668 fors_star_get_zeropoint_err, NULL,
00669 dzeropoint,
00670 fors_star_list_size(subset) >= 2 ? &red_chisq : NULL);
00671
00672 cpl_msg_debug(cpl_func, "Optimal zeropoint (%.2f mag, %.2f sigma rejection) = %f mag",
00673 5*cutoffE, cutoffk, zeropoint);
00674
00675 cuts.hi = zeropoint + cutoffE;
00676 cuts.lo = zeropoint - cutoffE;
00677 cuts.z = zeropoint;
00678 cuts.kappa = cutoffk;
00679
00680 {
00681 fors_star_list *tmp = fors_star_list_duplicate(subset, fors_star_duplicate);
00682 fors_star_list_delete(&subset, fors_star_delete);
00683
00684 subset = fors_star_list_extract(tmp, fors_star_duplicate,
00685 zeropoint_inside, &cuts);
00686
00687 if ( fors_star_list_size(subset) == 0 ) {
00688 cpl_msg_warning(cpl_func,
00689 "All stars rejected (%f mag, %f sigma). Cannot "
00690 "compute zeropoint", cutoffE, cutoffk);
00691 *n = 0;
00692 *dzeropoint = 0;
00693 cleanup;
00694 return 0;
00695 }
00696
00697 fors_star_list_delete(&tmp, fors_star_delete);
00698 }
00699
00700 zeropoint = fors_star_list_mean_optimal(subset,
00701 fors_star_get_zeropoint, NULL,
00702 fors_star_get_zeropoint_err, NULL,
00703 dzeropoint,
00704 fors_star_list_size(subset) >= 2 ? &red_chisq : NULL);
00705
00706 cpl_msg_info(cpl_func, "Optimal zeropoint (%.2f mag, %.2f sigma rejection) = %f mag",
00707 cutoffE, cutoffk, zeropoint);
00708
00709 *n = fors_star_list_size(subset);
00710 {
00711 int outliers =
00712 fors_star_list_size(identified) - fors_star_list_size(subset);
00713 cpl_msg_info(cpl_func, "%d outlier%s rejected, %d non-outliers",
00714 outliers, outliers == 1 ? "" : "s",
00715 *n);
00716 }
00717
00718 if ( *n == 0 ) {
00719 cpl_msg_warning(cpl_func,
00720 "All stars were rejected during zeropoint computation" );
00721 *dzeropoint = 0;
00722 cleanup;
00723 return 0;
00724 }
00725
00726
00727
00728
00729
00730
00731
00732
00733
00734
00735
00736
00737
00738
00739
00740
00741
00742
00743
00744
00745 cpl_matrix *covariance = cpl_matrix_new(*n,
00746 *n);
00747
00748
00749 fors_star_list *ident_dup = fors_star_list_duplicate(subset, fors_star_duplicate);
00750 {
00751
00752 fors_star *s, *t;
00753 int i, j;
00754 for (s = fors_star_list_first(subset), i = 0;
00755 s != NULL;
00756 s = fors_star_list_next(subset), i++) {
00757
00758 for (t = fors_star_list_first(ident_dup), j = 0;
00759 t != NULL;
00760 t = fors_star_list_next(ident_dup), j++) {
00761
00762 double cij;
00763
00764 if (fors_star_equal(s, t)) {
00765 cij = fors_star_get_zeropoint_err(s, NULL)*fors_star_get_zeropoint_err(s, NULL);
00766
00767 }
00768 else {
00769 cij = s->id->color * t->id->color * dcolor_term*dcolor_term
00770 + avg_airmass*avg_airmass*dext_coeff*dext_coeff;
00771 }
00772
00773 cpl_matrix_set(covariance, i, j, cij);
00774 }
00775 }
00776 }
00777
00778
00779
00780
00781
00782
00783
00784
00785
00786
00787
00788
00789 cpl_matrix *covariance_inverse = cpl_matrix_invert_create(covariance);
00790
00791 assure( !cpl_error_get_code(), return 0,
00792 "Could not invert zeropoints covariance matrix");
00793
00794
00795
00796
00797
00798 cpl_matrix_delete(covariance); covariance = NULL;
00799
00800 cpl_matrix *const_vector = cpl_matrix_new(*n, 1);
00801 cpl_matrix_fill(const_vector, 1.0);
00802
00803 cpl_matrix *weights = cpl_matrix_product_create(covariance_inverse, const_vector);
00804
00805 cpl_matrix_delete(const_vector); const_vector = NULL;
00806
00807
00808
00809
00810 {
00811 double wz = 0;
00812 double w = 0;
00813
00814 int i;
00815 fors_star *s;
00816 for (i = 0, s = fors_star_list_first(subset);
00817 s != NULL;
00818 s = fors_star_list_next(subset), i++) {
00819
00820 double weight = cpl_matrix_get(weights, i, 0);
00821
00822 cpl_msg_debug(cpl_func, "Weight_%d = %f", i, weight);
00823
00824 wz += weight * fors_star_get_zeropoint(s, NULL);
00825 w += weight;
00826
00827
00828 {
00829 fors_star *t;
00830
00831 for (t = fors_star_list_first(stars);
00832 t != NULL;
00833 t = fors_star_list_next(stars)) {
00834
00835 if (fors_star_equal(s, t)) {
00836 t->weight = weight;
00837 }
00838 }
00839 }
00840 }
00841
00842 cpl_matrix_delete(weights); weights = NULL;
00843
00844 cpl_msg_debug(cpl_func, "Sum of weights = %f", w);
00845
00846
00847
00848
00849
00850
00851
00852
00853
00854
00855
00856
00857
00858
00859 assure( w != 0, return 0, "Sum of optimal weights is zero!" );
00860 assure( sqrt(w) != 0, return 0, "Square root of sum of optimal weights is zero!" );
00861
00862 zeropoint = wz / w;
00863 *dzeropoint = 1 / sqrt(w);
00864 }
00865
00866
00867
00868
00869
00870
00871
00872
00873
00874
00875
00876 cpl_msg_info(cpl_func, "Optimal zeropoint = %f +- %f mag",
00877 zeropoint, *dzeropoint);
00878
00879 if (*n >= 2) {
00880 fors_star *s, *t;
00881 int i, j;
00882
00883 red_chisq = 0;
00884
00885 for (s = fors_star_list_first(subset), i = 0;
00886 s != NULL;
00887 s = fors_star_list_next(subset), i++) {
00888
00889 for (t = fors_star_list_first(ident_dup), j = 0;
00890 t != NULL;
00891 t = fors_star_list_next(ident_dup), j++) {
00892
00893 red_chisq +=
00894 (fors_star_get_zeropoint(s, NULL) - zeropoint) *
00895 (fors_star_get_zeropoint(t, NULL) - zeropoint) *
00896 cpl_matrix_get(covariance_inverse, i, j);
00897 }
00898 }
00899 red_chisq /= (*n - 1);
00900 }
00901
00902 cpl_matrix_delete(covariance_inverse); covariance_inverse = NULL;
00903
00904 fors_star_list_delete(&ident_dup, fors_star_delete);
00905
00906 cpl_msg_info(cpl_func, "Reduced chi square = %f", red_chisq);
00907
00908 cpl_msg_indent_less();
00909
00910 cleanup;
00911 return zeropoint;
00912 }
00913
00914