fors_pmos_science.c

00001 /* $Id: fors_pmos_science.c,v 1.34 2009/06/09 14:20:27 cizzo Exp $
00002  *
00003  * This file is part of the FORS Data Reduction Pipeline
00004  * Copyright (C) 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/09 14:20:27 $
00024  * $Revision: 1.34 $
00025  * $Name:  $
00026  */
00027 
00028 #ifdef HAVE_CONFIG_H
00029 #include <config.h>
00030 #endif
00031 
00032 #include <math.h>
00033 #include <string.h>
00034 #include <ctype.h>
00035 #include <assert.h>
00036 
00037 #include <cpl.h>
00038 #include <moses.h>
00039 #include <fors_dfs.h>
00040 #include <fors_utils.h>
00041 #include <fors_qc.h>
00042 
00043 static int fors_pmos_science_create(cpl_plugin *);
00044 static int fors_pmos_science_exec(cpl_plugin *);
00045 static int fors_pmos_science_destroy(cpl_plugin *);
00046 static int fors_pmos_science(cpl_parameterlist *, cpl_frameset *);
00047 
00048 static float * fors_check_angles(cpl_frameset *, int, const char *, int *);
00049 static int
00050 fors_find_angle_pos(float * angles, int nangles, float angle);
00051 
00052 static char fors_pmos_science_description[] =
00053 "This recipe is used to reduce scientific spectra using the extraction\n"
00054 "mask and the products created by the recipe fors_mpol_calib. The spectra\n"
00055 "are bias subtracted, flat fielded (if a normalised flat field is specified)\n"
00056 "and remapped eliminating the optical distortions. The wavelength calibration\n"
00057 "can be optionally upgraded using a number of sky lines: if no sky lines\n"
00058 "catalog of wavelengths is specified, an internal one is used instead.\n"
00059 "If the alignment to the sky lines is performed, the input dispersion\n"
00060 "coefficients table is upgraded and saved to disk, and a new CCD wavelengths\n"
00061 "map is created.\n"
00062 "This recipe accepts both FORS1 and FORS2 frames. A grism table (typically\n"
00063 "depending on the instrument mode, and in particular on the grism used)\n"
00064 "may also be specified: this table contains a default recipe parameter\n" 
00065 "setting to control the way spectra are extracted for a specific instrument\n"
00066 "mode, as it is used for automatic run of the pipeline on Paranal and in\n" 
00067 "Garching. If this table is specified, it will modify the default recipe\n" 
00068 "parameter setting, with the exception of those parameters which have been\n" 
00069 "explicitly modifyed on the command line. If a grism table is not specified,\n"
00070 "the input recipe parameters values will always be read from the command\n" 
00071 "line, or from an esorex configuration file if present, or from their\n" 
00072 "generic default values (that are rarely meaningful).\n" 
00073 "Either a scientific or a standard star exposure can be specified in input.\n"
00074 "The acronym SCI on products should be read STD in case of standard stars\n"
00075 "observations.\n\n"
00076 "Input files:\n\n"
00077 "  DO category:               Type:       Explanation:         Required:\n"
00078 "  SCIENCE_PMOS                  Raw         Scientific exposure     Y\n"
00079 "  or STANDARD_PMOS              Raw         Standard star exposure  Y\n"
00080 "  MASTER_BIAS                   Calib       Master bias             Y\n"
00081 "  GRISM_TABLE                   Calib       Grism table             .\n"
00082 "  MASTER_SKYLINECAT             Calib       Sky lines catalog       .\n"
00083 "\n"
00084 "  MASTER_NORM_FLAT_PMOS         Calib       Normalised flat field   .\n"
00085 "  DISP_COEFF_PMOS               Calib       Inverse dispersion      Y\n"
00086 "  CURV_COEFF_PMOS               Calib       Spectral curvature      Y\n"
00087 "  SLIT_LOCATION_PMOS            Calib       Slits positions table   Y\n"
00088 "  RETARDER_WAVEPLATE_CHROMATISM Calib       Chromatism correction   .\n"
00089 "\n"
00090 "Output files:\n\n"
00091 "  DO category:               Data type:  Explanation:\n"
00092 "  REDUCED_SCI_PMOS             FITS image  Extracted scientific spectra\n"
00093 "  REDUCED_SKY_SCI_PMOS         FITS image  Extracted sky spectra\n"
00094 "  REDUCED_ERROR_SCI_PMOS       FITS image  Errors on extracted spectra\n"
00095 "  REDUCED_X_SCI_PMOS           FITS image  X Stokes parameter (and L)\n"
00096 "  REDUCED_ERROR_X_SCI_PMOS     FITS image  Error on X Stokes parameter\n"
00097 "  REDUCED_NUL_X_SCI_PMOS       FITS image  Null parameter for X\n"
00098 "  REDUCED_ANGLE_SCI_PMOS       FITS image  Direction of linear polarization\n"
00099 "  REDUCED_ERROR_ANGLE_SCI_PMOS FITS image  Error on polarization direction\n"
00100 "  UNMAPPED_SCI_PMOS            FITS image  Sky subtracted scientific spectra\n"
00101 "  MAPPED_SCI_PMOS              FITS image  Rectified scientific spectra\n"
00102 "  MAPPED_ALL_SCI_PMOS          FITS image  Rectified science spectra with sky\n"
00103 "  MAPPED_SKY_SCI_PMOS          FITS image  Rectified sky spectra\n"
00104 "  UNMAPPED_SKY_SCI_PMOS        FITS image  Sky on CCD\n"
00105 "  OBJECT_TABLE_SCI_PMOS        FITS table  Positions of detected objects\n"
00106 "  OBJECT_TABLE_POL_SCI_PMOS    FITS table  Positions of real objects\n"
00107 "\n"
00108 "  Only if the sky-alignment of the wavelength solution is requested:\n"
00109 "  DISP_COEFF_SCI_PMOS          FITS table  Upgraded dispersion coefficients\n"
00110 "  WAVELENGTH_MAP_SCI_PMOS      FITS image  Upgraded wavelength map\n\n";
00111 
00112 #define fors_pmos_science_exit(message)            \
00113 {                                             \
00114 if (message) cpl_msg_error(recipe, message);  \
00115 cpl_free(instrume);                           \
00116 cpl_image_delete(dummy);                      \
00117 cpl_image_delete(mapped_sky);                 \
00118 cpl_image_delete(mapped_cleaned);             \
00119 cpl_image_delete(skymap);                     \
00120 cpl_image_delete(smapped);                    \
00121 cpl_table_delete(offsets);                    \
00122 cpl_table_delete(sky);                        \
00123 cpl_image_delete(bias);                       \
00124 cpl_image_delete(spectra);                    \
00125 cpl_image_delete(coordinate);                 \
00126 cpl_image_delete(norm_flat);                  \
00127 cpl_image_delete(rainbow);                    \
00128 cpl_image_delete(rectified);                  \
00129 cpl_image_delete(wavemap);                    \
00130 cpl_propertylist_delete(header);              \
00131 cpl_propertylist_delete(save_header);         \
00132 cpl_table_delete(grism_table);                \
00133 cpl_table_delete(idscoeff);                   \
00134 cpl_table_delete(maskslits);                  \
00135 cpl_table_delete(overscans);                  \
00136 cpl_table_delete(polytraces);                 \
00137 cpl_table_delete(wavelengths);                \
00138 cpl_vector_delete(lines);                     \
00139 cpl_msg_indent_less();                        \
00140 return -1;                                    \
00141 }
00142 
00143 
00144 #define fors_pmos_science_exit_memcheck(message)   \
00145 {                                             \
00146 if (message) cpl_msg_info(recipe, message);   \
00147 cpl_free(instrume);                           \
00148 cpl_image_delete(dummy);                      \
00149 cpl_image_delete(mapped_cleaned);             \
00150 cpl_image_delete(mapped_sky);                 \
00151 cpl_image_delete(skymap);                     \
00152 cpl_image_delete(smapped);                    \
00153 cpl_table_delete(offsets);                    \
00154 cpl_table_delete(sky);                        \
00155 cpl_image_delete(bias);                       \
00156 cpl_image_delete(spectra);                    \
00157 cpl_image_delete(coordinate);                 \
00158 cpl_image_delete(norm_flat);                  \
00159 cpl_image_delete(rainbow);                    \
00160 cpl_image_delete(rectified);                  \
00161 cpl_image_delete(wavemap);                    \
00162 cpl_propertylist_delete(header);              \
00163 cpl_propertylist_delete(save_header);         \
00164 cpl_table_delete(grism_table);                \
00165 cpl_table_delete(idscoeff);                   \
00166 cpl_table_delete(maskslits);                  \
00167 cpl_table_delete(overscans);                  \
00168 cpl_table_delete(polytraces);                 \
00169 cpl_table_delete(wavelengths);                \
00170 cpl_vector_delete(lines);                     \
00171 cpl_msg_indent_less();                        \
00172 return 0;                                     \
00173 }
00174 
00175 
00187 int cpl_plugin_get_info(cpl_pluginlist *list)
00188 {
00189     cpl_recipe *recipe = cpl_calloc(1, sizeof *recipe );
00190     cpl_plugin *plugin = &recipe->interface;
00191 
00192     cpl_plugin_init(plugin,
00193                     CPL_PLUGIN_API,
00194                     FORS_BINARY_VERSION,
00195                     CPL_PLUGIN_TYPE_RECIPE,
00196                     "fors_pmos_science",
00197                     "Extraction of scientific spectra",
00198                     fors_pmos_science_description,
00199                     "Carlo Izzo",
00200                     PACKAGE_BUGREPORT,
00201     "This file is currently part of the FORS Instrument Pipeline\n"
00202     "Copyright (C) 2002-2006 European Southern Observatory\n\n"
00203     "This program is free software; you can redistribute it and/or modify\n"
00204     "it under the terms of the GNU General Public License as published by\n"
00205     "the Free Software Foundation; either version 2 of the License, or\n"
00206     "(at your option) any later version.\n\n"
00207     "This program is distributed in the hope that it will be useful,\n"
00208     "but WITHOUT ANY WARRANTY; without even the implied warranty of\n"
00209     "MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the\n"
00210     "GNU General Public License for more details.\n\n"
00211     "You should have received a copy of the GNU General Public License\n"
00212     "along with this program; if not, write to the Free Software Foundation,\n"
00213     "Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA\n",
00214                     fors_pmos_science_create,
00215                     fors_pmos_science_exec,
00216                     fors_pmos_science_destroy);
00217 
00218     cpl_pluginlist_append(list, plugin);
00219     
00220     return 0;
00221 }
00222 
00223 
00234 static int fors_pmos_science_create(cpl_plugin *plugin)
00235 {
00236     cpl_recipe    *recipe;
00237     cpl_parameter *p;
00238 
00239 
00240     /* 
00241      * Check that the plugin is part of a valid recipe 
00242      */
00243 
00244     if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE) 
00245         recipe = (cpl_recipe *)plugin;
00246     else 
00247         return -1;
00248 
00249     /* 
00250      * Create the parameters list in the cpl_recipe object 
00251      */
00252 
00253     recipe->parameters = cpl_parameterlist_new(); 
00254 
00255 
00256     /*
00257      * Dispersion
00258      */
00259 
00260     p = cpl_parameter_new_value("fors.fors_pmos_science.dispersion",
00261                                 CPL_TYPE_DOUBLE,
00262                                 "Resampling step (Angstrom/pixel)",
00263                                 "fors.fors_pmos_science",
00264                                 0.0);
00265     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "dispersion");
00266     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00267     cpl_parameterlist_append(recipe->parameters, p);
00268 
00269     /*
00270      * Sky lines alignment
00271      */
00272 
00273     p = cpl_parameter_new_value("fors.fors_pmos_science.skyalign",
00274                                 CPL_TYPE_INT,
00275                                 "Polynomial order for sky lines alignment, "
00276                                 "or -1 to avoid alignment",
00277                                 "fors.fors_pmos_science",
00278                                 0);
00279     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "skyalign");
00280     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00281     cpl_parameterlist_append(recipe->parameters, p);
00282 
00283     /*
00284      * Line catalog table column containing the sky reference wavelengths
00285      */
00286 
00287     p = cpl_parameter_new_value("fors.fors_pmos_science.wcolumn",
00288                                 CPL_TYPE_STRING,
00289                                 "Name of sky line catalog table column "
00290                                 "with wavelengths",
00291                                 "fors.fors_pmos_science",
00292                                 "WLEN");
00293     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "wcolumn");
00294     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00295     cpl_parameterlist_append(recipe->parameters, p);
00296 
00297     /*
00298      * Start wavelength for spectral extraction
00299      */
00300 
00301     p = cpl_parameter_new_value("fors.fors_pmos_science.startwavelength",
00302                                 CPL_TYPE_DOUBLE,
00303                                 "Start wavelength in spectral extraction",
00304                                 "fors.fors_pmos_science",
00305                                 0.0);
00306     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "startwavelength");
00307     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00308     cpl_parameterlist_append(recipe->parameters, p);
00309 
00310     /*
00311      * End wavelength for spectral extraction
00312      */
00313 
00314     p = cpl_parameter_new_value("fors.fors_pmos_science.endwavelength",
00315                                 CPL_TYPE_DOUBLE,
00316                                 "End wavelength in spectral extraction",
00317                                 "fors.fors_pmos_science",
00318                                 0.0);
00319     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "endwavelength");
00320     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00321     cpl_parameterlist_append(recipe->parameters, p);
00322 
00323     /*
00324      * Flux conservation
00325      */
00326 
00327     p = cpl_parameter_new_value("fors.fors_pmos_science.flux",
00328                                 CPL_TYPE_BOOL,
00329                                 "Apply flux conservation",
00330                                 "fors.fors_pmos_science",
00331                                 TRUE);
00332     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "flux");
00333     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00334     cpl_parameterlist_append(recipe->parameters, p);
00335 
00336     /*
00337      * Apply flat field
00338      */
00339 
00340     p = cpl_parameter_new_value("fors.fors_pmos_science.flatfield",
00341                                 CPL_TYPE_BOOL,
00342                                 "Apply flat field",
00343                                 "fors.fors_pmos_science",
00344                                 TRUE);
00345     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "flatfield");
00346     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00347     cpl_parameterlist_append(recipe->parameters, p);
00348 
00349     p = cpl_parameter_new_value("fors.fors_pmos_science.skymedian",
00350                                 CPL_TYPE_BOOL,
00351                                 "Sky subtraction from extracted slit spectra",
00352                                 "fors.fors_pmos_science",
00353                                 FALSE);
00354     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "skymedian");
00355     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00356     cpl_parameterlist_append(recipe->parameters, p);
00357 
00358     /*
00359      * Local sky subtraction on CCD spectra
00360      */
00361 
00362     p = cpl_parameter_new_value("fors.fors_pmos_science.skylocal",
00363                                 CPL_TYPE_BOOL,
00364                                 "Sky subtraction from CCD slit spectra",
00365                                 "fors.fors_pmos_science",
00366                                 TRUE);
00367     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "skylocal");
00368     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00369     cpl_parameterlist_append(recipe->parameters, p);
00370 
00371     /*
00372      * Cosmic rays removal
00373      */
00374 
00375     p = cpl_parameter_new_value("fors.fors_pmos_science.cosmics",
00376                                 CPL_TYPE_BOOL,
00377                                 "Eliminate cosmic rays hits (only if local "
00378                                 "sky subtraction is also requested)",
00379                                 "fors.fors_pmos_science",
00380                                 FALSE);
00381     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "cosmics");
00382     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00383     cpl_parameterlist_append(recipe->parameters, p);
00384 
00385     /*
00386      * Slit margin
00387      */
00388 
00389     p = cpl_parameter_new_value("fors.fors_pmos_science.slit_margin",
00390                                 CPL_TYPE_INT,
00391                                 "Number of pixels to exclude at each slit "
00392                                 "in object detection and extraction",
00393                                 "fors.fors_pmos_science",
00394                                 3);
00395     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "slit_margin");
00396     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00397     cpl_parameterlist_append(recipe->parameters, p);
00398 
00399     /*
00400      * Extraction radius
00401      */
00402 
00403     p = cpl_parameter_new_value("fors.fors_pmos_science.ext_radius",
00404                                 CPL_TYPE_INT,
00405                                 "Maximum extraction radius for detected "
00406                                 "objects (pixel)",
00407                                 "fors.fors_pmos_science",
00408                                 6);
00409     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "ext_radius");
00410     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00411     cpl_parameterlist_append(recipe->parameters, p);
00412 
00413     /*
00414      * Contamination radius
00415      */
00416 
00417     p = cpl_parameter_new_value("fors.fors_pmos_science.cont_radius",
00418                                 CPL_TYPE_INT,
00419                                 "Minimum distance at which two objects "
00420                                 "of equal luminosity do not contaminate "
00421                                 "each other (pixel)",
00422                                 "fors.fors_pmos_science",
00423                                 0);
00424     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "cont_radius");
00425     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00426     cpl_parameterlist_append(recipe->parameters, p);
00427 
00428     /*
00429      * Object extraction method
00430      */
00431 
00432     p = cpl_parameter_new_value("fors.fors_pmos_science.ext_mode",
00433                                 CPL_TYPE_INT,
00434                                 "Object extraction method: 0 = aperture, "
00435                                 "1 = Horne optimal extraction",
00436                                 "fors.fors_pmos_science",
00437                                 1);
00438     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "ext_mode");
00439     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00440     cpl_parameterlist_append(recipe->parameters, p);
00441 
00442     /*
00443      * Normalise output by exposure time
00444      */
00445 
00446     p = cpl_parameter_new_value("fors.fors_pmos_science.time_normalise",
00447                                 CPL_TYPE_BOOL,
00448                                 "Normalise output spectra by the exposure time",
00449                                 "fors.fors_pmos_science",
00450                                 TRUE);
00451     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "time_normalise");
00452     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00453     cpl_parameterlist_append(recipe->parameters, p);
00454 
00455     /*
00456      * Apply chromatism correction to polarization angle
00457      */
00458 
00459     p = cpl_parameter_new_value("fors.fors_pmos_science.chromatism",
00460                                 CPL_TYPE_BOOL,
00461                                 "Chromatism correction to polarization angles",
00462                                 "fors.fors_pmos_science",
00463                                 TRUE);
00464     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "chromatism");
00465     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00466     cpl_parameterlist_append(recipe->parameters, p);
00467 
00468     /*
00469      * Create check products
00470      */
00471 
00472     p = cpl_parameter_new_value("fors.fors_pmos_science.check",
00473                                 CPL_TYPE_BOOL,
00474                                 "Create intermediate products",
00475                                 "fors.fors_pmos_science",
00476                                 FALSE);
00477     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "check");
00478     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00479     cpl_parameterlist_append(recipe->parameters, p);
00480 
00481     return 0;
00482 }
00483 
00484 
00493 static int fors_pmos_science_exec(cpl_plugin *plugin)
00494 {
00495     cpl_recipe *recipe;
00496     
00497     if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE) 
00498         recipe = (cpl_recipe *)plugin;
00499     else 
00500         return -1;
00501 
00502     return fors_pmos_science(recipe->parameters, recipe->frames);
00503 }
00504 
00505 
00514 static int fors_pmos_science_destroy(cpl_plugin *plugin)
00515 {
00516     cpl_recipe *recipe;
00517     
00518     if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE) 
00519         recipe = (cpl_recipe *)plugin;
00520     else 
00521         return -1;
00522 
00523     cpl_parameterlist_delete(recipe->parameters); 
00524 
00525     return 0;
00526 }
00527 
00528 
00538 static int fors_pmos_science(cpl_parameterlist *parlist, cpl_frameset *frameset)
00539 {
00540 
00541     const char *recipe = "fors_pmos_science";
00542 
00543 
00544     /*
00545      * Input parameters
00546      */
00547 
00548     double      dispersion;
00549     int         skyalign;
00550     const char *wcolumn;
00551     double      startwavelength;
00552     double      endwavelength;
00553     int         flux;
00554     int         flatfield;
00555     int         skylocal;
00556     int         skymedian;
00557     int         chromatism;
00558     int         cosmics;
00559     int         slit_margin;
00560     int         ext_radius;
00561     int         cont_radius;
00562     int         ext_mode;
00563     int         time_normalise;
00564     int         check;
00565 
00566     /*
00567      * CPL objects
00568      */
00569 
00570     cpl_image       **images;
00571 
00572     cpl_image       **reduceds       = NULL;
00573     cpl_image       **rerrors        = NULL;
00574     cpl_table       **slitss         = NULL;
00575     cpl_image       **mappeds        = NULL;
00576     cpl_image       **skylocalmaps   = NULL;
00577     
00578     int nobjects = 0;
00579 
00580     cpl_image        *bias           = NULL;
00581     cpl_image        *norm_flat      = NULL;
00582     cpl_image        *spectra        = NULL;
00583     cpl_image        *rectified      = NULL;
00584     cpl_image        *coordinate     = NULL;
00585     cpl_image        *rainbow        = NULL;
00586     cpl_image        *mapped         = NULL;
00587     cpl_image        *mapped_sky     = NULL;
00588     cpl_image        *mapped_cleaned = NULL;
00589     cpl_image        *smapped        = NULL;
00590     cpl_image        *wavemap        = NULL;
00591     cpl_image        *skymap         = NULL;
00592     cpl_image        *skylocalmap    = NULL;
00593     cpl_image        *dummy          = NULL;
00594 
00595     cpl_table        *grism_table    = NULL;
00596     cpl_table        *overscans      = NULL;
00597     cpl_table        *wavelengths    = NULL;
00598     cpl_table        *idscoeff       = NULL;
00599     cpl_table        *slits          = NULL;
00600     cpl_table        *origslits     = NULL;
00601     cpl_table        *maskslits      = NULL;
00602     cpl_table        *polytraces     = NULL;
00603     cpl_table        *offsets        = NULL;
00604     cpl_table        *sky            = NULL;
00605 
00606     cpl_vector       *lines          = NULL;
00607 
00608     cpl_propertylist *header         = NULL;
00609     cpl_propertylist *save_header    = NULL;
00610 
00611     /*
00612      * Auxiliary variables
00613      */
00614 
00615     char    version[80];
00616     char   *instrume = NULL;
00617     const char   *science_tag;
00618     const char   *master_norm_flat_tag;
00619     const char   *disp_coeff_tag;
00620     const char   *disp_coeff_sky_tag;
00621     const char   *wavelength_map_sky_tag;
00622     const char   *curv_coeff_tag;
00623     const char   *slit_location_tag;
00624     const char   *reduced_science_tag;
00625     const char   *reduced_sky_tag;
00626     const char   *reduced_error_tag;
00627     const char   *mapped_science_tag;
00628     const char   *unmapped_science_tag;
00629     const char   *mapped_science_sky_tag;
00630     const char   *mapped_sky_tag;
00631     const char   *unmapped_sky_tag;
00632     const char   *object_table_tag;
00633     const char   *object_table_pol_tag;
00634     const char   *skylines_offsets_tag;
00635     const char   *reduced_q_tag;
00636     const char   *reduced_u_tag;
00637     const char   *reduced_v_tag;
00638     const char   *reduced_l_tag;
00639     const char   *reduced_error_q_tag;
00640     const char   *reduced_error_u_tag;
00641     const char   *reduced_error_v_tag;
00642     const char   *reduced_error_l_tag;
00643     const char   *reduced_nul_q_tag;
00644     const char   *reduced_nul_u_tag;
00645     const char   *reduced_nul_v_tag;
00646     const char   *reduced_angle_tag;
00647     const char   *reduced_error_angle_tag;
00648     const char   *chrom_table_tag = "RETARDER_WAVEPLATE_CHROMATISM";
00649     float *angles = NULL;
00650     int     pmos, circ;
00651     int     nscience;
00652     double  alltime;
00653     double  mean_rms;
00654     int     nlines;
00655     int     rebin;
00656     double *line;
00657     int     nx = 0, ny;
00658     int     ccd_xsize, ccd_ysize;
00659     double  reference;
00660     double  gain;
00661     double  ron;
00662     int     standard;
00663     int     highres;
00664     int     i, j;
00665     cpl_error_code error;
00666 
00667     int * nobjs_per_slit;
00668 
00669     snprintf(version, 80, "%s-%s", PACKAGE, PACKAGE_VERSION);
00670 
00671     cpl_msg_set_indentation(2);
00672 
00673     if (dfs_files_dont_exist(frameset))
00674         fors_pmos_science_exit(NULL);
00675 
00676 
00677     /* 
00678      * Get configuration parameters
00679      */
00680 
00681     cpl_msg_info(recipe, "Recipe %s configuration parameters:", recipe);
00682     cpl_msg_indent_more();
00683 
00684     if (cpl_frameset_count_tags(frameset, "GRISM_TABLE") > 1)
00685         fors_pmos_science_exit("Too many in input: GRISM_TABLE");
00686 
00687     grism_table = dfs_load_table(frameset, "GRISM_TABLE", 1);
00688 
00689     dispersion = dfs_get_parameter_double(parlist, 
00690                     "fors.fors_pmos_science.dispersion", grism_table);
00691 
00692     if (dispersion <= 0.0)
00693         fors_pmos_science_exit("Invalid resampling step");
00694 
00695     skyalign = dfs_get_parameter_int(parlist, 
00696                     "fors.fors_pmos_science.skyalign", NULL);
00697 
00698     if (skyalign > 2)
00699         fors_pmos_science_exit("Max polynomial degree for sky alignment is 2");
00700 
00701     wcolumn = dfs_get_parameter_string(parlist, 
00702                     "fors.fors_pmos_science.wcolumn", NULL);
00703 
00704     startwavelength = dfs_get_parameter_double(parlist, 
00705                     "fors.fors_pmos_science.startwavelength", grism_table);
00706     if (startwavelength < 3000.0 || startwavelength > 13000.0)
00707         fors_pmos_science_exit("Invalid wavelength");
00708 
00709     endwavelength = dfs_get_parameter_double(parlist, 
00710                     "fors.fors_pmos_science.endwavelength", grism_table);
00711     if (endwavelength < 3000.0 || endwavelength > 13000.0)
00712         fors_pmos_science_exit("Invalid wavelength");
00713 
00714     if (endwavelength - startwavelength <= 0.0)
00715         fors_pmos_science_exit("Invalid wavelength interval");
00716 
00717     flux = dfs_get_parameter_bool(parlist, "fors.fors_pmos_science.flux", NULL);
00718 
00719     flatfield = dfs_get_parameter_bool(parlist, "fors.fors_pmos_science.flatfield", 
00720                                        NULL);
00721 
00722     skylocal  = dfs_get_parameter_bool(parlist, "fors.fors_pmos_science.skylocal", 
00723                                        NULL);
00724     skymedian = dfs_get_parameter_bool(parlist, "fors.fors_pmos_science.skymedian", 
00725                                        NULL);
00726 /* NSS
00727     skymedian = dfs_get_parameter_int(parlist, "fors.fors_pmos_science.skymedian", 
00728                                        NULL);
00729 */
00730     
00731     chromatism = 
00732         dfs_get_parameter_bool(parlist, "fors.fors_pmos_science.chromatism", 
00733                                NULL);
00734 
00735     if (skylocal && skymedian)
00736         fors_pmos_science_exit("Cannot apply sky subtraction both on extracted "
00737                           "and non-extracted spectra");
00738 
00739     cosmics = dfs_get_parameter_bool(parlist, 
00740                                      "fors.fors_pmos_science.cosmics", NULL);
00741 
00742     if (cosmics)
00743         if (!skylocal)
00744             fors_pmos_science_exit("Cosmic rays correction requires "
00745                                    "skylocal=true");
00746 
00747     slit_margin = dfs_get_parameter_int(parlist, 
00748                                         "fors.fors_pmos_science.slit_margin",
00749                                         NULL);
00750     if (slit_margin < 0)
00751         fors_pmos_science_exit("Value must be zero or positive");
00752 
00753     ext_radius = dfs_get_parameter_int(parlist, "fors.fors_pmos_science.ext_radius",
00754                                        NULL);
00755     if (ext_radius < 0)
00756         fors_pmos_science_exit("Value must be zero or positive");
00757 
00758     cont_radius = dfs_get_parameter_int(parlist, 
00759                                         "fors.fors_pmos_science.cont_radius",
00760                                         NULL);
00761     if (cont_radius < 0)
00762         fors_pmos_science_exit("Value must be zero or positive");
00763 
00764     ext_mode = dfs_get_parameter_int(parlist, "fors.fors_pmos_science.ext_mode",
00765                                        NULL);
00766     if (ext_mode < 0 || ext_mode > 1)
00767         fors_pmos_science_exit("Invalid object extraction mode");
00768 
00769     time_normalise = dfs_get_parameter_bool(parlist, 
00770                              "fors.fors_pmos_science.time_normalise", NULL);
00771 
00772     check = dfs_get_parameter_bool(parlist, "fors.fors_pmos_science.check", NULL);
00773     cpl_table_delete(grism_table); grism_table = NULL;
00774 
00775     if (cpl_error_get_code())
00776         fors_pmos_science_exit("Failure getting the configuration parameters");
00777 
00778     
00779     /* 
00780      * Check input set-of-frames
00781      */
00782 
00783     cpl_msg_indent_less();
00784     cpl_msg_info(recipe, "Check input set-of-frames:");
00785     cpl_msg_indent_more();
00786 
00787     if (!dfs_equal_keyword(frameset, "ESO INS GRIS1 ID"))
00788         fors_pmos_science_exit("Input frames are not from the same grism");
00789 
00790     if (!dfs_equal_keyword(frameset, "ESO INS FILT1 ID"))
00791         fors_pmos_science_exit("Input frames are not from the same filter");
00792 
00793     if (!dfs_equal_keyword(frameset, "ESO DET CHIP1 ID"))
00794         fors_pmos_science_exit("Input frames are not from the same chip");
00795 
00796     standard = 0;
00797     pmos = cpl_frameset_count_tags(frameset, "SCIENCE_PMOS");
00798 
00799     if (pmos == 0) {
00800         pmos = cpl_frameset_count_tags(frameset, "STANDARD_PMOS");
00801         standard = 1;
00802     }
00803 
00804     if (pmos == 0)
00805         fors_pmos_science_exit("Missing input scientific frame");
00806 
00807     angles = fors_check_angles(frameset, pmos, 
00808                                 standard ? "STANDARD_PMOS" : "SCIENCE_PMOS", 
00809                                 &circ);
00810     if (angles == NULL)
00811         fors_pmos_science_exit("Polarization angles could not be read");
00812 
00813     if (circ)
00814         chromatism = 0; /* Chromatism correction unrequired for 
00815                            circular polarimetry */
00816 
00817 
00818     nscience = pmos;
00819 
00820     reduceds = (cpl_image **)cpl_malloc(sizeof(cpl_image *) * nscience);
00821     rerrors  = (cpl_image **)cpl_malloc(sizeof(cpl_image *) * nscience);
00822     slitss   = (cpl_table **)cpl_malloc(sizeof(cpl_table *) * nscience);
00823     mappeds  = (cpl_image **)cpl_malloc(sizeof(cpl_image *) * nscience);
00824     skylocalmaps = (cpl_image **)cpl_malloc(sizeof(cpl_image *) * nscience);
00825 
00826     if (pmos) {
00827         cpl_msg_info(recipe, "PMOS data found");
00828         if (standard) {
00829             science_tag             = "STANDARD_PMOS";
00830             reduced_science_tag     = "REDUCED_STD_PMOS";
00831             unmapped_science_tag    = "UNMAPPED_STD_PMOS";
00832             mapped_science_tag      = "MAPPED_STD_PMOS";
00833             mapped_science_sky_tag  = "MAPPED_ALL_STD_PMOS";
00834             skylines_offsets_tag    = "SKY_SHIFTS_SLIT_STD_PMOS";
00835             wavelength_map_sky_tag  = "WAVELENGTH_MAP_STD_PMOS";
00836             disp_coeff_sky_tag      = "DISP_COEFF_STD_PMOS";
00837             mapped_sky_tag          = "MAPPED_SKY_STD_PMOS";
00838             unmapped_sky_tag        = "UNMAPPED_SKY_STD_PMOS";
00839             object_table_tag        = "OBJECT_TABLE_STD_PMOS";
00840             object_table_pol_tag    = "OBJECT_TABLE_POL_STD_PMOS";
00841             reduced_sky_tag         = "REDUCED_SKY_STD_PMOS";
00842             reduced_error_tag       = "REDUCED_ERROR_STD_PMOS";
00843             reduced_q_tag           = "REDUCED_Q_STD_PMOS";
00844             reduced_u_tag           = "REDUCED_U_STD_PMOS";
00845             reduced_v_tag           = "REDUCED_V_STD_PMOS";
00846             reduced_l_tag           = "REDUCED_L_STD_PMOS";
00847             reduced_error_q_tag     = "REDUCED_ERROR_Q_STD_PMOS";
00848             reduced_error_u_tag     = "REDUCED_ERROR_U_STD_PMOS";
00849             reduced_error_v_tag     = "REDUCED_ERROR_V_STD_PMOS";
00850             reduced_error_l_tag     = "REDUCED_ERROR_L_STD_PMOS";
00851             reduced_nul_q_tag       = "REDUCED_NUL_Q_STD_PMOS";
00852             reduced_nul_u_tag       = "REDUCED_NUL_U_STD_PMOS";
00853             reduced_nul_v_tag       = "REDUCED_NUL_V_STD_PMOS";
00854             reduced_angle_tag       = "REDUCED_ANGLE_STD_PMOS";
00855             reduced_error_angle_tag = "REDUCED_ERROR_ANGLE_STD_PMOS";
00856         }
00857         else {
00858             science_tag             = "SCIENCE_PMOS";
00859             reduced_science_tag     = "REDUCED_SCI_PMOS";
00860             unmapped_science_tag    = "UNMAPPED_SCI_PMOS";
00861             mapped_science_tag      = "MAPPED_SCI_PMOS";
00862             mapped_science_sky_tag  = "MAPPED_ALL_SCI_PMOS";
00863             skylines_offsets_tag    = "SKY_SHIFTS_SLIT_SCI_PMOS";
00864             wavelength_map_sky_tag  = "WAVELENGTH_MAP_SCI_PMOS";
00865             disp_coeff_sky_tag      = "DISP_COEFF_SCI_PMOS";
00866             mapped_sky_tag          = "MAPPED_SKY_SCI_PMOS";
00867             unmapped_sky_tag        = "UNMAPPED_SKY_SCI_PMOS";
00868             object_table_tag        = "OBJECT_TABLE_SCI_PMOS";
00869             object_table_pol_tag    = "OBJECT_TABLE_POL_SCI_PMOS";
00870             reduced_sky_tag         = "REDUCED_SKY_SCI_PMOS";
00871             reduced_error_tag       = "REDUCED_ERROR_SCI_PMOS";
00872             reduced_q_tag           = "REDUCED_Q_SCI_PMOS";
00873             reduced_u_tag           = "REDUCED_U_SCI_PMOS";
00874             reduced_v_tag           = "REDUCED_V_SCI_PMOS";
00875             reduced_l_tag           = "REDUCED_L_SCI_PMOS";
00876             reduced_error_q_tag     = "REDUCED_ERROR_Q_SCI_PMOS";
00877             reduced_error_u_tag     = "REDUCED_ERROR_U_SCI_PMOS";
00878             reduced_error_v_tag     = "REDUCED_ERROR_V_SCI_PMOS";
00879             reduced_error_l_tag     = "REDUCED_ERROR_L_SCI_PMOS";
00880             reduced_nul_q_tag       = "REDUCED_NUL_Q_SCI_PMOS";
00881             reduced_nul_u_tag       = "REDUCED_NUL_U_SCI_PMOS";
00882             reduced_nul_v_tag       = "REDUCED_NUL_V_SCI_PMOS";
00883             reduced_angle_tag       = "REDUCED_ANGLE_SCI_PMOS";
00884             reduced_error_angle_tag = "REDUCED_ERROR_ANGLE_SCI_PMOS";
00885         }
00886 
00887         master_norm_flat_tag    = "MASTER_NORM_FLAT_PMOS";
00888         disp_coeff_tag          = "DISP_COEFF_PMOS";
00889         curv_coeff_tag          = "CURV_COEFF_PMOS";
00890         slit_location_tag       = "SLIT_LOCATION_PMOS";
00891 
00892         if (!cpl_frameset_count_tags(frameset, master_norm_flat_tag)) {
00893             master_norm_flat_tag    = "MASTER_NORM_FLAT_LONG_PMOS";
00894             disp_coeff_tag          = "DISP_COEFF_LONG_PMOS";
00895             slit_location_tag       = "SLIT_LOCATION_LONG_PMOS";
00896         }
00897     }
00898 
00899     if (cpl_frameset_count_tags(frameset, "MASTER_BIAS") == 0)
00900         fors_pmos_science_exit("Missing required input: MASTER_BIAS");
00901 
00902     if (cpl_frameset_count_tags(frameset, "MASTER_BIAS") > 1)
00903         fors_pmos_science_exit("Too many in input: MASTER_BIAS");
00904 
00905     if (skyalign >= 0)
00906         if (cpl_frameset_count_tags(frameset, "MASTER_SKYLINECAT") > 1)
00907             fors_pmos_science_exit("Too many in input: MASTER_SKYLINECAT");
00908 
00909     if (cpl_frameset_count_tags(frameset, disp_coeff_tag) == 0) {
00910         cpl_msg_error(recipe, "Missing required input: %s", disp_coeff_tag);
00911         fors_pmos_science_exit(NULL);
00912     }
00913 
00914     if (cpl_frameset_count_tags(frameset, disp_coeff_tag) > 1) {
00915         cpl_msg_error(recipe, "Too many in input: %s", disp_coeff_tag);
00916         fors_pmos_science_exit(NULL);
00917     }
00918 
00919     if (cpl_frameset_count_tags(frameset, slit_location_tag) == 0) {
00920         cpl_msg_error(recipe, "Missing required input: %s",
00921                       slit_location_tag);
00922         fors_pmos_science_exit(NULL);
00923     }
00924 
00925     if (cpl_frameset_count_tags(frameset, slit_location_tag) > 1) {
00926         cpl_msg_error(recipe, "Too many in input: %s", slit_location_tag);
00927         fors_pmos_science_exit(NULL);
00928     }
00929 
00930     if (chromatism) {
00931         if (cpl_frameset_count_tags(frameset, chrom_table_tag) == 0) {
00932             cpl_msg_error(recipe, "Missing required input: %s",
00933                           chrom_table_tag);
00934             fors_pmos_science_exit(NULL);
00935         }
00936 
00937         if (cpl_frameset_count_tags(frameset, chrom_table_tag) > 1) {
00938             cpl_msg_error(recipe, "Too many in input: %s", chrom_table_tag);
00939             fors_pmos_science_exit(NULL);
00940         }
00941     }
00942 
00943     if (cpl_frameset_count_tags(frameset, master_norm_flat_tag) > 1) {
00944         if (flatfield) {
00945             cpl_msg_error(recipe, "Too many in input: %s", 
00946                           master_norm_flat_tag);
00947             fors_pmos_science_exit(NULL);
00948         }
00949         else {
00950             cpl_msg_warning(recipe, "%s in input are ignored, "
00951                             "since flat field correction was not requested", 
00952                             master_norm_flat_tag);
00953         }
00954     }
00955 
00956     if (cpl_frameset_count_tags(frameset, master_norm_flat_tag) == 1) {
00957         if (!flatfield) {
00958             cpl_msg_warning(recipe, "%s in input is ignored, "
00959                             "since flat field correction was not requested", 
00960                             master_norm_flat_tag);
00961         }
00962     }
00963 
00964     if (cpl_frameset_count_tags(frameset, master_norm_flat_tag) == 0) {
00965         if (flatfield) {
00966             cpl_msg_error(recipe, "Flat field correction was requested, "
00967                           "but no %s are found in input",
00968                           master_norm_flat_tag);
00969             fors_pmos_science_exit(NULL);
00970         }
00971     }
00972 
00973     cpl_msg_indent_less();
00974 
00975     /*
00976      * Get the reference wavelength and the rebin factor along the
00977      * dispersion direction from a scientific exposure
00978      */
00979 
00980     header = dfs_load_header(frameset, science_tag, 0);
00981 
00982     if (header == NULL)
00983         fors_pmos_science_exit("Cannot load scientific frame header");
00984 
00985     instrume = (char *)cpl_propertylist_get_string(header, "INSTRUME");
00986     if (instrume == NULL)
00987         fors_pmos_science_exit("Missing keyword INSTRUME in scientific header");
00988     instrume = cpl_strdup(instrume);
00989 
00990     if (instrume[4] == '1')
00991         snprintf(version, 80, "%s/%s", "fors1", VERSION);
00992     if (instrume[4] == '2')
00993         snprintf(version, 80, "%s/%s", "fors2", VERSION);
00994 
00995     reference = cpl_propertylist_get_double(header, "ESO INS GRIS1 WLEN");
00996 
00997     if (cpl_error_get_code() != CPL_ERROR_NONE)
00998         fors_pmos_science_exit("Missing keyword ESO INS GRIS1 WLEN in scientific "
00999                         "frame header");
01000 
01001     if (reference < 3000.0)   /* Perhaps in nanometers... */
01002         reference *= 10;
01003 
01004     if (reference < 3000.0 || reference > 13000.0) {
01005         cpl_msg_error(recipe, "Invalid central wavelength %.2f read from "
01006                       "keyword ESO INS GRIS1 WLEN in scientific frame header",
01007                       reference);
01008         fors_pmos_science_exit(NULL);
01009     }
01010 
01011     cpl_msg_info(recipe, "The central wavelength is: %.2f", reference);
01012 
01013     rebin = cpl_propertylist_get_int(header, "ESO DET WIN1 BINX");
01014 
01015     if (cpl_error_get_code() != CPL_ERROR_NONE)
01016         fors_pmos_science_exit("Missing keyword ESO DET WIN1 BINX in scientific "
01017                         "frame header");
01018 
01019     if (rebin != 1) {
01020         dispersion *= rebin;
01021         cpl_msg_warning(recipe, "The rebin factor is %d, and therefore the "
01022                         "resampling step used is %f A/pixel", rebin, 
01023                         dispersion);
01024     }
01025 
01026     gain = cpl_propertylist_get_double(header, "ESO DET OUT1 CONAD");
01027 
01028     if (cpl_error_get_code() != CPL_ERROR_NONE)
01029         fors_pmos_science_exit("Missing keyword ESO DET OUT1 CONAD in scientific "
01030                           "frame header");
01031 
01032     cpl_msg_info(recipe, "The gain factor is: %.2f e-/ADU", gain);
01033 
01034     ron = cpl_propertylist_get_double(header, "ESO DET OUT1 RON");
01035 
01036     if (cpl_error_get_code() != CPL_ERROR_NONE)
01037         fors_pmos_science_exit("Missing keyword ESO DET OUT1 RON in scientific "
01038                           "frame header");
01039 
01040     ron /= gain;     /* Convert from electrons to ADU */
01041 
01042     cpl_msg_info(recipe, "The read-out-noise is: %.2f ADU", ron);
01043 
01044     if (cpl_frameset_count_tags(frameset, curv_coeff_tag) == 0) {
01045         cpl_msg_error(recipe, "Missing required input: %s", curv_coeff_tag);
01046         fors_pmos_science_exit(NULL);
01047     }
01048 
01049     if (cpl_frameset_count_tags(frameset, curv_coeff_tag) > 1) {
01050         cpl_msg_error(recipe, "Too many in input: %s", curv_coeff_tag);
01051         fors_pmos_science_exit(NULL);
01052     }
01053 
01054     cpl_msg_info(recipe, "Load normalised flat field (if present)...");
01055     cpl_msg_indent_more();
01056 
01057     if (flatfield) {
01058         norm_flat = dfs_load_image(frameset, master_norm_flat_tag, 
01059                                    CPL_TYPE_FLOAT, 0, 1);
01060     }
01061 
01062     if (skyalign >= 0) {
01063 
01064         cpl_msg_indent_less();
01065         cpl_msg_info(recipe, "Load input sky line catalog...");
01066         cpl_msg_indent_more();
01067 
01068         wavelengths = dfs_load_table(frameset, "MASTER_SKYLINECAT", 1);
01069 
01070         if (wavelengths) {
01071             /*
01072              * Cast the wavelengths into a (double precision) CPL vector
01073              */
01074 
01075             nlines = cpl_table_get_nrow(wavelengths);
01076 
01077             if (nlines == 0)
01078                 fors_pmos_science_exit("Empty input sky line catalog");
01079 
01080             if (cpl_table_has_column(wavelengths, wcolumn) != 1) {
01081                 cpl_msg_error(recipe, "Missing column %s in input line "
01082                               "catalog table", wcolumn);
01083                 fors_pmos_science_exit(NULL);
01084             }
01085 
01086             line = cpl_malloc(nlines * sizeof(double));
01087     
01088             for (i = 0; i < nlines; i++)
01089                 line[i] = cpl_table_get(wavelengths, wcolumn, i, NULL);
01090 
01091             cpl_table_delete(wavelengths); wavelengths = NULL;
01092 
01093             lines = cpl_vector_wrap(nlines, line);
01094         }
01095         else {
01096             cpl_msg_info(recipe, "No sky line catalog found in input - fine!");
01097         }
01098     }
01099 
01100 
01101     cpl_propertylist_delete(header); header = NULL;
01102 
01103     /*
01104      * Load the wavelength calibration table
01105      */
01106 
01107     idscoeff = dfs_load_table(frameset, disp_coeff_tag, 1);
01108     if (idscoeff == NULL)
01109         fors_pmos_science_exit("Cannot load wavelength calibration table");
01110 
01111     for (j = 0; j < nscience; j++) {
01112         int k;
01113 
01114         cpl_msg_indent_less();
01115         cpl_msg_info(recipe, "Processing scientific exposure of angle %.2f "
01116                      "(%d out of %d) ...",
01117                      angles[j], j + 1, nscience);
01118         cpl_msg_indent_more();
01119 
01120         cpl_msg_info(recipe, "Load scientific exposure...");
01121         cpl_msg_indent_more();
01122 
01123         /*
01124          * FIXME: Horrible workaround to avoid the problem because of the
01125          * multiple encapsulation of cpl_frameset_find() in different 
01126          * loading functions
01127          */
01128 
01129         header = dfs_load_header(frameset, science_tag, 0);
01130 
01131         for (k = 0; k < j; k ++) {
01132             cpl_propertylist_delete(header);
01133             header = dfs_load_header(frameset, NULL, 0);
01134         }
01135 
01136         spectra = dfs_load_image(frameset, science_tag, CPL_TYPE_FLOAT, 0, 0);
01137 
01138         for (k = 0; k < j; k ++) {
01139             cpl_image_delete(spectra);
01140             spectra = dfs_load_image(frameset, NULL, CPL_TYPE_FLOAT, 0, 0);
01141         }
01142 
01143         if (spectra == NULL)
01144             fors_pmos_science_exit("Cannot load scientific frame");
01145             
01146         if (header == NULL)
01147             fors_pmos_science_exit("Cannot load scientific frame header");
01148 
01149         alltime = cpl_propertylist_get_double(header, "EXPTIME");
01150 
01151         if (cpl_error_get_code() != CPL_ERROR_NONE)
01152             fors_pmos_science_exit("Missing keyword EXPTIME in scientific "
01153                                    "frame header");
01154 
01155         /* Leave the header on for the next step... */
01156         //cpl_propertylist_delete(header); header = NULL;
01157 
01158         cpl_msg_info(recipe, "Scientific frame exposure time: %.2f s", 
01159                      alltime);
01160 
01161         cpl_msg_indent_less();
01162 
01163         /*
01164          * Remove the master bias
01165          */
01166 
01167         cpl_msg_info(recipe, "Remove the master bias...");
01168 
01169         bias = dfs_load_image(frameset, "MASTER_BIAS", CPL_TYPE_FLOAT, 0, 1);
01170 
01171         if (bias == NULL)
01172             fors_pmos_science_exit("Cannot load master bias");
01173 
01174         overscans = mos_load_overscans_vimos(header, 1);
01175 
01176         dummy = mos_remove_bias(spectra, bias, overscans);
01177         cpl_image_delete(spectra); spectra = dummy; dummy = NULL;
01178         cpl_image_delete(bias); bias = NULL;
01179         cpl_table_delete(overscans); overscans = NULL;
01180 
01181         if (spectra == NULL)
01182             fors_pmos_science_exit("Cannot remove bias from scientific frame");
01183 
01184         ccd_xsize = nx = cpl_image_get_size_x(spectra);
01185         ccd_ysize = ny = cpl_image_get_size_y(spectra);
01186 
01187         if (flatfield) {
01188 
01189             if (norm_flat) {
01190                 cpl_msg_info(recipe, "Apply flat field correction...");
01191                 if (cpl_image_divide(spectra, norm_flat) != CPL_ERROR_NONE) {
01192                     cpl_msg_error(recipe, "Failure of flat field correction: %s",
01193                                   cpl_error_get_message());
01194                     fors_pmos_science_exit(NULL);
01195                 }
01196             }
01197             else {
01198                 cpl_msg_error(recipe, "Cannot load input %s for flat field "
01199                               "correction", master_norm_flat_tag);
01200                 fors_pmos_science_exit(NULL);
01201             }
01202 
01203         }
01204 
01205         /*
01206          * Load the spectral curvature table
01207          */
01208 
01209         polytraces = dfs_load_table(frameset, curv_coeff_tag, 1);
01210         if (polytraces == NULL)
01211             fors_pmos_science_exit("Cannot load spectral curvature table");
01212 
01213         /*
01214          * Load the slit location table
01215          */
01216 
01217         slits = dfs_load_table(frameset, slit_location_tag, 1);
01218         if (slits == NULL)
01219             fors_pmos_science_exit("Cannot load slits location table");
01220 
01221         cpl_msg_info(recipe, "Processing scientific spectra...");
01222 
01223         /*
01224          * This one will also generate the spatial map from the spectral 
01225          * curvature table (in the case of multislit data)
01226          */
01227 
01228         coordinate = cpl_image_new(nx, ny, CPL_TYPE_FLOAT);
01229 
01230         smapped = mos_spatial_calibration(spectra, slits, polytraces, reference,
01231                                           startwavelength, endwavelength,
01232                                           dispersion, flux, coordinate);
01233 
01234         /*
01235          * Generate a rectified wavelength map from the wavelength calibration 
01236          * table
01237          */
01238 
01239         rainbow = mos_map_idscoeff(idscoeff, nx, reference, startwavelength, 
01240                                    endwavelength);
01241 
01242         if (dispersion > 1.0)
01243             highres = 0;
01244         else
01245             highres = 1;
01246 
01247         if (skyalign >= 0) {
01248             if (skyalign) {
01249                 cpl_msg_info(recipe, "Align wavelength solution to reference "
01250                              "skylines applying %d order residual fit...", skyalign);
01251             }
01252             else {
01253                 cpl_msg_info(recipe, "Align wavelength solution to reference "
01254                              "skylines applying median offset...");
01255             }
01256 
01257             if (!j) {
01258                 offsets = mos_wavelength_align(smapped, slits, reference, 
01259                                                startwavelength, endwavelength, 
01260                                                idscoeff, lines, highres, 
01261                                                skyalign, rainbow, 4);
01262                 if (offsets) {
01263                     if (standard)
01264                         cpl_msg_warning(recipe, "Alignment of the wavelength solution "
01265                                         "to reference sky lines may be unreliable in "
01266                                         "this case!");
01267 
01268                     if (dfs_save_table(frameset, offsets, skylines_offsets_tag,
01269                                        NULL, parlist, recipe, version)) {
01270                         fors_pmos_science_exit(NULL);
01271                     }
01272 
01273                 } else {
01274                     cpl_msg_warning(recipe, "Alignment of the wavelength solution "
01275                                     "to reference sky lines could not be done!");
01276                     skyalign = -1;
01277                 }
01278             }
01279 
01280 
01281         }
01282 
01283         wavemap = mos_map_wavelengths(coordinate, rainbow, slits, 
01284                                       polytraces, reference, 
01285                                       startwavelength, endwavelength,
01286                                       dispersion);
01287 
01288 
01289         cpl_image_delete(rainbow); rainbow = NULL;
01290         cpl_image_delete(coordinate); coordinate = NULL;
01291 
01292         /*
01293          * Here the wavelength calibrated slit spectra are created. This frame
01294          * contains sky_science.
01295          */
01296 
01297         mapped_sky = mos_wavelength_calibration(smapped, reference,
01298                                                 startwavelength, endwavelength,
01299                                                 dispersion, idscoeff, flux);
01300 
01301         if (!j) {
01302             cpl_msg_indent_less();
01303             cpl_msg_info(recipe, "Check applied wavelength against skylines...");
01304             cpl_msg_indent_more();
01305 
01306             mean_rms = mos_distortions_rms(mapped_sky, lines, startwavelength,
01307                                            dispersion, 6, highres);
01308 
01309 
01310             cpl_msg_info(recipe, "Mean residual: %f", mean_rms);
01311 
01312             mean_rms = cpl_table_get_column_mean(idscoeff, "error");
01313 
01314             cpl_msg_info(recipe, "Mean model accuracy: %f pixel (%f A)",
01315                          mean_rms, mean_rms * dispersion);
01316         }
01317 
01318         save_header = cpl_propertylist_duplicate(header);
01319 
01320         cpl_propertylist_update_double(header, "CRPIX1", 1.0);
01321         cpl_propertylist_update_double(header, "CRPIX2", 1.0);
01322         cpl_propertylist_update_double(header, "CRVAL1", 
01323                                        startwavelength + dispersion/2);
01324         cpl_propertylist_update_double(header, "CRVAL2", 1.0);
01325         /* cpl_propertylist_update_double(header, "CDELT1", dispersion);
01326            cpl_propertylist_update_double(header, "CDELT2", 1.0); */
01327         cpl_propertylist_update_double(header, "CD1_1", dispersion);
01328         cpl_propertylist_update_double(header, "CD1_2", 0.0);
01329         cpl_propertylist_update_double(header, "CD2_1", 0.0);
01330         cpl_propertylist_update_double(header, "CD2_2", 1.0);
01331         cpl_propertylist_update_string(header, "CTYPE1", "LINEAR");
01332         cpl_propertylist_update_string(header, "CTYPE2", "PIXEL");
01333 
01334         if (time_normalise) {
01335             dummy = cpl_image_divide_scalar_create(mapped_sky, alltime);
01336 
01337             if (!j) {
01338                 if(dfs_save_image_null(frameset, parlist,
01339                                        mapped_science_sky_tag,
01340                                        recipe, version)) {
01341                     fors_pmos_science_exit(NULL);
01342                 }
01343             }
01344 
01345             if (dfs_save_image_ext(dummy, mapped_science_sky_tag, header)) {
01346                 fors_pmos_science_exit(NULL);
01347             }
01348 
01349             cpl_image_delete(dummy); dummy = NULL;
01350         }
01351         else {
01352 
01353             if (!j) {
01354                 if(dfs_save_image_null(frameset, parlist,
01355                                        mapped_science_sky_tag,
01356                                        recipe, version)) {
01357                     fors_pmos_science_exit(NULL);
01358                 }
01359             }
01360 
01361             if (dfs_save_image_ext(mapped_sky,
01362                                    mapped_science_sky_tag, header)) {
01363                 fors_pmos_science_exit(NULL);
01364             }
01365 
01366         }
01367 
01368         if (skymedian == 0 && skylocal == 0) {
01369             cpl_image_delete(mapped_sky); mapped_sky = NULL;
01370         }
01371 
01372         if (skylocal) {
01373 
01374             cpl_msg_indent_less();
01375 
01376             cpl_msg_info(recipe, "Local sky determination...");
01377             cpl_msg_indent_more();
01378             skymap = mos_subtract_sky(spectra, slits, polytraces, reference,
01379                                   startwavelength, endwavelength, dispersion);
01380 
01381             if (skymap) {
01382                 if (time_normalise)
01383                     cpl_image_divide_scalar(skymap, alltime);
01384 
01385                 if (!j) {
01386                     if(dfs_save_image_null(frameset, parlist,
01387                                            unmapped_sky_tag,
01388                                            recipe, version)) {
01389                         fors_pmos_science_exit(NULL);
01390                     }
01391                 }
01392 
01393                 if (dfs_save_image_ext(skymap, unmapped_sky_tag,
01394                                        save_header)) {
01395                     fors_pmos_science_exit(NULL);
01396                 }
01397 
01398                 cpl_image_delete(skymap); skymap = NULL;
01399 
01400                 if (!j) {
01401                     if(dfs_save_image_null(frameset, parlist,
01402                                            unmapped_science_tag,
01403                                            recipe, version)) {
01404                         fors_pmos_science_exit(NULL);
01405                     }
01406                 }
01407 
01408                 if (dfs_save_image_ext(spectra, unmapped_science_tag,
01409                                        save_header)) {
01410                     fors_pmos_science_exit(NULL);
01411                 }
01412 
01413 //                cpl_propertylist_delete(save_header); save_header = NULL;
01414 
01415                 if (cosmics) {
01416                     cpl_msg_info(recipe, "Removing cosmic rays...");
01417                     mos_clean_cosmics(spectra, gain, -1., -1.);
01418                 }
01419 
01420                 /*
01421                  * The spatially rectified image, that contained the sky,
01422                  * is replaced by a sky-subtracted spatially rectified image:
01423                  */
01424 
01425                 cpl_image_delete(smapped); smapped = NULL;
01426 
01427                 smapped = mos_spatial_calibration(spectra, slits, polytraces, 
01428                                                   reference, startwavelength, 
01429                                                   endwavelength, dispersion, 
01430                                                   flux, NULL);
01431             }
01432             else {
01433                 cpl_msg_warning(recipe, "Sky subtraction failure");
01434                 if (cosmics)
01435                     cpl_msg_warning(recipe, "Cosmic rays removal not performed!");
01436                 cosmics = skylocal = 0;
01437             }
01438         }
01439 
01440         cpl_image_delete(spectra); spectra = NULL;
01441         cpl_table_delete(polytraces); polytraces = NULL;
01442 
01443         if (skyalign >= 0) {
01444             save_header = dfs_load_header(frameset, science_tag, 0);
01445 
01446             if (!j) {
01447                 if(dfs_save_image_null(frameset, parlist,
01448                                        wavelength_map_sky_tag,
01449                                        recipe, version)) {
01450                     fors_pmos_science_exit(NULL);
01451                 }
01452             }
01453 
01454             if (dfs_save_image_ext(wavemap, wavelength_map_sky_tag,
01455                                    save_header)) {
01456                 fors_pmos_science_exit(NULL);
01457             }
01458 
01459 //            cpl_propertylist_delete(save_header); save_header = NULL;
01460         }
01461 
01462         cpl_image_delete(wavemap); wavemap = NULL;
01463 
01464         mapped = mos_wavelength_calibration(smapped, reference,
01465                                             startwavelength, endwavelength,
01466                                             dispersion, idscoeff, flux);
01467 
01468         cpl_image_delete(smapped); smapped = NULL;
01469 
01470         if (skyalign >= 0) {
01471             if (!j) {
01472                 if (dfs_save_table(frameset, idscoeff, disp_coeff_sky_tag,
01473                                    NULL, parlist, recipe, version)) {
01474                     fors_pmos_science_exit(NULL);
01475                 }
01476             }
01477         }
01478 
01479         if (skymedian) {
01480             cpl_msg_indent_less();
01481             cpl_msg_info(recipe, "Local sky determination...");
01482             cpl_msg_indent_more();
01483        
01484             skylocalmap = mos_sky_local_old(mapped, slits);       
01485             cpl_image_subtract(mapped, skylocalmap);
01486             cpl_image_delete(skylocalmap); skylocalmap = NULL;
01487         }
01488 
01489         if (skymedian || skylocal) {
01490 
01491             skylocalmap = cpl_image_subtract_create(mapped_sky, mapped);
01492 
01493             cpl_image_delete(mapped_sky); mapped_sky = NULL;
01494 
01495             if (time_normalise) {
01496                 dummy = cpl_image_divide_scalar_create(skylocalmap, alltime);
01497 
01498                 if (!j) {
01499                     if(dfs_save_image_null(frameset, parlist,
01500                                            mapped_sky_tag,
01501                                            recipe, version)) {
01502                         fors_pmos_science_exit(NULL);
01503                     }
01504                 }
01505 
01506                 if (dfs_save_image_ext(dummy, mapped_sky_tag,
01507                                        header)) {
01508                     fors_pmos_science_exit(NULL);
01509                 }
01510 
01511                 cpl_image_delete(dummy); dummy = NULL;
01512             }
01513             else {
01514                 if (!j) {
01515                     if(dfs_save_image_null(frameset, parlist,
01516                                            mapped_sky_tag,
01517                                            recipe, version)) {
01518                         fors_pmos_science_exit(NULL);
01519                     }
01520                 }
01521 
01522                 if (dfs_save_image_ext(skylocalmap, mapped_sky_tag,
01523                                        header)) {
01524                     fors_pmos_science_exit(NULL);
01525                 }
01526             }
01527 
01528             skylocalmaps[j] = skylocalmap;
01529 
01530             cpl_msg_indent_less();
01531             cpl_msg_info(recipe, "Object detection...");
01532             cpl_msg_indent_more();
01533 
01534             if (!j)
01535                 origslits = cpl_table_duplicate(slits);
01536 
01537             if (cosmics || nscience > 1) {
01538                 dummy = mos_detect_objects(mapped, slits, slit_margin, ext_radius, 
01539                                            cont_radius);
01540             }
01541             else {
01542                 mapped_cleaned = cpl_image_duplicate(mapped);
01543                 mos_clean_cosmics(mapped_cleaned, gain, -1., -1.);
01544                 dummy = mos_detect_objects(mapped_cleaned, slits, slit_margin, 
01545                                            ext_radius, cont_radius);
01546 
01547                 cpl_image_delete(mapped_cleaned); mapped_cleaned = NULL;
01548             }
01549 
01550             cpl_image_delete(dummy); dummy = NULL;
01551 
01552             if (check) {
01553 
01554 /* Old saving:
01555 
01556                 if (!j) {
01557                     if (dfs_save_table(frameset, slits, object_table_tag,
01558                                        NULL, parlist, recipe, version)) {
01559                         fors_pmos_science_exit(NULL);
01560                     }
01561                 } else {
01562                     if (dfs_save_table_ext(slits, object_table_tag, NULL)) {
01563                         fors_pmos_science_exit(NULL);
01564                     }
01565                 }
01566 
01567 End old saving */
01568 
01569 /* OldNew saving:
01570                 if (!j) {
01571                     if(dfs_save_image_null(frameset, parlist,
01572                                            object_table_tag,
01573                                            recipe, version)) {
01574                         fors_pmos_science_exit(NULL);
01575                     }
01576                 }
01577 
01578                 if (dfs_save_table_ext(slits, object_table_tag, NULL)) {
01579                     fors_pmos_science_exit(NULL);
01580                 }
01581 End oldnew saving */
01582             }
01583         }
01584 
01585         slitss[j]  = slits;
01586         mappeds[j] = mapped;
01587 
01588         cpl_msg_indent_less();
01589 
01590         cpl_propertylist_delete(header); header = NULL;
01591         cpl_propertylist_delete(save_header); save_header = NULL;
01592     }
01593 
01594     cpl_table_delete(offsets); offsets = NULL;
01595     cpl_table_delete(idscoeff); idscoeff = NULL;
01596 
01597     cpl_image_delete(norm_flat); norm_flat = NULL;
01598     cpl_vector_delete(lines); lines = NULL;
01599 
01600         
01601     cpl_msg_indent_less();
01602     cpl_msg_info(recipe, 
01603                  "Check object detection in both beams for all angles...");
01604     cpl_msg_indent_more();
01605 
01606     /* House keeping - selection of objects for which information required 
01607        for Stokes parameters computation is present */
01608     error = mos_object_intersect(slitss, origslits, nscience);
01609     if (error == CPL_ERROR_DATA_NOT_FOUND) {
01610         cpl_msg_warning(recipe, "No objects found: no Stokes "
01611                        "parameters to compute!");
01612         for (j = 0; j < nscience; j++)
01613             cpl_table_delete(slitss[j]);
01614         cpl_free(slitss);
01615         cpl_table_delete(origslits);
01616         return 0;
01617     } else if (error) {
01618         fors_pmos_science_exit("Problem in polarimetric object selection");
01619     }
01620 
01621     if (dfs_save_table(frameset, origslits, object_table_pol_tag,
01622                        NULL, parlist, recipe, version)) {
01623         fors_pmos_science_exit(NULL);
01624     }
01625 
01626     /*
01627      * Save also object tables per angle after intersection
01628      */
01629 
01630     for (j = 0; j < nscience; j++) {
01631         if (!j) {
01632             if(dfs_save_image_null(frameset, parlist,
01633                                    object_table_tag,
01634                                    recipe, version)) {
01635                 fors_pmos_science_exit(NULL);
01636             }
01637         }
01638 
01639         if (dfs_save_table_ext(slitss[j], object_table_tag, NULL)) {
01640             fors_pmos_science_exit(NULL);
01641         }
01642     }
01643 
01644     nobjs_per_slit = fors_get_nobjs_perslit(origslits);
01645 
01646     cpl_msg_indent_less();
01647     cpl_msg_info(recipe, "Object extraction...");
01648     cpl_msg_indent_more();
01649 
01650     for (j = 0; j < nscience; j++) {
01651         int k;
01652 
01653         header = dfs_load_header(frameset, science_tag, 0);
01654 
01655         for (k = 0; k < j; k ++) {
01656             cpl_propertylist_delete(header);
01657             header = dfs_load_header(frameset, NULL, 0);
01658         }
01659 
01660         cpl_propertylist_update_double(header, "CRPIX1", 1.0);
01661         cpl_propertylist_update_double(header, "CRPIX2", 1.0);
01662         cpl_propertylist_update_double(header, "CRVAL1", 
01663                                        startwavelength + dispersion/2);
01664         cpl_propertylist_update_double(header, "CRVAL2", 1.0);
01665         /* cpl_propertylist_update_double(header, "CDELT1", dispersion);
01666            cpl_propertylist_update_double(header, "CDELT2", 1.0); */
01667         cpl_propertylist_update_double(header, "CD1_1", dispersion);
01668         cpl_propertylist_update_double(header, "CD1_2", 0.0);
01669         cpl_propertylist_update_double(header, "CD2_1", 0.0);
01670         cpl_propertylist_update_double(header, "CD2_2", 1.0);
01671         cpl_propertylist_update_string(header, "CTYPE1", "LINEAR");
01672         cpl_propertylist_update_string(header, "CTYPE2", "PIXEL");
01673 
01674         if (skymedian || skylocal) {
01675 
01676             cpl_msg_info(recipe, "Extracting at angle %.2f (%d out of %d) ...",
01677                          angles[j], j + 1, nscience);
01678 
01679             images = mos_extract_objects(mappeds[j], skylocalmaps[j],
01680                                          origslits, 
01681                                          ext_mode, ron, gain, 1);
01682 
01683             cpl_image_delete(skylocalmaps[j]); skylocalmaps[j] = NULL;
01684 
01685             if (images) {
01686                 if (time_normalise)
01687                     cpl_image_divide_scalar(images[0], alltime);
01688 
01689                 if (!j) {
01690                     if(dfs_save_image_null(frameset, parlist,
01691                                            reduced_science_tag,
01692                                            recipe, version)) {
01693                         fors_pmos_science_exit(NULL);
01694                     }
01695                 }
01696 
01697                 if (dfs_save_image_ext(images[0], reduced_science_tag,
01698                                        header)) {
01699                     fors_pmos_science_exit(NULL);
01700                 }
01701 
01702                 reduceds[j] = images[0];
01703 //                cpl_image_delete(images[0]);
01704     
01705                 if (time_normalise)
01706                     cpl_image_divide_scalar(images[1], alltime);
01707 
01708                 if (!j) {
01709                     if(dfs_save_image_null(frameset, parlist,
01710                                            reduced_sky_tag,
01711                                            recipe, version)) {
01712                         fors_pmos_science_exit(NULL);
01713                     }
01714                 }
01715 
01716                 if (dfs_save_image_ext(images[1], reduced_sky_tag,
01717                                        header)) {
01718                     fors_pmos_science_exit(NULL);
01719                 }
01720                 cpl_image_delete(images[1]);
01721     
01722                 if (time_normalise)
01723                     cpl_image_divide_scalar(images[2], alltime);
01724 
01725                 if (!j) {
01726                     if(dfs_save_image_null(frameset, parlist,
01727                                            reduced_error_tag,
01728                                            recipe, version)) {
01729                         fors_pmos_science_exit(NULL);
01730                     }
01731                 }
01732 
01733                 if (dfs_save_image_ext(images[2], reduced_error_tag,
01734                                        header)) {
01735                     fors_pmos_science_exit(NULL);
01736                 }
01737 
01738                 rerrors[j] = images[2];
01739 
01740                 cpl_free(images);
01741             }
01742             else {
01743                 cpl_msg_warning(recipe, "No objects found: the products "
01744                                 "%s, %s, and %s are not created", 
01745                                 reduced_science_tag, reduced_sky_tag, 
01746                                 reduced_error_tag);
01747             }
01748 
01749         }
01750 
01751         if (skymedian || skylocal) {
01752             if (time_normalise)
01753                 cpl_image_divide_scalar(mappeds[j], alltime);
01754 
01755             if (!j) {
01756                 if(dfs_save_image_null(frameset, parlist,
01757                                        mapped_science_tag,
01758                                        recipe, version)) {
01759                     fors_pmos_science_exit(NULL);
01760                 }
01761             }
01762 
01763             if (dfs_save_image_ext(mappeds[j], mapped_science_tag,
01764                                    header)) {
01765                 fors_pmos_science_exit(NULL);
01766             }
01767         }
01768 
01769         cpl_image_delete(mappeds[j]); mappeds[j] = NULL;
01770         cpl_propertylist_delete(header); header = NULL;
01771 
01772     }
01773 
01774     cpl_table_delete(origslits);
01775 
01776     /* Stokes computation */
01777 
01778     nobjects = cpl_image_get_size_y(reduceds[0]) / 2;
01779     nx       = cpl_image_get_size_x(reduceds[0]);
01780 
01781     header = cpl_propertylist_new();
01782     cpl_propertylist_update_double(header, "CRPIX1", 1.0);
01783     cpl_propertylist_update_double(header, "CRPIX2", 1.0);
01784     cpl_propertylist_update_double(header, "CRVAL1", 
01785                                    startwavelength + dispersion/2);
01786     cpl_propertylist_update_double(header, "CRVAL2", 1.0);
01787     /* cpl_propertylist_update_double(header, "CDELT1", dispersion);
01788        cpl_propertylist_update_double(header, "CDELT2", 1.0); */
01789     cpl_propertylist_update_double(header, "CD1_1", dispersion);
01790     cpl_propertylist_update_double(header, "CD1_2", 0.0);
01791     cpl_propertylist_update_double(header, "CD2_1", 0.0);
01792     cpl_propertylist_update_double(header, "CD2_2", 1.0);
01793     cpl_propertylist_update_string(header, "CTYPE1", "LINEAR");
01794     cpl_propertylist_update_string(header, "CTYPE2", "PIXEL");
01795     
01796     if (circ) {
01797 
01798         cpl_image        *pv_im          = NULL;
01799         cpl_image        *pvnull_im      = NULL;
01800         cpl_image        *perr_im        = NULL;
01801 
01802         double           *p_v            = NULL;
01803         double           *p_vnull        = NULL;
01804         double           *perr           = NULL;
01805 
01806         double            mean_vnull;
01807 
01808         int p = -1;
01809         int total = 0;
01810 
01811         pv_im     = cpl_image_new(nx, nobjects, CPL_TYPE_DOUBLE);
01812         perr_im   = cpl_image_new(nx, nobjects, CPL_TYPE_DOUBLE);
01813 
01814         p_v     = cpl_image_get_data_double(pv_im);
01815         perr    = cpl_image_get_data_double(perr_im);
01816 
01817         if (nscience / 2 > 1) {
01818             pvnull_im = cpl_image_new(nx, nobjects, CPL_TYPE_DOUBLE);
01819             p_vnull = cpl_image_get_data_double(pvnull_im);
01820         }
01821 
01822         for (j = 0; j < nobjects; j++) {
01823             int k, m;
01824 
01825             double * ip_v,
01826                    * ip_vnull, * iperr;
01827 
01828             float * data;
01829             float * iff,  * ierr;
01830 
01831             ip_v = p_v + (nobjects - 1 - j) * nx;
01832 
01833             if (nscience / 2 > 1)
01834                 ip_vnull = p_vnull + (nobjects - 1 - j) * nx;
01835 
01836             iperr = perr + (nobjects - 1 - j) * nx;
01837 
01838             while (total < j + 1) {
01839                 p++;
01840                 total += nobjs_per_slit[p];
01841             }
01842 
01843             for (k = 0; k < nscience / 2; k++) {
01844                 float * if_o,  * if_e,  * ifdelta_o, * ifdelta_e;
01845 
01846                 int pos   = fors_find_angle_pos(angles, nscience, 180 * k - 45);
01847                 int pos_d = fors_find_angle_pos(angles, nscience, 180 * k + 45);
01848 
01849                 data = cpl_image_get_data_float(reduceds[pos]);
01850 
01851                 if_o = data + (2 * (nobjects - total) + nobjs_per_slit[p] 
01852                                + (total - j - 1)) * nx;
01853 
01854                 if_e = data + (2 * (nobjects - total) 
01855                                + (total - j - 1)) * nx;
01856 
01857 //                if_o = data + (2 * (nobjects - 1 - j) + 1) * nx;
01858 //                if_e = data +  2 * (nobjects - 1 - j)      * nx;
01859 
01860                 data = cpl_image_get_data_float(reduceds[pos_d]);
01861 
01862                 ifdelta_o = data + (2 * (nobjects - total) + nobjs_per_slit[p] 
01863                                + (total - j - 1)) * nx;
01864 
01865                 ifdelta_e = data + (2 * (nobjects - total) 
01866                                + (total - j - 1)) * nx;
01867 
01868 //                ifdelta_o = data + (2 * (nobjects - 1 - j) + 1) * nx;
01869 //                ifdelta_e = data +  2 * (nobjects - 1 - j)      * nx;
01870 
01871                 for (m = 0; m < nx; m++) {
01872 
01873                     double quantity = if_o[m] + if_e[m] == 0.0 ? 0.0 :
01874                         (if_o[m]      - if_e[m]     ) /
01875                         (if_o[m]      + if_e[m]     ) -
01876                         (ifdelta_o[m] - ifdelta_e[m]) /
01877                         (ifdelta_o[m] + ifdelta_e[m]);
01878 
01879                     quantity = isfinite(quantity) ? quantity : 0.0;
01880 
01881                     /* PQ map computation */
01882                     ip_v[m] += quantity * 0.5 / (nscience / 2);
01883 
01884                     /* PQnull map computation */
01885                     if (nscience / 2 > 1) {
01886                         if (k % 2)
01887                             ip_vnull[m] += quantity * 0.5 / (nscience / 2);
01888                         else
01889                             ip_vnull[m] -= quantity * 0.5 / (nscience / 2);
01890                     }
01891                 }
01892             }
01893 
01894             /* Error map */
01895             data = cpl_image_get_data_float(reduceds[0]);
01896             iff  = data +  2 * (nobjects - 1 - j)      * nx;
01897 
01898             data = cpl_image_get_data_float(rerrors[0]);
01899             ierr = data +  2 * (nobjects - 1 - j)      * nx;
01900 
01901             for (m = 0; m < nx; m++)
01902                 iperr[m] = iff[m] <= 0.0 ? 
01903                     0.0 : ierr[m] / iff[m] * 0.5 / sqrt (nscience / 2);
01904 
01905             if (nscience / 2 > 1) {
01906                 float * weights;
01907                 float   max, sum, sum2, imean;
01908 
01909                 int k;
01910 
01911                 /* QC on U NULL */
01912                 weights = cpl_malloc(sizeof(float) * nx);
01913 
01914                 max = 0.0;
01915                 for (k = 0; k < nx; k++) {
01916                     if (max < iff[k]) max = iff[k];
01917                 }
01918             
01919                 for (k = 0; k < nx; k++) {
01920                     weights[k] = iff[k] < 0.0 ? 
01921                         0.0 : iff[k] * iff[k] / (max * max);
01922                 }
01923             
01924                 sum  = 0.0;
01925                 sum2 = 0.0;
01926                 for (k = 0; k < nx; k++) {
01927                     sum  += weights[k] * ip_vnull[k];
01928                     sum2 += weights[k];
01929                 }
01930 
01931                 cpl_free(weights);
01932 
01933                 imean = sum / sum2;
01934 
01935                 mean_vnull += (imean - mean_vnull) / (j + 1.0);
01936             }
01937         }
01938 
01939         if (dfs_save_image(frameset, pv_im, reduced_v_tag, header, 
01940                            parlist, recipe, version))
01941             fors_pmos_science_exit(NULL);
01942 
01943         if (nscience / 2 > 1) {
01944             char * pipefile, * keyname;
01945             cpl_propertylist * qheader = dfs_load_header(frameset, science_tag, 0);
01946 
01947             cpl_propertylist_update_double(qheader, "CRPIX1", 1.0);
01948             cpl_propertylist_update_double(qheader, "CRPIX2", 1.0);
01949             cpl_propertylist_update_double(qheader, "CRVAL1", 
01950                                                   startwavelength + dispersion/2);
01951             cpl_propertylist_update_double(qheader, "CRVAL2", 1.0);
01952             /* cpl_propertylist_update_double(header, "CDELT1", dispersion);
01953                cpl_propertylist_update_double(header, "CDELT2", 1.0); */
01954             cpl_propertylist_update_double(qheader, "CD1_1", dispersion);
01955             cpl_propertylist_update_double(qheader, "CD1_2", 0.0);
01956             cpl_propertylist_update_double(qheader, "CD2_1", 0.0);
01957             cpl_propertylist_update_double(qheader, "CD2_2", 1.0);
01958             cpl_propertylist_update_string(qheader, "CTYPE1", "LINEAR");
01959             cpl_propertylist_update_string(qheader, "CTYPE2", "PIXEL");
01960 
01961             fors_qc_start_group(qheader, "2.0", instrume);
01962 
01963             /*
01964              * QC1 group header
01965              */
01966 
01967             if (fors_qc_write_string("PRO.CATG", reduced_nul_v_tag,
01968                                      "Product category", instrume))
01969                 fors_pmos_science_exit("Cannot write product category to "
01970                                      "QC log file");
01971 
01972             if (fors_qc_keyword_to_paf(qheader, "ESO DPR TYPE", NULL,
01973                                        "DPR type", instrume))
01974                 fors_pmos_science_exit("Missing keyword DPR TYPE in arc "
01975                                      "lamp header");
01976 
01977             if (fors_qc_keyword_to_paf(qheader, "ESO TPL ID", NULL,
01978                                        "Template", instrume))
01979                 fors_pmos_science_exit("Missing keyword TPL ID in arc "
01980                                      "lamp header");
01981 
01982             if (fors_qc_keyword_to_paf(qheader, "ESO INS GRIS1 NAME", NULL,
01983                                        "Grism name", instrume))
01984                 fors_pmos_science_exit("Missing keyword INS GRIS1 NAME in arc "
01985                                      "lamp header");
01986 
01987             if (fors_qc_keyword_to_paf(qheader, "ESO INS GRIS1 ID", NULL,
01988                                        "Grism identifier", instrume))
01989                 fors_pmos_science_exit("Missing keyword INS GRIS1 ID in arc "
01990                                      "lamp header");
01991 
01992             if (cpl_propertylist_has(qheader, "ESO INS FILT1 NAME"))
01993                 fors_qc_keyword_to_paf(qheader, "ESO INS FILT1 NAME", NULL,
01994                                        "Filter name", instrume);
01995 
01996             if (fors_qc_keyword_to_paf(qheader, "ESO INS COLL NAME", NULL,
01997                                        "Collimator name", instrume))
01998                 fors_pmos_science_exit("Missing keyword INS COLL NAME in arc "
01999                                      "lamp header");
02000 
02001             if (fors_qc_keyword_to_paf(qheader, "ESO DET CHIP1 ID", NULL,
02002                                        "Chip identifier", instrume))
02003                 fors_pmos_science_exit("Missing keyword DET CHIP1 ID in arc "
02004                                      "lamp header");
02005 
02006             if (fors_qc_keyword_to_paf(qheader, "ARCFILE", NULL,
02007                                        "Archive name of input data", 
02008                                        instrume))
02009                 fors_pmos_science_exit("Missing keyword ARCFILE in arc "
02010                                      "lamp header");
02011 
02012             pipefile = dfs_generate_filename_tfits(reduced_nul_v_tag);
02013             if (fors_qc_write_string("PIPEFILE", pipefile,
02014                                      "Pipeline product name", instrume))
02015                 fors_pmos_science_exit("Cannot write PIPEFILE to QC log file");
02016             cpl_free(pipefile); pipefile = NULL;
02017 
02018 
02019             /*
02020              * QC1 parameters
02021              */
02022 
02023             keyname = "QC.NULL.V.MEAN";
02024                     
02025             if (fors_qc_write_qc_double(qheader, mean_vnull,
02026                                         keyname, NULL,
02027                                         "Mean V null parameter",
02028                                         instrume)) {
02029                 fors_pmos_science_exit("Cannot write mean Q null parameter "
02030                                      "to QC log file");
02031             }
02032 
02033             fors_qc_end_group();
02034 
02035             if (dfs_save_image(frameset, pvnull_im, reduced_nul_v_tag, qheader, 
02036                                parlist, recipe, version))
02037                 fors_pmos_science_exit(NULL);
02038 
02039             cpl_propertylist_delete(qheader);
02040         }
02041 
02042         if (dfs_save_image(frameset, perr_im, reduced_error_v_tag, header, 
02043                            parlist, recipe, version))
02044             fors_pmos_science_exit(NULL);
02045 
02046         cpl_image_delete(pv_im);
02047         cpl_image_delete(pvnull_im);
02048         cpl_image_delete(perr_im);
02049     } else {
02050         cpl_image        *pq_im          = NULL;
02051         cpl_image        *pu_im          = NULL;
02052         cpl_image        *pl_im          = NULL;
02053 
02054         cpl_image        *pqnull_im      = NULL;
02055         cpl_image        *punull_im      = NULL;
02056 
02057         cpl_image        *pqerr_im        = NULL;
02058         cpl_image        *puerr_im        = NULL;
02059         cpl_image        *plerr_im        = NULL;
02060 
02061         cpl_image        *pang_im        = NULL;
02062         cpl_image        *pangerr_im        = NULL;
02063 
02064         double           *p_q            = NULL;
02065         double           *p_u            = NULL;
02066         double           *p_l            = NULL;
02067 
02068         double           *p_qnull        = NULL;
02069         double           *p_unull        = NULL;
02070 
02071         double           *pqerr           = NULL;
02072         double           *puerr           = NULL;
02073         double           *plerr           = NULL;
02074 
02075         double           *pang           = NULL;
02076         double           *pangerr           = NULL;
02077 
02078         int k, m;
02079 
02080         cpl_image * correct_im = cpl_image_new(nx, 1, CPL_TYPE_DOUBLE);
02081         double    * correct    = cpl_image_get_data_double(correct_im);
02082 
02083         double mean_unull, mean_qnull;
02084 
02085         int p = -1;
02086         int total = 0;
02087             
02088         pq_im = cpl_image_new(nx, nobjects, CPL_TYPE_DOUBLE);
02089         pu_im = cpl_image_new(nx, nobjects, CPL_TYPE_DOUBLE);
02090         pl_im = cpl_image_new(nx, nobjects, CPL_TYPE_DOUBLE);
02091 
02092         pqerr_im   = cpl_image_new(nx, nobjects, CPL_TYPE_DOUBLE);
02093         puerr_im   = cpl_image_new(nx, nobjects, CPL_TYPE_DOUBLE);
02094         plerr_im   = cpl_image_new(nx, nobjects, CPL_TYPE_DOUBLE);
02095 
02096         pang_im   = cpl_image_new(nx, nobjects, CPL_TYPE_DOUBLE);
02097         pangerr_im   = cpl_image_new(nx, nobjects, CPL_TYPE_DOUBLE);
02098 
02099         p_q = cpl_image_get_data_double(pq_im);
02100         p_u = cpl_image_get_data_double(pu_im);
02101         p_l = cpl_image_get_data_double(pl_im);
02102 
02103         if (nscience / 4 > 1) {
02104             pqnull_im = cpl_image_new(nx, nobjects, CPL_TYPE_DOUBLE);
02105             punull_im = cpl_image_new(nx, nobjects, CPL_TYPE_DOUBLE);
02106 
02107             p_qnull = cpl_image_get_data_double(pqnull_im);
02108             p_unull = cpl_image_get_data_double(punull_im);
02109         } else {
02110             cpl_msg_warning(cpl_func, 
02111                             "Not enough pairs to compute null parameters");
02112         }
02113 
02114         pqerr = cpl_image_get_data_double(pqerr_im);
02115         puerr = cpl_image_get_data_double(puerr_im);
02116         plerr = cpl_image_get_data_double(plerr_im);
02117 
02118         pang = cpl_image_get_data_double(pang_im);
02119         pangerr = cpl_image_get_data_double(pangerr_im);
02120 
02121         if (chromatism) {
02122             cpl_table * chrotbl = 
02123                 dfs_load_table(frameset, chrom_table_tag, 1);
02124 
02125             int      nrow   = cpl_table_get_nrow(chrotbl);
02126             float  * lambda = cpl_table_get_data_float(chrotbl, "lambda");
02127             float  * theta  = cpl_table_get_data_float(chrotbl, "eps_theta");
02128 
02129             for (j = 0; j < nx; j++) {
02130                 double c_wave = 
02131                     startwavelength + dispersion / 2 + j * dispersion;
02132             
02133                 int found = 0;
02134 
02135                 for (k = 0; k < nrow - 1; k++) {
02136                     if (lambda[k] <= c_wave && c_wave < lambda[k + 1]) {
02137                         found = 1;
02138                         break;
02139                     }
02140                 }
02141 
02142                 if (found) {
02143                     correct[j] = (theta [k + 1] - theta [k]) /
02144                                  (lambda[k + 1] - lambda[k]) *
02145                                  (c_wave        - lambda[k])   + theta[k];
02146                     correct[j] *= M_PI / 180;   /* Radians */
02147                 }
02148                 else if (j)
02149                     correct[j] = correct[j-1];
02150                 else
02151                     correct[j] = 0.0;
02152             }
02153 
02154             cpl_table_delete(chrotbl);
02155         }
02156 
02157         for (j = 0; j < nobjects; j++) {
02158             double * ip_q,     * ip_u, * ip_l, 
02159                 * ip_qnull, * ip_unull, * ipqerr, * ipuerr, * iplerr,
02160                 * ipang, * ipangerr;
02161 
02162             float * data;
02163             float * iffq,  * ierrq, * iffu, * ierru;
02164 
02165             int pos, pos_d;
02166 
02167             ip_q = p_q + (nobjects - 1 - j) * nx;
02168             ip_u = p_u + (nobjects - 1 - j) * nx;
02169             ip_l = p_l + (nobjects - 1 - j) * nx;
02170 
02171             if (nscience / 4 > 1) {
02172                 ip_qnull = p_qnull + (nobjects - 1 - j) * nx;
02173                 ip_unull = p_unull + (nobjects - 1 - j) * nx;
02174             }
02175 
02176             ipqerr = pqerr + (nobjects - 1 - j) * nx;
02177             ipuerr = puerr + (nobjects - 1 - j) * nx;
02178             iplerr = plerr + (nobjects - 1 - j) * nx;
02179 
02180             ipang = pang + (nobjects - 1 - j) * nx;
02181             ipangerr = pangerr + (nobjects - 1 - j) * nx;
02182 
02183             while (total < j + 1) {
02184                 p++;
02185                 total += nobjs_per_slit[p];
02186             }
02187 
02188             for (k = 0; k < nscience / 4; k++) {
02189                 float * if_o,  * if_e,  * ifdelta_o, * ifdelta_e;
02190 
02191                 /* First P_Q */
02192 
02193                 pos   = fors_find_angle_pos(angles, nscience, 90 * k);
02194                 pos_d = fors_find_angle_pos(angles, nscience, 90 * k + 45);
02195 
02196                 data = cpl_image_get_data_float(reduceds[pos]);
02197 
02198                 if_o = data + (2 * (nobjects - total) + nobjs_per_slit[p] 
02199                                + (total - j - 1)) * nx;
02200 
02201                 if_e = data + (2 * (nobjects - total) 
02202                                + (total - j - 1)) * nx;
02203 
02204 //                if_o = data + (2 * (nobjects - 1 - j) + 1) * nx;
02205 //                if_e = data +  2 * (nobjects - 1 - j)      * nx;
02206 
02207                 data = cpl_image_get_data_float(reduceds[pos_d]);
02208 
02209                 ifdelta_o = data + (2 * (nobjects - total) + nobjs_per_slit[p] 
02210                                + (total - j - 1)) * nx;
02211 
02212                 ifdelta_e = data + (2 * (nobjects - total) 
02213                                + (total - j - 1)) * nx;
02214 
02215 //                ifdelta_o = data + (2 * (nobjects - 1 - j) + 1) * nx;
02216 //                ifdelta_e = data +  2 * (nobjects - 1 - j)      * nx;
02217 
02218 
02219                 for (m = 0; m < nx; m++) {
02220 
02221                     double quantity = fabs(if_o[m] + if_e[m]) < FLT_MIN ? 0.0 :
02222                         (if_o[m]      - if_e[m]     ) /
02223                         (if_o[m]      + if_e[m]     ) -
02224                         (ifdelta_o[m] - ifdelta_e[m]) /
02225                         (ifdelta_o[m] + ifdelta_e[m]);
02226 
02227                     quantity = isfinite(quantity) ? quantity : 0.0;
02228 
02229                     /* PQ map computation */
02230                     ip_q[m] += quantity * 0.5 / (nscience / 4);
02231 
02232                     /* PQnull map computation */
02233                     if (nscience / 4 > 1) {
02234                         if (k % 2)
02235                             ip_qnull[m] += quantity * 0.5 / (nscience / 4);
02236                         else
02237                             ip_qnull[m] -= quantity * 0.5 / (nscience / 4);
02238                     }
02239                 }
02240 
02241                 /* Now P_U */
02242 
02243                 pos   = fors_find_angle_pos(angles, nscience, 90 * k + 22.5);
02244                 pos_d = fors_find_angle_pos(angles, nscience, 90 * k + 67.5);
02245 
02246                 data = cpl_image_get_data_float(reduceds[pos]);
02247 
02248                 if_o = data + (2 * (nobjects - total) + nobjs_per_slit[p] 
02249                                + (total - j - 1)) * nx;
02250 
02251                 if_e = data + (2 * (nobjects - total) 
02252                                + (total - j - 1)) * nx;
02253 
02254 //                if_o = data + (2 * (nobjects - 1 - j) + 1) * nx;
02255 //                if_e = data +  2 * (nobjects - 1 - j)      * nx;
02256 
02257                 data = cpl_image_get_data_float(reduceds[pos_d]);
02258 
02259                 ifdelta_o = data + (2 * (nobjects - total) + nobjs_per_slit[p] 
02260                                + (total - j - 1)) * nx;
02261 
02262                 ifdelta_e = data + (2 * (nobjects - total) 
02263                                + (total - j - 1)) * nx;
02264 
02265 //                ifdelta_o = data + (2 * (nobjects - 1 - j) + 1) * nx;
02266 //                ifdelta_e = data +  2 * (nobjects - 1 - j)      * nx;
02267 
02268                 for (m = 0; m < nx; m++) {
02269 
02270                     double quantity = fabs(if_o[m] + if_e[m]) < FLT_MIN ? 0.0 :
02271                         (if_o[m]      - if_e[m]     ) /
02272                         (if_o[m]      + if_e[m]     ) -
02273                         (ifdelta_o[m] - ifdelta_e[m]) /
02274                         (ifdelta_o[m] + ifdelta_e[m]);
02275 
02276                     quantity = isfinite(quantity) ? quantity : 0.0;
02277 
02278                     /* PU map computation */
02279                     ip_u[m] += quantity * 0.5 / (nscience / 4);
02280 
02281                     /* PUnull map computation */
02282                     if (nscience / 4 > 1) {
02283                         if (k % 2)
02284                             ip_unull[m] += quantity * 0.5 / (nscience / 4);
02285                         else
02286                             ip_unull[m] -= quantity * 0.5 / (nscience / 4);
02287                     }
02288                 }
02289             }
02290 
02291             /* Error map */
02292 
02293             pos   = fors_find_angle_pos(angles, nscience, 0.0);
02294 
02295             data = cpl_image_get_data_float(reduceds[pos]);
02296             iffq  = data +  2 * (nobjects - 1 - j)      * nx;
02297 
02298             data = cpl_image_get_data_float(rerrors[pos]);
02299             ierrq = data +  2 * (nobjects - 1 - j)      * nx;
02300             
02301             pos   = fors_find_angle_pos(angles, nscience, 22.5);
02302 
02303             data = cpl_image_get_data_float(reduceds[pos]);
02304             iffu  = data +  2 * (nobjects - 1 - j)      * nx;
02305 
02306             data = cpl_image_get_data_float(rerrors[pos]);
02307             ierru = data +  2 * (nobjects - 1 - j)      * nx;
02308 
02309             for (m = 0; m < nx; m++) {
02310 
02311                 double radicand; 
02312 
02313                 ipqerr[m] = iffq[m] <= 0.0 ? 
02314                     0.0 : ierrq[m] / iffq[m] * 0.5 / sqrt (nscience / 4);
02315 
02316                 ipuerr[m] = iffu[m] <= 0.0 ? 
02317                     0.0 : ierru[m] / iffu[m] * 0.5 / sqrt (nscience / 4);
02318 
02319                 iplerr[m] = CPL_MATH_SQRT1_2 * 0.5 * (ipqerr[m] + ipuerr[m]);
02320 
02321                 /* PL computation */
02322                 ip_l[m] = sqrt(ip_u[m] * ip_u[m] + ip_q[m] * ip_q[m]);
02323 
02324                 /* P angle computation */
02325 /*** Lander
02326                 ipang[m] = (ip_q[m] == 0.0 ?
02327                     (ip_u[m] > 0.0 ? 45.0 : 135.0)
02328                     : 0.5 * (atan2(ip_u[m], ip_q[m]) * 180 / M_PI + 
02329                              ((atan2(ip_u[m], ip_q[m]) > 0.0 ? 0.0 : 360.0))));
02330 ***/
02331 
02332 /*** Carlo */
02333                 if (fabs(ip_q[m]) < 0.00001) {
02334                     if (ip_u[m] > 0.0) {
02335                         ipang[m] = 45.0;
02336                     }
02337                     else {
02338                         ipang[m] = 135.0;
02339                     }
02340                 }
02341                 else {
02342                     ipang[m] = 0.5 * atan(ip_u[m] / ip_q[m]) * 180 / M_PI;
02343                     if (ip_q[m] > 0.0) {
02344                         if (ip_u[m] < 0.0) {
02345                             ipang[m] += 180.;
02346                         }
02347                     }
02348                     else {
02349                         ipang[m] += 90.;
02350                     }
02351                 }
02352 
02353 /***/
02354 
02355 
02356                 /* Error on the angle computation */
02357                 radicand = ip_q[m] * ip_q[m] * ipuerr[m] * ipuerr[m] + 
02358                            ip_u[m] * ip_u[m] * ipqerr[m] * ipqerr[m];
02359   
02360                 ipangerr[m] = ip_l[m] == 0.0 ? 0.0 :
02361                      sqrt(radicand) * 0.5 / (ip_l[m] * ip_l[m]) * 180 / M_PI;
02362 
02363                 /*
02364                  * This is a quick and dirty patch for FORS2 had the
02365                  * polariser mounted +90 with respect to FORS1. I must
02366                  * hardcode it, because there is no such info in the 
02367                  * header.
02368                  */
02369 
02370                 if (instrume[4] == '2') {
02371 
02372                     double w_rotation = - M_PI / 2; /* Rotation of polariser */
02373 
02374                     ipang[m] -= w_rotation * 180 / M_PI;
02375 
02376                     ip_q[m] = ip_q[m] * cos(2 * w_rotation)
02377                             + ip_u[m] * sin(2 * w_rotation);
02378 
02379                     ip_u[m] = ip_u[m] * cos(2 * w_rotation)
02380                             - ip_q[m] * sin(2 * w_rotation);
02381                 }
02382 
02383                 if (chromatism) {
02384                     ipang[m] -= correct[m] * 180 / M_PI;
02385 
02386                     ip_q[m] = ip_q[m] * cos(2 * correct[m])
02387                             + ip_u[m] * sin(2 * correct[m]);
02388     
02389                     ip_u[m] = ip_u[m] * cos(2 * correct[m])
02390                             - ip_q[m] * sin(2 * correct[m]);
02391                 }
02392             }
02393 
02394             if (nscience / 4 > 1) {
02395                 float * weights;
02396                 float   max, sum, sum2, imean;
02397 
02398                 int k;
02399 
02400                 /* QC on U NULL */
02401                 weights = cpl_malloc(sizeof(float) * nx);
02402 
02403                 max = 0.0;
02404                 for (k = 0; k < nx; k++) {
02405                     if (max < iffq[k]) max = iffq[k];
02406                 }
02407             
02408                 for (k = 0; k < nx; k++) {
02409                     weights[k] = iffq[k] < 0.0 ? 
02410                         0.0 : iffq[k] * iffq[k] / (max * max);
02411                 }
02412             
02413                 sum  = 0.0;
02414                 sum2 = 0.0;
02415                 for (k = 0; k < nx; k++) {
02416                     sum  += weights[k] * ip_qnull[k];
02417                     sum2 += weights[k];
02418                 }
02419 
02420                 cpl_free(weights);
02421 
02422                 imean = sum / sum2;
02423 
02424                 mean_qnull += (imean - mean_qnull) / (j + 1.0);
02425                   
02426                 /* QC on U NULL */
02427                 weights = cpl_malloc(sizeof(float) * nx);
02428             
02429                 max = 0.0;
02430                 for (k = 0; k < nx; k++) {
02431                     if (max < iffu[k]) max = iffu[k];
02432                 }
02433             
02434                 for (k = 0; k < nx; k++) {
02435                     weights[k] = iffu[k] < 0.0 ? 
02436                         0.0 : iffu[k] * iffu[k] / (max * max);
02437                 }
02438             
02439                 sum  = 0.0;
02440                 sum2 = 0.0;
02441                 for (k = 0; k < nx; k++) {
02442                     sum  += weights[k] * ip_unull[k];
02443                     sum2 += weights[k];
02444                 }
02445 
02446                 cpl_free(weights);
02447 
02448                 imean = sum / sum2;
02449 
02450                 mean_unull += (imean - mean_unull) / (j + 1.0);
02451             }
02452         }
02453 
02454         cpl_image_delete(correct_im);
02455 
02456         if (dfs_save_image(frameset, pq_im, reduced_q_tag, header, 
02457                            parlist, recipe, version))
02458             fors_pmos_science_exit(NULL);
02459 
02460         if (dfs_save_image(frameset, pu_im, reduced_u_tag, header, 
02461                            parlist, recipe, version))
02462             fors_pmos_science_exit(NULL);
02463 
02464         if (dfs_save_image(frameset, pl_im, reduced_l_tag, header, 
02465                            parlist, recipe, version))
02466             fors_pmos_science_exit(NULL);
02467 
02468         if (nscience / 4 > 1) {
02469             char *pipefile; 
02470             char *keyname;
02471             cpl_propertylist *qheader = dfs_load_header(frameset, 
02472                                                         science_tag, 0);
02473 
02474             cpl_propertylist_update_double(qheader, "CRPIX1", 1.0);
02475             cpl_propertylist_update_double(qheader, "CRPIX2", 1.0);
02476             cpl_propertylist_update_double(qheader, "CRVAL1", 
02477                                                   startwavelength + dispersion/2);
02478             cpl_propertylist_update_double(qheader, "CRVAL2", 1.0);
02479             /* cpl_propertylist_update_double(header, "CDELT1", dispersion);
02480                cpl_propertylist_update_double(header, "CDELT2", 1.0); */
02481             cpl_propertylist_update_double(qheader, "CD1_1", dispersion);
02482             cpl_propertylist_update_double(qheader, "CD1_2", 0.0);
02483             cpl_propertylist_update_double(qheader, "CD2_1", 0.0);
02484             cpl_propertylist_update_double(qheader, "CD2_2", 1.0);
02485             cpl_propertylist_update_string(qheader, "CTYPE1", "LINEAR");
02486             cpl_propertylist_update_string(qheader, "CTYPE2", "PIXEL");
02487 
02488             fors_qc_start_group(qheader, "2.0", instrume);
02489 
02490             /*
02491              * QC1 group header
02492              */
02493 
02494             if (fors_qc_write_string("PRO.CATG", reduced_nul_q_tag,
02495                                      "Product category", instrume))
02496                 fors_pmos_science_exit("Cannot write product category to "
02497                                      "QC log file");
02498 
02499             if (fors_qc_keyword_to_paf(qheader, "ESO DPR TYPE", NULL,
02500                                        "DPR type", instrume))
02501                 fors_pmos_science_exit("Missing keyword DPR TYPE in arc "
02502                                      "lamp header");
02503 
02504             if (fors_qc_keyword_to_paf(qheader, "ESO TPL ID", NULL,
02505                                        "Template", instrume))
02506                 fors_pmos_science_exit("Missing keyword TPL ID in arc "
02507                                      "lamp header");
02508 
02509             if (fors_qc_keyword_to_paf(qheader, "ESO INS GRIS1 NAME", NULL,
02510                                        "Grism name", instrume))
02511                 fors_pmos_science_exit("Missing keyword INS GRIS1 NAME in arc "
02512                                      "lamp header");
02513 
02514             if (fors_qc_keyword_to_paf(qheader, "ESO INS GRIS1 ID", NULL,
02515                                        "Grism identifier", instrume))
02516                 fors_pmos_science_exit("Missing keyword INS GRIS1 ID in arc "
02517                                      "lamp header");
02518 
02519             if (cpl_propertylist_has(qheader, "ESO INS FILT1 NAME"))
02520                 fors_qc_keyword_to_paf(qheader, "ESO INS FILT1 NAME", NULL,
02521                                        "Filter name", instrume);
02522 
02523             if (fors_qc_keyword_to_paf(qheader, "ESO INS COLL NAME", NULL,
02524                                        "Collimator name", instrume))
02525                 fors_pmos_science_exit("Missing keyword INS COLL NAME in arc "
02526                                      "lamp header");
02527 
02528             if (fors_qc_keyword_to_paf(qheader, "ESO DET CHIP1 ID", NULL,
02529                                        "Chip identifier", instrume))
02530                 fors_pmos_science_exit("Missing keyword DET CHIP1 ID in arc "
02531                                      "lamp header");
02532 
02533             if (fors_qc_keyword_to_paf(qheader, "ARCFILE", NULL,
02534                                        "Archive name of input data", 
02535                                        instrume))
02536                 fors_pmos_science_exit("Missing keyword ARCFILE in arc "
02537                                      "lamp header");
02538 
02539             pipefile = dfs_generate_filename_tfits(reduced_nul_q_tag);
02540             if (fors_qc_write_string("PIPEFILE", pipefile,
02541                                      "Pipeline product name", instrume))
02542                 fors_pmos_science_exit("Cannot write PIPEFILE to QC log file");
02543             cpl_free(pipefile); pipefile = NULL;
02544 
02545 
02546             /*
02547              * QC1 parameters
02548              */
02549 
02550             keyname = "QC.NULL.Q.MEAN";
02551                     
02552             if (fors_qc_write_qc_double(qheader, mean_qnull,
02553                                         keyname, NULL,
02554                                         "Mean Q null parameter",
02555                                         instrume)) {
02556                 fors_pmos_science_exit("Cannot write mean Q null parameter "
02557                                      "to QC log file");
02558             }
02559 
02560             fors_qc_end_group();
02561 
02562             if (dfs_save_image(frameset, pqnull_im, reduced_nul_q_tag, qheader, 
02563                                parlist, recipe, version))
02564                 fors_pmos_science_exit(NULL);
02565 
02566             cpl_propertylist_delete(qheader);
02567 
02568             qheader = dfs_load_header(frameset, science_tag, 0);
02569 
02570             cpl_propertylist_update_double(qheader, "CRPIX1", 1.0);
02571             cpl_propertylist_update_double(qheader, "CRPIX2", 1.0);
02572             cpl_propertylist_update_double(qheader, "CRVAL1", 
02573                                                   startwavelength + dispersion/2);
02574             cpl_propertylist_update_double(qheader, "CRVAL2", 1.0);
02575             /* cpl_propertylist_update_double(header, "CDELT1", dispersion);
02576                cpl_propertylist_update_double(header, "CDELT2", 1.0); */
02577             cpl_propertylist_update_double(qheader, "CD1_1", dispersion);
02578             cpl_propertylist_update_double(qheader, "CD1_2", 0.0);
02579             cpl_propertylist_update_double(qheader, "CD2_1", 0.0);
02580             cpl_propertylist_update_double(qheader, "CD2_2", 1.0);
02581             cpl_propertylist_update_string(qheader, "CTYPE1", "LINEAR");
02582             cpl_propertylist_update_string(qheader, "CTYPE2", "PIXEL");
02583 
02584             fors_qc_start_group(qheader, "2.0", instrume);
02585 
02586             /*
02587              * QC1 group header
02588              */
02589 
02590             if (fors_qc_write_string("PRO.CATG", reduced_nul_u_tag,
02591                                      "Product category", instrume))
02592                 fors_pmos_science_exit("Cannot write product category to "
02593                                      "QC log file");
02594 
02595             if (fors_qc_keyword_to_paf(qheader, "ESO DPR TYPE", NULL,
02596                                        "DPR type", instrume))
02597                 fors_pmos_science_exit("Missing keyword DPR TYPE in arc "
02598                                      "lamp header");
02599 
02600             if (fors_qc_keyword_to_paf(qheader, "ESO TPL ID", NULL,
02601                                        "Template", instrume))
02602                 fors_pmos_science_exit("Missing keyword TPL ID in arc "
02603                                      "lamp header");
02604 
02605             if (fors_qc_keyword_to_paf(qheader, "ESO INS GRIS1 NAME", NULL,
02606                                        "Grism name", instrume))
02607                 fors_pmos_science_exit("Missing keyword INS GRIS1 NAME in arc "
02608                                      "lamp header");
02609 
02610             if (fors_qc_keyword_to_paf(qheader, "ESO INS GRIS1 ID", NULL,
02611                                        "Grism identifier", instrume))
02612                 fors_pmos_science_exit("Missing keyword INS GRIS1 ID in arc "
02613                                      "lamp header");
02614 
02615             if (cpl_propertylist_has(qheader, "ESO INS FILT1 NAME"))
02616                 fors_qc_keyword_to_paf(qheader, "ESO INS FILT1 NAME", NULL,
02617                                        "Filter name", instrume);
02618 
02619             if (fors_qc_keyword_to_paf(qheader, "ESO INS COLL NAME", NULL,
02620                                        "Collimator name", instrume))
02621                 fors_pmos_science_exit("Missing keyword INS COLL NAME in arc "
02622                                      "lamp header");
02623 
02624             if (fors_qc_keyword_to_paf(qheader, "ESO DET CHIP1 ID", NULL,
02625                                        "Chip identifier", instrume))
02626                 fors_pmos_science_exit("Missing keyword DET CHIP1 ID in arc "
02627                                      "lamp header");
02628 
02629             if (fors_qc_keyword_to_paf(qheader, "ARCFILE", NULL,
02630                                        "Archive name of input data", 
02631                                        instrume))
02632                 fors_pmos_science_exit("Missing keyword ARCFILE in arc "
02633                                      "lamp header");
02634 
02635             pipefile = dfs_generate_filename_tfits(reduced_nul_u_tag);
02636             if (fors_qc_write_string("PIPEFILE", pipefile,
02637                                      "Pipeline product name", instrume))
02638                 fors_pmos_science_exit("Cannot write PIPEFILE to QC log file");
02639             cpl_free(pipefile); pipefile = NULL;
02640 
02641 
02642             /*
02643              * QC1 parameters
02644              */
02645 
02646             keyname = "QC.NULL.U.MEAN";
02647                     
02648             if (fors_qc_write_qc_double(qheader, mean_unull,
02649                                         keyname, NULL,
02650                                         "Mean U null parameter",
02651                                         instrume)) {
02652                 fors_pmos_science_exit("Cannot write mean U null parameter "
02653                                      "to QC log file");
02654             }
02655 
02656             fors_qc_end_group();
02657 
02658             if (dfs_save_image(frameset, punull_im, reduced_nul_u_tag, qheader, 
02659                                parlist, recipe, version))
02660                 fors_pmos_science_exit(NULL);
02661 
02662             cpl_propertylist_delete(qheader);
02663         }
02664 
02665         if (dfs_save_image(frameset, pqerr_im, reduced_error_q_tag, header, 
02666                            parlist, recipe, version))
02667             fors_pmos_science_exit(NULL);
02668 
02669         if (dfs_save_image(frameset, puerr_im, reduced_error_u_tag, header, 
02670                            parlist, recipe, version))
02671             fors_pmos_science_exit(NULL);
02672 
02673         if (dfs_save_image(frameset, plerr_im, reduced_error_l_tag, header, 
02674                            parlist, recipe, version))
02675             fors_pmos_science_exit(NULL);
02676 
02677         if (dfs_save_image(frameset, pang_im, reduced_angle_tag, header, 
02678                            parlist, recipe, version))
02679             fors_pmos_science_exit(NULL);
02680 
02681         if (dfs_save_image(frameset, pangerr_im, reduced_error_angle_tag, 
02682                            header, parlist, recipe, version))
02683             fors_pmos_science_exit(NULL);
02684 
02685         cpl_image_delete(pq_im);
02686         cpl_image_delete(pu_im);
02687         cpl_image_delete(pl_im);
02688 
02689         cpl_image_delete(pqnull_im);
02690         cpl_image_delete(punull_im);
02691 
02692         cpl_image_delete(pqerr_im);
02693         cpl_image_delete(puerr_im);
02694         cpl_image_delete(plerr_im);
02695         cpl_image_delete(pang_im);
02696         cpl_image_delete(pangerr_im);
02697     }
02698 
02699     cpl_propertylist_delete(header);
02700 
02701     /* End of Stokes computation */
02702 
02703     for (j = 0; j < nscience; j++) {
02704         cpl_image_delete(reduceds[j]);
02705         cpl_image_delete(rerrors[j]);
02706         cpl_table_delete(slitss[j]);
02707         cpl_image_delete(mappeds[j]);
02708     }
02709 
02710     cpl_free(reduceds);
02711     cpl_free(rerrors);
02712     cpl_free(slitss);
02713     cpl_free(mappeds);
02714 
02715     cpl_free(instrume); instrume = NULL;
02716 
02717     cpl_free(skylocalmaps);
02718 
02719     cpl_free(nobjs_per_slit);
02720 
02721     if (cpl_error_get_code()) {
02722         cpl_msg_error(cpl_error_get_where(), cpl_error_get_message());
02723         fors_pmos_science_exit(NULL);
02724     }
02725     else 
02726         return 0;
02727 }
02728 
02729 /*----------------------------------------------------------------------------*/
02740 /*----------------------------------------------------------------------------*/
02741 static float * fors_check_angles(cpl_frameset * frameset,
02742                                  int pmos, const char *tag, int * circ)
02743 {
02744     float     *angles  = NULL;
02745     cpl_frame *c_frame = NULL;
02746     char      *ret_id  = NULL;
02747 
02748     int i = 0;
02749 
02750     angles = cpl_malloc(sizeof(float) * pmos);
02751 
02752     for (c_frame  = cpl_frameset_find(frameset, tag);
02753          c_frame != NULL; c_frame = cpl_frameset_find(frameset, NULL)) {
02754 
02755         cpl_propertylist * header =
02756             cpl_propertylist_load(cpl_frame_get_filename(c_frame), 0);
02757         
02758         if (!ret_id) {
02759             ret_id = cpl_strdup(cpl_propertylist_get_string(header, 
02760                                                         "ESO INS OPTI4 ID"));
02761 
02762             if (ret_id[1] != '5' && ret_id[1] != '4') {
02763                 cpl_msg_error(cpl_func, 
02764                               "Unknown retarder plate id: %s", ret_id);
02765                 return NULL;
02766             }
02767         } else {
02768             char * c_ret_id = (char *)
02769                 cpl_propertylist_get_string(header, "ESO INS OPTI4 ID");
02770             if (ret_id[1] != c_ret_id[1]) {
02771                 cpl_msg_error(cpl_func, "Input frames are not from the same "
02772                               "retarder plate");
02773                 return NULL;
02774             }
02775         }
02776         
02777         if (ret_id[1] == '5') {  /* Linear polarimetry */
02778             angles[i] = (float)
02779                 cpl_propertylist_get_double(header, "ESO INS RETA2 ROT");
02780             *circ = 0;
02781         } else {                 /* Circular polarimetry */
02782             angles[i] = (float)
02783                 cpl_propertylist_get_double(header, "ESO INS RETA4 ROT");
02784             *circ = 1;
02785         }
02786 
02787         cpl_propertylist_delete(header);
02788         i++;
02789     }
02790 
02791     cpl_free(ret_id);
02792 
02793     if (*circ) {
02794         if (pmos != 2 && pmos != 4) {
02795             cpl_msg_error(cpl_func, "Wrong angle configuration: %d angles "
02796                           "found, but either 2 or 4 are required for "
02797                           "circular polarization measurements!", pmos);
02798             return NULL;
02799         }
02800     } else {
02801         if (pmos != 4 && pmos != 8 && pmos != 16) {
02802             cpl_msg_error(cpl_func, "Wrong angle configuration: %d angles "
02803                           "found, but either 4, 8, or 16 are required for "
02804                           "linear polarization measurements!", pmos);
02805             return NULL;
02806         }
02807     }
02808     
02809     /* Check completeness */
02810 
02811     if (*circ) {
02812         for (i = 0; i < pmos; i++) {
02813             if (fors_find_angle_pos(angles, pmos, 90.0 * i - 45.0) < 0) {
02814                 const char *cangles;
02815                 switch (pmos) {
02816                 case 2: cangles  = "-45.0, 45.0"; break;
02817                 case 4: cangles  = "-45.0, 45.0, 135.0, 225.0"; break;
02818                 default: assert(0);
02819                 }  
02820 
02821                 cpl_msg_error(cpl_func, "Wrong angle configuration: missing "
02822                               "angle %.2f. All angles %s must be provided.",
02823                               angles[i], cangles);
02824                 return NULL;
02825             }
02826         }
02827     }
02828     else {
02829         for (i = 0; i < pmos; i++) {
02830             if (fors_find_angle_pos(angles, pmos, 22.5 * i) < 0) {
02831                 const char *cangles;
02832                 switch (pmos) {
02833                 case 4: cangles  = "0.0, 22.5, 45.0, 67.5"; break;
02834                 case 8: cangles  = "0.0, 22.5, 45.0, 67.5, "
02835                                    "90.0, 112.5, 135.0, 157.5"; break;
02836                 case 16: cangles = "0.0, 22.5, 45.0, 67.5, "
02837                                    "90.0, 112.5, 135.0, 157.5, "
02838                                    "180.0, 202.5, 225.0, 247.5, "
02839                                    "270.0, 292.5, 315.0, 337.5"; break;
02840                 default: assert(0);
02841                 }  
02842 
02843                 cpl_msg_error(cpl_func, "Wrong angle configuration: missing "
02844                               "angle %.2f. All angles %s must be provided.",
02845                               angles[i], cangles);
02846                 return NULL;
02847             }
02848         }
02849     }
02850 
02851     return angles;
02852 }
02853 
02854 /*----------------------------------------------------------------------------*/
02862 /*----------------------------------------------------------------------------*/
02863 static int
02864 fors_find_angle_pos(float * angles, int nangles, float angle)
02865 {
02866     int i, match = 0;
02867 
02868     for (i = 0; i < nangles; i++) {
02869         if (fabs(angles[i]         - angle) < 1.0 || 
02870             fabs(angles[i] - 360.0 - angle) < 1.0) {
02871             match = 1;
02872             break;
02873         }
02874     }
02875 
02876     return match ? i : -1;
02877 }
02878 

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