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_utils.h>
00040 #include <fors_double.h>
00041
00042 #include <cpl.h>
00043 #include <math.h>
00044 #include <assert.h>
00045
00052 const char *const fors_photometry_name = "fors_photometry";
00053 const char *const fors_photometry_description_short = "Compute corrected flatfield";
00054 const char *const fors_photometry_author = "Jonas M. Larsen";
00055 const char *const fors_photometry_email = PACKAGE_BUGREPORT;
00056 const char *const fors_photometry_description =
00057 "Input files:\n"
00058 " DO category: Type: Explanation: Number:\n"
00059 " ALIGNED_PHOT FITS table Photometry 1+\n"
00060 " MASTER_SKY_FLAT_IMG FITS image Master flat field 1\n"
00061 "\n"
00062 "Output files:\n"
00063 " DO category: Data type: Explanation:\n"
00064 " CORRECTION_MAP FITS image Correction map (magnitude)\n"
00065 " CORRECTION_FACTOR FITS image Correction map (flux)\n"
00066 " MASTER_FLAT_IMG FITS image Corrected master flat field\n";
00067
00068 typedef struct entry
00069 {
00070 int frame_number;
00071 int star_number;
00072
00073
00074 double airmass, gain, exptime;
00075 fors_star *star;
00076 bool fit_mag;
00077 } entry;
00078
00079
00080 #undef LIST_ELEM
00081 #define LIST_ELEM entry
00082 #undef LIST_DEFINE
00083 #include <list.h>
00084 #define LIST_DEFINE
00085 #include <list.h>
00086
00087 const double arcsec_tol = 5.0;
00088
00089 static entry_list *
00090 read_input(const cpl_frameset *aligned_phot_frames,
00091 const fors_setting *setting,
00092 bool override_fit_m,
00093 int *n,
00094 int *nfit,
00095 int *m,
00096 bool **fit_mag,
00097 fors_std_star_list **to_be_destroyed);
00098
00099 static void
00100 build_equations(const entry_list *obs,
00101 bool fit_z, bool fit_c, bool fit_x,
00102 int degreef1, int degreef2, int degreep,
00103 int n, int m,
00104 int n_parameters,
00105 const bool *fit_mag,
00106 double ext_coeff, double dext_coeff,
00107 double dcolor_coeff,
00108 cpl_matrix **A,
00109 cpl_matrix **M,
00110 cpl_matrix **covarianceM);
00111
00117 entry *entry_new(int frame_number, int star_number,
00118 double airmass, double gain, double exptime,
00119 fors_star *star,
00120 bool fit_mag)
00121 {
00122 entry *e = cpl_malloc(sizeof(*e));
00123
00124 e->frame_number = frame_number;
00125 e->star_number = star_number;
00126 e->airmass = airmass;
00127 e->gain = gain;
00128 e->exptime = exptime;
00129 e->star = star;
00130 e->fit_mag = fit_mag;
00131
00132 return e;
00133 }
00134
00139 void
00140 entry_delete(entry **e)
00141 {
00142 if (e && *e) {
00143 fors_star_delete(&(*e)->star);
00144 cpl_free(*e); *e = NULL;
00145 }
00146 return;
00147 }
00148
00153 static void
00154 entry_delete_but_standard(entry **e)
00155 {
00156 if (e && *e) {
00157 fors_star_delete_but_standard(&(*e)->star);
00158 cpl_free(*e); *e = NULL;
00159 }
00160 return;
00161 }
00162
00163
00169 void
00170 entry_list_print(const entry_list *l,
00171 cpl_msg_severity level)
00172 {
00173 const entry *e;
00174 for (e = entry_list_first_const(l);
00175 e != NULL;
00176 e = entry_list_next_const(l)) {
00177
00178 fors_msg(level, "%d, %d: airmass = %f, gain = %f, exptime = %f s",
00179 e->frame_number, e->star_number,
00180 e->airmass, e->gain, e->exptime);
00181
00182 fors_star_print(level, e->star);
00183 }
00184
00185 return;
00186 }
00187
00188
00194 void fors_photometry_define_parameters(cpl_parameterlist *parameters)
00195 {
00196 const char *context = cpl_sprintf("fors.%s", fors_photometry_name);
00197
00198 const char *name, *full_name;
00199 cpl_parameter *p;
00200
00201 name = "fitz";
00202 full_name = cpl_sprintf("%s.%s", context, name);
00203 p = cpl_parameter_new_value(full_name,
00204 CPL_TYPE_BOOL,
00205 "Fit individual frame zeropoints",
00206 context,
00207 true);
00208 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
00209 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00210 cpl_parameterlist_append(parameters, p);
00211 cpl_free((void *)full_name); full_name = NULL;
00212
00213 name = "fitm";
00214 full_name = cpl_sprintf("%s.%s", context, name);
00215 p = cpl_parameter_new_value(full_name,
00216 CPL_TYPE_BOOL,
00217 "Always fit star magnitudes",
00218 context,
00219 true);
00220 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
00221 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00222 cpl_parameterlist_append(parameters, p);
00223 cpl_free((void *)full_name); full_name = NULL;
00224
00225 name = "fitx";
00226 full_name = cpl_sprintf("%s.%s", context, name);
00227 p = cpl_parameter_new_value(full_name,
00228 CPL_TYPE_BOOL,
00229 "Fit extinction coefficient",
00230 context,
00231 false);
00232 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
00233 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00234 cpl_parameterlist_append(parameters, p);
00235 cpl_free((void *)full_name); full_name = NULL;
00236
00237 name = "fitc";
00238 full_name = cpl_sprintf("%s.%s", context, name);
00239 p = cpl_parameter_new_value(full_name,
00240 CPL_TYPE_BOOL,
00241 "Fit color term",
00242 context,
00243 false);
00244 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
00245 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00246 cpl_parameterlist_append(parameters, p);
00247 cpl_free((void *)full_name); full_name = NULL;
00248
00249
00250
00251 name = "degreef1";
00252 full_name = cpl_sprintf("%s.%s", context, name);
00253 p = cpl_parameter_new_value(full_name,
00254 CPL_TYPE_INT,
00255 "correction map polynomial degree (x)",
00256 context,
00257 2);
00258 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
00259 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00260 cpl_parameterlist_append(parameters, p);
00261 cpl_free((void *)full_name); full_name = NULL;
00262
00263
00264 name = "degreef2";
00265 full_name = cpl_sprintf("%s.%s", context, name);
00266 p = cpl_parameter_new_value(full_name,
00267 CPL_TYPE_INT,
00268 "correction map polynomial degree (y), "
00269 "or negative for "
00270 "triangular coefficient matrix",
00271 context,
00272 -1);
00273 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
00274 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00275 cpl_parameterlist_append(parameters, p);
00276 cpl_free((void *)full_name); full_name = NULL;
00277
00278 name = "degreep";
00279 full_name = cpl_sprintf("%s.%s", context, name);
00280 p = cpl_parameter_new_value(full_name,
00281 CPL_TYPE_INT,
00282 "extinction/color coupling degree",
00283 context,
00284 0);
00285 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
00286 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00287 cpl_parameterlist_append(parameters, p);
00288 cpl_free((void *)full_name); full_name = NULL;
00289
00290 cpl_free((void *)context);
00291
00292 return;
00293 }
00294
00295 static void
00296 myprintf(const char *format, ...)
00297 {
00298 va_list al;
00299
00300 va_start(al, format);
00301
00302 va_end(al);
00303
00304 return;
00305 }
00306
00307
00308
00309 static cpl_matrix * matrix_product_normal_create(const cpl_matrix * self)
00310 {
00311 double sum;
00312 cpl_matrix * product;
00313 const double * ai = cpl_matrix_get_data_const(self);
00314 const double * aj;
00315 double * bwrite;
00316 const int m = cpl_matrix_get_nrow(self);
00317 const int n = cpl_matrix_get_ncol(self);
00318 int i, j, k;
00319
00320
00321 cpl_ensure(self != NULL, CPL_ERROR_NULL_INPUT, NULL);
00322
00323 #if 0
00324
00325
00326
00327 product = cpl_matrix_new(m, m);
00328 bwrite = cpl_matrix_get_data(product);
00329 #else
00330 bwrite = (double *) cpl_malloc(m * m * sizeof(double));
00331 product = cpl_matrix_wrap(m, m, bwrite);
00332 #endif
00333
00334
00335 for (i = 0; i < m; i++, bwrite += m, ai += n) {
00336 aj = ai;
00337 for (j = i; j < m; j++, aj += n) {
00338 sum = 0.0;
00339 for (k = 0; k < n; k++) {
00340 sum += ai[k] * aj[k];
00341 }
00342 bwrite[j] = sum;
00343 }
00344 }
00345
00346 return product;
00347 }
00348
00349 #undef cleanup
00350 #define cleanup
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372
00373
00374
00375
00376
00377
00378
00379
00380
00381
00382
00383
00384
00385
00386 static cpl_matrix *
00387 solve_normal(const cpl_matrix *coeff,
00388 const cpl_matrix *rhs,
00389 const cpl_matrix *cov_rhs,
00390 double *red_chisq)
00391 {
00392 cpl_matrix * solution;
00393 cpl_matrix * cov1;
00394 cpl_matrix * At;
00395 cpl_matrix * AtC;
00396 cpl_matrix * AtCA;
00397 cpl_error_code error;
00398
00399 cpl_ensure(coeff != NULL, CPL_ERROR_NULL_INPUT, NULL);
00400 cpl_ensure(rhs != NULL, CPL_ERROR_NULL_INPUT, NULL);
00401 cpl_ensure(cov_rhs != NULL, CPL_ERROR_NULL_INPUT, NULL);
00402
00403
00404
00405 cov1 = cpl_matrix_invert_create(cov_rhs);
00406
00407 assure( cov1 != NULL, return NULL,
00408 "Could not invert covariance matrix. Make sure that provided "
00409 "errors are positive");
00410
00411
00412
00413
00414
00415
00416
00417 At = cpl_matrix_transpose_create(coeff);
00418
00419 AtC = cpl_matrix_product_create(At, cov1);
00420
00421 cpl_matrix_delete(At);
00422
00423 AtCA = cpl_matrix_product_create(AtC, coeff);
00424
00425 solution = cpl_matrix_product_create(AtC, rhs);
00426 cpl_matrix_set_size(solution,
00427 cpl_matrix_get_nrow(solution),
00428 1 + cpl_matrix_get_nrow(solution));
00429 {
00430 int i = 0;
00431 for (i = 0; i < cpl_matrix_get_nrow(solution); i++) {
00432 cpl_matrix_set(solution, i, i+1, 1);
00433 }
00434 }
00435
00436 cpl_matrix_delete(AtC);
00437
00438
00439
00440
00441 error = cpl_matrix_decomp_chol(AtCA);
00442 if (!error) {
00443 error = cpl_matrix_solve_chol(AtCA, solution);
00444 }
00445
00446 cpl_matrix_delete(AtCA);
00447
00448 if (error) {
00449 cpl_matrix_delete(solution);
00450 cpl_ensure(0, error, NULL);
00451 }
00452
00453
00454 if (red_chisq != NULL) {
00455
00456
00457 cpl_matrix *p = cpl_matrix_duplicate(solution);
00458 cpl_matrix_set_size(p, cpl_matrix_get_nrow(p), 1);
00459
00460 cpl_matrix *Ap = cpl_matrix_product_create(coeff, p);
00461
00462 cpl_matrix *bAp = cpl_matrix_duplicate(rhs);
00463 cpl_matrix_subtract(bAp, Ap);
00464
00465 cpl_matrix *bApt = cpl_matrix_transpose_create(bAp);
00466
00467 cpl_matrix *C1bAp = cpl_matrix_product_create(cov1, bAp);
00468
00469 cpl_matrix *chi2 = cpl_matrix_product_create(bApt, C1bAp);
00470
00471 assure( cpl_matrix_get_nrow(chi2) == 1 &&
00472 cpl_matrix_get_ncol(chi2) == 1, return NULL,
00473 "Something is wrong" );
00474
00475 *red_chisq = cpl_matrix_get(chi2, 0, 0) /
00476 (cpl_matrix_get_nrow(rhs) - cpl_matrix_get_nrow(p));
00477
00478 cpl_matrix_delete(p);
00479 cpl_matrix_delete(Ap);
00480 cpl_matrix_delete(bAp);
00481 cpl_matrix_delete(bApt);
00482 cpl_matrix_delete(C1bAp);
00483 cpl_matrix_delete(chi2);
00484 }
00485
00486 cpl_matrix_delete(cov1);
00487
00488 return solution;
00489 }
00490
00491 #undef cleanup
00492 #define cleanup \
00493 do { \
00494 cpl_frameset_delete((cpl_frameset *)aligned_phot_frames); \
00495 cpl_frameset_delete((cpl_frameset *)master_flat_frame); \
00496 cpl_frameset_delete((cpl_frameset *)phot_table); \
00497 fors_setting_delete(&setting); \
00498 fors_image_delete(&master_flat); \
00499 fors_image_delete(&correction); \
00500 } while (0)
00501
00511 void fors_photometry(cpl_frameset *frames, const cpl_parameterlist *parameters)
00512 {
00513
00514 const cpl_frameset *aligned_phot_frames = NULL;
00515 fors_image *master_flat = NULL;
00516 const cpl_frameset *master_flat_frame = NULL;
00517 const cpl_frameset *phot_table = NULL;
00518
00519
00520 fors_image *correction = NULL;
00521
00522
00523 const char *context = cpl_sprintf("fors.%s", fors_photometry_name);
00524 fors_setting *setting = NULL;
00525
00526
00527
00528 cpl_propertylist *qc = NULL;
00529 double qc_zeropoint = -1.0;
00530 double qc_zeropoint_err = -1.0;
00531 double qc_extinction = -1.0;
00532 double qc_extinction_err = -1.0;
00533 double qc_colorterm = -1.0;
00534 double qc_colorterm_err = -1.0;
00535
00536
00537 aligned_phot_frames = fors_frameset_extract(frames, ALIGNED_PHOT);
00538 assure( cpl_frameset_get_size(aligned_phot_frames) > 0, return,
00539 "No %s provided", ALIGNED_PHOT);
00540
00541 master_flat_frame = fors_frameset_extract(frames, MASTER_SKY_FLAT_IMG);
00542 assure( cpl_frameset_get_size(master_flat_frame) > 0, return,
00543 "No %s provided", MASTER_SKY_FLAT_IMG);
00544
00545 phot_table = fors_frameset_extract(frames, PHOT_TABLE);
00546 assure( cpl_frameset_get_size(phot_table) == 1, return,
00547 "One %s required. %d found",
00548 PHOT_TABLE, cpl_frameset_get_size(phot_table));
00549
00550
00551 bool fit_z, override_fit_m, fit_x, fit_c;
00552 int degreef1, degreef2, degreep;
00553
00554 cpl_msg_indent_more();
00555 const char *name = cpl_sprintf("%s.%s", context, "degreef1");
00556 degreef1 = dfs_get_parameter_int_const(parameters,
00557 name);
00558 cpl_free((void *)name); name = NULL;
00559 cpl_msg_indent_less();
00560 assure( !cpl_error_get_code(), return, NULL );
00561
00562
00563 cpl_msg_indent_more();
00564 name = cpl_sprintf("%s.%s", context, "degreef2");
00565 degreef2 = dfs_get_parameter_int_const(parameters,
00566 name);
00567 cpl_free((void *)name); name = NULL;
00568 cpl_msg_indent_less();
00569 assure( !cpl_error_get_code(), return, NULL );
00570
00571
00572
00573 cpl_msg_indent_more();
00574 name = cpl_sprintf("%s.%s", context, "degreep");
00575 degreep = dfs_get_parameter_int_const(parameters,
00576 name);
00577 cpl_free((void *)name); name = NULL;
00578 cpl_msg_indent_less();
00579 assure( !cpl_error_get_code(), return, NULL );
00580
00581
00582 cpl_msg_indent_more();
00583 name = cpl_sprintf("%s.%s", context, "fitz");
00584 fit_z = dfs_get_parameter_bool_const(parameters,
00585 name);
00586 cpl_free((void *)name); name = NULL;
00587 cpl_msg_indent_less();
00588 assure( !cpl_error_get_code(), return, NULL );
00589
00590
00591 cpl_msg_indent_more();
00592 name = cpl_sprintf("%s.%s", context, "fitm");
00593 override_fit_m = dfs_get_parameter_bool_const(parameters,
00594 name);
00595 cpl_free((void *)name); name = NULL;
00596 cpl_msg_indent_less();
00597 assure( !cpl_error_get_code(), return, NULL );
00598
00599
00600 cpl_msg_indent_more();
00601 name = cpl_sprintf("%s.%s", context, "fitx");
00602 fit_x = dfs_get_parameter_bool_const(parameters,
00603 name);
00604 cpl_free((void *)name); name = NULL;
00605 cpl_msg_indent_less();
00606 assure( !cpl_error_get_code(), return, NULL );
00607
00608
00609 cpl_msg_indent_more();
00610 name = cpl_sprintf("%s.%s", context, "fitc");
00611 fit_c = dfs_get_parameter_bool_const(parameters,
00612 name);
00613 cpl_free((void *)name); name = NULL;
00614 cpl_msg_indent_less();
00615 assure( !cpl_error_get_code(), return, NULL );
00616
00617 cpl_free((void *)context); context = NULL;
00618
00619
00620 setting = fors_setting_new(cpl_frameset_get_first_const(aligned_phot_frames));
00621 assure( !cpl_error_get_code(), return, "Could not get instrument setting" );
00622
00623
00624 double color_coeff, dcolor_coeff;
00625 double ext_coeff, dext_coeff;
00626 fors_phot_table_load(cpl_frameset_get_first_const(phot_table), setting,
00627 &color_coeff, &dcolor_coeff,
00628 &ext_coeff, &dext_coeff,
00629 NULL, NULL);
00630 assure( !cpl_error_get_code(), return, "Could not load photometry table" );
00631
00632 if (fit_x && fit_z) {
00633 cpl_msg_warning(cpl_func,
00634 "Fitting extinction coefficient "
00635 "and frame zeropoints simultaneously will probably fail");
00636
00637
00638
00639
00640 }
00641
00642
00643 int n, nfit, m;
00644 bool *fit_mag;
00645 fors_std_star_list *to_be_destroyed;
00646 entry_list *obs = read_input(aligned_phot_frames, setting, override_fit_m,
00647 &n, &nfit, &m, &fit_mag, &to_be_destroyed);
00648 assure( !cpl_error_get_code(), return, "Could not read input tables" );
00649
00650 entry_list_print(obs, CPL_MSG_DEBUG);
00651
00652 {
00653 int ntot = entry_list_size(obs);
00654 cpl_msg_info(cpl_func,
00655 "Found %d table%s, %d star%s, %d unique star%s, "
00656 "%d magnitude%s to determine",
00657 m, m != 1 ? "s" : "",
00658 ntot, ntot != 1 ? "s" : "",
00659 n, n != 1 ? "s" : "",
00660 nfit, nfit != 1 ? "s" : "");
00661 }
00662
00663
00664 int n_parameters = 0;
00665 {
00666 int k, l;
00667 for (k = 0; k <= degreef1; k++) {
00668 for (l = 0;
00669 (degreef2 >= 0) ? (l <= degreef2) : (k + l <= degreef1);
00670 l++) if (k+l > 0) {
00671 myprintf("%5s%d,%d ", "f", k, l);
00672 n_parameters++;
00673 }
00674 }
00675
00676
00677 for (k = 0; k < n; k++) if (fit_mag[k]) {
00678 myprintf("M%d ", k);
00679 n_parameters++;
00680 }
00681
00682 for (k = 0; k < m; k++) if (fit_z || k == 0) {
00683 myprintf("z%d ", k);
00684 n_parameters++;
00685 }
00686
00687 if (fit_x) {
00688 myprintf("x ");
00689 n_parameters++;
00690 }
00691
00692 if (fit_c) {
00693 myprintf("c ");
00694 n_parameters++;
00695 }
00696
00697 for (k = 0; k <= degreep; k++) {
00698 for (l = 0; k + l <= degreep; l++) if (k+l > 1) {
00699 myprintf("%5s%d,%d ", "p", k, l);
00700 n_parameters++;
00701 }
00702 }
00703
00704
00705 myprintf("rhs");
00706 myprintf("\n");
00707 }
00708
00709
00710 {
00711 int z_pars = fit_z ? m : 1;
00712 int x_pars = fit_x ? 1 : 0;
00713 int c_pars = fit_c ? 1 : 0;
00714 int p_pars = degreep >= 2 ? ((degreep+1)*(degreep+2))/2-3 : 0;
00715
00716 assure( (degreef2 >= 0 && n_parameters ==
00717 (degreef1+1)*(degreef2+1)-1 + nfit + z_pars + x_pars + c_pars + p_pars)
00718 ||
00719 (degreef2 < 0 && n_parameters ==
00720 ((degreef1+1)*(degreef1+2))/2-1 + nfit + z_pars + x_pars + c_pars + p_pars),
00721 return, "Something is wrong");
00722 }
00723
00724
00725 cpl_matrix *A;
00726 cpl_matrix *M;
00727 cpl_matrix *covarianceM;
00728 build_equations(obs,
00729 fit_z, fit_c, fit_x,
00730 degreef1, degreef2, degreep,
00731 n, m, n_parameters, fit_mag,
00732 ext_coeff, dext_coeff, dcolor_coeff,
00733 &A, &M, &covarianceM);
00734
00735 assure( !cpl_error_get_code(), return, "Could not setup equation system" );
00736
00737
00738
00739
00740
00741
00742 if (0){
00743 int i, j;
00744 for (i = 0; i < cpl_matrix_get_nrow(covarianceM); i++) {
00745 for (j = 0; j < cpl_matrix_get_nrow(covarianceM); j++) {
00746 cpl_matrix_set(covarianceM, i, j, (i == j) ? 1 : 0);
00747 }
00748 }
00749 }
00750
00751 {
00752 int Arows = cpl_matrix_get_nrow(A);
00753 cpl_msg_info(cpl_func, "Solving %d equation%s for %d parameter%s",
00754 Arows, Arows != 1 ? "s" : "",
00755 n_parameters, n_parameters != 1 ? "s" : "");
00756 }
00757
00758
00759
00760
00761
00762
00763
00764 double red_chisq;
00765 cpl_matrix *par = solve_normal(A, M, covarianceM, &red_chisq);
00766
00767 cpl_matrix_delete(A); A = NULL;
00768 cpl_matrix_delete(M); M = NULL;
00769 cpl_matrix_delete(covarianceM); covarianceM = NULL;
00770
00771 assure( !cpl_error_get_code(), return, "Could not solve equation system");
00772
00773
00774
00775 {
00776 cpl_polynomial *f = cpl_polynomial_new(2);
00777 int k, l;
00778 int i = 0;
00779
00780 for (k = 0; k <= degreef1; k++) {
00781 for (l = 0;
00782 (degreef2 >= 0) ? (l <= degreef2) : (k + l <= degreef1);
00783 l++) if (k+l > 0) {
00784 double dcoeff = sqrt(cpl_matrix_get(par, i, i+1));
00785 double coeff = cpl_matrix_get(par, i, 0);
00786 i++;
00787
00788 cpl_msg_info(cpl_func, "%s%d,%d = %f +- %f", "f", k, l,
00789 coeff, dcoeff);
00790
00791 int powers[2];
00792 powers[0] = k;
00793 powers[1] = l;
00794 cpl_polynomial_set_coeff(f, powers, coeff);
00795 }
00796 }
00797
00798
00799 for (k = 0; k < n; k++) {
00800
00801 double input_mag = 9999999;
00802 const char *object = "???";
00803
00804 entry *e;
00805 for (e = entry_list_first(obs);
00806 e != NULL;
00807 e = entry_list_next(obs)) {
00808
00809 if (e->star_number == k) {
00810
00811 input_mag = e->star->id->magnitude;
00812 object = e->star->id->name;
00813 break;
00814 }
00815 }
00816
00817 if (!fit_mag[k]) {
00818 double dmag = 0;
00819 cpl_msg_info(cpl_func, "M%d = %f +- %f mag (fixed, %s)",
00820 k, input_mag, dmag, object);
00821 }
00822 else {
00823 cpl_msg_info(cpl_func, "M%d = %f +- %f mag (free, %s: input: %f mag)",
00824 k,
00825 cpl_matrix_get(par, i, 0),
00826 sqrt(cpl_matrix_get(par, i, i+1)),
00827 object,
00828 input_mag);
00829 i++;
00830 }
00831 }
00832
00833
00834 fors_std_star_list_delete(&to_be_destroyed, fors_std_star_delete);
00835 to_be_destroyed = NULL;
00836 entry_list_delete(&obs, entry_delete_but_standard); obs = NULL;
00837 cpl_free(fit_mag); fit_mag = NULL;
00838
00839
00840
00841 for (k = 0; k < m; k++) if (fit_z || k == 0) {
00842 if (k == 0 && !fit_z) {
00843 qc_zeropoint = cpl_matrix_get(par, i, 0);
00844 qc_zeropoint_err = sqrt(cpl_matrix_get(par, i, i+1));
00845 }
00846 cpl_msg_info(cpl_func, "z%d = %f +- %f mag", k,
00847 cpl_matrix_get(par, i, 0),
00848 sqrt(cpl_matrix_get(par, i, i+1)));
00849 i++;
00850 }
00851
00852
00853 if (fit_x) {
00854 qc_extinction = cpl_matrix_get(par, i, 0);
00855 qc_extinction_err = sqrt(cpl_matrix_get(par, i, i+1));
00856 cpl_msg_info(cpl_func, "x = %f +- %f mag/airmass",
00857 cpl_matrix_get(par, i, 0),
00858 sqrt(cpl_matrix_get(par, i, i+1)));
00859 i++;
00860 }
00861
00862
00863 if (fit_c) {
00864
00865
00866
00867
00868 qc_colorterm = cpl_matrix_get(par, i, 0);
00869 qc_colorterm_err = sqrt(cpl_matrix_get(par, i, i+1));
00870 cpl_msg_info(cpl_func, "c = %f +- %f",
00871 cpl_matrix_get(par, i, 0),
00872 sqrt(cpl_matrix_get(par, i, i+1)));
00873 i++;
00874 }
00875
00876
00877 for (k = 0; k <= degreep; k++) {
00878 for (l = 0; k + l <= degreep; l++) if (k+l > 1) {
00879 double dcoeff = sqrt(cpl_matrix_get(par, i, i+1));
00880 double coeff = cpl_matrix_get(par, i, 0);
00881 i++;
00882
00883 cpl_msg_info(cpl_func, "%s%d,%d = %f +- %f", "p", k, l,
00884 coeff, dcoeff);
00885 }
00886 }
00887
00888 cpl_msg_info(cpl_func, "Reduced chi square = %f", red_chisq);
00889
00890
00891
00892
00893 cpl_polynomial *f_variance = cpl_polynomial_new(2);
00894
00895 i = 0;
00896 for (k = 0; k <= degreef1; k++) {
00897 for (l = 0;
00898 (degreef2 >= 0) ? (l <= degreef2) : (k + l <= degreef1);
00899 l++) if (k+l > 0) {
00900
00901 int j = 0;
00902 int kk, ll;
00903 for (kk = 0; kk <= degreef1; kk++) {
00904 for (ll = 0;
00905 (degreef2 >= 0) ? (ll <= degreef2) : (kk + ll <= degreef1);
00906 ll++) if (kk+ll > 0) {
00907
00908
00909
00910
00911
00912
00913
00914
00915
00916 int powers[2];
00917 powers[0] = k + kk;
00918 powers[1] = l + ll;
00919
00920 double coeff = cpl_polynomial_get_coeff(f_variance, powers);
00921 cpl_polynomial_set_coeff(f_variance, powers,
00922 coeff +
00923 cpl_matrix_get(par, i, 1+j));
00924 j++;
00925 }
00926 }
00927 i++;
00928 }
00929 }
00930
00931
00932
00933
00934
00935
00936 master_flat = fors_image_load(cpl_frameset_get_first_const(master_flat_frame),
00937 NULL, setting, NULL);
00938 assure( !cpl_error_get_code(), return, "Could not load master flat");
00939
00940 int nx = fors_image_get_size_x(master_flat);
00941 int ny = fors_image_get_size_y(master_flat);
00942 cpl_image *correction_map = cpl_image_new(nx, ny, CPL_TYPE_FLOAT);
00943 cpl_image *correction_map_v = cpl_image_new(nx, ny, CPL_TYPE_FLOAT);
00944
00945 cpl_msg_info(cpl_func, "Creating correction map (magnitude)");
00946 cpl_image_fill_polynomial(correction_map, f,
00947 1.0, 1.0, 1.0, 1.0);
00948 cpl_image_fill_polynomial(correction_map_v, f_variance,
00949 1.0, 1.0, 1.0, 1.0);
00950
00951 correction = fors_image_new(correction_map, correction_map_v);
00952
00953 cpl_polynomial_delete(f); f = NULL;
00954 cpl_polynomial_delete(f_variance); f_variance = NULL;
00955 }
00956
00957 cpl_matrix_delete(par); par = NULL;
00958
00959 if (qc_zeropoint_err > 0.0 ||
00960 qc_extinction_err > 0.0 ||
00961 qc_colorterm_err > 0.0) {
00962
00963
00964
00965
00966
00967 qc = cpl_propertylist_new();
00968
00969 fors_qc_start_group(qc, fors_qc_dic_version, setting->instrument);
00970
00971
00972
00973
00974
00975
00976
00977
00978
00979 if (qc_zeropoint_err > 0.0) {
00980 fors_qc_write_qc_double(qc,
00981 qc_zeropoint,
00982 "QC.INSTRUMENT.ZEROPOINT",
00983 "mag",
00984 "Instrument zeropoint",
00985 setting->instrument);
00986
00987 fors_qc_write_qc_double(qc,
00988 qc_zeropoint_err,
00989 "QC.INSTRUMENT.ZEROPOINT.ERROR",
00990 "mag",
00991 "Instrument zeropoint error",
00992 setting->instrument);
00993 }
00994
00995 if (qc_extinction_err > 0.0) {
00996 fors_qc_write_qc_double(qc,
00997 qc_extinction,
00998 "QC.ATMOSPHERIC.EXTINCTION",
00999 "mag/airmass",
01000 "Atmospheric extinction",
01001 setting->instrument);
01002
01003 fors_qc_write_qc_double(qc,
01004 qc_extinction_err,
01005 "QC.ATMOSPHERIC.EXTINCTION.ERROR",
01006 "mag/airmass",
01007 "Atmospheric extinction error",
01008 setting->instrument);
01009 }
01010
01011 if (qc_colorterm_err > 0.0) {
01012 fors_qc_write_qc_double(qc,
01013 qc_colorterm,
01014 "QC.COLOR.CORRECTION",
01015 NULL,
01016 "Linear color correction term",
01017 setting->instrument);
01018
01019 fors_qc_write_qc_double(qc,
01020 qc_colorterm_err,
01021 "QC.COLOR.CORRECTION.ERROR",
01022 NULL,
01023 "Linear color correction term error",
01024 setting->instrument);
01025 }
01026
01027 fors_qc_end_group();
01028
01029
01030
01031
01032 }
01033
01034 fors_dfs_save_image(frames, correction, CORRECTION_MAP,
01035 qc, parameters, fors_photometry_name,
01036 cpl_frameset_get_first_const(master_flat_frame));
01037
01038 cpl_propertylist_delete(qc);
01039
01040 assure( !cpl_error_get_code(), return, "Saving %s failed",
01041 CORRECTION_MAP);
01042
01043
01044
01045
01046 cpl_msg_info(cpl_func, "Creating correction map (flux)");
01047
01048 fors_image_multiply_scalar(correction, -0.4, -1);
01049 fors_image_exponential(correction, 10, -1);
01050
01051
01052 fors_image_divide_scalar(correction,
01053 fors_image_get_median(correction, NULL), -1.0);
01054
01055 fors_dfs_save_image(frames, correction, CORRECTION_FACTOR,
01056 NULL, parameters, fors_photometry_name,
01057 cpl_frameset_get_first_const(master_flat_frame));
01058
01059 assure( !cpl_error_get_code(), return, "Saving %s failed",
01060 CORRECTION_FACTOR);
01061
01062
01063 cpl_msg_info(cpl_func, "Creating corrected master flat");
01064 fors_image_multiply(master_flat, correction);
01065
01066 fors_dfs_save_image(frames, master_flat, MASTER_FLAT_IMG,
01067 NULL, parameters, fors_photometry_name,
01068 cpl_frameset_get_first_const(master_flat_frame));
01069 assure( !cpl_error_get_code(), return, "Saving %s failed",
01070 MASTER_FLAT_IMG);
01071
01072 cleanup;
01073 return;
01074 }
01075
01076
01077 #undef cleanup
01078 #define cleanup
01079
01080
01081
01082
01083
01084
01085
01086
01087
01088
01089
01090
01091
01092
01093
01094
01095
01096
01097
01098
01099
01100
01101 static entry_list *
01102 read_input(const cpl_frameset *aligned_phot_frames,
01103 const fors_setting *setting,
01104 bool override_fit_m,
01105 int *n,
01106 int *nfit,
01107 int *m,
01108 bool **fit_mag,
01109 fors_std_star_list **to_be_destroyed)
01110 {
01111 entry_list *obs = entry_list_new();
01112
01113 fors_std_star_list *stars = fors_std_star_list_new();
01114
01115
01116 *to_be_destroyed = NULL;
01117 *fit_mag = NULL;
01118
01119
01120
01121
01122
01123
01124 const cpl_frame *f;
01125 for (f = cpl_frameset_get_first_const(aligned_phot_frames), *m = 0;
01126 f != NULL;
01127 f = cpl_frameset_get_next_const(aligned_phot_frames), (*m)++) {
01128
01129 const char *filename = cpl_frame_get_filename(f);
01130
01131 cpl_msg_info(cpl_func, "Loading %s", filename);
01132
01133 cpl_table *aligned_phot = cpl_table_load(filename, 1, 1);
01134
01135 assure( !cpl_error_get_code(), return NULL, "Could not load %s",
01136 filename );
01137
01138 cpl_propertylist *header = cpl_propertylist_load(filename, 0);
01139
01140 assure( !cpl_error_get_code(), return NULL, "Could not load %s header",
01141 filename );
01142
01143
01144
01145
01146
01147 fors_setting *setting_f;
01148 fors_setting_verify(setting, f, &setting_f);
01149
01150 double airmass = cpl_propertylist_get_double(header, "AIRMASS");
01151
01152 assure( !cpl_error_get_code(), return NULL, "%s: Could not read %s",
01153 filename, "AIRMASS" );
01154
01155
01156
01157
01158
01159
01160
01161
01162 struct {
01163 const char *name;
01164 cpl_type type;
01165 } col[] =
01166 {{"RA", CPL_TYPE_DOUBLE},
01167 {"DEC", CPL_TYPE_DOUBLE},
01168 {"X", CPL_TYPE_DOUBLE},
01169 {"Y", CPL_TYPE_DOUBLE},
01170 {"FWHM", CPL_TYPE_DOUBLE},
01171 {"A", CPL_TYPE_DOUBLE},
01172 {"B", CPL_TYPE_DOUBLE},
01173 {"THETA", CPL_TYPE_DOUBLE},
01174 {"INSTR_MAG", CPL_TYPE_DOUBLE},
01175 {"DINSTR_MAG", CPL_TYPE_DOUBLE},
01176 {"CAT_MAG", CPL_TYPE_DOUBLE},
01177 {"DCAT_MAG", CPL_TYPE_DOUBLE},
01178 {"MAG", CPL_TYPE_DOUBLE},
01179 {"DMAG", CPL_TYPE_DOUBLE},
01180 {"CLASS_STAR", CPL_TYPE_DOUBLE},
01181 {"USE_CAT", CPL_TYPE_INT},
01182 {"OBJECT", CPL_TYPE_STRING}};
01183 {
01184 unsigned i = 0;
01185 for (i = 0; i < sizeof(col) / sizeof(*col); i++) {
01186 assure( cpl_table_has_column(aligned_phot, col[i].name),
01187 return NULL,
01188 "%s: Missing column %s", filename, col[i].name);
01189 assure( cpl_table_get_column_type(aligned_phot, col[i].name)
01190 == col[i].type,
01191 return NULL,
01192 "%s: column %s: Type is %s, %s expected ",
01193 filename,
01194 col[i].name,
01195 fors_type_get_string(
01196 cpl_table_get_column_type(aligned_phot,
01197 col[i].name)
01198 ),
01199 fors_type_get_string(col[i].type));
01200 }
01201 }
01202
01203
01204
01205
01206
01207
01208 int i;
01209 for (i = 0; i < cpl_table_get_nrow(aligned_phot); i++) {
01210 int null;
01211 double ra = cpl_table_get_double(aligned_phot,
01212 "RA", i, &null);
01213 double dec = cpl_table_get_double(aligned_phot,
01214 "DEC", i, &null);
01215 double x = cpl_table_get_double(aligned_phot,
01216 "X", i, &null);
01217 double y = cpl_table_get_double(aligned_phot,
01218 "Y", i, &null);
01219 double fwhm = cpl_table_get_double(aligned_phot,
01220 "FWHM", i, &null);
01221 double smaj = cpl_table_get_double(aligned_phot,
01222 "A", i, &null);
01223 double smin = cpl_table_get_double(aligned_phot,
01224 "B", i, &null);
01225 double ornt = cpl_table_get_double(aligned_phot,
01226 "THETA", i, &null);
01227 double mag = cpl_table_get_double(aligned_phot,
01228 "INSTR_MAG", i, &null);
01229 double dmag = cpl_table_get_double(aligned_phot,
01230 "DINSTR_MAG", i, &null);
01231 double indx = cpl_table_get_double(aligned_phot,
01232 "CLASS_STAR", i, &null);
01233 double cat_mag = cpl_table_get_double(aligned_phot,
01234 "CAT_MAG", i, &null);
01235 double dcat_mag = cpl_table_get_double(aligned_phot,
01236 "DCAT_MAG", i, &null);
01237 double cmag = cpl_table_get_double(aligned_phot,
01238 "MAG", i, &null);
01239 double dcmag = cpl_table_get_double(aligned_phot,
01240 "DMAG", i, &null);
01241 double color = cpl_table_get_double(aligned_phot,
01242 "COLOR", i, &null);
01243
01244 bool fit_m = override_fit_m ||
01245 (cpl_table_get_int(aligned_phot,
01246 "USE_CAT", i, &null) == 0);
01247
01248 const char *object = cpl_table_get_string(aligned_phot,
01249 "OBJECT", i);
01250
01251 cpl_msg_debug(cpl_func, "%s: %f, %f",
01252 object != NULL ? object : "-", ra, dec);
01253
01254 assure( dmag > 0, return NULL,
01255 "Non-positive magnitude error: %f mag at row %d",
01256 dmag, i+1 );
01257
01258 if (object != NULL) {
01259
01260
01261
01262
01263
01264
01265
01266 fors_std_star *std_star =
01267 fors_std_star_new(ra, dec,
01268 cmag, dcmag,
01269 cat_mag, dcat_mag,
01270 color, object);
01271
01272 fors_star *star = fors_star_new(x, y,
01273 fwhm,
01274 smaj, smin,
01275 ornt,
01276 mag, dmag,
01277 indx);
01278
01279 entry_list_insert(obs,
01280 entry_new(*m,
01281 -1,
01282 airmass,
01283 setting_f->average_gain,
01284 setting_f->exposure_time,
01285 star,
01286 fit_m));
01287
01288
01289
01290
01291
01292
01293
01294
01295
01296
01297 star->id = std_star;
01298
01299
01300
01301
01302
01303
01304 bool is_new = false;
01305
01306 if (fors_std_star_list_size(stars) == 0) {
01307 is_new = true;
01308 }
01309 else {
01310
01311
01312
01313
01314
01315
01316
01317 fors_std_star *nearest =
01318 fors_std_star_list_kth_val(
01319 stars, 1,
01320 (fors_std_star_list_func_eval)
01321 fors_std_star_dist_arcsec,
01322 std_star);
01323
01324 cpl_msg_debug(cpl_func, "dist = %f arcseconds",
01325 fors_std_star_dist_arcsec(nearest, std_star));
01326
01327 if (fors_std_star_dist_arcsec(nearest, std_star)
01328 > arcsec_tol) {
01329 is_new = true;
01330 }
01331 else {
01332 star->id = nearest;
01333 }
01334 }
01335
01336 if (is_new) {
01337
01338
01339
01340
01341
01342 fors_std_star_list_insert(stars, std_star);
01343
01344 } else {
01345
01346
01347
01348
01349
01350 fors_std_star_delete(&std_star);
01351 }
01352
01353 }
01354
01355 }
01356
01357 cpl_propertylist_delete(header); header = NULL;
01358 cpl_table_delete(aligned_phot); aligned_phot = NULL;
01359 fors_setting_delete(&setting_f); setting_f = NULL;
01360
01361 }
01362
01363
01364 *to_be_destroyed = stars;
01365 *n = fors_std_star_list_size(stars);
01366
01367 assure( entry_list_size(obs) > 0, return NULL, "No stars found");
01368
01369 *fit_mag = cpl_calloc(*n, sizeof(**fit_mag));
01370
01371
01372
01373
01374
01375 {
01376 fors_std_star *ref;
01377 int i;
01378
01379 for (ref = fors_std_star_list_first(stars), i = 0;
01380 ref != NULL;
01381 ref = fors_std_star_list_next(stars), i++) {
01382
01383 bool is_first = true;
01384 entry *e;
01385 for (e = entry_list_first(obs);
01386 e != NULL;
01387 e = entry_list_next(obs)) {
01388
01389 if (e->star->id == ref) {
01390 e->star_number = i;
01391
01392 if (is_first) {
01393 is_first = false;
01394
01395 if (i == 0) {
01396
01397
01398 e->fit_mag = false;
01399 }
01400
01401 (*fit_mag)[i] = e->fit_mag;
01402 }
01403 else {
01404
01405
01406
01407 e->fit_mag = (*fit_mag)[i];
01408 }
01409 }
01410 }
01411 }
01412
01413 assure( i == *n, return NULL, "Something is wrong" );
01414 }
01415
01416 {
01417 *nfit = 0;
01418
01419 int k;
01420 for (k = 0; k < *n; k++) {
01421 if ((*fit_mag)[k]) (*nfit)++;
01422 }
01423 }
01424
01425 return obs;
01426 }
01427
01428
01429
01430
01431
01432
01433
01434
01435
01436
01437
01438
01439
01440
01441
01442
01443
01444
01445
01446
01447
01448 static void
01449 build_equations(const entry_list *obs,
01450 bool fit_z, bool fit_c, bool fit_x,
01451 int degreef1, int degreef2, int degreep,
01452 int n, int m,
01453 int n_parameters,
01454 const bool *fit_mag,
01455 double ext_coeff, double dext_coeff,
01456 double dcolor_coeff,
01457 cpl_matrix **A,
01458 cpl_matrix **M,
01459 cpl_matrix **covarianceM)
01460 {
01461 int Arows = 0;
01462 int Arows_alloc = 1;
01463
01464 *A = cpl_matrix_new(Arows_alloc, n_parameters);
01465 *M = cpl_matrix_new(Arows_alloc, 1);
01466
01467 *covarianceM = cpl_matrix_new(entry_list_size(obs),
01468 entry_list_size(obs));
01469
01470 {
01471 entry_list *obs_dup = entry_list_duplicate(obs, NULL);
01472
01473 const entry *e;
01474 for (e = entry_list_first_const(obs);
01475 e != NULL;
01476 e = entry_list_next_const(obs)) {
01477
01478 assure( e->star_number >= 0, return,
01479 "Star not identified, should not happen" );
01480
01481
01482 if (Arows_alloc == Arows) {
01483
01484 Arows_alloc *= 2;
01485 cpl_matrix_set_size(*A,
01486 Arows_alloc,
01487 cpl_matrix_get_ncol(*A));
01488 cpl_matrix_set_size(*M,
01489 Arows_alloc,
01490 1);
01491 }
01492
01493
01494 double rhs = e->star->magnitude;
01495
01496
01497
01498
01499 if (e->fit_mag) {
01500
01501 }
01502 else if (fit_c) {
01503
01504 rhs -= e->star->id->cat_magnitude;
01505 }
01506 else {
01507
01508 rhs -= e->star->id->magnitude;
01509 }
01510
01511
01512 rhs += 2.5*log10(e->gain);
01513 rhs += 2.5*log10(e->exptime);
01514
01515 if (!fit_x) {
01516 rhs -= e->airmass * ext_coeff;
01517 }
01518
01519
01520
01521
01522
01523
01524
01525
01526
01527
01528
01529
01530
01531
01532 {
01533 const entry *f;
01534 int j = 0;
01535 for (f = entry_list_first_const(obs_dup), j = 0;
01536 f != NULL;
01537 f = entry_list_next_const(obs_dup), j++) {
01538
01539 double cov = 0;
01540
01541 if (f == e) {
01542
01543
01544
01545 cov += e->star->dmagnitude * e->star->dmagnitude;
01546
01547
01548
01549
01550
01551 assure( Arows == j, return, "Something is wrong");
01552 }
01553
01554 if (f->star_number == e->star_number) {
01555 if (!fit_mag[f->star_number]) {
01556 if (fit_c) {
01557
01558
01559
01560 cov +=
01561 f->star->id->dcat_magnitude *
01562 f->star->id->dcat_magnitude;
01563 }
01564 else {
01565
01566
01567
01568
01569 cov +=
01570 f->star->id->dmagnitude *
01571 f->star->id->dmagnitude;
01572 }
01573 }
01574 }
01575
01576 if (!fit_mag[e->star_number] && !fit_mag[f->star_number] && !fit_c) {
01577
01578 if (f->star_number != e->star_number) {
01579 cov +=
01580 dcolor_coeff*dcolor_coeff *
01581 e->star->id->color *
01582 f->star->id->color;
01583 }
01584 }
01585
01586
01587 if (!fit_x) {
01588 cov +=
01589 dext_coeff*dext_coeff *
01590 e->airmass *
01591 f->airmass;
01592 }
01593
01594 cpl_matrix_set(*covarianceM, Arows, j, cov);
01595 }
01596 }
01597
01598
01599 int par_no = 0;
01600 int k, l;
01601 for (k = 0; k <= degreef1; k++) {
01602 for (l = 0;
01603 (degreef2 >= 0) ? (l <= degreef2) : (k + l <= degreef1);
01604 l++) if (k+l > 0) {
01605 double val =
01606 pow(e->star->pixel->x, k)*
01607 pow(e->star->pixel->y, l);
01608 myprintf("%.2e ", val);
01609 cpl_matrix_set(*A, Arows, par_no++, val);
01610 }
01611 }
01612
01613
01614 for (k = 0; k < n; k++) if (fit_mag[k]) {
01615
01616 double val = (k == e->star_number) ? 1 : 0;
01617
01618 myprintf("%f ", val);
01619 cpl_matrix_set(*A, Arows, par_no++, val);
01620 }
01621
01622
01623 for (k = 0; k < m; k++) if (fit_z || k == 0) {
01624
01625 double val;
01626 if (fit_z) {
01627 val = (k == e->frame_number) ? -1 : 0;
01628 }
01629 else {
01630
01631 val = -1;
01632 }
01633
01634 myprintf("%f ", val);
01635 cpl_matrix_set(*A, Arows, par_no++, val);
01636 }
01637
01638
01639 if (fit_x) {
01640 double val = e->airmass;
01641
01642 myprintf("%f ", val);
01643 cpl_matrix_set(*A, Arows, par_no++, val);
01644 }
01645
01646
01647 if (fit_c) {
01648 double val = e->star->id->color;
01649
01650 myprintf("%f ", val);
01651 cpl_matrix_set(*A, Arows, par_no++, val);
01652 }
01653
01654
01655 for (k = 0; k <= degreep; k++) {
01656 for (l = 0; k + l <= degreep; l++) if (k+l > 1) {
01657 double val =
01658 pow(e->airmass, k)*
01659 pow(e->star->id->color, l);
01660 myprintf("%.2e ", val);
01661 cpl_matrix_set(*A, Arows, par_no++, val);
01662 }
01663 }
01664
01665 assure( !cpl_error_get_code(), return, NULL);
01666
01667 myprintf("%g", rhs);
01668 cpl_matrix_set(*M, Arows, 0, rhs);
01669
01670 myprintf("\n");
01671 Arows++;
01672
01673 cpl_msg_debug(cpl_func,
01674 "j = %d; i = %d; (x, y) = (%.2f, %.2f); "
01675 "m = %.2f +- %.2f; airmass = %f; "
01676 "exptime = %.2f s; %s",
01677 e->frame_number, e->star_number,
01678 e->star->pixel->x,
01679 e->star->pixel->y,
01680 e->star->magnitude, e->star->dmagnitude,
01681 e->airmass,
01682 e->exptime,
01683 e->star->id->name);
01684 }
01685 entry_list_delete(&obs_dup, NULL); obs_dup = NULL;
01686 }
01687
01688 Arows_alloc = Arows;
01689 cpl_matrix_set_size(*A,
01690 Arows_alloc,
01691 cpl_matrix_get_ncol(*A));
01692 cpl_matrix_set_size(*M,
01693 Arows_alloc,
01694 1);
01695 return;
01696 }
01697