fors_photometry_impl.c

00001 /* $Id: fors_photometry_impl.c,v 1.25 2008/08/07 10:04:21 cizzo Exp $
00002  *
00003  * This file is part of the FORS Library
00004  * Copyright (C) 2002-2006 European Southern Observatory
00005  *
00006  * This program is free software; you can redistribute it and/or modify
00007  * it under the terms of the GNU General Public License as published by
00008  * the Free Software Foundation; either version 2 of the License, or
00009  * (at your option) any later version.
00010  *
00011  * This program is distributed in the hope that it will be useful,
00012  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014  * GNU General Public License for more details.
00015  *
00016  * You should have received a copy of the GNU General Public License
00017  * along with this program; if not, write to the Free Software
00018  * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301 USA
00019  */
00020 
00021 /*
00022  * $Author: cizzo $
00023  * $Date: 2008/08/07 10:04:21 $
00024  * $Revision: 1.25 $
00025  * $Name:  $
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;  /* Counting from zero */
00071     int star_number;   /* Star identification number,
00072               counting from zero */
00073 
00074     double airmass, gain, exptime;
00075     fors_star *star;
00076     bool fit_mag;      /* Fit the magnitude of this star? */
00077 } entry;
00078 
00079 /* Declare and define entry_list */
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     //vprintf(format, al);
00302     va_end(al);
00303 
00304     return;
00305 }
00306 
00307 
00308 /* Internal CPL function, duplicated */
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     /* Initialize all values to zero.
00325        This is done to avoid access of uninitilized memory,  in case
00326        someone passes the matrix to for example cpl_matrix_dump(). */
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     /* The result at (i,j) is the dot-product of i'th and j'th row */
00335     for (i = 0; i < m; i++, bwrite += m, ai += n) {
00336         aj = ai; /* aj points to first entry in j'th row */
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    @brief Linear correlated weighted least squares fit
00353    @param coeff      design matrix (A)
00354    @param rhs        right hand side (b)
00355    @param cov_rhs    covariance of b (C)
00356    @param red_chisq  (output) reduced chi squared, or NULL
00357 
00358    Similar to cpl_matrix_solve_normal, except the output matrix is not
00359    a m x 1 matrix
00360 
00361    x0
00362    x1
00363    x2
00364    :
00365 
00366    but an m x (m+1) matrix
00367 
00368    x0  C00 C01 C02 ...
00369    x1  C10 C11 
00370    x2  C20     . 
00371    :   :        .
00372 
00373    where the first column is the least chi squared solution to the
00374    overdetermined equation system Ax = b, given the covariance, C, of b.
00375 
00376    and the last m columns is the covariance matrix of the solution
00377    (A^t C^-1 A)^-1
00378 
00379    The reduced chi squared is given by
00380    
00381      (b-Ap)^t C^-1 (b-Ap) / (degrees of freedom)
00382 
00383    where  degrees of freedom = |b| - |p|
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;  /* C^-1 */
00394     cpl_matrix * At;    /* A^t */
00395     cpl_matrix * AtC;   /* A^t C^-1 */
00396     cpl_matrix * AtCA;  /* A^t C^-1 A */
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     /* The overall time is probably dominated by this
00404        matrix inversion which is O(n^3) if C is nxn */
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     /* The covariance matrix is singular if one (or more) eigenvalue is
00411        zero, i.e. if there exist a unitary transformation (rotation of
00412        coordinates) that makes the variance of one of the new coordinates
00413        zero (and therefore a failure at this place is probably because the
00414        user has provided some non-positive error).
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     //cpl_matrix_dump(AtCA, stdout);
00439     //cpl_matrix_dump(solution, stdout);
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     /* Get first column vector, p, of solution */
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     /* Input */
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     /* Products */
00520     fors_image *correction        = NULL;
00521 
00522     /* Other */
00523     const char *context = cpl_sprintf("fors.%s", fors_photometry_name);
00524     fors_setting *setting = NULL;
00525 
00526     /* Photometric parameters */
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     /* Find input */
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     /* Get parameteres */
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     /* Get setting */
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     /* Load filter coefficients */
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     /* because the airmass term
00637          x * airmass_j 
00638        can be absorbed in z_j
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     /* Count fitting parameters, print matrix header */
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     /* magnitudes */
00677     for (k = 0; k < n; k++) if (fit_mag[k]) {
00678         myprintf("M%d ", k);
00679         n_parameters++;
00680     }
00681     /* frame zeropoints */
00682     for (k = 0; k < m; k++) if (fit_z || k == 0) {
00683       myprintf("z%d ", k);
00684       n_parameters++;
00685     }
00686     /* extinction */
00687     if (fit_x) {
00688         myprintf("x ");
00689         n_parameters++;
00690     }
00691     /* color coeff. */
00692     if (fit_c) {
00693         myprintf("c ");
00694         n_parameters++;
00695     }
00696     /* color/extinction coupling */
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     /* Consistency check */
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     /* We have all information now, build the matrices */
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     //cpl_matrix_dump(A, stdout);
00738     //cpl_matrix_dump(M, stdout);
00739     //cpl_matrix_dump(covarianceM, stdout);
00740     
00741     /* Code for unweighted fit */
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     /* Solve this overdetermined set of equations in the least chi squared
00759        sense using Cholesky-decomposition, output matrix
00760        is the solution vector (1st column) and the covariance matrix
00761        in the remaining columns.
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     //cpl_matrix_dump(par, stdout);
00773 
00774     /* Print solution, convert to CPL polynomial */
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     /* magnitudes */
00799     for (k = 0; k < n; k++) {
00800 
00801         double input_mag = 9999999; /* color corrected */
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 //printf("BEFORE %p\n", obs);
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 //printf("AFTER %p\n", obs);
00839 
00840     /* frame zeropoints */
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     /* extinction */
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     /* color */
00863     if (fit_c) {
00864         /* Note different sign convention. The values
00865            provided in PHOT_TABLEs are actually -c.
00866            As in fors_std_cat_load() 
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     /* coupling */
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     /* Compute pointwise variance of f, 
00891        which happens to also be a polynomial */
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             /* Add cov_ij to the proper coeff:
00909                cov_ij * dp/d(ai) * dp/d(aj) =
00910                cov_ij * (x^degx[i] * y^degy[i]) * (x^degx[i] * y^degy[i]) =
00911                cov_ij * x^(degx[i]+degx[j]) * y^(degy[i] + degy[j]),
00912                
00913                i.e. add cov_ij to coeff (degx[i]+degx[j], degy[i]+degy[j]) 
00914             */
00915             
00916             int powers[2];
00917             powers[0] = k + kk; /* sum of x degrees */
00918             powers[1] = l + ll; /* sum of y degrees */
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     //cpl_polynomial_dump(f, stdout);
00932     //cpl_polynomial_dump(f_variance, stdout);
00933 
00934 
00935     /* Compute correction image, apply to master flat */
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          * Write QCs
00965          */
00966 
00967         qc = cpl_propertylist_new();
00968 
00969         fors_qc_start_group(qc, fors_qc_dic_version, setting->instrument);
00970 
00971 /*
00972     fors_qc_write_group_heading(cpl_frameset_get_first_const(master_flat_frame),
00973                                 CORRECTION_MAP,
00974                                 setting->instrument);
00975     assure( !cpl_error_get_code(), return, "Could not write %s QC parameters",
00976             CORRECTION_MAP);
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          * End write QCs
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     /* Convert from magnitude to flux.
01044        F = 10^(-0.4 m)
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     /* Normalize to median = 1 */
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   @brief  read all input data from tables
01082   @param  aligned_phot_frames input frames
01083   @param  setting             instrument setting
01084   @param  override_fit_m      whether to always fit magnitudes (true)
01085   @param  n                   (output) number of unique stars
01086   @param  nfit                (output) number of magnitudes to fit
01087   @param  m                   (output) number of frames
01088   @param  fit_mag             (output) array of length n, fit_mag[k] = true   
01089                               iff the magnitude of star number k is to be fitted
01090   @param  to_be_destroyed     This is added as a workaround to the
01091                               impossibility of destroying the returned 
01092                               entry_list, which points to stars in the
01093                               to_be_destroyed standard stars list
01094                               (C.Izzo, 29.01.08).
01095 
01096   @return set of observations
01097 
01098   We match by (RA, DEC) and tolerance, not
01099   by OBJECT string which is not unique across different catalogs.
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     /* Set of unique stars  */
01115    
01116     *to_be_destroyed = NULL;
01117     *fit_mag = NULL;
01118 
01119     /*
01120      * Loop on all aligned photometric tables in input, and count the
01121      * found frames in m.
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          * Load and verify setting for this frame 
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         /* FIXME:
01156          * This whole section is for verifying the input table. The
01157          * structure of this table is probably defined in some other
01158          * place too... Likely this is a dependency between modules
01159          * that should be eliminated. Please check... (C.Izzo, 28.01.08)
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     } /* Table check done */
01202     
01203 
01204     /* 
01205          * Get IDed stars in this table: 
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              /* Use catalog magnitude xor fit the magnitude: */
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                  * Uncomment, to test a mixture of fixed/free magnitudes:
01262                  *
01263                if (i % 2 == 0) fit_m = true;
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, /* ID, cannot identify yet */
01282                         airmass, 
01283                         setting_f->average_gain,
01284                         setting_f->exposure_time,
01285                         star,
01286                         fit_m));
01287 
01288                 /*
01289                  * The star and the standard star are the same star. 
01290                  * The star is detector-oriented, the standard star
01291                  * is its identification and includes physical propeties.
01292                  * Here the information is linked together, and the
01293                  * standard star is "owned" by the star object - so
01294                  * it should not be destroyed.  (C.Izzo, 28.01.08)
01295                  */
01296 
01297         star->id = std_star;
01298 
01299                 /*
01300                  * The standard star may be already in the list of
01301                  * standard stars (stars). We should avoid duplications:
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                      * Only if the nearest standard star is farther 
01313                      * away than 5 arcsecs, insert this standard star 
01314                      * in the standard star list.
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                      * Add the newly found standard star to the list
01340                      */
01341 
01342             fors_std_star_list_insert(stars, std_star);
01343 
01344         } else {
01345 
01346                     /*
01347                      * Destroy it, we have it already in the list
01348                      */
01349 
01350             fors_std_star_delete(&std_star);
01351         }
01352 
01353         } /* if identified */
01354 
01355     }/* for each star */
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     }/* For each table */
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      * For each (unique) standard star, label entries with a running ID 
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                 /* Never fit the magnitude of the first
01397                    star, only differences can be determined */
01398                 e->fit_mag = false;
01399             }
01400             
01401             (*fit_mag)[i] = e->fit_mag;
01402             }
01403             else {
01404             /* The USE_CAT values for the same star in different
01405                ALIGNED_PHOT tables could be inconsistent, 
01406                make sure to synchronize */
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   @brief   build matrix equation
01430   @param   obs          list of star observations
01431   @param   fit_z        fit individual frame zeropoints?
01432   @param   fit_c        fit color coeff.?
01433   @param   fit_x        fit extinction coeff.?
01434   @param   degreef1     polynomial degree
01435   @param   degreef2     polynomial degree
01436   @param   degreep      polynomial degree
01437   @param   n            number of unique stars
01438   @param   m            number of frames
01439   @param   n_parameters number of fit parameters
01440   @param   fit_mag      whether to fit star magnitudes
01441   @param   ext_coeff    extinction coefficient
01442   @param   dext_coeff   ext_coeff error
01443   @param   dcolor_coeff color coefficient error
01444   @param   A            (output) design matrix
01445   @param   M            (output) right hand side
01446   @param   covarianceM  (output) covariance of M
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; /* We cannot predict the number of rows
01463                 in the matrix, resize when necessary */
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         /* Add matrix row */
01482         if (Arows_alloc == Arows) {
01483         /* Amortized constant cost */
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         /* Equation right hand side */
01494         double rhs = e->star->magnitude;  /* instrumental magnitude */
01495         
01496         /* Handle this term:
01497            Ui + c(Ui-Bi) 
01498         */
01499         if (e->fit_mag) {
01500         /* Stay on left hand side */
01501         }
01502         else if (fit_c) {
01503         /* Move just Ui to rhs */
01504         rhs -= e->star->id->cat_magnitude;
01505         }
01506         else {
01507         /* Move color corrected catalog magnitude to rhs */
01508         rhs -= e->star->id->magnitude;
01509         }
01510         
01511         /* Always correct for gain, exptime */
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         /* Compute covariance:
01520 
01521               C(rhs_ij, rhs_kl)   =   
01522              1)     Variance(u_ij) * delta_ik*delta_jk    always
01523              2)   + Variance(Ui)   * delta_ik             if !fit_mag[i==k]
01524              3)   + Variance(Ui-Bi)* c^2 * delta_ik       if !fit_mag[i==k]
01525                   + Covariance(Ui, Ui-Bi) * c * delta_ik
01526              4)   + Variance(c)    * (Ui-Bi)*(Uk-Bk)      if !fit_mag[i] &&
01527                                                           !fit_mag[k]
01528              5)   + Variance(x)    * airmass_j*airmass_l  if !fit_x
01529 
01530         */
01531 
01532         {
01533         const entry *f;
01534         int j = 0;  /* called k in the comment above */
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             /* Diagonal terms */
01543 
01544             /* 1 */
01545             cov += e->star->dmagnitude * e->star->dmagnitude;
01546             
01547             /* This code depends on the loop order through obs
01548                and obs_dup being the same, make this assumption
01549                explicit 
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                 /* 2 
01558                    variance of  Ui
01559                 */
01560                 cov += 
01561                     f->star->id->dcat_magnitude * 
01562                     f->star->id->dcat_magnitude;
01563                 }
01564                 else {
01565                 /* 2, 3, 4
01566                    full variance of   Ui + c(Ui-Bi)
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             /* 4, the diagonal case is handled above */
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             /* 5 */
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         /* Loop through A columns, one column per free parameter */
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         /* magnitudes */
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         /* frame zeropoints */
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             /* Always z0 independent on actual frame number */
01631             val = -1;
01632         }
01633         
01634         myprintf("%f  ", val);
01635         cpl_matrix_set(*A, Arows, par_no++, val);
01636         }
01637 
01638         /* extinction */
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         /* color coeff */
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         /* coupling */
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 

Generated on Wed Sep 10 07:31:52 2008 for FORS Pipeline Reference Manual by  doxygen 1.4.6