fors_zeropoint_impl.c

00001 /* $Id: fors_zeropoint_impl.c,v 1.68 2008/08/27 12:55:29 cizzo Exp $
00002  *
00003  * This file is part of the FORS Library
00004  * Copyright (C) 2002-2006 European Southern Observatory
00005  *
00006  * This program is free software; you can redistribute it and/or modify
00007  * it under the terms of the GNU General Public License as published by
00008  * the Free Software Foundation; either version 2 of the License, or
00009  * (at your option) any later version.
00010  *
00011  * This program is distributed in the hope that it will be useful,
00012  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014  * GNU General Public License for more details.
00015  *
00016  * You should have received a copy of the GNU General Public License
00017  * along with this program; if not, write to the Free Software
00018  * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301 USA
00019  */
00020 
00021 /*
00022  * $Author: cizzo $
00023  * $Date: 2008/08/27 12:55:29 $
00024  * $Revision: 1.68 $
00025  * $Name:  $
00026  */
00027 
00028 #ifdef HAVE_CONFIG_H
00029 #include <config.h>
00030 #endif
00031 
00032 #include <fors_zeropoint_impl.h>
00033 
00034 #include <fors_extract.h>
00035 #include <fors_identify.h>
00036 #include <fors_tools.h>
00037 #include <fors_star.h>
00038 #include <fors_std_star.h>
00039 #include <fors_image.h>
00040 #include <fors_data.h>
00041 #include <fors_setting.h>
00042 #include <fors_dfs.h>
00043 #include <fors_qc.h>
00044 #include <fors_pfits.h>
00045 #include <fors_utils.h>
00046 
00047 #include <cpl.h>
00048 
00049 #include <math.h>
00050 
00057 const char *const fors_zeropoint_name = "fors_zeropoint";
00058 const char *const fors_zeropoint_description_short = "Compute zeropoint";
00059 const char *const fors_zeropoint_author = "Jonas M. Larsen";
00060 const char *const fors_zeropoint_email = PACKAGE_BUGREPORT;
00061 const char *const fors_zeropoint_description = 
00062 "Input files:\n"
00063 "  DO category:               Type:       Explanation:              Number:\n"
00064 "  STANDARD_IMG               FITS image  Phot. standard field        1\n"
00065 "  MASTER_BIAS                FITS image  Master bias                 1\n"
00066 "  MASTER_SKY_FLAT_IMG        FITS image  Master sky flatfield        1\n"
00067 "  FLX_STD_IMG                FITS table  Standard star catalog       1+\n"
00068 "  PHOT_TABLE                 FITS table  Filter ext. coeff, color    1\n"
00069 "\n"
00070 "Output files:\n"
00071 "  DO category:               Data type:  Explanation:\n"
00072 "  SOURCES_STD_IMG            FITS image  Unfiltered SExtractor output\n"
00073 "  ALIGNED_PHOT               FITS table\n"
00074 "  PHOT_BACKGROUND_STD_IMG    FITS image  Reduced science image background\n"
00075 "  STANDARD_REDUCED_IMG       FITS image  Reduced std image\n";
00076 
00077 static double
00078 get_zeropoint(fors_star_list *stars, 
00079               double cutoffE,
00080               double cutoffk,
00081           double dext_coeff,
00082           double dcolor_term,
00083           double avg_airmass,
00084               double *dzeropoint,
00085               int *n);
00086 
00091 void fors_zeropoint_define_parameters(cpl_parameterlist *parameters)
00092 {
00093     const char *context = cpl_sprintf("fors.%s", fors_zeropoint_name);
00094     
00095     fors_extract_define_parameters(parameters, context);
00096 
00097     fors_identify_define_parameters(parameters, context);
00098     
00099     const char *name, *full_name;
00100     cpl_parameter *p;
00101 
00102     name = "magcutE";
00103     full_name = cpl_sprintf("%s.%s", context, name);
00104     p = cpl_parameter_new_value(full_name,
00105                                 CPL_TYPE_DOUBLE,
00106                                 "Zeropoint absolute cutoff (magnitude)",
00107                                 context,
00108                                 1.0);
00109     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
00110     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00111     cpl_parameterlist_append(parameters, p);
00112     cpl_free((void *)full_name); full_name = NULL;
00113 
00114     name = "magcutk";
00115     full_name = cpl_sprintf("%s.%s", context, name);
00116     p = cpl_parameter_new_value(full_name,
00117                                 CPL_TYPE_DOUBLE,
00118                                 "Zeropoint kappa rejection parameter",
00119                                 context,
00120                                 5.0);
00121     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, name);
00122     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00123     cpl_parameterlist_append(parameters, p);
00124     cpl_free((void *)full_name); full_name = NULL;
00125 
00126 
00127     cpl_free((void *)context);
00128 
00129     return;
00130 }
00131 
00132 #undef cleanup
00133 #define cleanup \
00134 do { \
00135     cpl_frameset_delete(std_frame); \
00136     cpl_frameset_delete(master_bias_frame); \
00137     cpl_frameset_delete(master_flat_frame); \
00138     cpl_frameset_delete(std_cat); \
00139     cpl_frameset_delete(phot_table); \
00140     fors_image_delete(&std); \
00141     fors_image_delete_const(&master_bias); \
00142     fors_image_delete(&master_flat); \
00143     cpl_table_delete(aligned_phot); \
00144     cpl_image_delete(background); \
00145     cpl_table_delete(sources); \
00146     fors_extract_method_delete(&em); \
00147     fors_identify_method_delete(&im); \
00148     fors_std_star_list_delete(&cat, fors_std_star_delete); \
00149     fors_star_list_delete(&stars, fors_star_delete); \
00150     cpl_free((void *)context); \
00151     fors_setting_delete(&setting); \
00152     cpl_propertylist_delete(qc); \
00153     cpl_propertylist_delete(product_header); \
00154 } while (0)
00155 
00164 void fors_zeropoint(cpl_frameset *frames, const cpl_parameterlist *parameters)
00165 {
00166     /* Raw */
00167     cpl_frameset *std_frame      = NULL;
00168     fors_image *std              = NULL;
00169 
00170     /* Calibration */
00171     cpl_frameset *master_bias_frame = NULL;
00172     const fors_image *master_bias   = NULL; 
00173 
00174     cpl_frameset *master_flat_frame = NULL;
00175     fors_image *master_flat         = NULL; 
00176 
00177     cpl_frameset *std_cat           = NULL;
00178     fors_std_star_list *cat         = NULL;
00179 
00180     cpl_frameset *phot_table        = NULL;
00181     double color_term, dcolor_term;
00182     double ext_coeff, dext_coeff;
00183     double expected_zeropoint, dexpected_zeropoint;
00184 
00185     /* Products */
00186     int nzeropoint = -1;             /* Suppress warning */
00187     cpl_table *aligned_phot = NULL;
00188     cpl_propertylist *qc = cpl_propertylist_new();
00189     double zeropoint, dzeropoint;
00190     fors_extract_sky_stats sky_stats;
00191     cpl_image *background = NULL;
00192     cpl_table *sources    = NULL;
00193     cpl_propertylist *product_header = cpl_propertylist_new();
00194     
00195     /* Parameters */
00196     extract_method  *em = NULL;
00197     identify_method *im = NULL;
00198     double cutoffE, cutoffk;
00199 
00200     /* Other */
00201     const char *context   = cpl_sprintf("fors.%s", fors_zeropoint_name);
00202     fors_star_list *stars = NULL;
00203     fors_setting *setting = NULL;
00204     double avg_airmass;
00205 
00206     /* Get parameters */
00207     em = fors_extract_method_new(parameters, context);
00208     assure( !cpl_error_get_code(), return, 
00209             "Could not get extraction parameters" );
00210     
00211     im = fors_identify_method_new(parameters, context);
00212     assure( !cpl_error_get_code(), return, 
00213             "Could not get identification parameters" );
00214 
00215 
00216     cpl_msg_indent_more();
00217     const char *name = cpl_sprintf("%s.%s", context, "magcutE");
00218     cutoffE = dfs_get_parameter_double_const(parameters, 
00219                                             name);
00220     cpl_free((void *)name); name = NULL;
00221     cpl_msg_indent_less();
00222     assure( !cpl_error_get_code(), return, NULL );
00223         
00224     cpl_msg_indent_more();
00225     name = cpl_sprintf("%s.%s", context, "magcutk");
00226     cutoffk = dfs_get_parameter_double_const(parameters, 
00227                                              name);
00228     cpl_free((void *)name); name = NULL;
00229     cpl_msg_indent_less();
00230     assure( !cpl_error_get_code(), return, NULL );
00231         
00232     /* Find raw */
00233     std_frame = fors_frameset_extract(frames, STANDARD_IMG);
00234     assure( cpl_frameset_get_size(std_frame) == 1, return, 
00235             "Exactly 1 %s required. %d found", 
00236             STANDARD_IMG, cpl_frameset_get_size(std_frame) );
00237 
00238     /* Find calibration */
00239     master_bias_frame = fors_frameset_extract(frames, MASTER_BIAS);
00240     assure( cpl_frameset_get_size(master_bias_frame) == 1, return, 
00241             "One %s required. %d found", 
00242             MASTER_BIAS, cpl_frameset_get_size(master_bias_frame) );
00243 
00244     master_flat_frame = fors_frameset_extract(frames, MASTER_SKY_FLAT_IMG);
00245     assure( cpl_frameset_get_size(master_flat_frame) == 1, return, 
00246             "One %s required. %d found", 
00247             MASTER_SKY_FLAT_IMG, cpl_frameset_get_size(master_flat_frame) );
00248 
00249     std_cat = fors_frameset_extract(frames, FLX_STD_IMG);
00250     assure( cpl_frameset_get_size(std_cat) >= 1, return, 
00251             "One or more %s required. %d found",
00252             FLX_STD_IMG, cpl_frameset_get_size(std_cat));
00253 
00254     phot_table = fors_frameset_extract(frames, PHOT_TABLE);
00255     assure( cpl_frameset_get_size(phot_table) == 1, return, 
00256             "One %s required. %d found",
00257             PHOT_TABLE, cpl_frameset_get_size(phot_table));
00258 
00259     /* Done finding frames */
00260 
00261     /* Get setting */
00262     setting = fors_setting_new(cpl_frameset_get_first(std_frame));
00263     assure( !cpl_error_get_code(), return, "Could not get instrument setting" );
00264     
00265     master_bias = fors_image_load(cpl_frameset_get_first(master_bias_frame), 
00266                                   NULL, setting, NULL);
00267     assure( !cpl_error_get_code(), return, 
00268             "Could not load master bias");
00269 
00270     /* Load raw frames, subtract bias */
00271     std = fors_image_load(cpl_frameset_get_first(std_frame), 
00272                           master_bias, setting, NULL);
00273     assure( !cpl_error_get_code(), return, "Could not load standard image");
00274     fors_image_delete_const(&master_bias);
00275 
00276     /* Load master flat */
00277     master_flat = fors_image_load(cpl_frameset_get_first(master_flat_frame), 
00278                                   NULL, setting, NULL);
00279     assure( !cpl_error_get_code(), return, "Could not load master flat");
00280     
00281     /* Divide by flat */
00282     fors_image_divide_scalar(master_flat,
00283                              fors_image_get_median(master_flat, NULL), -1.0);
00284     fors_image_divide(std, master_flat);
00285     assure( !cpl_error_get_code(), return, "Could not divide by master flat");
00286     fors_image_delete(&master_flat);
00287 
00288     /* Extract sources */
00289     stars = fors_extract(std, setting, em, 
00290              &sky_stats, &background, &sources);
00291     assure( !cpl_error_get_code(), return, "Could not extract objects");
00292 
00293     if (setting->filter_name != NULL) {
00294         /* Load filter coefficients */
00295         fors_phot_table_load(cpl_frameset_get_first(phot_table), setting,
00296                              &color_term, &dcolor_term,
00297                  &ext_coeff, &dext_coeff,
00298                  &expected_zeropoint, &dexpected_zeropoint);
00299         assure( !cpl_error_get_code(), return, 
00300                 "Could not load photometry table" );
00301 
00302         /* Load std catalogue(s) */
00303         cat = fors_std_cat_load(std_cat,
00304                                 cpl_frameset_get_first(std_frame),
00305                                 setting,
00306                                 color_term, dcolor_term);
00307 
00308         assure( !cpl_error_get_code(), return, 
00309                 "Failed to create catalogue patterns" );
00310     
00311 
00312         /* Identify */
00313         fors_identify(stars, cat, im);
00314         assure( !cpl_error_get_code(), return, "Failed to identify sources");
00315     
00316         /* Correct for atmospheric extinction */
00317         avg_airmass = fors_star_ext_corr(stars, setting, ext_coeff, dext_coeff,
00318                                          cpl_frameset_get_first(std_frame));
00319         assure( !cpl_error_get_code(), return, 
00320                 "Extinction correction failed");
00321 
00322         /* Get zeropoint */
00323         zeropoint = get_zeropoint(stars, cutoffE, cutoffk, 
00324                       dext_coeff, dcolor_term,
00325                       avg_airmass,
00326                       &dzeropoint, &nzeropoint);
00327         assure( !cpl_error_get_code(), return, 
00328                 "Failed to get zeropoint" );
00329     }
00330     else {
00331        cpl_msg_warning(cpl_func, "Zeropoint computation is not supported "
00332                        "for non-standard filters");
00333        color_term = 0.0;
00334        dcolor_term = 9999.0;
00335        ext_coeff = 0.0;
00336        dext_coeff = 9999.0;
00337        expected_zeropoint = 0.0;
00338        dexpected_zeropoint = 9999.0;
00339        zeropoint = 0.0;
00340        dzeropoint = 0.0;
00341        nzeropoint = 0;
00342     }
00343 
00344     /* QC */
00345     fors_qc_start_group(qc, fors_qc_dic_version, setting->instrument);
00346     
00347     fors_qc_write_group_heading(cpl_frameset_get_first(std_frame),
00348                                 ALIGNED_PHOT,
00349                                 setting->instrument);
00350     assure( !cpl_error_get_code(), return, "Could not write %s QC parameters", 
00351             ALIGNED_PHOT);
00352 
00353 
00354     fors_qc_write_qc_double(qc,
00355                             zeropoint,
00356                             "QC.ZPOINT",
00357                             "mag",
00358                             "Frame zeropoint",
00359                             setting->instrument);
00360     fors_qc_write_qc_double(qc,
00361                             dzeropoint,
00362                             "QC.ZPOINTRMS",
00363                             "mag",
00364                             "Uncertainty of frame zeropoint",
00365                             setting->instrument);
00366     fors_qc_write_qc_int(qc,
00367                          nzeropoint,
00368                          "QC.ZPOINT.NSTARS",
00369                          NULL,
00370                          "Number of stars used for zeropoint computation",
00371                          setting->instrument);
00372     
00373     double derived_ext_coeff, derived_ext_coeff_err;
00374 
00375     if (setting->filter_name != NULL) {
00376         cpl_msg_info(cpl_func,
00377              "Computing extinction "
00378              "(assuming zeropoint = %.3f +- %.3f mag)",
00379                      expected_zeropoint,
00380              dexpected_zeropoint);
00381 
00382         cpl_msg_indent_more();
00383 
00384         if (nzeropoint > 0) {
00385             derived_ext_coeff = ext_coeff +
00386                 (expected_zeropoint - zeropoint) / avg_airmass;
00387         /* Things are very correlated here
00388            (e.g. ext_coeff was used to compute zeropoint).
00389        
00390            The final error on the ext.coeff. depends only
00391            on the reference and computed zeropoints' errors. 
00392            The airmass is assumed errorless.
00393 
00394            We assume the 2 zeropoints' errors are not correlated and
00395            add in quadrature.
00396 
00397            The derived extinction's error does not suffer from the part
00398            of dzeropoint which was due to the error of the assumed 
00399                extinction.
00400         */
00401         derived_ext_coeff_err = sqrt(
00402             dexpected_zeropoint*dexpected_zeropoint +
00403             dzeropoint*dzeropoint) / avg_airmass - dext_coeff*dext_coeff;
00404     
00405             cpl_msg_info(cpl_func, "Atmospheric extinction = "
00406                          "%f +- %f mag/airmass", derived_ext_coeff,
00407                  derived_ext_coeff_err);
00408         }
00409         else {
00410             cpl_msg_warning(cpl_func, "Too few stars available, "
00411                             "setting extinction to zero");
00412             derived_ext_coeff = 0;
00413             derived_ext_coeff_err = 9999;
00414         }
00415     }
00416     else {
00417         derived_ext_coeff = 0;
00418         derived_ext_coeff_err = 9999;
00419     }
00420 
00421     fors_qc_write_qc_double(qc,
00422                             derived_ext_coeff,
00423                             "QC.EXTCOEFF",
00424                             "mag/airmass",
00425                             "Atmospheric extinction (mag/airmass)",
00426                             setting->instrument);
00427 
00428     fors_qc_write_qc_double(qc,
00429                             derived_ext_coeff_err,
00430                             "QC.EXTCOEFFERR",
00431                             "mag/airmass",
00432                             "Error on atmospheric extinction (mag/airmass)",
00433                             setting->instrument);
00434     cpl_msg_indent_less();
00435 
00436     fors_qc_end_group();
00437 
00438 
00439     /* Convert to CPL table */
00440     aligned_phot = fors_create_sources_table(stars);
00441     assure( !cpl_error_get_code(), return,
00442             "Failed to create aligned photometry table");
00443 
00444     /* Save products */
00445     fors_dfs_save_table(frames, sources, SOURCES_STD,
00446                         NULL, parameters, fors_zeropoint_name, 
00447                         cpl_frameset_get_first(std_frame));
00448     assure( !cpl_error_get_code(), return, "Saving %s failed",
00449             SOURCES_STD);
00450 
00451     /* Add keywords necessary for fors_photometry 
00452        (just reuse/overload to the qc propertylist variable)
00453     */
00454     cpl_propertylist_update_double(qc, FORS_PFITS_EXPOSURE_TIME, 
00455                    setting->exposure_time);
00456     cpl_propertylist_update_double(qc, "AIRMASS", avg_airmass);
00457 
00458     fors_dfs_save_table(frames, aligned_phot, ALIGNED_PHOT,
00459                         qc, parameters, fors_zeropoint_name, 
00460                         cpl_frameset_get_first(std_frame));
00461     assure( !cpl_error_get_code(), return, "Saving %s failed",
00462             ALIGNED_PHOT);
00463     
00464     cpl_propertylist_delete(qc); qc = NULL;
00465 
00466 
00467     fors_dfs_add_wcs(product_header,
00468              cpl_frameset_get_first(std_frame), setting);
00469     fors_dfs_add_exptime(product_header, 
00470                          cpl_frameset_get_first(std_frame), 0.);
00471     
00472     fors_dfs_save_image(frames, std, STANDARD_REDUCED_IMG,
00473             product_header, parameters, fors_zeropoint_name, 
00474             cpl_frameset_get_first(std_frame));
00475     assure( !cpl_error_get_code(), return, "Saving %s failed",
00476         STANDARD_REDUCED_IMG);
00477     
00478     dfs_save_image(frames, background, PHOT_BACKGROUND_STD_IMG,
00479            product_header, parameters, fors_zeropoint_name, 
00480            setting->version);
00481     assure( !cpl_error_get_code(), return, "Saving %s failed",
00482         PHOT_BACKGROUND_STD_IMG);
00483     
00484     cpl_image_delete(background); background = NULL;
00485     
00486     if (setting->filter_name == NULL) {
00487 
00488         /* No debug image can be created */
00489 
00490         cleanup;
00491         return;
00492     }
00493 
00494     /* Create debugging image
00495 
00496        Legend:
00497 
00498         ------  (horiz. line):   detected source
00499 
00500           |
00501           |    (vert. line):   catalog position
00502           |
00503 
00504           _
00505          / \     (circle) :   identified
00506          \_/
00507     */
00508 
00509     double color = fors_image_get_min(std);
00510     {
00511         fors_star *s;
00512         for (s = fors_star_list_first(stars);
00513              s != NULL; 
00514              s = fors_star_list_next(stars)) {
00515             fors_image_draw(std, 0,
00516                             s->pixel->x,
00517                             s->pixel->y,
00518                             10, color);
00519 
00520             if (s->id != NULL) {
00521                 fors_image_draw(std, 2,
00522                                 s->pixel->x,
00523                                 s->pixel->y,
00524                                 10, color);
00525             }
00526         }
00527     }
00528     {        
00529         fors_std_star *s;
00530         for (s = fors_std_star_list_first(cat);
00531              s != NULL; 
00532              s = fors_std_star_list_next(cat)) {
00533             fors_image_draw(std, 1,
00534                             s->pixel->x,
00535                             s->pixel->y,
00536                             10, color);
00537         }
00538         fors_dfs_save_image(frames, std, "DEBUG",
00539                             product_header, parameters, fors_zeropoint_name, 
00540                             cpl_frameset_get_first(std_frame));
00541         assure( !cpl_error_get_code(), return, "Saving %s failed",
00542                 "DEBUG");
00543     }
00544 
00545     cleanup;
00546     return;
00547 }
00548 
00557 static bool
00558 zeropoint_inside(const fors_star *s,
00559                  void *data) 
00560 {
00561     struct {
00562         double hi, lo;   /* magnitude */
00563         double z, kappa; /* avg zeropoint, kappa */
00564     } *cuts = data;
00565     
00566     double z  = fors_star_get_zeropoint(s, NULL);
00567     double dz = fors_star_get_zeropoint_err(s, NULL);
00568 
00569     return
00570         (cuts->lo                 <= z && z <= cuts->hi) ||
00571         (cuts->z - cuts->kappa*dz <= z && z <= cuts->z + cuts->kappa*dz);
00572 }
00573 
00574 #undef cleanup
00575 #define cleanup \
00576 do { \
00577     fors_star_list_delete(&subset, fors_star_delete); \
00578     fors_star_list_delete(&identified, fors_star_delete); \
00579 } while(0)
00580 
00592 static double
00593 get_zeropoint(fors_star_list *stars, 
00594               double cutoffE,
00595               double cutoffk,
00596           double dext_coeff,
00597           double dcolor_term,
00598           double avg_airmass,
00599               double *dzeropoint,
00600               int *n)
00601 {
00602     fors_star_list *subset = NULL;
00603     fors_star_list *identified = 
00604     fors_star_list_extract(stars, fors_star_duplicate,
00605                    fors_star_is_identified, NULL);
00606 
00607     assure( stars != NULL, return 0, NULL );
00608     assure( dzeropoint != NULL, return 0, NULL );
00609     assure( n != NULL, return 0, NULL );
00610 
00611     if ( fors_star_list_size(identified) == 0 ) {
00612         cpl_msg_warning(cpl_func, 
00613                         "No identified stars for zeropoint computation");
00614         *n = 0;
00615         *dzeropoint = 0;
00616         cleanup;
00617         return 0;
00618     }
00619     
00620     cpl_msg_info(cpl_func, "Computing zeropoint (assuming extinction)");
00621     cpl_msg_indent_more();
00622 
00623     double zeropoint;
00624     double red_chisq = -1.0;
00625 
00626     /* This method does not take into account that the error bars are
00627        correlated, and therefore computes an unrealistically low
00628        dzeropoint */
00629     zeropoint = fors_star_list_mean_optimal(identified,
00630                         fors_star_get_zeropoint, NULL,
00631                         fors_star_get_zeropoint_err, NULL,
00632                         dzeropoint,
00633                         fors_star_list_size(identified) >= 2 ? &red_chisq : NULL);
00634 
00635     cpl_msg_info(cpl_func, "Optimal zeropoint (no rejection, %d stars) = %f mag",
00636                  fors_star_list_size(identified), zeropoint);  
00637     
00638     /* Reject stars that are absolute (0.3 mag) outliers and
00639        kappa sigma outliers. For robustness (against error in the initial
00640        zeropoint estimates) apply the absolute cut in two steps
00641        and update the estimated zeropoint after the first step.
00642     */
00643     struct {
00644         double hi, lo;   /* magnitude */
00645         double z, kappa; /* avg zeropoint, kappa */
00646     } cuts;
00647     cuts.hi = zeropoint + 5*cutoffE;
00648     cuts.lo = zeropoint - 5*cutoffE;
00649     cuts.z = zeropoint;
00650     cuts.kappa = cutoffk;
00651     
00652     subset = fors_star_list_extract(identified, fors_star_duplicate,
00653                                     zeropoint_inside, &cuts);
00654 
00655 
00656     if ( fors_star_list_size(subset) == 0 ) {
00657         cpl_msg_warning(cpl_func, 
00658                         "All stars rejected (%f mag). Cannot "
00659             "compute zeropoint", 5*cutoffE);
00660         *n = 0;
00661         *dzeropoint = 0;
00662         cleanup;
00663         return 0;
00664     }
00665 
00666     zeropoint = fors_star_list_mean_optimal(subset,
00667                         fors_star_get_zeropoint, NULL,
00668                         fors_star_get_zeropoint_err, NULL,
00669                         dzeropoint,
00670                         fors_star_list_size(subset) >= 2 ? &red_chisq : NULL);
00671 
00672     cpl_msg_debug(cpl_func, "Optimal zeropoint (%.2f mag, %.2f sigma rejection) = %f mag",
00673           5*cutoffE, cutoffk, zeropoint);  
00674     
00675     cuts.hi = zeropoint + cutoffE;
00676     cuts.lo = zeropoint - cutoffE;
00677     cuts.z = zeropoint;
00678     cuts.kappa = cutoffk;
00679     
00680     {
00681     fors_star_list *tmp = fors_star_list_duplicate(subset, fors_star_duplicate);
00682     fors_star_list_delete(&subset, fors_star_delete);
00683 
00684     subset = fors_star_list_extract(tmp, fors_star_duplicate,
00685                     zeropoint_inside, &cuts);
00686 
00687     if ( fors_star_list_size(subset) == 0 ) {
00688         cpl_msg_warning(cpl_func, 
00689                 "All stars rejected (%f mag, %f sigma). Cannot "
00690                 "compute zeropoint", cutoffE, cutoffk);
00691         *n = 0;
00692         *dzeropoint = 0;
00693         cleanup;
00694         return 0;
00695     }
00696 
00697     fors_star_list_delete(&tmp, fors_star_delete);
00698     }
00699 
00700     zeropoint = fors_star_list_mean_optimal(subset,
00701                         fors_star_get_zeropoint, NULL,
00702                         fors_star_get_zeropoint_err, NULL,
00703                         dzeropoint,
00704                         fors_star_list_size(subset) >= 2 ? &red_chisq : NULL);
00705     
00706     cpl_msg_info(cpl_func, "Optimal zeropoint (%.2f mag, %.2f sigma rejection) = %f mag",
00707          cutoffE, cutoffk, zeropoint);  
00708     
00709     *n = fors_star_list_size(subset);
00710     {
00711         int outliers = 
00712             fors_star_list_size(identified) - fors_star_list_size(subset);
00713         cpl_msg_info(cpl_func, "%d outlier%s rejected, %d non-outliers",
00714                      outliers, outliers == 1 ? "" : "s",
00715                      *n);
00716     }
00717 
00718     if ( *n == 0 ) {
00719         cpl_msg_warning(cpl_func, 
00720                         "All stars were rejected during zeropoint computation" );
00721         *dzeropoint = 0;
00722         cleanup;
00723         return 0;
00724     }
00725 
00726 
00727 
00728 
00729 
00730 
00731 
00732     /*
00733       Build zeropoint covariance matrix.
00734       We have already the variances from fors_star_get_zeropoint_err().
00735       Non-diagonal terms are
00736          Cij = airmass^2 * Variance(ext.coeff) + color_i * color_j * Variance(color.coeff)
00737 
00738       It was considered and tried to subtract the term
00739                airmass^2 * Variance(ext.coeff)
00740       from every Cij. This has no effect on the relative weights, only the weight's overall
00741       normalization. Since we use the normalization for computing the zeropoint error this
00742       term is kept.
00743 
00744     */
00745     cpl_matrix *covariance = cpl_matrix_new(*n,
00746                         *n);
00747 
00748     /* Duplicate the list to allow simultaneous iterations */
00749     fors_star_list *ident_dup = fors_star_list_duplicate(subset, fors_star_duplicate);
00750     {
00751       
00752       fors_star *s, *t;
00753       int i, j;
00754       for (s = fors_star_list_first(subset), i = 0;
00755        s != NULL;
00756        s = fors_star_list_next(subset), i++) {
00757 
00758     for (t = fors_star_list_first(ident_dup), j = 0;
00759          t != NULL;
00760          t = fors_star_list_next(ident_dup), j++) {
00761       
00762       double cij;
00763 
00764       if (fors_star_equal(s, t)) {
00765           cij = fors_star_get_zeropoint_err(s, NULL)*fors_star_get_zeropoint_err(s, NULL);
00766           /*  -avg_airmass*avg_airmass*dext_coeff*dext_coeff */
00767       }
00768       else {
00769           cij = s->id->color * t->id->color * dcolor_term*dcolor_term
00770           + avg_airmass*avg_airmass*dext_coeff*dext_coeff;
00771       }
00772       
00773       cpl_matrix_set(covariance, i, j, cij);
00774     }
00775       }
00776     }
00777     /* cpl_matrix_dump(covariance, stdout); */
00778     /* cpl_matrix_dump(cpl_matrix_invert_create(covariance), stdout); */
00779 
00780     /*
00781       Compute optimal weights, w, as
00782       
00783       w = C^-1 * const
00784 
00785       C    : nxn covariance matrix
00786       const: nx1 constant vector with all elements equal to 1
00787      */
00788 
00789     cpl_matrix *covariance_inverse = cpl_matrix_invert_create(covariance);
00790 
00791     assure( !cpl_error_get_code(), return 0,
00792         "Could not invert zeropoints covariance matrix");
00793 
00794     /*  cpl_matrix_dump(cpl_matrix_product_create(covariance_inverse, covariance), stdout); */
00795 
00796     /* fprintf(stderr, "is_identity = %d\n", cpl_matrix_is_identity(cpl_matrix_product_create(covariance_inverse, covariance),1e-10)); */
00797 
00798     cpl_matrix_delete(covariance); covariance = NULL;
00799 
00800     cpl_matrix *const_vector = cpl_matrix_new(*n, 1);
00801     cpl_matrix_fill(const_vector, 1.0);
00802     
00803     cpl_matrix *weights = cpl_matrix_product_create(covariance_inverse, const_vector);
00804 
00805     cpl_matrix_delete(const_vector); const_vector = NULL;
00806 
00807     /* cpl_matrix_dump(weights, stdout); */
00808     
00809     
00810     {
00811     double wz = 0;
00812     double w = 0;
00813     
00814     int i;
00815     fors_star *s;
00816     for (i = 0, s = fors_star_list_first(subset);
00817          s != NULL;
00818          s = fors_star_list_next(subset), i++) {
00819         
00820         double weight = cpl_matrix_get(weights, i, 0);
00821 
00822         cpl_msg_debug(cpl_func, "Weight_%d = %f", i, weight);
00823         
00824         wz += weight * fors_star_get_zeropoint(s, NULL);
00825         w += weight;
00826 
00827         /* Loop through original list to record the weight of this star */
00828         {
00829         fors_star *t;
00830         
00831         for (t = fors_star_list_first(stars);
00832              t != NULL;
00833              t = fors_star_list_next(stars)) {
00834             
00835             if (fors_star_equal(s, t)) {
00836             t->weight = weight;
00837             }
00838         }
00839         }
00840     }
00841 
00842         cpl_matrix_delete(weights); weights = NULL;
00843 
00844     cpl_msg_debug(cpl_func, "Sum of weights = %f", w);
00845 
00846     /* 
00847        C is positive definite (because all eigenvalues are positive). Therefore
00848        C^-1 is also positive definite, a property of positive definite matrices.
00849        
00850        Positive definite matrices always have the property that
00851                  z* C z > 0     for any non-zero (complex) vector z
00852 
00853        The sum of the weights is just
00854                 const.* w = const.* C^-1 const.
00855            where const. is our constant vector filled with 1. Therefore the sum of
00856            the weights should always be positive. But make the paranoia check anyway:
00857      */
00858 
00859     assure( w != 0, return 0, "Sum of optimal weights is zero!" );
00860     assure( sqrt(w) != 0, return 0, "Square root of sum of optimal weights is zero!" );
00861     
00862     zeropoint = wz / w;
00863     *dzeropoint = 1 / sqrt(w);
00864     }
00865 
00866     /* Previous code: weighted average, uncorrelated errors
00867 
00868     zeropoint = fors_star_list_mean_optimal(
00869         subset, 
00870         fors_star_get_zeropoint, NULL,
00871         fors_star_get_zeropoint_err, NULL,
00872         dzeropoint,
00873         *n >= 2 ? &red_chisq : NULL);
00874     */
00875     
00876     cpl_msg_info(cpl_func, "Optimal zeropoint = %f +- %f mag", 
00877          zeropoint, *dzeropoint);
00878 
00879     if (*n >= 2) {
00880     fors_star *s, *t;
00881     int i, j;
00882 
00883     red_chisq = 0;
00884 
00885     for (s = fors_star_list_first(subset), i = 0;
00886          s != NULL;
00887          s = fors_star_list_next(subset), i++) {
00888         
00889         for (t = fors_star_list_first(ident_dup), j = 0;
00890          t != NULL;
00891          t = fors_star_list_next(ident_dup), j++) {
00892         
00893         red_chisq += 
00894             (fors_star_get_zeropoint(s, NULL) - zeropoint) *
00895             (fors_star_get_zeropoint(t, NULL) - zeropoint) *
00896             cpl_matrix_get(covariance_inverse, i, j);
00897         }
00898     }
00899     red_chisq /= (*n - 1);
00900     }
00901 
00902     cpl_matrix_delete(covariance_inverse); covariance_inverse = NULL;
00903 
00904     fors_star_list_delete(&ident_dup, fors_star_delete);
00905     
00906     cpl_msg_info(cpl_func, "Reduced chi square = %f", red_chisq);
00907     
00908     cpl_msg_indent_less();
00909 
00910     cleanup;
00911     return zeropoint;
00912 }
00913 
00914 

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