fors_photometry_impl.c

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

Generated on Wed Jun 24 14:11:17 2009 for FORS Pipeline Reference Manual by  doxygen 1.4.7