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_photometry_impl.h>
00033
00034 #include <fors_data.h>
00035 #include <fors_qc.h>
00036 #include <fors_dfs.h>
00037 #include <fors_tools.h>
00038 #include <fors_pfits.h>
00039 #include <fors_polynomial.h>
00040 #include <fors_utils.h>
00041 #include <fors_double.h>
00042 #include <fors_extract.h>
00043
00044 #include <cpl.h>
00045 #include <math.h>
00046 #include <assert.h>
00047 #include <string.h>
00048 #include <stdlib.h>
00049
00056 const char *const fors_photometry_name = "fors_photometry";
00057 const char *const fors_photometry_description_short = "Compute corrected flatfield";
00058 const char *const fors_photometry_author = "Jonas M. Larsen";
00059 const char *const fors_photometry_email = PACKAGE_BUGREPORT;
00060 const char *const fors_photometry_description =
00061 "Input files:\n"
00062 " DO category: Type: Explanation: Number:\n"
00063 " PHOT_TABLE FITS table Expected extinction params 1\n"
00064 " ALIGNED_PHOT FITS table Photometry 1+\n"
00065 " MASTER_SKY_FLAT_IMG FITS image Master flat field 1\n"
00066 "\n"
00067 "Output files:\n"
00068 " DO category: Data type: Explanation:\n"
00069 " PHOT_COEFF_TABLE FITS image Observed extinction coefficients\n"
00070 " CORRECTION_MAP FITS image Correction map (magnitude)\n"
00071 " CORRECTION_FACTOR FITS image Correction map (flux)\n"
00072 " MASTER_FLAT_IMG FITS image Corrected master flat field\n";
00073
00074 typedef enum fors_fit_ncoeff {
00075 FORS_FIT_NCOEFF_NO = 0,
00076 FORS_FIT_NCOEFF_ONE,
00077 FORS_FIT_NCOEFF_PERFRAME,
00078 FORS_FIT_NCOEFF_PERNIGHT
00079 } fors_fit_ncoeff;
00080
00081 const struct fors_fit_ncoeff_paropts {
00082 const char *no,
00083 *one,
00084 *perframe,
00085 *pernight;
00086 } fors_fit_ncoeff_paropts = { "no",
00087 "one",
00088 "perframe",
00089 "pernight"};
00090
00091 typedef struct entry
00092 {
00093 int frame_index,
00094 star_index,
00095 atm_ext_index;
00096 int atm_ext_identifier;
00097 double airmass,
00098 gain,
00099 exptime;
00100 fors_star *star;
00101 } entry;
00102
00103
00104 #undef LIST_ELEM
00105 #define LIST_ELEM entry
00106 #undef LIST_DEFINE
00107 #include <list.h>
00108 #define LIST_DEFINE
00109 #include <list.h>
00110
00111 const double arcsec_tol = 5.0;
00112
00113 entry *entry_new( int frame_index,
00114 int star_index,
00115 double airmass,
00116 double gain,
00117 double exptime,
00118 int atm_ext_identifier,
00119 fors_star *star);
00120
00121 void
00122 entry_delete( entry **e);
00123
00124 static void
00125 entry_delete_but_standard( entry **e);
00126
00127 static double
00128 entry_get_powers_x_y( const entry *e,
00129 const cpl_array *powers);
00130
00131 void
00132 entry_list_print( const entry_list *l,
00133 cpl_msg_severity level);
00134
00135 static double
00136 entry_get_powers_airmass_color( const entry *e,
00137 const cpl_array *powers);
00138
00139 static entry_list *
00140 fors_photometry_read_input( const cpl_frameset *alphot_frames,
00141 const fors_setting *setting,
00142 int (*get_atm_ext_id_function)(
00143 const cpl_propertylist *header),
00144 bool import_unknown,
00145 int *n_frames,
00146 fors_std_star_list **std_star_list,
00147 int filter
00148 );
00149
00150 static fors_std_star*
00151 fors_photometry_read_input_listinsert_star_if_new(
00152 fors_std_star_list *std_list,
00153 fors_std_star *std,
00154 double mind_arcsec);
00155
00156 static void
00157 fors_delete_star_lists( entry_list **el,
00158 fors_std_star_list **sl);
00159
00160 static bool
00161 fors_fits_compare_string( const char *s1,
00162 const char *s2);
00163
00164 static void
00165 fors_matrix_null( cpl_matrix **m);
00166
00167 static void
00168 fors_matrix_append_delete( cpl_matrix **m1,
00169 cpl_matrix **m2);
00170
00171 static double
00172 fors_property_get_num( const cpl_property *prop);
00173
00174 static double
00175 fors_photometry_parameter_get_num( const cpl_parameterlist *parameters,
00176 const char *name,
00177 cpl_type type);
00178
00179 static const char*
00180 fors_photometry_parameter_get_string( const cpl_parameterlist *parameters,
00181 const char *name);
00182
00183 static fors_fit_ncoeff
00184 fors_photometry_parameter_get_ncoeff( const cpl_parameterlist *parameters,
00185 const char *name);
00186
00187 void fors_photometry_define_parameters( cpl_parameterlist *parameters);
00188
00189 static cpl_polynomial*
00190 fors_photometry_define_polyf( int degreef1,
00191 int degreef2);
00192
00193 static cpl_polynomial*
00194 fors_photometry_define_polyp( int degreep);
00195
00196 static cpl_error_code
00197 fors_photometry_poly_new_from_coefficients( const cpl_polynomial *p_def,
00198 const cpl_matrix *coeffs,
00199 const cpl_matrix *cov_coeffs,
00200 cpl_polynomial **poly,
00201 cpl_polynomial **var_poly);
00202
00203 int
00204 fors_photometry_get_timezone_observer( const cpl_propertylist *header);
00205
00206 int
00207 fors_photometry_get_night_id( const cpl_propertylist *header);
00208
00209 static int
00210 fors_photometry_atm_ext_create_index_by_identifier(
00211 entry_list *obs_list);
00212
00213 static int
00214 fors_photometry_atm_ext_create_indices( entry_list *obsl,
00215 fors_fit_ncoeff fit_e);
00216
00217
00218 static cpl_table *
00219 fors_photometry_atm_ext_print_index_by_framename(
00220 const entry_list *obs_list,
00221 const cpl_frameset *frames);
00222
00223 static cpl_error_code
00224 fors_photometry_check_input_value( double value,
00225 double value_error,
00226 const char *value_name,
00227 const char *input_name,
00228 double min_limit,
00229 double max_limit,
00230 double max_error);
00231
00232 static cpl_error_code
00233 fors_photometry_check_fitparam_atm_ext( entry_list *obsl,
00234 fors_fit_ncoeff fit_e,
00235 bool fit_z);
00236
00237 static cpl_error_code
00238 fors_photometry_adjust_fit_mag_flags( fors_std_star_list *stdl,
00239 entry_list *obsl,
00240 bool override_fit_m,
00241 int *n_mag_fits);
00242
00243 static cpl_error_code
00244 fors_photometry_remove_unnecessary( fors_std_star_list *std_list,
00245 entry_list *obs_list,
00246 int *n_mag_fits);
00247
00248 static cpl_array*
00249 fors_photometry_count_observations( fors_std_star_list *std_list,
00250 entry_list *obs_list);
00251
00252 static cpl_matrix*
00253 build_equations_lhs_matrix_from_parameters( const entry_list *obs_list,
00254 const fors_std_star_list *std_list,
00255 bool fit_z,
00256 bool fit_c,
00257 int *n_fit_e_cols);
00258
00259 static cpl_matrix*
00260 build_equations_lhs_matrix_from_poly( const entry_list *obs_list,
00261 const cpl_polynomial *poly,
00262 const char *pname,
00263 double (*powerfunc)(
00264 const entry*,
00265 const cpl_array*));
00266
00267 static cpl_error_code
00268 build_equations_rhs_cov( const entry_list *obs_list,
00269 const fors_std_star_list *std_list,
00270 bool fit_z,
00271 bool fit_c,
00272 bool fit_e,
00273 double color_coeff,
00274 double dcolor_coeff,
00275 double ext_coeff,
00276 double dext_coeff,
00277 double zpoint,
00278 double dzpoint,
00279 cpl_matrix **rhs,
00280 cpl_matrix **rhs_cov);
00281
00282
00288
00289 entry *entry_new( int frame_index,
00290 int star_index,
00291 double airmass,
00292 double gain,
00293 double exptime,
00294 int atm_ext_identifier,
00295 fors_star *star)
00296 {
00297 entry *e = cpl_malloc(sizeof(*e));
00298
00299 e->frame_index = frame_index;
00300 e->star_index = star_index;
00301 e->atm_ext_index = -1;
00302 e->airmass = airmass;
00303 e->gain = gain;
00304 e->exptime = exptime;
00305
00306 e->atm_ext_identifier = atm_ext_identifier;
00307
00308 e->star = star;
00309
00310 return e;
00311 }
00312
00313
00318
00319 void
00320 entry_delete( entry **e)
00321 {
00322 if (e && *e) {
00323 fors_star_delete(&(*e)->star);
00324 cpl_free(*e); *e = NULL;
00325 }
00326 return;
00327 }
00328
00329
00334
00335 static void
00336 entry_delete_but_standard( entry **e)
00337 {
00338 if (e && *e) {
00339 fors_star_delete_but_standard(&(*e)->star);
00340 cpl_free(*e); *e = NULL;
00341 }
00342 return;
00343 }
00344
00345
00346
00352
00353 void
00354 entry_list_print( const entry_list *l,
00355 cpl_msg_severity level)
00356 {
00357 const entry *e;
00358
00359 fors_msg(level, "Observation list:");
00360 cpl_msg_indent_more();
00361 for (e = entry_list_first_const(l);
00362 e != NULL;
00363 e = entry_list_next_const(l)) {
00364
00365 fors_msg(level,
00366 "frame %d, star %d: airmass = %f, gain = %f, exptime = %f s",
00367 e->frame_index, e->star_index,
00368 e->airmass, e->gain, e->exptime);
00369
00370 fors_star_print(level, e->star);
00371 }
00372 cpl_msg_indent_less();
00373
00374 return;
00375 }
00376
00377
00378 #undef cleanup
00379 #define cleanup
00380
00381 static double
00382 entry_get_powers_x_y( const entry *e,
00383 const cpl_array *powers)
00384 {
00385 passure(powers != NULL && e != NULL, return sqrt(-1));
00386 passure(cpl_array_get_size(powers) == 2, return sqrt(-1));
00387 return pow(e->star->pixel->x, cpl_array_get(powers, 0, NULL))
00388 * pow(e->star->pixel->y, cpl_array_get(powers, 1, NULL));
00389 }
00390
00391
00392 #undef cleanup
00393 #define cleanup
00394
00395 static double
00396 entry_get_powers_airmass_color( const entry *e,
00397 const cpl_array *powers)
00398 {
00399 passure(powers != NULL && e != NULL, return sqrt(-1));
00400 passure(cpl_array_get_size(powers) == 2, return sqrt(-1));
00401 return pow(e->airmass, cpl_array_get(powers, 0, NULL))
00402 * pow(e->star->id->color, cpl_array_get(powers, 1, NULL));
00403 }
00404
00405
00412
00413 static void
00414 fors_delete_star_lists( entry_list **el,
00415 fors_std_star_list **sl)
00416 {
00417 entry *e;
00418 if (el != NULL && *el != NULL)
00419 {
00420 for ( e = entry_list_first(*el);
00421 e != NULL;
00422 e = entry_list_next(*el))
00423 {
00424 e->star->id = NULL;
00425 }
00426 }
00427
00428 fors_std_star_list_delete(sl, fors_std_star_delete);
00429 entry_list_delete(el, entry_delete);
00430 }
00431
00432
00441
00442 static bool
00443 fors_fits_compare_string( const char *s1,
00444 const char *s2)
00445 {
00446 const char *m1 = "",
00447 *m2 = "";
00448 int len1,
00449 len2;
00450
00451 if (s1 != NULL)
00452 m1 = s1;
00453 if (s2 != NULL)
00454 m2 = s2;
00455
00456 len1 = strlen(m1);
00457 len2 = strlen(m2);
00458
00459 while (len1 > 0 && m1[len1-1] == ' ') len1--;
00460 while (len2 > 0 && m2[len2-1] == ' ') len2--;
00461
00462 if (len1 != len2)
00463 return false;
00464
00465 if (strncmp(m1, m2, len1) != 0)
00466 return false;
00467
00468 return true;
00469 }
00470
00471
00472
00478
00479 static void
00480 fors_matrix_null( cpl_matrix **m)
00481 {
00482 if (m != NULL)
00483 {
00484 cpl_matrix_delete(*m);
00485 *m = NULL;
00486 }
00487 }
00488
00489
00490
00491 static void
00492 fors_matrix_append_delete( cpl_matrix **m1,
00493 cpl_matrix **m2)
00494 {
00495 if (m1 == NULL || m2 == NULL)
00496 return;
00497
00498 if (*m2 != NULL)
00499 {
00500 if (*m1 == NULL)
00501 {
00502 *m1 = *m2;
00503 *m2 = NULL;
00504 }
00505 else
00506 {
00507 cpl_matrix_append(*m1, *m2, 0);
00508 fors_matrix_null(m2);
00509 }
00510 }
00511
00512 }
00513
00514
00515
00516 static double
00517 fors_property_get_num( const cpl_property *prop)
00518 {
00519 double retval = 0;
00520 cpl_type type;
00521
00522 cassure_automsg( prop != NULL,
00523 CPL_ERROR_NULL_INPUT,
00524 return 0);
00525
00526 type = cpl_property_get_type(prop);
00527
00528 switch (type)
00529 {
00530 case CPL_TYPE_BOOL:
00531 retval = cpl_property_get_bool(prop);
00532 break;
00533 case CPL_TYPE_INT:
00534 retval = cpl_property_get_int(prop);
00535 break;
00536 case CPL_TYPE_FLOAT:
00537 retval = cpl_property_get_float(prop);
00538 break;
00539 case CPL_TYPE_DOUBLE:
00540 retval = cpl_property_get_double(prop);
00541 break;
00542 default:
00543 cpl_error_set_message( cpl_func, CPL_ERROR_INVALID_TYPE,
00544 "type must be bool, int, float or double");
00545 }
00546
00547 switch (type)
00548 {
00549 case CPL_TYPE_BOOL:
00550 return (fabs(retval) > 0.5) ? 1 : 0;
00551 case CPL_TYPE_INT:
00552 return round(retval);
00553 default:
00554 return retval;
00555 }
00556 }
00557
00558
00559
00560 static double
00561 fors_photometry_parameter_get_num( const cpl_parameterlist *parameters,
00562 const char *name,
00563 cpl_type type)
00564 {
00565 char *descriptor;
00566 double retval = -1;
00567
00568 cpl_msg_indent_more();
00569 descriptor = cpl_sprintf("fors.%s.%s", fors_photometry_name, name);
00570 switch (type)
00571 {
00572 case CPL_TYPE_BOOL:
00573 retval = dfs_get_parameter_bool_const(parameters, descriptor);
00574 break;
00575 case CPL_TYPE_INT:
00576 retval = dfs_get_parameter_int_const(parameters, descriptor);
00577 break;
00578
00579
00580
00581 case CPL_TYPE_DOUBLE:
00582 retval = dfs_get_parameter_double_const(parameters, descriptor);
00583 break;
00584 default:
00585 cpl_error_set_message( cpl_func, CPL_ERROR_INVALID_TYPE,
00586 "type must be bool, int"
00587
00588 " or double");
00589 }
00590 cpl_free(descriptor);
00591 cpl_msg_indent_less();
00592
00593 switch (type)
00594 {
00595 case CPL_TYPE_BOOL:
00596 return (fabs(retval) > 0.5) ? true : false;
00597 case CPL_TYPE_INT:
00598 return round(retval);
00599 default:
00600 return retval;
00601 }
00602 }
00603
00604
00605
00606 static const char*
00607 fors_photometry_parameter_get_string( const cpl_parameterlist *parameters,
00608 const char *name)
00609 {
00610 char *descriptor;
00611 const char *retval = NULL;
00612
00613 cpl_msg_indent_more();
00614 descriptor = cpl_sprintf("fors.%s.%s", fors_photometry_name, name);
00615
00616 retval = dfs_get_parameter_string_const(parameters, descriptor);
00617
00618 cpl_free(descriptor);
00619 cpl_msg_indent_less();
00620
00621 return retval;
00622 }
00623
00624
00625
00626 static fors_fit_ncoeff
00627 fors_photometry_parameter_get_ncoeff( const cpl_parameterlist *parameters,
00628 const char *name)
00629 {
00630 const char *fit_n_str;
00631 fit_n_str = fors_photometry_parameter_get_string(
00632 parameters,
00633 name);
00634 if (fit_n_str == NULL)
00635 {
00636 cpl_error_set_message( cpl_func,
00637 CPL_ERROR_ILLEGAL_INPUT,
00638 "parameter %s not found",
00639 name);
00640 return -1;
00641 }
00642
00643 if (strcmp(fit_n_str, fors_fit_ncoeff_paropts.no) == 0)
00644 return FORS_FIT_NCOEFF_NO;
00645 else if (strcmp(fit_n_str, fors_fit_ncoeff_paropts.one) == 0)
00646 return FORS_FIT_NCOEFF_ONE;
00647 else if (strcmp(fit_n_str, fors_fit_ncoeff_paropts.perframe)== 0)
00648 return FORS_FIT_NCOEFF_PERFRAME;
00649 else if (strcmp(fit_n_str, fors_fit_ncoeff_paropts.pernight)== 0)
00650 return FORS_FIT_NCOEFF_PERNIGHT;
00651 else
00652 {
00653 cpl_error_set_message( cpl_func,
00654 CPL_ERROR_ILLEGAL_INPUT,
00655 "unknown parameter value \"%s\" "
00656 "for %s",
00657 fit_n_str,
00658 name);
00659 return -1;
00660 }
00661 }
00662
00663
00668
00669 void fors_photometry_define_parameters( cpl_parameterlist *parameters)
00670 {
00671 const char *context = cpl_sprintf("fors.%s", fors_photometry_name);
00672
00673 const char *name, *full_name;
00674 cpl_parameter *p;
00675
00676 name = "fitz";
00677 full_name = cpl_sprintf("%s.%s", context, name);
00678 p = cpl_parameter_new_value(full_name,
00679 CPL_TYPE_BOOL,
00680 "Fit zeropoint",
00681 context,
00682 true);
00683 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
00684 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00685 cpl_parameterlist_append(parameters, p);
00686 cpl_free((void *)full_name); full_name = NULL;
00687
00688 name = "fit_all_mag";
00689 full_name = cpl_sprintf("%s.%s", context, name);
00690 p = cpl_parameter_new_value(full_name,
00691 CPL_TYPE_BOOL,
00692 "Always fit star magnitudes",
00693 context,
00694 false);
00695 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
00696 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00697 cpl_parameterlist_append(parameters, p);
00698 cpl_free((void *)full_name); full_name = NULL;
00699
00700 name = "fite";
00701 full_name = cpl_sprintf("%s.%s", context, name);
00702 p = cpl_parameter_new_enum( full_name,
00703 CPL_TYPE_STRING,
00704 "Fit atmospheric extinctions",
00705 context,
00706 fors_fit_ncoeff_paropts.pernight,
00707 4,
00708 fors_fit_ncoeff_paropts.no,
00709 fors_fit_ncoeff_paropts.one,
00710 fors_fit_ncoeff_paropts.perframe,
00711 fors_fit_ncoeff_paropts.pernight);
00712 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
00713 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00714 cpl_parameterlist_append(parameters, p);
00715 cpl_free((void *)full_name); full_name = NULL;
00716
00717 name = "fitc";
00718 full_name = cpl_sprintf("%s.%s", context, name);
00719 p = cpl_parameter_new_value(full_name,
00720 CPL_TYPE_BOOL,
00721 "Fit color correction term",
00722 context,
00723 false);
00724 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
00725 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00726 cpl_parameterlist_append(parameters, p);
00727 cpl_free((void *)full_name); full_name = NULL;
00728
00729 name = "use_all_stars";
00730 full_name = cpl_sprintf("%s.%s", context, name);
00731 p = cpl_parameter_new_value(full_name,
00732 CPL_TYPE_BOOL,
00733 "Use also non-standard stars to fit "
00734 "polynomial f",
00735 context,
00736 false);
00737 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
00738 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00739 cpl_parameterlist_append(parameters, p);
00740 cpl_free((void *)full_name); full_name = NULL;
00741
00742 name = "degreef1";
00743 full_name = cpl_sprintf("%s.%s", context, name);
00744 p = cpl_parameter_new_value(full_name,
00745 CPL_TYPE_INT,
00746 "FLatfield correction map polynomial degree "
00747 "(x)",
00748 context,
00749 0);
00750 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
00751 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00752 cpl_parameterlist_append(parameters, p);
00753 cpl_free((void *)full_name); full_name = NULL;
00754
00755 name = "degreef2";
00756 full_name = cpl_sprintf("%s.%s", context, name);
00757 p = cpl_parameter_new_value(full_name,
00758 CPL_TYPE_INT,
00759 "Flatfield correction map polynomial degree "
00760 "(y), or negative for "
00761 "triangular coefficient matrix",
00762 context,
00763 -1);
00764 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
00765 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00766 cpl_parameterlist_append(parameters, p);
00767 cpl_free((void *)full_name); full_name = NULL;
00768
00769 name = "degreep";
00770 full_name = cpl_sprintf("%s.%s", context, name);
00771 p = cpl_parameter_new_value(full_name,
00772 CPL_TYPE_INT,
00773 "Extinction/color coupling degree",
00774 context,
00775 0);
00776 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
00777 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00778 cpl_parameterlist_append(parameters, p);
00779 cpl_free((void *)full_name); full_name = NULL;
00780
00781 cpl_free((void *)context);
00782
00783 return;
00784 }
00785
00786
00787 #undef cleanup
00788 #define cleanup \
00789 do { \
00790 cpl_polynomial_delete(polyf); polyf = NULL; \
00791 } while (0)
00792
00803
00804 static cpl_polynomial*
00805 fors_photometry_define_polyf( int degreef1,
00806 int degreef2)
00807 {
00808 int xpow,
00809 ypow;
00810 cpl_polynomial *polyf = NULL;
00811
00812
00813 cleanup;
00814
00815
00816 assure(!cpl_error_get_code(), return NULL, "Previous error caught.");
00817
00818 if (degreef1 < 0)
00819 {
00820 cpl_error_set_message( cpl_func, CPL_ERROR_ILLEGAL_INPUT,
00821 "!(degreef1 >= 0)");
00822 return NULL;
00823 }
00824
00825
00826 polyf = cpl_polynomial_new(2);
00827 for (xpow = 0; xpow <= degreef1; xpow++)
00828 {
00829 for (ypow = 0;
00830 (degreef2 >= 0) ? (ypow <= degreef2) : (xpow + ypow <= degreef1);
00831 ypow++)
00832 {
00833 if (xpow+ypow > 0)
00834 {
00835 int pows[2];
00836 pows[0] = xpow;
00837 pows[1] = ypow;
00838 cpl_polynomial_set_coeff(polyf, pows, 1.0);
00839 }
00840 }
00841 }
00842
00843 fors_polynomial_dump(polyf, "polynomial definition f", CPL_MSG_DEBUG, NULL);
00844
00845
00846 if (degreef2 >= 0)
00847 {
00848 cassure(fors_polynomial_count_coeff(polyf)
00849 == (degreef1+1)*(degreef2+1)-1,
00850 CPL_ERROR_UNSPECIFIED,
00851 return polyf,
00852 "Consistency check for rectangular f polynomial failed");
00853 }
00854 else
00855 {
00856 cassure(fors_polynomial_count_coeff(polyf)
00857 == ((degreef1+1)*(degreef1+2))/2-1,
00858 CPL_ERROR_UNSPECIFIED,
00859 return polyf,
00860 "Consistency check for triangular f polynomial failed");
00861 }
00862
00863 return polyf;
00864 }
00865
00866
00867 #undef cleanup
00868 #define cleanup \
00869 do { \
00870 cpl_polynomial_delete(polyp); polyp = NULL; \
00871 } while (0)
00872
00883
00884 static cpl_polynomial*
00885 fors_photometry_define_polyp( int degreep)
00886 {
00887 int k,
00888 l;
00889 cpl_polynomial *polyp = NULL;
00890
00891
00892 cleanup;
00893
00894
00895 assure(!cpl_error_get_code(), return NULL, "Previous error caught.");
00896
00897 cassure( degreep >= 0,
00898 CPL_ERROR_ILLEGAL_INPUT,
00899 return polyp,
00900 "!(degreep >= 0)");
00901
00902
00903 polyp = cpl_polynomial_new(2);
00904
00905 for (k = 0; k <= degreep; k++) {
00906 for (l = 0; k + l <= degreep; l++)
00907 if (k+l > 1) {
00908 int pows[2];
00909 pows[0] = k;
00910 pows[1] = l;
00911 cpl_polynomial_set_coeff(polyp, pows, 1.0);
00912 }
00913 }
00914
00915 fors_polynomial_dump(polyp, "polynomial definition p", CPL_MSG_DEBUG, NULL);
00916
00917
00918 {
00919 int npars = degreep >= 2 ? ((degreep+1)*(degreep+2))/2-3 : 0;
00920 cassure(fors_polynomial_count_coeff(polyp) == npars,
00921 CPL_ERROR_UNSPECIFIED,
00922 return polyp,
00923 "Consistency check for triangular p polynomial failed");
00924 }
00925
00926 return polyp;
00927 }
00928
00929
00930 #undef cleanup
00931 #define cleanup \
00932 do { \
00933 if (poly != NULL) { cpl_polynomial_delete(*poly); *poly = NULL; } \
00934 if (var_poly != NULL) {cpl_polynomial_delete(*var_poly); *var_poly = NULL;}\
00935 } while (0)
00936
00950
00951 static cpl_error_code
00952 fors_photometry_poly_new_from_coefficients( const cpl_polynomial *p_def,
00953 const cpl_matrix *coeffs,
00954 const cpl_matrix *cov_coeffs,
00955 cpl_polynomial **poly,
00956 cpl_polynomial **var_poly)
00957 {
00958 int n_coeffs;
00959 cpl_errorstate errstat = cpl_errorstate_get();
00960
00961
00962 cleanup;
00963
00964 cassure_automsg( p_def != NULL,
00965 CPL_ERROR_NULL_INPUT,
00966 return cpl_error_get_code());
00967 cassure_automsg( poly != NULL,
00968 CPL_ERROR_NULL_INPUT,
00969 return cpl_error_get_code());
00970
00971 n_coeffs = fors_polynomial_count_coeff( p_def);
00972 cassure_automsg( n_coeffs == 0 || coeffs != NULL,
00973 CPL_ERROR_NULL_INPUT,
00974 return cpl_error_get_code());
00975 cassure_automsg( n_coeffs == 0
00976 || var_poly == NULL
00977 || cov_coeffs != NULL,
00978 CPL_ERROR_NULL_INPUT,
00979 return cpl_error_get_code());
00980 if (n_coeffs > 0)
00981 {
00982
00983 cassure_automsg( cpl_matrix_get_ncol(coeffs) == 1,
00984 CPL_ERROR_ILLEGAL_INPUT,
00985 return cpl_error_get_code());
00986 cassure_automsg( cpl_matrix_get_nrow(coeffs)
00987 == n_coeffs,
00988 CPL_ERROR_INCOMPATIBLE_INPUT,
00989 return cpl_error_get_code());
00990 if (var_poly != NULL)
00991 {
00992 cassure_automsg( cpl_matrix_get_nrow(cov_coeffs)
00993 == n_coeffs,
00994 CPL_ERROR_INCOMPATIBLE_INPUT,
00995 return cpl_error_get_code());
00996 cassure( cpl_matrix_get_nrow(cov_coeffs)
00997 == cpl_matrix_get_ncol(cov_coeffs),
00998 CPL_ERROR_INCOMPATIBLE_INPUT,
00999 return cpl_error_get_code(),
01000 "cov_coeffs is not square");
01001 }
01002 *poly = cpl_polynomial_duplicate( p_def);
01003 fors_polynomial_set_existing_coeff( *poly,
01004 cpl_matrix_get_data_const(coeffs),
01005 n_coeffs);
01006 passure(cpl_errorstate_is_equal(errstat), return cpl_error_get_code());
01007
01008 if (var_poly != NULL)
01009 *var_poly = fors_polynomial_create_variance_polynomial(
01010 p_def,
01011 cov_coeffs);
01012 passure(cpl_errorstate_is_equal(errstat), return cpl_error_get_code());
01013 }
01014 else
01015 {
01016 *poly = cpl_polynomial_new( cpl_polynomial_get_dimension(
01017 p_def));
01018 if (var_poly != NULL)
01019 *var_poly = cpl_polynomial_new( cpl_polynomial_get_dimension(
01020 p_def));
01021 }
01022 passure(cpl_errorstate_is_equal(errstat), return cpl_error_get_code());
01023
01024 return (cpl_errorstate_is_equal(errstat) ?
01025 CPL_ERROR_NONE :
01026 cpl_error_get_code());
01027 }
01028
01029
01030 #undef cleanup
01031 #define cleanup
01032
01037
01038 int
01039 fors_photometry_get_timezone_observer( const cpl_propertylist *header)
01040 {
01041 const cpl_property *prop;
01042
01043 cassure_automsg( header != NULL,
01044 CPL_ERROR_NULL_INPUT,
01045 return 0);
01046
01047 do {
01048 const char *origin;
01049
01050 prop = cpl_propertylist_get_property_const(header, "ORIGIN");
01051 if (prop == NULL)
01052 {
01053 cpl_error_set_message( cpl_func,
01054 CPL_ERROR_DATA_NOT_FOUND,
01055 "Couldn't find the keyword ORIGIN");
01056 return 0;
01057 }
01058
01059 if (cpl_property_get_type(prop) != CPL_TYPE_STRING)
01060 break;
01061
01062 if ((origin = cpl_property_get_string(prop)) == NULL)
01063 break;
01064
01065 if (!fors_fits_compare_string(origin, "ESO"))
01066 break;
01067
01068
01069 return -3;
01070
01071 } while (0);
01072
01073 cpl_error_set_message( cpl_func,
01074 CPL_ERROR_ILLEGAL_INPUT,
01075 "Don't know the originator of the "
01076 "frame specified in ORIGIN");
01077 return 0;
01078 }
01079
01080
01081 #undef cleanup
01082 #define cleanup
01083
01090
01091 int
01092 fors_photometry_get_night_id( const cpl_propertylist *header)
01093 {
01094 const cpl_property *prop;
01095 cpl_errorstate errstat = cpl_errorstate_get();
01096
01097 cassure_automsg( header != NULL,
01098 CPL_ERROR_NULL_INPUT,
01099 return 0);
01100
01101
01102 prop = cpl_propertylist_get_property_const(header, "MJD-OBS");
01103 if (prop != NULL)
01104 {
01105 double mjd,
01106 jd,
01107 timezone;
01108 int localstartday;
01109 mjd = fors_property_get_num(prop);
01110 assure( cpl_errorstate_is_equal(errstat),
01111 return 0,
01112 "Could not interprete Modified "
01113 "Julian Date keyword MJD-OBS");
01114
01115
01116
01117
01118
01119
01120
01121
01122 jd = mjd + 2400000.5;
01123
01124
01125
01126
01127
01128 timezone = fors_photometry_get_timezone_observer(header);
01129
01130
01131
01132 jd += (double)timezone / 24.0;
01133
01134
01135
01136 localstartday = floor(jd);
01137 cpl_msg_debug( cpl_func,
01138 "Julian day no. of observation "
01139 "night: %d",
01140 localstartday);
01141
01142 return localstartday;
01143 }
01144
01145
01146
01147 cpl_error_set_message( cpl_func,
01148 CPL_ERROR_DATA_NOT_FOUND,
01149 "Couldn't find the keyword "
01150 "MJD-OBS");
01151 return 0;
01152 }
01153
01154
01155 #undef cleanup
01156 #define cleanup \
01157 do { \
01158 cpl_free(ident_array); ident_array = NULL; \
01159 } while (0)
01160
01168
01169 static int
01170 fors_photometry_atm_ext_create_index_by_identifier(
01171 entry_list *obs_list)
01172 {
01173 entry *e;
01174 int *ident_array;
01175 int n_entries,
01176 n_idents = 0;
01177 cpl_errorstate errstat = cpl_errorstate_get();
01178
01179 cassure_automsg( obs_list != NULL,
01180 CPL_ERROR_NULL_INPUT,
01181 return 0);
01182
01183 n_entries = entry_list_size(obs_list);
01184 ident_array = cpl_malloc(n_entries * sizeof(*ident_array));
01185
01186 for ( e = entry_list_first(obs_list);
01187 e != NULL;
01188 e = entry_list_next(obs_list))
01189 {
01190 int i;
01191 bool found = false;
01192 for (i = 0; i < n_idents && !found; i++)
01193 {
01194 if (e->atm_ext_identifier == ident_array[i])
01195 {
01196 found = true;
01197 e->atm_ext_index = i;
01198 }
01199 }
01200 if (!found)
01201 {
01202 ident_array[n_idents] = e->atm_ext_identifier;
01203 e->atm_ext_index = n_idents;
01204
01205 cpl_msg_debug( cpl_func,
01206 "Creating atm. extinction index "
01207 "%2d for identifier %d",
01208 n_idents,
01209 ident_array[n_idents]);
01210
01211 n_idents++;
01212 }
01213 }
01214
01215 passure( cpl_errorstate_is_equal(errstat),
01216 return 0);
01217
01218 cpl_free(ident_array);
01219 return n_idents;
01220 }
01221
01222
01223 #undef cleanup
01224 #define cleanup
01225
01241
01242 static int
01243 fors_photometry_atm_ext_create_indices( entry_list *obsl,
01244 fors_fit_ncoeff fit_e)
01245 {
01246 entry *e;
01247 int n_atm_ext_indices = 0;
01248 cpl_errorstate errstat = cpl_errorstate_get();
01249
01250 cassure_automsg( obsl != NULL,
01251 CPL_ERROR_NULL_INPUT,
01252 return -1);
01253
01254 if (fit_e != FORS_FIT_NCOEFF_NO
01255 && fit_e != FORS_FIT_NCOEFF_ONE
01256 && fit_e != FORS_FIT_NCOEFF_PERFRAME)
01257 {
01258
01259
01260 n_atm_ext_indices =
01261 fors_photometry_atm_ext_create_index_by_identifier(obsl);
01262 }
01263 else
01264 {
01265 if (fit_e == FORS_FIT_NCOEFF_NO)
01266 {
01267 for (e = entry_list_first(obsl); e!=NULL; e = entry_list_next(obsl))
01268 e->atm_ext_index = -1;
01269 n_atm_ext_indices = 0;
01270 }
01271 else if (fit_e == FORS_FIT_NCOEFF_ONE)
01272 {
01273 for (e = entry_list_first(obsl); e!=NULL; e = entry_list_next(obsl))
01274 e->atm_ext_index = 0;
01275 n_atm_ext_indices = 1;
01276 }
01277 else if (fit_e == FORS_FIT_NCOEFF_PERFRAME)
01278 {
01279 for (e = entry_list_first(obsl); e!=NULL; e = entry_list_next(obsl))
01280 {
01281 e->atm_ext_index = e->frame_index;
01282 if (e->frame_index >= n_atm_ext_indices)
01283 n_atm_ext_indices = e->frame_index + 1;
01284 }
01285 }
01286 }
01287
01288 if (!cpl_errorstate_is_equal(errstat))
01289 cpl_error_set_where(cpl_func);
01290
01291 return (cpl_errorstate_is_equal(errstat) ? n_atm_ext_indices : -1);
01292 }
01293
01294
01295 #undef cleanup
01296 #define cleanup \
01297 do { \
01298 if (frame_printed != NULL) \
01299 { \
01300 cpl_free(frame_printed); \
01301 frame_printed = NULL; \
01302 } \
01303 } while (0)
01304
01310
01311
01312 static cpl_table *
01313 fors_photometry_atm_ext_print_index_by_framename(
01314 const entry_list *obs_list,
01315 const cpl_frameset *frames)
01316 {
01317 const entry *e;
01318 bool *frame_printed = NULL;
01319 int n_frames,
01320 ext_index,
01321 max_ext_index,
01322 n;
01323 int row = 0;
01324 cpl_table *summary;
01325 cpl_errorstate errstat = cpl_errorstate_get();
01326
01327 cassure_automsg( obs_list != NULL,
01328 CPL_ERROR_NULL_INPUT,
01329 return NULL);
01330 cassure_automsg( frames != NULL,
01331 CPL_ERROR_NULL_INPUT,
01332 return NULL);
01333
01334 n_frames = cpl_frameset_get_size(frames);
01335 frame_printed = cpl_malloc(n_frames * sizeof(*frame_printed));
01336
01337 summary = cpl_table_new(n_frames);
01338 cpl_table_new_column(summary, "filename", CPL_TYPE_STRING);
01339 cpl_table_new_column(summary, "index", CPL_TYPE_INT);
01340
01341 max_ext_index = -1;
01342 for ( e = entry_list_first_const(obs_list);
01343 e != NULL;
01344 e = entry_list_next_const(obs_list))
01345 {
01346 if (e->atm_ext_index > max_ext_index)
01347 max_ext_index = e->atm_ext_index;
01348 }
01349
01350 if (max_ext_index >= 0)
01351 cpl_msg_info( cpl_func,
01352 "Assignment of atmospheric "
01353 "extinction indices:");
01354
01355 for (ext_index = 0; ext_index <= max_ext_index; ext_index++)
01356 {
01357 bool first_file = true;
01358 char estr[15];
01359 for (n = 0; n < n_frames; n++)
01360 frame_printed[n] = false;
01361
01362 cpl_msg_indent_more();
01363 sprintf(estr, "E_%d: ", ext_index);
01364 estr[9] = '\0';
01365
01366 for ( e = entry_list_first_const(obs_list);
01367 e != NULL;
01368 e = entry_list_next_const(obs_list))
01369 {
01370 if (e->atm_ext_index == ext_index)
01371 {
01372 cassure_automsg( e->frame_index >= 0
01373 && e->frame_index < n_frames,
01374 CPL_ERROR_ILLEGAL_INPUT,
01375 {
01376 cpl_msg_indent_less();
01377 return NULL;
01378 });
01379
01380 if (!frame_printed[e->frame_index])
01381 {
01382 const cpl_frame *f;
01383
01384 f = cpl_frameset_get_frame_const(frames, e->frame_index);
01385 if (first_file)
01386 {
01387 cpl_msg_info( cpl_func,
01388 "%s%s",
01389 estr,
01390 cpl_frame_get_filename(f));
01391 }
01392 else
01393 {
01394 cpl_msg_info( cpl_func,
01395 " %s",
01396 cpl_frame_get_filename(f));
01397 }
01398 frame_printed[e->frame_index] = true;
01399 first_file = false;
01400
01401
01402
01403
01404
01405
01406 cpl_table_set_string(summary, "filename", row,
01407 cpl_frame_get_filename(f));
01408 cpl_table_set_int(summary, "index", row, ext_index);
01409 row++;
01410
01411 }
01412 }
01413
01414
01415 }
01416 cpl_msg_indent_less();
01417 }
01418
01419
01420
01421 cleanup;
01422
01423 passure( cpl_errorstate_is_equal(errstat),
01424 return NULL);
01425
01426 return summary;
01427 }
01428
01429
01430 #undef cleanup
01431 #define cleanup
01432
01445
01446 static cpl_error_code
01447 fors_photometry_check_input_value( double value,
01448 double value_error,
01449 const char *value_name,
01450 const char *input_name,
01451 double min_limit,
01452 double max_limit,
01453 double max_error)
01454 {
01455 cpl_errorstate errstat = cpl_errorstate_get();
01456
01457 if (value < min_limit || value > max_limit)
01458 {
01459 cpl_error_set_message( cpl_func,
01460 CPL_ERROR_ILLEGAL_INPUT,
01461 "invalid %s (%f)"
01462 "%s%s"
01463 ", either correct input, or try to "
01464 "re-run this recipe with fitting "
01465 "%s enabled",
01466 value_name,
01467 value,
01468 (input_name != NULL)?" read from ":"",
01469 (input_name != NULL)?input_name:"",
01470 value_name);
01471 }
01472 else if (value_error > max_error || value_error < 0)
01473 {
01474 char exceed_max_err[30];
01475 sprintf(exceed_max_err, "> %f", max_error);
01476 cpl_error_set_message( cpl_func,
01477 CPL_ERROR_ILLEGAL_INPUT,
01478 "unreliable %s ((error = %g) %s)"
01479 "%s%s"
01480 ", either recompute input, or try "
01481 "to re-run this recipe with "
01482 "fitting %s enabled",
01483 value_name,
01484 value_error,
01485 (value_error < 0) ?
01486 "< 0" : exceed_max_err,
01487 (input_name != NULL)?" read from ":"",
01488 (input_name != NULL)?input_name:"",
01489 value_name);
01490 }
01491 else
01492 {
01493 cpl_msg_info( cpl_func,
01494 "Using input value%s%s: "
01495 "%s = %f +- %f",
01496 (input_name != NULL)?" from ":"",
01497 (input_name != NULL)?input_name:"",
01498 value_name,
01499 value,
01500 value_error);
01501 }
01502
01503 return (cpl_errorstate_is_equal(errstat) ?
01504 CPL_ERROR_NONE :
01505 cpl_error_get_code());
01506 }
01507
01508
01509 #undef cleanup
01510 #define cleanup \
01511 do { \
01512 cpl_array_delete(airmasses); airmasses = NULL; \
01513 } while (0)
01514
01520
01521 static cpl_error_code
01522 fors_photometry_check_fitparam_atm_ext( entry_list *obsl,
01523 fors_fit_ncoeff fit_e,
01524 bool fit_z)
01525 {
01526 entry *e;
01527 int n_atm_ext_indices = 0;
01528 cpl_array *airmasses = NULL;
01529 cpl_errorstate errstat = cpl_errorstate_get();
01530
01531 cassure_automsg( obsl != NULL,
01532 CPL_ERROR_NULL_INPUT,
01533 return cpl_error_get_code());
01534
01535
01536
01537
01538 if (!fit_z || fit_e == FORS_FIT_NCOEFF_NO)
01539 {
01540 return CPL_ERROR_NONE;
01541 }
01542
01543
01544 for (e = entry_list_first(obsl); e!=NULL; e = entry_list_next(obsl))
01545 {
01546 if (e->atm_ext_index >= n_atm_ext_indices)
01547 n_atm_ext_indices = e->atm_ext_index + 1;
01548 }
01549
01550
01551
01552
01553
01554
01555
01556
01557 if (n_atm_ext_indices > 0)
01558 {
01559 bool multiple_found = false;
01560
01561 airmasses = cpl_array_new(n_atm_ext_indices, CPL_TYPE_DOUBLE);
01562
01563 for ( e = entry_list_first(obsl);
01564 e != NULL;
01565 e = entry_list_next(obsl))
01566 {
01567 double first_airmass;
01568 int is_set;
01569 first_airmass = cpl_array_get_double(
01570 airmasses,
01571 e->atm_ext_index,
01572 &is_set);
01573 passure( cpl_errorstate_is_equal(errstat),
01574 return cpl_error_get_code());
01575 is_set = (is_set == 0);
01576 if (!is_set)
01577 cpl_array_set_double(airmasses, e->atm_ext_index, e->airmass);
01578 else
01579 {
01580
01581 if (fabs(e->airmass - first_airmass) > 10*DBL_EPSILON)
01582 {
01583 multiple_found = true;
01584 break;
01585 }
01586
01587
01588
01589 }
01590 }
01591
01592 if (!multiple_found)
01593 {
01594 if (n_atm_ext_indices > 1)
01595 cpl_msg_error( cpl_func,
01596 "No atmospheric extinction was "
01597 "observed at different airmasses.");
01598 else
01599 cpl_msg_error( cpl_func,
01600 "Atmospheric extinction was not "
01601 "observed at different airmasses.");
01602 }
01603 cassure( multiple_found,
01604 CPL_ERROR_ILLEGAL_INPUT,
01605 return cpl_error_get_code(),
01606 "For fitting the zeropoint and "
01607 "atmospheric extinction, "
01608 "there must be >= 2 different "
01609 "airmasses for at least 1 "
01610 "atmospheric extinction");
01611 }
01612
01613 cpl_array_delete(airmasses);
01614
01615 return (cpl_errorstate_is_equal(errstat) ?
01616 CPL_ERROR_NONE :
01617 cpl_error_get_code());
01618 }
01619
01620
01621
01622 static void
01623 myprintf(const char *format, ...)
01624 {
01625 va_list al;
01626
01627 va_start(al, format);
01628
01629 va_end(al);
01630
01631 return;
01632 }
01633
01634
01635
01636
01637
01638 static cpl_matrix * matrix_product_normal_create(const cpl_matrix * self)
01639 {
01640 double sum;
01641 cpl_matrix * product;
01642 const double * ai = cpl_matrix_get_data_const(self);
01643 const double * aj;
01644 double * bwrite;
01645 const int m = cpl_matrix_get_nrow(self);
01646 const int n = cpl_matrix_get_ncol(self);
01647 int i, j, k;
01648
01649
01650 cpl_ensure(self != NULL, CPL_ERROR_NULL_INPUT, NULL);
01651
01652 #if 0
01653
01654
01655
01656 product = cpl_matrix_new(m, m);
01657 bwrite = cpl_matrix_get_data(product);
01658 #else
01659 bwrite = (double *) cpl_malloc(m * m * sizeof(double));
01660 product = cpl_matrix_wrap(m, m, bwrite);
01661 #endif
01662
01663
01664 for (i = 0; i < m; i++, bwrite += m, ai += n) {
01665 aj = ai;
01666 for (j = i; j < m; j++, aj += n) {
01667 sum = 0.0;
01668 for (k = 0; k < n; k++) {
01669 sum += ai[k] * aj[k];
01670 }
01671 bwrite[j] = sum;
01672 }
01673 }
01674
01675 return product;
01676 }
01677
01678
01679 #undef cleanup
01680 #define cleanup \
01681 do { \
01682 fors_matrix_null(&solution); \
01683 fors_matrix_null(&cov1); \
01684 fors_matrix_null(&At); \
01685 fors_matrix_null(&AtC); \
01686 fors_matrix_null(&AtCA); \
01687 fors_matrix_null(&p); \
01688 fors_matrix_null(&Ap); \
01689 fors_matrix_null(&bAp); \
01690 fors_matrix_null(&bApt); \
01691 fors_matrix_null(&C1bAp); \
01692 fors_matrix_null(&chi2); \
01693 } while (0)
01694
01695
01696
01697
01698
01699
01700
01701
01702
01703
01704
01705
01706
01707
01708
01709
01710
01711
01712
01713
01714
01715
01716
01717
01718
01719
01720
01721
01722
01723
01724
01725
01726
01727
01728
01729
01730 static cpl_matrix *
01731 solve_normal(const cpl_matrix *coeff,
01732 const cpl_matrix *rhs,
01733 const cpl_matrix *cov_rhs,
01734 double *red_chisq)
01735 {
01736 cpl_matrix *solution = NULL;
01737 cpl_matrix *cov1 = NULL;
01738 cpl_matrix *At = NULL;
01739 cpl_matrix *AtC = NULL;
01740 cpl_matrix *AtCA = NULL;
01741 cpl_matrix *p = NULL;
01742 cpl_matrix *Ap = NULL;
01743 cpl_matrix *bAp = NULL;
01744 cpl_matrix *bApt = NULL;
01745 cpl_matrix *C1bAp = NULL;
01746 cpl_matrix *chi2 = NULL;
01747 cpl_error_code error;
01748
01749 cpl_ensure(coeff != NULL, CPL_ERROR_NULL_INPUT, NULL);
01750 cpl_ensure(rhs != NULL, CPL_ERROR_NULL_INPUT, NULL);
01751 cpl_ensure(cov_rhs != NULL, CPL_ERROR_NULL_INPUT, NULL);
01752
01753
01754
01755 cov1 = cpl_matrix_invert_create(cov_rhs);
01756
01757 assure( cov1 != NULL, return NULL,
01758 "Could not invert covariance matrix. Make sure that provided "
01759 "errors are positive");
01760
01761
01762
01763
01764
01765
01766
01767 At = cpl_matrix_transpose_create(coeff);
01768
01769 AtC = cpl_matrix_product_create(At, cov1);
01770
01771 fors_matrix_null(&At);
01772
01773 AtCA = cpl_matrix_product_create(AtC, coeff);
01774
01775 solution = cpl_matrix_product_create(AtC, rhs);
01776 cpl_matrix_set_size(solution,
01777 cpl_matrix_get_nrow(solution),
01778 1 + cpl_matrix_get_nrow(solution));
01779 {
01780 int i = 0;
01781 for (i = 0; i < cpl_matrix_get_nrow(solution); i++) {
01782 cpl_matrix_set(solution, i, i+1, 1);
01783 }
01784 }
01785
01786 fors_matrix_null(&AtC);
01787
01788
01789
01790
01791 error = cpl_matrix_decomp_chol(AtCA);
01792 if (!error) {
01793 error = cpl_matrix_solve_chol(AtCA, solution);
01794 }
01795
01796 fors_matrix_null(&AtCA);
01797
01798 if (error) {
01799 cleanup;
01800 cpl_ensure(0, error, NULL);
01801 }
01802
01803
01804 if (red_chisq != NULL) {
01805
01806
01807 p = cpl_matrix_duplicate(solution);
01808 cpl_matrix_set_size(p, cpl_matrix_get_nrow(p), 1);
01809
01810 Ap = cpl_matrix_product_create(coeff, p);
01811
01812 bAp = cpl_matrix_duplicate(rhs);
01813 cpl_matrix_subtract(bAp, Ap);
01814
01815 bApt = cpl_matrix_transpose_create(bAp);
01816
01817 C1bAp = cpl_matrix_product_create(cov1, bAp);
01818
01819 chi2 = cpl_matrix_product_create(bApt, C1bAp);
01820
01821 passure(cpl_matrix_get_nrow(chi2) == 1 &&
01822 cpl_matrix_get_ncol(chi2) == 1, return NULL);
01823
01824 *red_chisq = cpl_matrix_get(chi2, 0, 0) /
01825 (cpl_matrix_get_nrow(rhs) - cpl_matrix_get_nrow(p));
01826
01827 }
01828
01829 cpl_matrix *return_solution = solution; solution = NULL;
01830 cleanup;
01831
01832 return return_solution;
01833 }
01834
01835
01836 #undef cleanup
01837 #define cleanup \
01838 do { \
01839 cpl_frameset_delete((cpl_frameset *)aligned_phot_frames); \
01840 cpl_frameset_delete((cpl_frameset *)master_flat_frame); \
01841 cpl_frameset_delete((cpl_frameset *)phot_table); \
01842 fors_setting_delete(&setting); \
01843 fors_image_delete(&master_flat); \
01844 fors_image_delete(&correction); \
01845 cpl_table_delete(phot_coeff); \
01846 cpl_table_delete(summary); \
01847 fors_delete_star_lists(&obs, &std_star_list); \
01848 cpl_array_delete(n_std_star_obs); n_std_star_obs = NULL; \
01849 fors_matrix_null(&eqn_lhs); \
01850 fors_matrix_null(&eqn_rhs); \
01851 fors_matrix_null(&eqn_cov_rhs); \
01852 fors_matrix_null(&eqn_result); \
01853 fors_matrix_null(&tmp_mat); \
01854 fors_matrix_null(&result_polyf); \
01855 fors_matrix_null(&result_cov_polyf); \
01856 fors_matrix_null(&result_params); \
01857 fors_matrix_null(&result_cov_params); \
01858 fors_matrix_null(&result_polyp); \
01859 fors_matrix_null(&result_cov_polyp); \
01860 cpl_polynomial_delete(polyf); polyf = NULL; \
01861 cpl_polynomial_delete(polyf_definition); polyf_definition = NULL; \
01862 cpl_polynomial_delete(polyf_variance); polyf_variance = NULL; \
01863 cpl_polynomial_delete(polyp); polyp = NULL; \
01864 cpl_polynomial_delete(polyp_definition); polyp_definition = NULL; \
01865 } while (0)
01866
01875
01876 void fors_photometry(cpl_frameset *frames, const cpl_parameterlist *parameters)
01877 {
01878
01879 const cpl_frameset *aligned_phot_frames = NULL,
01880 *master_flat_frame = NULL,
01881 *phot_table = NULL;
01882 fors_image *master_flat = NULL;
01883
01884
01885 fors_image *correction = NULL;
01886 cpl_table *phot_coeff = NULL;
01887 cpl_table *summary = NULL;
01888
01889
01890 fors_std_star_list *std_star_list = NULL;
01891 entry_list *obs = NULL;
01892 cpl_array *n_std_star_obs = NULL;
01893
01894
01895 cpl_matrix *eqn_lhs = NULL,
01896 *eqn_rhs = NULL,
01897 *eqn_cov_rhs = NULL,
01898 *eqn_result = NULL,
01899 *result_polyf = NULL,
01900 *result_cov_polyf = NULL,
01901 *result_params = NULL,
01902 *result_cov_params = NULL,
01903 *result_polyp = NULL,
01904 *result_cov_polyp = NULL,
01905 *tmp_mat = NULL;
01906
01907
01908 cpl_polynomial *polyf = NULL,
01909 *polyf_definition = NULL,
01910 *polyf_variance = NULL,
01911 *polyp = NULL,
01912 *polyp_definition = NULL;
01913
01914
01915 fors_setting *setting = NULL;
01916 const char *tag = NULL;
01917 int row;
01918
01919
01920 cpl_propertylist *qc = NULL;
01921 double qc_zeropoint = -1.0;
01922 double qc_zeropoint_err = -1.0;
01923 double qc_extinction = -1.0;
01924 double qc_extinction_err = -1.0;
01925 double qc_colorterm = -1.0;
01926 double qc_colorterm_err = -1.0;
01927
01928
01929
01930 aligned_phot_frames = fors_frameset_extract(frames, ALIGNED_PHOT);
01931 cassure(cpl_frameset_get_size(aligned_phot_frames) > 0,
01932 CPL_ERROR_DATA_NOT_FOUND,
01933 return,
01934 "No %s provided", ALIGNED_PHOT);
01935
01936 master_flat_frame = fors_frameset_extract(frames, MASTER_SKY_FLAT_IMG);
01937 cassure(cpl_frameset_get_size(master_flat_frame) > 0,
01938 CPL_ERROR_DATA_NOT_FOUND,
01939 return,
01940 "No %s provided", MASTER_SKY_FLAT_IMG);
01941
01942 phot_table = fors_frameset_extract(frames, PHOT_TABLE);
01943 cassure(cpl_frameset_get_size(phot_table) == 1,
01944 CPL_ERROR_DATA_NOT_FOUND,
01945 return,
01946 "One %s required. %d found",
01947 PHOT_TABLE, cpl_frameset_get_size(phot_table));
01948
01949
01950
01951 bool fit_z,
01952 override_fit_m,
01953 fit_c,
01954 use_all;
01955 int degreef1,
01956 degreef2,
01957 degreep;
01958 fors_fit_ncoeff fit_e;
01959
01960 degreef1 = fors_photometry_parameter_get_num(
01961 parameters,
01962 "degreef1",
01963 CPL_TYPE_INT);
01964 degreef2 = fors_photometry_parameter_get_num(
01965 parameters,
01966 "degreef2",
01967 CPL_TYPE_INT);
01968 degreep = fors_photometry_parameter_get_num(
01969 parameters,
01970 "degreep",
01971 CPL_TYPE_INT);
01972 fit_z = fors_photometry_parameter_get_num(
01973 parameters,
01974 "fitz",
01975 CPL_TYPE_BOOL);
01976 override_fit_m = fors_photometry_parameter_get_num(
01977 parameters,
01978 "fit_all_mag",
01979 CPL_TYPE_BOOL);
01980 fit_c = fors_photometry_parameter_get_num(
01981 parameters,
01982 "fitc",
01983 CPL_TYPE_BOOL);
01984 use_all = fors_photometry_parameter_get_num(
01985 parameters,
01986 "use_all_stars",
01987 CPL_TYPE_BOOL);
01988 fit_e = fors_photometry_parameter_get_ncoeff(
01989 parameters,
01990 "fite");
01991 assure( !cpl_error_get_code(), return, NULL );
01992
01993 if (fit_e == FORS_FIT_NCOEFF_PERFRAME
01994 && fit_z == true)
01995 {
01996 cpl_error_set_message( cpl_func,
01997 CPL_ERROR_ILLEGAL_INPUT,
01998 "Fitting one atmospheric "
01999 "extinction per frame and the "
02000 "zeropoint is ambiguous and "
02001 "therefore not possible");
02002 return;
02003 }
02004
02005
02006
02007
02008 setting = fors_setting_new( cpl_frameset_get_first_const(
02009 aligned_phot_frames));
02010 assure( !cpl_error_get_code(),
02011 return,
02012 "Could not get instrument setting");
02013
02014
02015 struct phot_input {
02016 double color_coeff,
02017 dcolor_coeff,
02018 ext,
02019 dext,
02020 zpoint,
02021 dzpoint;
02022 } phot_input;
02023 cpl_msg_info(cpl_func, "Loading photometry table");
02024 fors_phot_table_load( cpl_frameset_get_first_const(
02025 phot_table),
02026 setting,
02027 &phot_input.color_coeff,
02028 &phot_input.dcolor_coeff,
02029 &phot_input.ext,
02030 &phot_input.dext,
02031 &phot_input.zpoint,
02032 &phot_input.dzpoint);
02033 assure( !cpl_error_get_code(),
02034 return,
02035 "Could not load photometry table");
02036
02037
02038 cpl_msg_indent_more();
02039 if (!fit_c)
02040 {
02041 fors_photometry_check_input_value( phot_input.color_coeff,
02042 phot_input.dcolor_coeff,
02043 "color correction term",
02044 NULL,
02045 -10,
02046 10,
02047 1);
02048 }
02049 if (fit_e == FORS_FIT_NCOEFF_NO)
02050 {
02051 fors_photometry_check_input_value( phot_input.ext,
02052 phot_input.dext,
02053 "atmospheric extinction",
02054 NULL,
02055 0,
02056 5,
02057 1);
02058 }
02059 if (!fit_z)
02060 {
02061 fors_photometry_check_input_value( phot_input.zpoint,
02062 phot_input.dzpoint,
02063 "zeropoint",
02064 NULL,
02065 0,
02066 50,
02067 1);
02068 }
02069 cpl_msg_indent_less();
02070 if (cpl_error_get_code() != CPL_ERROR_NONE)
02071 return;
02072
02073
02074
02075 int n_mag_fits,
02076 n_frames;
02077 int (*get_atm_ext_id_func)(const cpl_propertylist*);
02078
02079 cpl_msg_info( cpl_func,
02080 "Importing %s tables:",
02081 ALIGNED_PHOT);
02082 cpl_msg_indent_more();
02083
02084 switch (fit_e)
02085 {
02086 case FORS_FIT_NCOEFF_PERNIGHT:
02087 get_atm_ext_id_func = fors_photometry_get_night_id;
02088 break;
02089 default:
02090 get_atm_ext_id_func = NULL;
02091 }
02092
02093 obs = fors_photometry_read_input( aligned_phot_frames,
02094 setting,
02095 get_atm_ext_id_func,
02096 use_all,
02097 &n_frames,
02098 &std_star_list,
02099 degreef1 < 1);
02100 assure(!cpl_error_get_code(), return, NULL);
02101
02102 fors_photometry_adjust_fit_mag_flags( std_star_list,
02103 obs,
02104 override_fit_m,
02105 &n_mag_fits);
02106 assure(!cpl_error_get_code(), return, NULL);
02107
02108 fors_photometry_atm_ext_create_indices( obs,
02109 fit_e);
02110 passure(!cpl_error_get_code(), return);
02111
02112 if (fit_e != FORS_FIT_NCOEFF_NO && fit_e != FORS_FIT_NCOEFF_ONE)
02113 {
02114 summary = fors_photometry_atm_ext_print_index_by_framename(
02115 obs,
02116 aligned_phot_frames);
02117 assure(!cpl_error_get_code(), return, "Internal error" );
02118 }
02119
02120 fors_photometry_remove_unnecessary( std_star_list,
02121 obs,
02122 &n_mag_fits);
02123 assure(!cpl_error_get_code(), return, NULL);
02124
02125 fors_photometry_check_fitparam_atm_ext(obs, fit_e, fit_z);
02126 assure(!cpl_error_get_code(), return, NULL);
02127
02128 entry_list_print(obs, CPL_MSG_DEBUG);
02129
02130 {
02131 int ntot,
02132 n_std_stars;
02133 ntot = entry_list_size(obs);
02134 n_std_stars = fors_std_star_list_size(std_star_list);
02135
02136 cpl_msg_info(cpl_func,
02137 "Found %d table%s, %d star%s, %d unique star%s, "
02138 "%d magnitude%s to determine",
02139 n_frames, n_frames != 1 ? "s" : "",
02140 ntot, ntot != 1 ? "s" : "",
02141 n_std_stars, n_std_stars != 1 ? "s" : "",
02142 n_mag_fits, n_mag_fits != 1 ? "s" : "");
02143 }
02144 cpl_msg_indent_less();
02145
02146
02147
02148 int n_coeff_polyf = 0,
02149 n_coeff_polyp = 0,
02150 n_coeff_params = 0,
02151 n_coeff_params_ext = 0;
02152 fors_matrix_null(&eqn_lhs);
02153
02154 polyf_definition = fors_photometry_define_polyf(
02155 degreef1,
02156 degreef2);
02157 polyp_definition = fors_photometry_define_polyp(
02158 degreep);
02159 assure(!cpl_error_get_code(), return, NULL);
02160
02161
02162 tmp_mat = build_equations_lhs_matrix_from_poly(
02163 obs,
02164 polyf_definition,
02165 "f",
02166 &entry_get_powers_x_y);
02167 assure(!cpl_error_get_code(), return, NULL);
02168 if (tmp_mat != NULL) n_coeff_polyf = cpl_matrix_get_ncol(tmp_mat);
02169 fors_matrix_append_delete(&eqn_lhs, &tmp_mat);
02170
02171 tmp_mat = build_equations_lhs_matrix_from_parameters(
02172 obs,
02173 std_star_list,
02174 fit_z,
02175 fit_c,
02176 &n_coeff_params_ext);
02177 assure(!cpl_error_get_code(), return, NULL);
02178 if (tmp_mat != NULL) n_coeff_params = cpl_matrix_get_ncol(tmp_mat);
02179 fors_matrix_append_delete(&eqn_lhs, &tmp_mat);
02180
02181 tmp_mat = build_equations_lhs_matrix_from_poly(
02182 obs,
02183 polyp_definition,
02184 "p",
02185 &entry_get_powers_airmass_color);
02186 assure(!cpl_error_get_code(), return, NULL);
02187 if (tmp_mat != NULL) n_coeff_polyp = cpl_matrix_get_ncol(tmp_mat);
02188 fors_matrix_append_delete(&eqn_lhs, &tmp_mat);
02189
02190
02191 build_equations_rhs_cov( obs,
02192 std_star_list,
02193 fit_z,
02194 fit_c,
02195 (n_coeff_params_ext > 0),
02196 phot_input.color_coeff,
02197 phot_input.dcolor_coeff,
02198 phot_input.ext,
02199 phot_input.dext,
02200 phot_input.zpoint,
02201 phot_input.dzpoint,
02202 &eqn_rhs,
02203 &eqn_cov_rhs);
02204
02205
02206
02207
02208
02209
02210
02211 assure( !cpl_error_get_code(), return, "Could not setup rhs equations" );
02212
02213 int n_parameters = cpl_matrix_get_ncol(eqn_lhs);
02214 fors_matrix_null(&tmp_mat);
02215
02216
02217
02218 {
02219 int eqn_nrows = cpl_matrix_get_nrow(eqn_lhs);
02220 cpl_msg_info(cpl_func, "Solving %d equation%s for %d parameter%s",
02221 eqn_nrows, eqn_nrows != 1 ? "s" : "",
02222 n_parameters, n_parameters != 1 ? "s" : "");
02223 }
02224
02225
02226
02227
02228
02229
02230 double red_chisq;
02231 eqn_result = solve_normal( eqn_lhs,
02232 eqn_rhs,
02233 eqn_cov_rhs,
02234 &red_chisq);
02235 fors_matrix_null(&eqn_lhs);
02236 fors_matrix_null(&eqn_rhs);
02237 fors_matrix_null(&eqn_cov_rhs);
02238
02239 assure( !cpl_error_get_code(), return, "Could not solve equation system");
02240
02241
02242
02243
02244
02245 int offset = 0;
02246 {
02247 int size = n_coeff_polyf;
02248 if (size > 0)
02249 {
02250 result_polyf = cpl_matrix_extract(
02251 eqn_result,
02252 offset, 0,
02253 1, 1,
02254 size, 1);
02255 result_cov_polyf = cpl_matrix_extract(
02256 eqn_result,
02257 offset, 1+offset,
02258 1, 1,
02259 size, size);
02260 offset += size;
02261 }
02262 }
02263 {
02264 int size = n_coeff_params;
02265 if (size > 0)
02266 {
02267 result_params = cpl_matrix_extract(
02268 eqn_result,
02269 offset, 0,
02270 1, 1,
02271 size, 1);
02272 result_cov_params = cpl_matrix_extract(
02273 eqn_result,
02274 offset, 1+offset,
02275 1, 1,
02276 size, size);
02277 }
02278 offset += size;
02279 }
02280 {
02281 int size = n_coeff_polyp;
02282 if (size > 0)
02283 {
02284 result_polyp = cpl_matrix_extract(
02285 eqn_result,
02286 offset, 0,
02287 1, 1,
02288 size, 1);
02289 result_cov_polyp = cpl_matrix_extract(eqn_result,
02290 offset, 1+offset,
02291 1, 1,
02292 size, size);
02293 }
02294 offset += size;
02295 }
02296 passure(!cpl_error_get_code(), return);
02297
02298 fors_photometry_poly_new_from_coefficients(
02299 polyf_definition,
02300 result_polyf,
02301 result_cov_polyf,
02302 &polyf,
02303 &polyf_variance);
02304 passure(!cpl_error_get_code(), return);
02305 fors_photometry_poly_new_from_coefficients(
02306 polyp_definition,
02307 result_polyp,
02308 NULL,
02309 &polyp,
02310 NULL);
02311 passure(!cpl_error_get_code(), return);
02312
02313 cpl_msg_indent_more();
02314 fors_polynomial_dump(polyf, "f", CPL_MSG_INFO, polyf_definition);
02315 fors_polynomial_dump(polyp, "p", CPL_MSG_INFO, polyp_definition);
02316 cpl_msg_indent_less();
02317
02318 {
02319
02320 int k,
02321 i = 0;
02322 const fors_std_star *std;
02323
02324
02325 n_std_star_obs = fors_photometry_count_observations(
02326 std_star_list,
02327 obs);
02328 passure(!cpl_error_get_code(), return);
02329
02330 cpl_msg_indent_more();
02331 for ( std = fors_std_star_list_first_const(std_star_list), k = 0;
02332 std != NULL;
02333 std = fors_std_star_list_next_const(std_star_list), k++)
02334 {
02335 if (std->trusted)
02336 {
02337 cpl_msg_info(cpl_func, "M%d = %f +- %f mag (fixed, %s, %d obs)",
02338 k,
02339 std->magnitude,
02340 std->dmagnitude,
02341 std->name,
02342 cpl_array_get_int(n_std_star_obs, k, NULL));
02343 }
02344 else {
02345 cpl_msg_info(cpl_func, "M%d = %f +- %f mag (free, %s, %d obs)",
02346 k,
02347 cpl_matrix_get(result_params, i, 0),
02348 sqrt(cpl_matrix_get(result_cov_params, i, i)),
02349 std->name,
02350 cpl_array_get_int(n_std_star_obs, k, NULL));
02351 i++;
02352 }
02353 }
02354 cpl_msg_indent_less();
02355 passure(!cpl_error_get_code(), return);
02356
02357 cpl_array_delete(n_std_star_obs); n_std_star_obs = NULL;
02358 fors_std_star_list_delete(&std_star_list, fors_std_star_delete);
02359 entry_list_delete(&obs, entry_delete_but_standard); obs = NULL;
02360
02361
02362 cpl_msg_indent_more();
02363 if (fit_z)
02364 {
02365 qc_zeropoint = cpl_matrix_get(result_params, i, 0);
02366 qc_zeropoint_err = sqrt(cpl_matrix_get(result_cov_params, i, i));
02367 cpl_msg_info(cpl_func, "Z = %f +- %f mag",
02368 qc_zeropoint,
02369 qc_zeropoint_err);
02370 i++;
02371 }
02372 cpl_msg_indent_less();
02373 passure(!cpl_error_get_code(), return);
02374
02375
02376 cpl_msg_indent_more();
02377 if (n_coeff_params_ext > 0)
02378 {
02379 qc_extinction = cpl_matrix_get(result_params, i, 0);
02380 qc_extinction_err = sqrt(cpl_matrix_get(result_cov_params, i, i));
02381 if (n_coeff_params_ext == 1)
02382 {
02383 cpl_msg_info(cpl_func, "E = %f +- %f mag/airmass",
02384 qc_extinction,
02385 qc_extinction_err);
02386 i++;
02387 }
02388 else
02389 {
02390 if (summary) {
02391 cpl_table_new_column(summary, "EXT", CPL_TYPE_DOUBLE);
02392 cpl_table_new_column(summary, "DEXT", CPL_TYPE_DOUBLE);
02393 if (fit_e == FORS_FIT_NCOEFF_PERFRAME) {
02394 cpl_table_new_column(summary, "MJD-OBS",
02395 CPL_TYPE_DOUBLE);
02396 }
02397 if (fit_e == FORS_FIT_NCOEFF_PERNIGHT) {
02398 cpl_table_new_column(summary, "MJD-NIGHT",
02399 CPL_TYPE_INT);
02400 }
02401
02402 }
02403
02404 for (k = 0; k < n_coeff_params_ext; k++)
02405 {
02406 double ext = cpl_matrix_get(result_params, i, 0);
02407 double dext = sqrt(cpl_matrix_get(result_cov_params, i, i));
02408
02409 cpl_msg_info(cpl_func, "E_%d = %f +- %f mag/airmass",
02410 k, ext, dext);
02411
02412 if (summary) {
02413 cpl_table_select_all(summary);
02414 cpl_table_and_selected_int(summary, "index",
02415 CPL_EQUAL_TO, k);
02416 for (row = 0; row < n_frames; row++) {
02417 if (cpl_table_is_selected(summary, row)) {
02418 const char *filename =
02419 cpl_table_get_string(summary, "filename", row);
02420 cpl_propertylist *plist =
02421 cpl_propertylist_load_regexp(filename, 0,
02422 "MJD-OBS|ORIGIN",
02423 0);
02424
02425 cpl_table_set_double(summary,
02426 "EXT", row, ext);
02427 cpl_table_set_double(summary,
02428 "DEXT", row, dext);
02429 if (fit_e == FORS_FIT_NCOEFF_PERFRAME) {
02430 cpl_table_set_double(summary,
02431 "MJD-OBS", row,
02432 cpl_propertylist_get_double(plist,
02433 "MJD-OBS"));
02434 }
02435 if (fit_e == FORS_FIT_NCOEFF_PERNIGHT) {
02436 cpl_table_set_int(summary, "MJD-NIGHT",
02437 row,
02438 fors_photometry_get_night_id(plist));
02439 }
02440
02441 cpl_propertylist_delete(plist);
02442 }
02443 }
02444 }
02445
02446 i++;
02447 }
02448 }
02449 }
02450 cpl_msg_indent_less();
02451 passure(!cpl_error_get_code(), return);
02452
02453 if (summary) {
02454 cpl_table_select_all(summary);
02455 cpl_table_erase_column(summary, "index");
02456 }
02457
02458
02459 cpl_msg_indent_more();
02460 if (fit_c) {
02461
02462
02463
02464
02465 qc_colorterm = cpl_matrix_get(result_params, i, 0);
02466 qc_colorterm = -qc_colorterm;
02467 qc_colorterm_err = sqrt(cpl_matrix_get(result_cov_params, i, i));
02468 cpl_msg_info(cpl_func, "C_correction = %f +- %f",
02469 qc_colorterm, qc_colorterm_err);
02470 i++;
02471 }
02472 cpl_msg_indent_less();
02473 passure(!cpl_error_get_code(), return);
02474
02475
02476 if (qc_zeropoint_err > 1.0) {
02477 cpl_error_set_message(cpl_func, CPL_ERROR_INCOMPATIBLE_INPUT,
02478 "Unreliable zeropoint!");
02479 cleanup;
02480 return;
02481 }
02482 if (qc_extinction_err > 1.0) {
02483 cpl_error_set_message(cpl_func, CPL_ERROR_INCOMPATIBLE_INPUT,
02484 "Unreliable atmospheric extinction!");
02485 cleanup;
02486 return;
02487 }
02488 if (qc_colorterm_err > 1.0) {
02489 cpl_error_set_message(cpl_func, CPL_ERROR_INCOMPATIBLE_INPUT,
02490 "Unreliable color correction term!");
02491 cleanup;
02492 return;
02493 }
02494 if (qc_extinction_err > 0.0 && qc_extinction <= 0.0) {
02495 cpl_error_set_message(cpl_func, CPL_ERROR_INCOMPATIBLE_INPUT,
02496 "Impossible atmospheric extinction!");
02497 cleanup;
02498 return;
02499 }
02500 passure(!cpl_error_get_code(), return);
02501
02502 cpl_msg_indent_more();
02503 cpl_msg_info(cpl_func, "Reduced chi square = %f", red_chisq);
02504 cpl_msg_indent_less();
02505
02506 }
02507 passure(!cpl_error_get_code(), return);
02508
02509
02510 passure(polyf != NULL, return);
02511 passure(polyp != NULL, return);
02512
02513
02514
02515 {
02516 int nx,
02517 ny;
02518
02519 master_flat = fors_image_load( cpl_frameset_get_first_const(
02520 master_flat_frame),
02521 NULL, setting, NULL);
02522 assure( !cpl_error_get_code(), return, "Could not load master flat");
02523
02524 nx = fors_image_get_size_x(master_flat);
02525 ny = fors_image_get_size_y(master_flat);
02526 cpl_image *correction_map = cpl_image_new(nx, ny, CPL_TYPE_FLOAT);
02527 cpl_image *correction_map_v = cpl_image_new(nx, ny, CPL_TYPE_FLOAT);
02528
02529 cpl_msg_info(cpl_func, "Creating correction map (magnitude)");
02530 cpl_image_fill_polynomial(correction_map, polyf,
02531 1.0, 1.0, 1.0, 1.0);
02532 passure(!cpl_error_get_code(), return);
02533 cpl_image_fill_polynomial(correction_map_v, polyf_variance,
02534 1.0, 1.0, 1.0, 1.0);
02535 passure(!cpl_error_get_code(), return);
02536
02537 correction = fors_image_new(correction_map, correction_map_v);
02538 }
02539 passure(!cpl_error_get_code(), return);
02540
02541 cpl_polynomial_delete(polyf); polyf = NULL;
02542 cpl_polynomial_delete(polyf_variance); polyf_variance = NULL;
02543 fors_matrix_null(&eqn_result);
02544
02545 if (qc_zeropoint_err > 0.0 ||
02546 qc_extinction_err > 0.0 ||
02547 qc_colorterm_err > 0.0) {
02548
02549 phot_coeff = fors_phot_coeff_create(setting,
02550 qc_colorterm,
02551 qc_colorterm_err,
02552 qc_extinction,
02553 qc_extinction_err,
02554 qc_zeropoint,
02555 qc_zeropoint_err);
02556
02557
02558
02559
02560
02561 qc = cpl_propertylist_new();
02562
02563 fors_qc_start_group(qc, fors_qc_dic_version, setting->instrument);
02564
02565
02566
02567
02568
02569
02570
02571
02572
02573 if (qc_zeropoint_err > 0.0) {
02574 fors_qc_write_qc_double(qc,
02575 qc_zeropoint,
02576 "QC.INSTRUMENT.ZEROPOINT",
02577 "mag",
02578 "Instrument zeropoint",
02579 setting->instrument);
02580
02581 fors_qc_write_qc_double(qc,
02582 qc_zeropoint_err,
02583 "QC.INSTRUMENT.ZEROPOINT.ERROR",
02584 "mag",
02585 "Instrument zeropoint error",
02586 setting->instrument);
02587 }
02588
02589 if (qc_extinction_err > 0.0 && summary == NULL) {
02590 fors_qc_write_qc_double(qc,
02591 qc_extinction,
02592 "QC.ATMOSPHERIC.EXTINCTION",
02593 "mag/airmass",
02594 "Atmospheric extinction",
02595 setting->instrument);
02596
02597 fors_qc_write_qc_double(qc,
02598 qc_extinction_err,
02599 "QC.ATMOSPHERIC.EXTINCTION.ERROR",
02600 "mag/airmass",
02601 "Atmospheric extinction error",
02602 setting->instrument);
02603 }
02604
02605 if (qc_colorterm_err > 0.0) {
02606 fors_qc_write_qc_double(qc,
02607 qc_colorterm,
02608 "QC.COLOR.CORRECTION",
02609 NULL,
02610 "Linear color correction term",
02611 setting->instrument);
02612
02613 fors_qc_write_qc_double(qc,
02614 qc_colorterm_err,
02615 "QC.COLOR.CORRECTION.ERROR",
02616 NULL,
02617 "Linear color correction term error",
02618 setting->instrument);
02619 }
02620
02621 fors_qc_end_group();
02622
02623
02624
02625
02626 }
02627 passure(!cpl_error_get_code(), return);
02628
02629 if (summary && phot_coeff) {
02630 cpl_table_erase_column(phot_coeff, "EXT");
02631 cpl_table_erase_column(phot_coeff, "DEXT");
02632 if (1 == cpl_table_get_ncol(phot_coeff)) {
02633 cpl_table_delete(phot_coeff);
02634 phot_coeff = NULL;
02635 }
02636 }
02637
02638 if (phot_coeff) {
02639 fors_dfs_save_table(frames, phot_coeff, PHOT_COEFF_TABLE,
02640 qc, parameters, fors_photometry_name,
02641 cpl_frameset_get_first_const(master_flat_frame));
02642 cpl_propertylist_delete(qc); qc = NULL;
02643 }
02644
02645 assure( !cpl_error_get_code(), return, "Saving %s failed",
02646 PHOT_COEFF_TABLE);
02647
02648 if (summary) {
02649 if (fit_e == FORS_FIT_NCOEFF_PERFRAME) {
02650 tag = EXTINCTION_PER_FRAME;
02651 }
02652 if (fit_e == FORS_FIT_NCOEFF_PERNIGHT) {
02653 tag = EXTINCTION_PER_NIGHT;
02654 }
02655 fors_dfs_save_table(frames, summary, tag, NULL, parameters,
02656 fors_photometry_name,
02657 cpl_frameset_get_first_const(master_flat_frame));
02658 }
02659
02660 assure( !cpl_error_get_code(), return, "Saving %s failed", tag);
02661
02662 if (degreef1 > 0) {
02663
02664 fors_dfs_save_image(frames, correction, CORRECTION_MAP,
02665 qc, parameters, fors_photometry_name,
02666 cpl_frameset_get_first_const(master_flat_frame));
02667
02668 }
02669
02670 assure( !cpl_error_get_code(), return, "Saving %s failed",
02671 CORRECTION_MAP);
02672
02673 cpl_propertylist_delete(qc);
02674
02675
02676
02677
02678 if (degreef1 > 0) {
02679 cpl_msg_info(cpl_func, "Creating correction map (flux)");
02680
02681 fors_image_multiply_scalar(correction, -0.4, -1);
02682 fors_image_exponential(correction, 10, -1);
02683
02684
02685 fors_image_divide_scalar(correction,
02686 fors_image_get_median(correction, NULL), -1.0);
02687
02688 fors_dfs_save_image(frames, correction, CORRECTION_FACTOR,
02689 NULL, parameters, fors_photometry_name,
02690 cpl_frameset_get_first_const(master_flat_frame));
02691
02692 assure( !cpl_error_get_code(), return, "Saving %s failed",
02693 CORRECTION_FACTOR);
02694 }
02695
02696 if (degreef1 > 0) {
02697 cpl_msg_info(cpl_func, "Creating corrected master flat");
02698 fors_image_multiply(master_flat, correction);
02699
02700 fors_dfs_save_image(frames, master_flat, MASTER_FLAT_IMG,
02701 NULL, parameters, fors_photometry_name,
02702 cpl_frameset_get_first_const(master_flat_frame));
02703 assure( !cpl_error_get_code(), return, "Saving %s failed",
02704 MASTER_FLAT_IMG);
02705 }
02706
02707 cleanup;
02708 return;
02709 }
02710
02711
02712 #undef cleanup
02713 #define cleanup \
02714 do { \
02715 fors_std_star_delete(&std_star); \
02716 cpl_propertylist_delete(header); header = NULL; \
02717 cpl_table_delete(aligned_phot); aligned_phot = NULL; \
02718 fors_setting_delete(&setting_f); \
02719 fors_delete_star_lists(&obs, std_star_list); \
02720 fors_star_delete_but_standard(&obs_star); \
02721 fors_std_star_delete(&std_star); \
02722 \
02723 *n_frames = 0; \
02724 } while (0)
02725
02740
02741 static entry_list *
02742 fors_photometry_read_input( const cpl_frameset *alphot_frames,
02743 const fors_setting *setting,
02744 int (*get_atm_ext_id_function)(
02745 const cpl_propertylist *header),
02746 bool import_unknown,
02747 int *n_frames,
02748 fors_std_star_list **std_star_list,
02749 int filter)
02750 {
02751 entry_list *obs = NULL;
02752 fors_setting *setting_f = NULL;
02753 fors_star *obs_star = NULL;
02754 fors_std_star *std_star = NULL;
02755 cpl_propertylist *header = NULL;
02756 cpl_table *aligned_phot = NULL;
02757 const cpl_frame *f;
02758 int iframe,
02759 inonstd = 0;
02760 cpl_errorstate errstat = cpl_errorstate_get();
02761
02762
02763 cleanup;
02764
02765
02766 obs = entry_list_new();
02767 *std_star_list = fors_std_star_list_new();
02768
02769
02770
02771
02772
02773 for (f = cpl_frameset_get_first_const(alphot_frames), iframe = 0;
02774 f != NULL;
02775 f = cpl_frameset_get_next_const(alphot_frames), iframe++)
02776 {
02777 const char *filename;
02778 int atm_ext_id = 0;
02779 double airmass;
02780
02781 filename = cpl_frame_get_filename(f);
02782 cpl_msg_info(cpl_func, "Loading %s", filename);
02783 cpl_msg_indent_more();
02784
02785 aligned_phot = cpl_table_load(filename, 1, 1);
02786 assure( cpl_errorstate_is_equal(errstat),
02787 return NULL,
02788 "Could not load %s",
02789 filename);
02790
02791 if (filter && cpl_table_has_column(aligned_phot, "WEIGHT")) {
02792 cpl_table_and_selected_double(aligned_phot,
02793 "WEIGHT", CPL_LESS_THAN, 1.0);
02794 cpl_table_and_selected_double(aligned_phot,
02795 "WEIGHT", CPL_GREATER_THAN, -1.0);
02796 cpl_table_erase_selected(aligned_phot);
02797 }
02798
02799
02800 header = cpl_propertylist_load(filename, 0);
02801 assure( cpl_errorstate_is_equal(errstat),
02802 return NULL,
02803 "Could not load %s header",
02804 filename);
02805
02806
02807
02808
02809 fors_setting_verify(setting, f, &setting_f);
02810
02811 airmass = cpl_propertylist_get_double(header, "AIRMASS");
02812 assure( cpl_errorstate_is_equal(errstat),
02813 return NULL,
02814 "%s: Could not read %s",
02815 filename, "AIRMASS");
02816
02817
02818
02819
02820
02821
02822
02823 struct {
02824 const char *name;
02825 cpl_type type;
02826 bool optional;
02827 } col[] =
02828 {{"RA", CPL_TYPE_DOUBLE, false},
02829 {"DEC", CPL_TYPE_DOUBLE, false},
02830 {"X", CPL_TYPE_DOUBLE, false},
02831 {"Y", CPL_TYPE_DOUBLE, false},
02832 {"FWHM", CPL_TYPE_DOUBLE, false},
02833 {"A", CPL_TYPE_DOUBLE, false},
02834 {"B", CPL_TYPE_DOUBLE, false},
02835 {"THETA", CPL_TYPE_DOUBLE, false},
02836 {"INSTR_MAG", CPL_TYPE_DOUBLE, false},
02837 {"DINSTR_MAG", CPL_TYPE_DOUBLE, false},
02838 {"CAT_MAG", CPL_TYPE_DOUBLE, false},
02839 {"DCAT_MAG", CPL_TYPE_DOUBLE, false},
02840 {"MAG", CPL_TYPE_DOUBLE, false},
02841 {"DMAG", CPL_TYPE_DOUBLE, false},
02842 {"COLOR", CPL_TYPE_DOUBLE, false},
02843 {"DCOLOR", CPL_TYPE_DOUBLE, true},
02844 {"COV_CATM_COL", CPL_TYPE_DOUBLE, true},
02845 {"CLASS_STAR", CPL_TYPE_DOUBLE, false},
02846 {"USE_CAT", CPL_TYPE_INT, false},
02847 {"OBJECT", CPL_TYPE_STRING, false}};
02848 {
02849 unsigned i = 0;
02850 for (i = 0; i < sizeof(col) / sizeof(*col); i++)
02851 {
02852 bool exists;
02853 exists = cpl_table_has_column(aligned_phot, col[i].name);
02854 cassure(exists || col[i].optional,
02855 CPL_ERROR_DATA_NOT_FOUND,
02856 return NULL,
02857 "%s: Missing column %s", filename, col[i].name);
02858 cassure((!exists) ||
02859 cpl_table_get_column_type(aligned_phot, col[i].name)
02860 == col[i].type,
02861 CPL_ERROR_INVALID_TYPE,
02862 return NULL,
02863 "%s: column %s: Type is %s, %s expected ",
02864 filename,
02865 col[i].name,
02866 fors_type_get_string(
02867 cpl_table_get_column_type(aligned_phot,
02868 col[i].name)
02869 ),
02870 fors_type_get_string(col[i].type));
02871 }
02872 }
02873
02874 if (get_atm_ext_id_function != NULL)
02875 {
02876 atm_ext_id = get_atm_ext_id_function(header);
02877 assure( cpl_errorstate_is_equal(errstat),
02878 return NULL,
02879 "Getting atmospheric extinction "
02880 "identifier failed.");
02881 }
02882
02883
02884
02885
02886 int i;
02887 for (i = 0; i < cpl_table_get_nrow(aligned_phot); i++)
02888 {
02889 obs_star = fors_star_new_from_table(
02890 aligned_phot,
02891 i,
02892 "X", "Y",
02893 "FWHM",
02894 "A", "B",
02895 "THETA",
02896 "INSTR_MAG", "DINSTR_MAG",
02897 "CLASS_STAR");
02898 std_star = fors_std_star_new_from_table(
02899 aligned_phot,
02900 i,
02901 "RA", "DEC",
02902 "MAG", "DMAG",
02903 "CAT_MAG", "DCAT_MAG",
02904 "COLOR", NULL,
02905 NULL,
02906 NULL, NULL,
02907 "OBJECT");
02908
02909 if (cpl_table_has_column(aligned_phot, "DCOLOR")
02910 && cpl_table_has_column(aligned_phot, "COV_CATM_COL"))
02911 {
02912 std_star->dcolor = cpl_table_get(
02913 aligned_phot,
02914 "DCOLOR",
02915 i,
02916 NULL);
02917 std_star->cov_catm_color = cpl_table_get(
02918 aligned_phot,
02919 "COV_CATM_COL",
02920 i,
02921 NULL);
02922 }
02923
02924
02925
02926
02927
02928
02929
02930
02931 obs_star->id = std_star;
02932
02933
02934 std_star->trusted = (0 != cpl_table_get_int(
02935 aligned_phot, "USE_CAT", i, NULL));
02936
02937 assure( cpl_errorstate_is_equal(errstat)
02938 && obs_star != NULL
02939 && std_star != NULL,
02940 return NULL,
02941 "Reading from aligned photometry "
02942 "table failed.");
02943
02944 cassure( obs_star->dmagnitude > 0,
02945 CPL_ERROR_ILLEGAL_INPUT,
02946 return NULL,
02947 "Non-positive magnitude error: "
02948 "%f mag at row %d",
02949 obs_star->dmagnitude,
02950 i+1);
02951
02952
02953
02954
02955 if (! fors_extract_check_sex_star(obs_star, NULL))
02956 {
02957
02958 cpl_msg_warning(cpl_func, "Rejecting object no. %d from "
02959 "frame %d, "
02960 "which should have been caught "
02961 "by recent fors_zeropoint recipe. "
02962 "Consider reprocessing input data.",
02963 i+1, iframe+1);
02964 fors_std_star_delete(&std_star);
02965 }
02966 else if (std_star->name != NULL && (std_star->name)[0] != '\0')
02967 {
02968
02969
02970 std_star = fors_photometry_read_input_listinsert_star_if_new(
02971 *std_star_list,
02972 std_star,
02973 arcsec_tol);
02974 }
02975 else if (import_unknown)
02976 {
02977 fors_std_star *s_in_list;
02978
02979 std_star->trusted = false;
02980
02981 s_in_list = fors_photometry_read_input_listinsert_star_if_new(
02982 *std_star_list,
02983 std_star,
02984 arcsec_tol);
02985
02986 if (s_in_list == std_star)
02987 {
02988 char name[16];
02989 sprintf(name, "non-std-%d", inonstd);
02990 fors_std_star_set_name(std_star, name);
02991 inonstd++;
02992 }
02993
02994 std_star = s_in_list;
02995 }
02996 else
02997 {
02998 fors_std_star_delete(&std_star);
02999 }
03000
03001 if (std_star != NULL)
03002 {
03003 entry *e;
03004 obs_star->id = std_star;
03005
03006
03007 e = entry_new( iframe,
03008 -1,
03009 airmass,
03010 setting_f->average_gain,
03011 setting_f->exposure_time,
03012 atm_ext_id,
03013 obs_star);
03014 entry_list_insert( obs, e);
03015 }
03016 else
03017 {
03018 fors_star_delete_but_standard(&obs_star);
03019 }
03020
03021 }
03022
03023 cpl_propertylist_delete(header); header = NULL;
03024 cpl_table_delete(aligned_phot); aligned_phot = NULL;
03025 fors_setting_delete(&setting_f);
03026
03027 cpl_msg_indent_less();
03028 }
03029
03030 *n_frames = iframe;
03031
03032
03033
03034
03035
03036
03037
03038 entry_list_reverse(obs);
03039 fors_std_star_list_reverse(*std_star_list);
03040
03041
03042 {
03043 fors_std_star *ref;
03044 entry *e;
03045
03046 for (e = entry_list_first(obs); e != NULL; e = entry_list_next(obs))
03047 {
03048 int i;
03049 for ( ref = fors_std_star_list_first(*std_star_list), i = 0;
03050 ref != NULL;
03051 ref = fors_std_star_list_next(*std_star_list), i++)
03052 {
03053 if (e->star->id == ref)
03054 {
03055 e->star_index = i;
03056 break;
03057 }
03058 }
03059 }
03060 }
03061
03062 cassure( entry_list_size(obs) > 0,
03063 CPL_ERROR_DATA_NOT_FOUND,
03064 return NULL,
03065 "No stars found");
03066
03067 return obs;
03068 }
03069
03070
03071 #undef cleanup
03072 #define cleanup
03073
03086
03087 static fors_std_star*
03088 fors_photometry_read_input_listinsert_star_if_new(
03089 fors_std_star_list *std_list,
03090 fors_std_star *std,
03091 double mind_arcsec)
03092 {
03093 cpl_errorstate errstat = cpl_errorstate_get();
03094 bool found = false;
03095
03096 cassure_automsg( std_list != NULL,
03097 CPL_ERROR_NULL_INPUT,
03098 return std);
03099 cassure_automsg( std != NULL,
03100 CPL_ERROR_NULL_INPUT,
03101 return std);
03102 cassure_automsg( mind_arcsec > 0,
03103 CPL_ERROR_ILLEGAL_INPUT,
03104 return std);
03105
03106 if (fors_std_star_list_size(std_list) > 0)
03107 {
03108 fors_std_star *nearest;
03109 double dist_arcsec;
03110
03111
03112
03113 nearest = fors_std_star_list_kth_val(
03114 std_list, 1,
03115 (fors_std_star_list_func_eval)
03116 fors_std_star_dist_arcsec,
03117 std);
03118 passure(cpl_errorstate_is_equal(errstat), return std);
03119
03120 dist_arcsec = fors_std_star_dist_arcsec(nearest, std);
03121 passure(cpl_errorstate_is_equal(errstat), return std);
03122
03123 cpl_msg_debug( cpl_func,
03124 "dist = %f arcseconds",
03125 dist_arcsec);
03126
03127 if (dist_arcsec < mind_arcsec)
03128 {
03129
03130 nearest->trusted &= std->trusted;
03131
03132 fors_std_star_delete(&std);
03133 std = nearest;
03134 found = true;
03135 }
03136 }
03137
03138 if (!found)
03139 {
03140 fors_std_star_list_insert(std_list, std);
03141 }
03142
03143 return std;
03144 }
03145
03146
03147 #undef cleanup
03148 #define cleanup
03149
03150
03151
03152
03153
03154
03155
03156
03157
03158 static cpl_error_code
03159 fors_photometry_adjust_fit_mag_flags( fors_std_star_list *stdl,
03160 entry_list *obsl,
03161 bool override_fit_m,
03162 int *n_mag_fits)
03163 {
03164 fors_std_star *std;
03165 cpl_errorstate errstat = cpl_errorstate_get();
03166
03167 cassure_automsg( stdl != NULL,
03168 CPL_ERROR_NULL_INPUT,
03169 return cpl_error_get_code());
03170 cassure_automsg( obsl != NULL,
03171 CPL_ERROR_NULL_INPUT,
03172 return cpl_error_get_code());
03173 cassure_automsg( n_mag_fits != NULL,
03174 CPL_ERROR_NULL_INPUT,
03175 return cpl_error_get_code());
03176
03177 *n_mag_fits = 0;
03178
03179
03180
03181
03182
03183 if (override_fit_m)
03184 {
03185 std = fors_std_star_list_first(stdl);
03186
03187 while (std != NULL && std->trusted == false)
03188 {
03189 std = fors_std_star_list_next(stdl);
03190 (*n_mag_fits)++;
03191 }
03192
03193 if (std != NULL)
03194 std = fors_std_star_list_next(stdl);
03195
03196 for ( ; std != NULL; std = fors_std_star_list_next(stdl))
03197 {
03198 std->trusted = false;
03199 (*n_mag_fits)++;
03200 }
03201 }
03202 else
03203 {
03204 for ( std = fors_std_star_list_first(stdl);
03205 std != NULL;
03206 std = fors_std_star_list_next(stdl))
03207 {
03208 if(! std->trusted)
03209 (*n_mag_fits)++;
03210 }
03211 }
03212
03213 return (cpl_errorstate_is_equal(errstat) ?
03214 CPL_ERROR_NONE :
03215 cpl_error_get_code());
03216 }
03217
03218
03219 #undef cleanup
03220 #define cleanup \
03221 do { \
03222 cpl_array_delete(n_obs_per_std); n_obs_per_std = NULL; \
03223 fors_std_star_list_delete(&std_list_copy, NULL); \
03224 entry_list_delete(&obs_list_copy, NULL); \
03225 } while (0)
03226
03233
03234 static cpl_error_code
03235 fors_photometry_remove_unnecessary( fors_std_star_list *std_list,
03236 entry_list *obs_list,
03237 int *n_mag_fits)
03238 {
03239 cpl_array *n_obs_per_std = NULL;
03240 fors_std_star_list *std_list_copy = NULL;
03241 entry_list *obs_list_copy = NULL;
03242 fors_std_star *std;
03243 entry *obs;
03244 int n_std_stars,
03245 n_removed = 0,
03246 n;
03247 cpl_errorstate errstat = cpl_errorstate_get();
03248
03249 cassure_automsg( std_list != NULL,
03250 CPL_ERROR_NULL_INPUT,
03251 return cpl_error_get_code());
03252 cassure_automsg( obs_list != NULL,
03253 CPL_ERROR_NULL_INPUT,
03254 return cpl_error_get_code());
03255 cassure_automsg( n_mag_fits != NULL,
03256 CPL_ERROR_NULL_INPUT,
03257 return cpl_error_get_code());
03258
03259 *n_mag_fits = 0;
03260
03261 n_std_stars = fors_std_star_list_size(std_list);
03262
03263 n_obs_per_std = fors_photometry_count_observations(std_list, obs_list);
03264 passure(cpl_errorstate_is_equal(errstat), return cpl_error_get_code());
03265
03266
03267
03268
03269 std_list_copy = fors_std_star_list_duplicate(std_list, NULL);
03270 obs_list_copy = entry_list_duplicate(obs_list, NULL);
03271
03272
03273 for ( obs = entry_list_first(obs_list_copy);
03274 obs != NULL;
03275 obs = entry_list_next(obs_list_copy))
03276 {
03277 int n_obs;
03278 bool remove = false;
03279
03280 remove |= (obs->star_index < 0 || obs->star_index >= n_std_stars);
03281
03282 n_obs = cpl_array_get_int(n_obs_per_std, obs->star_index, NULL);
03283 assure( cpl_errorstate_is_equal(errstat),
03284 return cpl_error_get_code(),
03285 NULL);
03286 remove |= (n_obs < 2 && !obs->star->id->trusted);
03287
03288 if (remove)
03289 {
03290 entry_list_remove(obs_list, obs);
03291 entry_delete_but_standard(&obs);
03292 }
03293 }
03294
03295 for ( std = fors_std_star_list_first(std_list_copy), n = 0;
03296 std != NULL;
03297 std = fors_std_star_list_next(std_list_copy), n++)
03298 {
03299 int n_obs;
03300
03301 n_obs = cpl_array_get_int(n_obs_per_std, n, NULL);
03302 assure( cpl_errorstate_is_equal(errstat),
03303 return cpl_error_get_code(),
03304 NULL);
03305
03306 if (n_obs < 2 && !std->trusted)
03307 {
03308 fors_std_star_list_remove(std_list, std);
03309 fors_std_star_delete(&std);
03310 n_removed++;
03311 }
03312 else if(!std->trusted)
03313 {
03314 (*n_mag_fits)++;
03315 }
03316 }
03317
03318 cleanup;
03319
03320
03321 for ( obs = entry_list_first(obs_list);
03322 obs != NULL;
03323 obs = entry_list_next(obs_list))
03324 {
03325 for ( std = fors_std_star_list_first(std_list), n = 0;
03326 std != NULL;
03327 std = fors_std_star_list_next(std_list), n++)
03328 {
03329 if (obs->star->id == std)
03330 {
03331 obs->star_index = n;
03332 break;
03333 }
03334 }
03335 }
03336
03337 if (n_removed > 0)
03338 cpl_msg_info( cpl_func,
03339 "Discarded %d untrusted/fitted objects which were "
03340 "observed only once (and therefore don't contribute).",
03341 n_removed);
03342
03343 cleanup;
03344 return (cpl_errorstate_is_equal(errstat) ?
03345 CPL_ERROR_NONE :
03346 cpl_error_get_code());
03347 }
03348
03349
03350 #undef cleanup
03351 #define cleanup \
03352 do { \
03353 cpl_array_unwrap(n_obs_a); n_obs_a = NULL; \
03354 cpl_free(n_obs); n_obs = NULL; \
03355 } while (0)
03356
03363
03364 static cpl_array*
03365 fors_photometry_count_observations( fors_std_star_list *std_list,
03366 entry_list *obs_list)
03367 {
03368 int *n_obs = NULL;
03369 cpl_array *n_obs_a = NULL;
03370 entry *e;
03371 int n_std_stars;
03372
03373 cassure_automsg( std_list != NULL,
03374 CPL_ERROR_NULL_INPUT,
03375 return n_obs_a);
03376 cassure_automsg( obs_list != NULL,
03377 CPL_ERROR_NULL_INPUT,
03378 return n_obs_a);
03379
03380 n_std_stars = fors_std_star_list_size(std_list);
03381 n_obs = cpl_calloc(n_std_stars, sizeof(*n_obs));
03382
03383 for ( e = entry_list_first(obs_list);
03384 e != NULL;
03385 e = entry_list_next(obs_list))
03386 {
03387 ppassure( e->star_index >= 0
03388 && e->star_index < n_std_stars,
03389 CPL_ERROR_UNSPECIFIED,
03390 return n_obs_a);
03391 n_obs[e->star_index]++;
03392 }
03393
03394 n_obs_a = cpl_array_wrap_int(n_obs, n_std_stars);
03395 return n_obs_a;
03396 }
03397
03398
03399 #undef cleanup
03400 #define cleanup \
03401 do { \
03402 fors_matrix_null(&lhs); \
03403 if (n_fit_e_cols != NULL) \
03404 *n_fit_e_cols = 0; \
03405 } while (0)
03406
03430
03431 static cpl_matrix*
03432 build_equations_lhs_matrix_from_parameters( const entry_list *obs_list,
03433 const fors_std_star_list *std_list,
03434 bool fit_z,
03435 bool fit_c,
03436 int *n_fit_e_cols)
03437 {
03438 int n_std_stars,
03439 n_obs,
03440 n_col,
03441 n_fit_std_mag = 0,
03442 n_frames,
03443 n_atm_ext,
03444 row;
03445 const entry *e;
03446 const fors_std_star
03447 *std;
03448 cpl_matrix *lhs = NULL;
03449 bool printed = false;
03450
03451
03452 cleanup;
03453
03454
03455 assure(!cpl_error_get_code(), return lhs, "Previous error caught.");
03456
03457 ppassure( obs_list != NULL
03458 && std_list != NULL,
03459 CPL_ERROR_NULL_INPUT,
03460 return lhs);
03461
03462 n_std_stars = fors_std_star_list_size(std_list);
03463 n_obs = entry_list_size(obs_list);
03464
03465 cassure( n_std_stars > 0 && n_obs > 0,
03466 CPL_ERROR_DATA_NOT_FOUND,
03467 return lhs,
03468 "Empty input list");
03469
03470
03471 n_obs = entry_list_size(obs_list);
03472 for ( std = fors_std_star_list_first_const(std_list);
03473 std != NULL;
03474 std = fors_std_star_list_next_const(std_list))
03475 {
03476 n_fit_std_mag += !(std->trusted);
03477 }
03478
03479 n_frames = 0;
03480 n_atm_ext = 0;
03481 for ( e = entry_list_first_const(obs_list);
03482 e != NULL;
03483 e = entry_list_next_const(obs_list))
03484 {
03485
03486 if (e->frame_index + 1 > n_frames)
03487 n_frames = e->frame_index + 1;
03488 if (e->atm_ext_index + 1 > n_atm_ext)
03489 n_atm_ext = e->atm_ext_index + 1;
03490 }
03491 if (n_atm_ext < 0) n_atm_ext = 0;
03492 passure(!cpl_error_get_code(), return lhs);
03493
03494 n_col = n_fit_std_mag
03495 + ((fit_z) ? 1 : 0)
03496 + n_atm_ext
03497 + ((fit_c) ? 1 : 0);
03498
03499 if (n_col == 0)
03500 {
03501 cleanup;
03502 return lhs;
03503 }
03504
03505 lhs = cpl_matrix_new(n_obs, n_col);
03506 passure(!cpl_error_get_code(), return lhs);
03507
03508
03509
03510 for (e = entry_list_first_const(obs_list), row = 0;
03511 e != NULL;
03512 e = entry_list_next_const(obs_list), row++)
03513 {
03514 int col = 0,
03515 k;
03516
03517
03518 ppassure( e->star_index >= 0,
03519 CPL_ERROR_ILLEGAL_INPUT,
03520 return lhs);
03521
03522 if (n_fit_std_mag > 0)
03523 {
03524 for ( std = fors_std_star_list_first_const(std_list), k = 0;
03525 std != NULL;
03526 std = fors_std_star_list_next_const(std_list), k++)
03527 {
03528 if (!(std->trusted))
03529 {
03530 if (!printed)
03531 cpl_msg_debug( cpl_func,
03532 "Creating column for mag(M%d)",
03533 k);
03534 if (e->star->id == std)
03535 cpl_matrix_set(lhs, row, col, 1);
03536 col++;
03537 }
03538 }
03539 }
03540
03541 if (fit_z)
03542 {
03543 if (!printed)
03544 cpl_msg_debug(cpl_func, "Creating column for Z");
03545 cpl_matrix_set(lhs, row, col++, -1);
03546 }
03547
03548 if (n_atm_ext > 0)
03549 {
03550 for (k = 0; k < n_atm_ext; k++)
03551 {
03552 if (!printed)
03553 cpl_msg_debug(cpl_func, "Creating column for E_%d", k);
03554 double val = (k == e->atm_ext_index) ? e->airmass : 0;
03555 cpl_matrix_set(lhs, row, col++, val);
03556 }
03557 }
03558
03559 if (fit_c)
03560 {
03561 double c;
03562 if (!printed)
03563 cpl_msg_debug(cpl_func, "Creating column for color "
03564 "correction term");
03565 c = e->star->id->color;
03566
03567
03568
03569 if (!(e->star->id->trusted))
03570 c = 0;
03571 cpl_matrix_set(lhs, row, col++, c);
03572 }
03573 printed = true;
03574 }
03575 passure(!cpl_error_get_code(), return lhs);
03576
03577 if (n_fit_e_cols != NULL)
03578 *n_fit_e_cols = n_atm_ext;
03579
03580 return lhs;
03581 }
03582
03583
03584 #undef cleanup
03585 #define cleanup \
03586 do { \
03587 fors_matrix_null(&mat); \
03588 cpl_array_delete(Apowers); Apowers = NULL; \
03589 } while (0)
03590
03608
03609 static cpl_matrix*
03610 build_equations_lhs_matrix_from_poly( const entry_list *obs_list,
03611 const cpl_polynomial *poly,
03612 const char *pname,
03613 double (*powerfunc)(
03614 const entry*,
03615 const cpl_array*))
03616 {
03617 int n_obs,
03618 n_coeff,
03619 n_dims,
03620 row;
03621 cpl_matrix *mat = NULL;
03622 cpl_array *Apowers = NULL;
03623 int *ipowers;
03624 const entry *e;
03625 cpl_error_code errc;
03626 bool printed = false;
03627
03628 assure(!(errc=cpl_error_get_code()), return NULL, "Previous error caught.");
03629
03630
03631 ppassure( poly != NULL && obs_list != NULL,
03632 CPL_ERROR_NULL_INPUT,
03633 return NULL);
03634
03635
03636 n_obs = entry_list_size(obs_list);
03637 n_coeff = fors_polynomial_count_coeff(poly);
03638 passure(!cpl_error_get_code(), return NULL);
03639
03640 if (n_coeff == 0)
03641 return NULL;
03642
03643 mat = cpl_matrix_new(n_obs, n_coeff);
03644
03645
03646 n_dims = cpl_polynomial_get_dimension(poly);
03647 Apowers = cpl_array_new(n_dims, CPL_TYPE_INT);
03648 cpl_array_fill_window_int(Apowers, 0, n_dims, 0);
03649 ipowers = cpl_array_get_data_int(Apowers);
03650 passure(!cpl_error_get_code(), return NULL);
03651
03652
03653
03654 for (e = entry_list_first_const(obs_list), row = 0;
03655 e != NULL;
03656 e = entry_list_next_const(obs_list), row++)
03657 {
03658 int col = 0;
03659 bool overflow;
03660
03661 overflow = fors_polynomial_powers_find_first_coeff(poly, ipowers);
03662 while (!overflow)
03663 {
03664 if (!printed)
03665 {
03666 char *cn = fors_polynomial_sprint_coeff(poly, ipowers, pname);
03667 if (cn != NULL)
03668 {
03669 cpl_msg_debug(cpl_func, "Creating column for %s", cn);
03670 cpl_free(cn);
03671 }
03672 }
03673 cpl_matrix_set(mat, row, col++, (*powerfunc)(e, Apowers));
03674 passure(!cpl_error_get_code(), return NULL);
03675
03676 overflow = fors_polynomial_powers_find_next_coeff(poly, ipowers);
03677 }
03678
03679 printed = true;
03680 }
03681
03682 cpl_array_delete(Apowers);
03683
03684 return mat;
03685 }
03686
03687
03688 #undef cleanup
03689 #define cleanup \
03690 do { \
03691 fors_matrix_null(&rhs_input_cov); \
03692 fors_matrix_null(&rhs_jacobian); \
03693 fors_matrix_null(&rhs_input); \
03694 fors_matrix_null(&rhs_jacobian_T); \
03695 fors_matrix_null(&tmp_matrix); \
03696 fors_matrix_null(rhs); \
03697 fors_matrix_null(rhs_cov); \
03698 } while (0)
03699
03714
03715 static cpl_error_code
03716 build_equations_rhs_cov( const entry_list *obs_list,
03717 const fors_std_star_list *std_list,
03718 bool fit_z,
03719 bool fit_c,
03720 bool fit_e,
03721 double color_coeff,
03722 double dcolor_coeff,
03723 double ext_coeff,
03724 double dext_coeff,
03725 double zpoint,
03726 double dzpoint,
03727 cpl_matrix **rhs,
03728 cpl_matrix **rhs_cov)
03729 {
03730
03731
03732
03733
03734
03735
03736
03737
03738
03739
03740
03741
03742
03743
03744
03745 cpl_matrix *rhs_input_cov = NULL,
03746 *rhs_jacobian = NULL,
03747 *rhs_input = NULL,
03748 *rhs_jacobian_T = NULL,
03749 *tmp_matrix = NULL;
03750 int n_std_stars,
03751 n_obs,
03752 n_col,
03753 r,
03754 c;
03755 const fors_std_star
03756 *std;
03757 const entry *obs;
03758 bool printed_debug = false;
03759 cpl_errorstate errstat = cpl_errorstate_get();
03760
03761
03762 cleanup;
03763
03764
03765 cassure_automsg( obs_list != NULL,
03766 CPL_ERROR_NULL_INPUT,
03767 return cpl_error_get_code());
03768 cassure_automsg( std_list != NULL,
03769 CPL_ERROR_NULL_INPUT,
03770 return cpl_error_get_code());
03771 cassure_automsg( rhs != NULL,
03772 CPL_ERROR_NULL_INPUT,
03773 return cpl_error_get_code());
03774 cassure_automsg( rhs_cov != NULL,
03775 CPL_ERROR_NULL_INPUT,
03776 return cpl_error_get_code());
03777
03778 n_std_stars = fors_std_star_list_size(std_list);
03779 n_obs = entry_list_size(obs_list);
03780
03781 cassure( n_std_stars > 0 && n_obs > 0,
03782 CPL_ERROR_DATA_NOT_FOUND,
03783 return cpl_error_get_code(),
03784 "Empty input list");
03785
03786
03787 n_col = n_std_stars*2 + 3;
03788
03789
03790
03791
03792
03793
03794
03795
03796
03797
03798
03799
03800
03801
03802 rhs_input_cov = cpl_matrix_new(n_col, n_col);
03803 for ( std = fors_std_star_list_first_const(std_list), c = 0;
03804 std != NULL;
03805 std = fors_std_star_list_next_const(std_list), c += 2)
03806 {
03807 double dcatm = std->dcat_magnitude,
03808 dcolor = std->dcolor,
03809 cov_catmag_color = std->cov_catm_color;
03810
03811
03812
03813 if (!(dcolor > 0) || isnan(cov_catmag_color))
03814 {
03815
03816
03817
03818
03819
03820
03821
03822
03823
03824
03825
03826
03827
03828
03829 cov_catmag_color = 0;
03830 dcolor = 0;
03831 if (std->trusted)
03832
03833
03834 {
03835 if (!fit_c)
03836 {
03837 cassure( dcatm > 0,
03838 CPL_ERROR_ILLEGAL_INPUT,
03839 return cpl_error_get_code(),
03840 "Could not determine color "
03841 "corrected magnitude with error "
03842 "estimate of object %s",
03843 (std->name != NULL) ?
03844 std->name : "unknown");
03845 dcatm = std->dmagnitude;
03846 if (!printed_debug)
03847 {
03848 cpl_msg_debug( cpl_func,
03849 "Having old fors_zeropoint data. "
03850 "Using color corrected magnitudes "
03851 "instead of catalogue magnitude "
03852 "and color separately.");
03853 printed_debug = true;
03854 }
03855 }
03856 }
03857
03858
03859
03860 }
03861
03862 cpl_matrix_set(rhs_input_cov, c, c, dcatm*dcatm);
03863 cpl_matrix_set(rhs_input_cov, c+1, c+1, dcolor*dcolor);
03864 cpl_matrix_set(rhs_input_cov, c, c+1, cov_catmag_color);
03865 cpl_matrix_set(rhs_input_cov, c+1, c, cov_catmag_color);
03866 }
03867 cpl_matrix_set(rhs_input_cov, c, c, dcolor_coeff*dcolor_coeff);
03868 cpl_matrix_set(rhs_input_cov, c+1, c+1, dext_coeff*dext_coeff);
03869 cpl_matrix_set(rhs_input_cov, c+1, c+1, dzpoint*dzpoint);
03870
03871 passure(cpl_errorstate_is_equal(errstat), return cpl_error_get_code());
03872
03873
03874
03875
03876
03877
03878
03879
03880
03881
03882
03883
03884
03885
03886
03887
03888
03889
03890
03891
03892
03893
03894 rhs_jacobian = cpl_matrix_new(n_obs, n_col);
03895 for ( obs = entry_list_first_const(obs_list), r = 0;
03896 obs != NULL;
03897 obs = entry_list_next_const(obs_list), r++)
03898 {
03899 bool fit_mag,
03900 compensate_color;
03901
03902 c = obs->star_index * 2;
03903 fit_mag = !(obs->star->id->trusted);
03904
03905
03906 compensate_color = (!fit_c) && (!fit_mag);
03907
03908 cpl_matrix_set(rhs_jacobian, r, c, -(!fit_mag)*1.0);
03909 cpl_matrix_set(rhs_jacobian, r, c+1, +(compensate_color)*color_coeff);
03910
03911 cpl_matrix_set(rhs_jacobian, r, n_col-3, + (compensate_color)
03912 * obs->star->id->color);
03913 cpl_matrix_set(rhs_jacobian, r, n_col-2, -(!fit_e)*obs->airmass);
03914 cpl_matrix_set(rhs_jacobian, r, n_col-1, +(!fit_z)*1.0);
03915 }
03916
03917 passure(cpl_errorstate_is_equal(errstat), return cpl_error_get_code());
03918
03919
03920
03921
03922
03923
03924
03925
03926
03927
03928 rhs_input = cpl_matrix_new(n_col, 1);
03929 for ( std = fors_std_star_list_first_const(std_list), r = 0;
03930 std != NULL;
03931 std = fors_std_star_list_next_const(std_list), r += 2)
03932 {
03933 double catm = std->cat_magnitude,
03934 color = std->color;
03935
03936
03937
03938 if (!(std->dcolor > 0) || isnan(std->cov_catm_color))
03939 {
03940
03941 color = 0;
03942 if (std->trusted)
03943 {
03944 if (!fit_c)
03945 catm = std->magnitude;
03946 }
03947 }
03948
03949 cpl_matrix_set(rhs_input, r, 0, catm);
03950 cpl_matrix_set(rhs_input, r+1, 0, color);
03951 }
03952
03953
03954 cpl_matrix_set(rhs_input, r+1, 0, ext_coeff);
03955 cpl_matrix_set(rhs_input, r+2, 0, zpoint);
03956
03957 passure(cpl_errorstate_is_equal(errstat), return cpl_error_get_code());
03958
03959
03960
03961 *rhs = cpl_matrix_product_create(rhs_jacobian, rhs_input);
03962 passure(cpl_errorstate_is_equal(errstat), return cpl_error_get_code());
03963
03964 rhs_jacobian_T = cpl_matrix_transpose_create(rhs_jacobian);
03965 tmp_matrix = cpl_matrix_product_create(rhs_input_cov, rhs_jacobian_T);
03966 passure(cpl_errorstate_is_equal(errstat), return cpl_error_get_code());
03967 *rhs_cov = cpl_matrix_product_create(rhs_jacobian, tmp_matrix);
03968 passure(cpl_errorstate_is_equal(errstat), return cpl_error_get_code());
03969
03970
03971
03972
03973
03974
03975
03976 for ( obs = entry_list_first_const(obs_list), r = 0;
03977 obs != NULL;
03978 obs = entry_list_next_const(obs_list), r++)
03979 {
03980 cpl_matrix_set(*rhs, r, 0, cpl_matrix_get(*rhs, r, 0)
03981 + obs->star->magnitude
03982 + 2.5*log10(obs->gain)
03983 + 2.5*log10(obs->exptime));
03984 cpl_matrix_set(*rhs_cov, r, r, cpl_matrix_get(*rhs_cov, r, r)
03985 + obs->star->dmagnitude
03986 * obs->star->dmagnitude);
03987 }
03988
03989 fors_matrix_null(&rhs_input_cov);
03990 fors_matrix_null(&rhs_jacobian);
03991 fors_matrix_null(&rhs_input);
03992 fors_matrix_null(&rhs_jacobian_T);
03993 fors_matrix_null(&tmp_matrix);
03994
03995 return ( cpl_errorstate_is_equal(errstat) ?
03996 CPL_ERROR_NONE :
03997 cpl_error_get_code());
03998 }
03999
04000