fors_tools.c

00001 /* $Id: fors_tools.c,v 1.15 2008/02/27 15:10:19 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/02/27 15:10:19 $
00024  * $Revision: 1.15 $
00025  * $Name:  $
00026  */
00027 
00028 #ifdef HAVE_CONFIG_H
00029 #include <config.h>
00030 #endif
00031 
00032 #include <fors_tools.h>
00033 
00034 #include <fors_pfits.h>
00035 #include <fors_utils.h>
00036 
00037 #include <cpl.h>
00038 #include <math.h>
00039 #include <stdbool.h>
00040 
00041 /*----------------------------------------------------------------------------*/
00045 /*----------------------------------------------------------------------------*/
00046 
00049 #undef cleanup
00050 #define cleanup \
00051 do { \
00052     cpl_propertylist_delete(header); \
00053 } while(0)
00054 
00063 double
00064 fors_star_ext_corr(fors_star_list *stars, 
00065                    const fors_setting *setting,
00066                    double ext_coeff,
00067            double dext_coeff,
00068                    const cpl_frame *raw_frame)
00069 {
00070     cpl_propertylist *header = NULL;
00071     
00072     
00073     cpl_msg_info(cpl_func, "Extinction correction");
00074     
00075     assure( cpl_frame_get_filename(raw_frame) != NULL, return -1, NULL );
00076     
00077     header = cpl_propertylist_load(cpl_frame_get_filename(raw_frame), 0);
00078     assure( !cpl_error_get_code(), return -1,
00079             "Failed to load %s primary header",
00080             cpl_frame_get_filename(raw_frame));
00081 
00082 
00083     double avg_airmass = fors_get_airmass(header);
00084     assure( !cpl_error_get_code(), return -1,
00085         "%s: Could not read airmass",
00086         cpl_frame_get_filename(raw_frame));
00087 
00088     cpl_msg_indent_more();
00089     cpl_msg_info(cpl_func, "Exposure time = %f s", setting->exposure_time);
00090     cpl_msg_info(cpl_func, "Gain          = %f ADU/e-", setting->average_gain);
00091     cpl_msg_info(cpl_func, "Ext. coeff.   = %f +- %f mag/airmass", 
00092          ext_coeff, dext_coeff);
00093     cpl_msg_info(cpl_func, "Avg. airmass  = %f airmass", avg_airmass);  
00094     /* The quantity and the unit are both 'airmass' */
00095 
00096     cpl_msg_indent_less();
00097     
00098     {
00099         fors_star *star;
00100         
00101         for (star = fors_star_list_first(stars);
00102              star != NULL;
00103              star = fors_star_list_next(stars)) {
00104             star->magnitude_corr = star->magnitude 
00105                 + 2.5*log(setting->average_gain)/log(10)
00106                 + 2.5*log(setting->exposure_time)/log(10)
00107                 - ext_coeff * avg_airmass;
00108 
00109         /* Propagate error from ext.coeff.
00110            gain, exptime and airmass
00111            have zero error */
00112         star->dmagnitude_corr = sqrt(star->dmagnitude * star->dmagnitude
00113                      + dext_coeff*dext_coeff * avg_airmass*avg_airmass);
00114     }
00115     }
00116 
00117     cleanup;
00118     return avg_airmass;
00119 }
00120 
00128 cpl_table *
00129 fors_create_sources_table(fors_star_list *sources)
00130 {
00131     cpl_table *t = NULL;
00132     
00133     t = cpl_table_new(fors_star_list_size(sources));
00134     cpl_table_new_column(t, "X", CPL_TYPE_DOUBLE);
00135     cpl_table_new_column(t, "Y", CPL_TYPE_DOUBLE);
00136     cpl_table_new_column(t, "FWHM", CPL_TYPE_DOUBLE);
00137     cpl_table_new_column(t, "A", CPL_TYPE_DOUBLE);
00138     cpl_table_new_column(t, "B", CPL_TYPE_DOUBLE);
00139     cpl_table_new_column(t, "THETA", CPL_TYPE_DOUBLE);
00140     cpl_table_new_column(t, "ELL", CPL_TYPE_DOUBLE);
00141     cpl_table_new_column(t, "INSTR_MAG", CPL_TYPE_DOUBLE);
00142     cpl_table_new_column(t, "DINSTR_MAG", CPL_TYPE_DOUBLE);
00143     cpl_table_new_column(t, "INSTR_CMAG", CPL_TYPE_DOUBLE);
00144     cpl_table_new_column(t, "DINSTR_CMAG", CPL_TYPE_DOUBLE);
00145     cpl_table_new_column(t, "CLASS_STAR", CPL_TYPE_DOUBLE);
00146 
00147     cpl_table_new_column(t, "OBJECT", CPL_TYPE_STRING);
00148     cpl_table_new_column(t, "RA", CPL_TYPE_DOUBLE);
00149     cpl_table_new_column(t, "DEC", CPL_TYPE_DOUBLE);
00150     cpl_table_new_column(t, "MAG", CPL_TYPE_DOUBLE);
00151     cpl_table_new_column(t, "DMAG", CPL_TYPE_DOUBLE);
00152     cpl_table_new_column(t, "CAT_MAG", CPL_TYPE_DOUBLE);
00153     cpl_table_new_column(t, "DCAT_MAG", CPL_TYPE_DOUBLE);
00154     cpl_table_new_column(t, "COLOR", CPL_TYPE_DOUBLE);    
00155     cpl_table_new_column(t, "USE_CAT", CPL_TYPE_INT);
00156     /* Shift in x and y between initial guess FITS header WCS position
00157        and measured position */
00158     cpl_table_new_column(t, "SHIFT_X", CPL_TYPE_DOUBLE); 
00159     cpl_table_new_column(t, "SHIFT_Y", CPL_TYPE_DOUBLE);
00160     cpl_table_new_column(t, "ZEROPOINT", CPL_TYPE_DOUBLE);
00161     cpl_table_new_column(t, "DZEROPOINT", CPL_TYPE_DOUBLE);
00162     cpl_table_new_column(t, "WEIGHT", CPL_TYPE_DOUBLE);
00163 
00164     {
00165         fors_star *s;
00166         int i;
00167         for (s = fors_star_list_first(sources), i = 0;
00168              s != NULL;
00169              s = fors_star_list_next(sources), i++) {
00170 
00171             const fors_std_star *id = s->id;
00172             
00173             cpl_table_set_double(t, "X", i, s->pixel->x);
00174             cpl_table_set_double(t, "Y", i, s->pixel->y);
00175             cpl_table_set_double(t, "FWHM", i, s->fwhm);
00176             cpl_table_set_double(t, "A", i, s->semi_major);
00177             cpl_table_set_double(t, "B", i, s->semi_minor);
00178             cpl_table_set_double(t, "THETA", i, s->orientation);
00179             cpl_table_set_double(t, "ELL", i, fors_star_ellipticity(s, NULL));
00180             cpl_table_set_double(t, "INSTR_MAG", i, s->magnitude);
00181             cpl_table_set_double(t, "DINSTR_MAG", i, s->dmagnitude);
00182             cpl_table_set_double(t, "INSTR_CMAG", i, s->magnitude_corr);
00183             cpl_table_set_double(t, "DINSTR_CMAG", i, s->dmagnitude_corr);
00184             cpl_table_set_double(t, "CLASS_STAR", i, s->stellarity_index);
00185         cpl_table_set_double(t, "WEIGHT", i, s->weight);
00186         
00187         if (id != NULL) {
00188                 cpl_table_set_string(t, "OBJECT" , i, id->name); /* possibly NULL */
00189                 cpl_table_set_double(t, "RA"     , i, id->ra);
00190                 cpl_table_set_double(t, "DEC"    , i, id->dec);
00191                 cpl_table_set_double(t, "MAG"    , i, id->magnitude);
00192                 cpl_table_set_double(t, "DMAG"   , i, id->dmagnitude);
00193                 cpl_table_set_double(t, "CAT_MAG", i, id->cat_magnitude);
00194                 cpl_table_set_double(t, "DCAT_MAG", i, id->dcat_magnitude);
00195                 cpl_table_set_double(t, "COLOR"  , i, id->color);
00196                 cpl_table_set_double(t, "SHIFT_X", i, s->pixel->x - id->pixel->x);
00197                 cpl_table_set_double(t, "SHIFT_Y", i, s->pixel->y - id->pixel->y);
00198                 cpl_table_set_double(t, "ZEROPOINT", i, 
00199                                      fors_star_get_zeropoint(s, NULL));
00200                 cpl_table_set_double(t, "DZEROPOINT", i, 
00201                                      fors_star_get_zeropoint_err(s, NULL));
00202         cpl_table_set_int   (t, "USE_CAT", i, 1); /* fit this magnitude
00203                                  in fors_photometry?
00204                                  (1 = true) */
00205             }
00206             else {
00207                 cpl_table_set_invalid(t, "RA" , i);
00208                 cpl_table_set_invalid(t, "DEC", i);
00209                 cpl_table_set_invalid(t, "MAG", i);
00210                 cpl_table_set_invalid(t, "DMAG", i);
00211                 cpl_table_set_invalid(t, "SHIFT_X", i);
00212                 cpl_table_set_invalid(t, "SHIFT_Y", i);
00213                 cpl_table_set_invalid(t, "ZEROPOINT", i);
00214                 cpl_table_set_invalid(t, "DZEROPOINT", i);
00215             }
00216         }
00217     }
00218 
00219     return t;
00220 }
00221 
00222 #undef cleanup
00223 #define cleanup \
00224 do { \
00225     fors_image_delete(&image); \
00226     fors_image_delete(&image2); \
00227 } while(0)
00228 
00235 double
00236 fors_fixed_pattern_noise(const fors_image *master,
00237                          double convert_ADU,
00238                          double master_noise)
00239 {
00240     double master_fixed_pattern_noise;
00241     fors_image *image = NULL;
00242     fors_image *image2 = NULL;
00243 
00244     assure( master != NULL, return -1, NULL );
00245 
00246     /* Use central 101x101 window 
00247        and 101x101 window shifted (10, 10) from center
00248     */
00249     if (fors_image_get_size_x(master) >= 121 &&
00250         fors_image_get_size_y(master) >= 121) {
00251         
00252         int mid_x = (fors_image_get_size_x(master) + 1) / 2;
00253         int mid_y = (fors_image_get_size_y(master) + 1) / 2;
00254         
00255         image = fors_image_duplicate(master);       
00256         fors_image_crop(image, 
00257                         mid_x - 50, mid_y - 50,
00258                         mid_x + 50, mid_y + 50);
00259         
00260         image2 = fors_image_duplicate(master);       
00261         fors_image_crop(image2, 
00262                         mid_x + 10 - 50, mid_y + 10 - 50,
00263                         mid_x + 10 + 50, mid_y + 10 + 50);
00264 
00265         fors_image_subtract(image, image2);
00266 
00267         master_fixed_pattern_noise = 
00268             fors_image_get_stdev(image, NULL) / sqrt(2);
00269 
00270         /* Convert to ADU */
00271         master_fixed_pattern_noise *= convert_ADU;
00272 
00273         /* Subtract photon noise */
00274         if (master_fixed_pattern_noise >= master_noise) {
00275             
00276             master_fixed_pattern_noise = sqrt(master_fixed_pattern_noise*
00277                                               master_fixed_pattern_noise
00278                                               -
00279                                               master_noise*
00280                                               master_noise);
00281         }
00282         else {
00283             cpl_msg_warning(cpl_func,
00284                             "Zero-shift noise (%f ADU) is greater than "
00285                             "accumulated zero-shift and fixed pattern noise (%f ADU), "
00286                             "setting fixed pattern noise to zero",
00287                             master_noise,
00288                             master_fixed_pattern_noise);
00289             master_fixed_pattern_noise = 0;
00290         }
00291     }
00292     else {
00293         cpl_msg_warning(cpl_func,
00294                         "Master flat too small (%dx%d), "
00295                         "need size 121x121 to compute master flat "
00296                         "fixed pattern noise",
00297                         fors_image_get_size_x(master),
00298                         fors_image_get_size_y(master));
00299         master_fixed_pattern_noise = -1;
00300     }
00301 
00302     cleanup;
00303     return master_fixed_pattern_noise;
00304 }
00305 
00306 
00307 #undef cleanup
00308 #define cleanup \
00309 do { \
00310     fors_image_delete(&image); \
00311     fors_image_delete(&image2); \
00312 } while(0)
00313 
00320 double
00321 fors_fixed_pattern_noise_bias(const fors_image *first_raw,
00322                               const fors_image *second_raw,
00323                               double ron)
00324 {
00325     double bias_fixed_pattern_noise;
00326     fors_image *image = NULL;
00327     fors_image *image2 = NULL;
00328     int nx, ny;
00329 
00330     assure( first_raw != NULL, return -1, NULL );
00331     assure( second_raw != NULL, return -1, NULL );
00332 
00333     /* 
00334      * Extract the largest possible two windows shifted (10, 10) 
00335      */
00336 
00337     nx = fors_image_get_size_x(first_raw);
00338     ny = fors_image_get_size_y(first_raw);
00339 
00340     image = fors_image_duplicate(first_raw);       
00341     fors_image_crop(image, 
00342                     1, 1,
00343                     nx - 10, ny - 10);
00344         
00345     image2 = fors_image_duplicate(second_raw);       
00346     fors_image_crop(image2, 
00347                     11, 11,
00348                     nx, ny);
00349 
00350     fors_image_subtract(image, image2);
00351 
00352     bias_fixed_pattern_noise = fors_image_get_stdev_robust(image, 50, NULL) 
00353                              / sqrt(2);
00354 
00355     /* 
00356      * Subtract ron quadratically
00357      */
00358 
00359     if (bias_fixed_pattern_noise > ron) {
00360             
00361         bias_fixed_pattern_noise = sqrt(bias_fixed_pattern_noise *
00362                                         bias_fixed_pattern_noise
00363                                         -
00364                                         ron * ron);
00365     }
00366     else {
00367         cpl_msg_warning(cpl_func,
00368                         "Zero-shift noise (%f ADU) is greater than "
00369                         "accumulated zero-shift and fixed pattern "
00370                         "noise (%f ADU), "
00371                         "setting fixed pattern noise to zero",
00372                         ron,
00373                         bias_fixed_pattern_noise);
00374         bias_fixed_pattern_noise = 0;
00375     }
00376 
00377     cleanup;
00378     return bias_fixed_pattern_noise;
00379 }
00380 
00381 
00382 #undef cleanup
00383 #define cleanup
00384 
00389 double
00390 fors_get_airmass(const cpl_propertylist *header)
00391 {
00392   double airmass_start, airmass_end;
00393   airmass_start = cpl_propertylist_get_double(header, FORS_PFITS_AIRMASS_START);
00394   assure( !cpl_error_get_code(), return -1, 
00395       "Could not read %s from header", 
00396       FORS_PFITS_AIRMASS_START);
00397   
00398   airmass_end = cpl_propertylist_get_double(header, FORS_PFITS_AIRMASS_END);
00399   assure( !cpl_error_get_code(), return -1,
00400       "Could not read %s from header", 
00401       FORS_PFITS_AIRMASS_END);
00402   
00403   return 0.5 * (airmass_start + airmass_end);
00404 }
00405 

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