vimos_science.c

00001 /* $Id: vimos_science.c,v 1.6 2008/08/05 14:21:35 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: 2008/08/05 14:21:35 $
00024  * $Revision: 1.6 $
00025  * $Name:  $
00026  */
00027 
00028 #ifdef HAVE_CONFIG_H
00029 #include <config.h>
00030 #endif
00031 
00032 #include <math.h>
00033 #include <cpl.h>
00034 #include <moses.h>
00035 #include <fors_dfs.h>
00036 
00037 static int vimos_science_create(cpl_plugin *);
00038 static int vimos_science_exec(cpl_plugin *);
00039 static int vimos_science_destroy(cpl_plugin *);
00040 static int vimos_science(cpl_parameterlist *, cpl_frameset *);
00041 
00042 static char vimos_science_description[] =
00043 "This recipe is used to reduce scientific spectra using the extraction\n"
00044 "mask and the products created by the recipe vimos_calib. The spectra are\n"
00045 "bias subtracted, flat fielded (if a normalised flat field is specified)\n"
00046 "and remapped eliminating the optical distortions. The wavelength calibration\n"
00047 "can be optionally upgraded using a number of sky lines: if no sky lines\n"
00048 "catalog of wavelengths is specified, an internal one is used instead.\n"
00049 "If the alignment to the sky lines is performed, the input dispersion\n"
00050 "coefficients table is upgraded and saved to disk, and a new CCD wavelengths\n"
00051 "map is created. A grism table (typically depending on the grism used)\n"
00052 "may also be specified: this table contains a default recipe parameter\n" 
00053 "setting to control the way spectra are extracted for a specific instrument\n"
00054 "mode, as it is used for automatic run of the pipeline on Paranal and in\n" 
00055 "Garching. If this table is specified, it will modify the default recipe\n" 
00056 "parameter setting, with the exception of those parameters which have been\n" 
00057 "explicitly modifyed on the command line. If a grism table is not specified,\n"
00058 "the input recipe parameters values will always be read from the command\n" 
00059 "line, or from an esorex configuration file if present, or from their\n" 
00060 "generic default values (that are rarely meaningful).\n" 
00061 "MOS_SLIT_LOCATION and MOS_CURV_COEFF tables are not (yet) expected for\n"
00062 "long-slit-like data.\n\n"
00063 "Input files:\n\n"
00064 "  DO category:                Type:       Explanation:         Required:\n"
00065 "  MOS_SCIENCE                 Raw         Scientific exposure     Y\n"
00066 "  or MOS_STANDARD             Raw         Standard star exposure  Y\n"
00067 "\n"
00068 "  MASTER_BIAS                 Calib       Master bias             Y\n"
00069 "  SKY_LINE_CATALOG            Calib       Sky lines catalog       .\n"
00070 "  MOS_MASTER_SCREEN_FLAT      Calib       Normalised flat field   .\n"
00071 "  MOS_DISP_COEFF              Calib       Inverse dispersion      Y\n"
00072 "  MOS_CURV_COEFF              Calib       Spectral curvature      Y\n"
00073 "  MOS_SLIT_LOCATION           Calib       Slits positions table   Y\n"
00074 "  GRISM_TABLE                 Calib       Grism table             .\n\n"
00075 "Output files:\n\n"
00076 "  DO category:                Data type:  Explanation:\n"
00077 "  MOS_SCIENCE_REDUCED         FITS image  Extracted scientific spectra\n"
00078 "  or MOS_STANDARD_REDUCED     FITS image  Extracted standard star spectrum\n"
00079 "  MOS_SKY_REDUCED             FITS image  Extracted sky spectra\n"
00080 "  MOS_ERROR_REDUCED           FITS image  Errors on extracted spectra\n"
00081 "\n"
00082 "  MOS_UNMAPPED_SCIENCE        FITS image  Sky subtracted scientific spectra\n"
00083 "  or MOS_UNMAPPED_STANDARD    FITS image  Sky subtracted standard spectrum\n"
00084 "\n"
00085 "  MOS_SCIENCE_EXTRACTED       FITS image  Rectified scientific spectra\n"
00086 "  or MOS_STANDARD_EXTRACTED   FITS image  Rectified standard star spectrum\n"
00087 "\n"
00088 "  MOS_SCIENCE_SKY_EXTRACTED   FITS image  Rectified science spectra with sky\n"
00089 "or MOS_STANDARD_SKY_EXTRACTED FITS image  Rectified std spectrum with sky\n"
00090 "\n"
00091 "  MOS_SCIENCE_SKY             FITS image  Rectified sky spectra\n"
00092 "  MOS_UNMAPPED_SKY            FITS image  Sky on CCD\n"
00093 "  MOS_GLOBAL_SKY_SPECTRUM     FITS table  Global sky spectrum\n"
00094 "  OBJECT_TABLE                FITS table  Positions of detected objects\n"
00095 "\n"
00096 "  Only if the sky-alignment of the wavelength solution is requested:\n"
00097 "  MOS_SKYLINES_OFFSETS_LONG   FITS table  Sky lines offsets (LSS-like data)\n"
00098 "or MOS_SKYLINES_OFFSETS_SLIT  FITS table  Sky lines offsets (MOS-like data)\n"
00099 "  MOS_DISP_COEFF_SKY          FITS table  Upgraded dispersion coefficients\n"
00100 "  MOS_WAVELENGTH_MAP_SKY      FITS image  Upgraded wavelength map\n\n";
00101 
00102 #define vimos_science_exit(message)           \
00103 {                                             \
00104 if (message) cpl_msg_error(recipe, message);  \
00105 cpl_free(exptime);                            \
00106 cpl_image_delete(dummy);                      \
00107 cpl_image_delete(mapped);                     \
00108 cpl_image_delete(mapped_sky);                 \
00109 cpl_image_delete(mapped_cleaned);             \
00110 cpl_image_delete(skylocalmap);                \
00111 cpl_image_delete(skymap);                     \
00112 cpl_image_delete(smapped);                    \
00113 cpl_table_delete(offsets);                    \
00114 cpl_table_delete(sky);                        \
00115 cpl_image_delete(bias);                       \
00116 cpl_image_delete(spectra);                    \
00117 cpl_image_delete(coordinate);                 \
00118 cpl_image_delete(norm_flat);                  \
00119 cpl_image_delete(rainbow);                    \
00120 cpl_image_delete(rectified);                  \
00121 cpl_image_delete(wavemap);                    \
00122 cpl_propertylist_delete(header);              \
00123 cpl_propertylist_delete(save_header);         \
00124 cpl_table_delete(grism_table);                \
00125 cpl_table_delete(idscoeff);                   \
00126 cpl_table_delete(maskslits);                  \
00127 cpl_table_delete(overscans);                  \
00128 cpl_table_delete(polytraces);                 \
00129 cpl_table_delete(slits);                      \
00130 cpl_table_delete(wavelengths);                \
00131 cpl_vector_delete(lines);                     \
00132 cpl_msg_indent_less();                        \
00133 return -1;                                    \
00134 }
00135 
00136 
00137 #define vimos_science_exit_memcheck(message)   \
00138 {                                             \
00139 if (message) cpl_msg_info(recipe, message);   \
00140 cpl_free(exptime);                            \
00141 cpl_image_delete(dummy);                      \
00142 cpl_image_delete(mapped);                     \
00143 cpl_image_delete(mapped_cleaned);             \
00144 cpl_image_delete(mapped_sky);                 \
00145 cpl_image_delete(skylocalmap);                \
00146 cpl_image_delete(skymap);                     \
00147 cpl_image_delete(smapped);                    \
00148 cpl_table_delete(offsets);                    \
00149 cpl_table_delete(sky);                        \
00150 cpl_image_delete(bias);                       \
00151 cpl_image_delete(spectra);                    \
00152 cpl_image_delete(coordinate);                 \
00153 cpl_image_delete(norm_flat);                  \
00154 cpl_image_delete(rainbow);                    \
00155 cpl_image_delete(rectified);                  \
00156 cpl_image_delete(wavemap);                    \
00157 cpl_propertylist_delete(header);              \
00158 cpl_propertylist_delete(save_header);         \
00159 cpl_table_delete(grism_table);                \
00160 cpl_table_delete(idscoeff);                   \
00161 cpl_table_delete(maskslits);                  \
00162 cpl_table_delete(overscans);                  \
00163 cpl_table_delete(polytraces);                 \
00164 cpl_table_delete(slits);                      \
00165 cpl_table_delete(wavelengths);                \
00166 cpl_vector_delete(lines);                     \
00167 cpl_msg_indent_less();                        \
00168 return 0;                                     \
00169 }
00170 
00171 
00183 int cpl_plugin_get_info(cpl_pluginlist *list)
00184 {
00185     cpl_recipe *recipe = cpl_calloc(1, sizeof *recipe );
00186     cpl_plugin *plugin = &recipe->interface;
00187 
00188     cpl_plugin_init(plugin,
00189                     CPL_PLUGIN_API,
00190                     FORS_BINARY_VERSION,
00191                     CPL_PLUGIN_TYPE_RECIPE,
00192                     "vimos_science",
00193                     "Extraction of scientific spectra",
00194                     vimos_science_description,
00195                     "Carlo Izzo",
00196                     PACKAGE_BUGREPORT,
00197     "This file is currently part of the FORS Instrument Pipeline\n"
00198     "Copyright (C) 2002-2006 European Southern Observatory\n\n"
00199     "This program is free software; you can redistribute it and/or modify\n"
00200     "it under the terms of the GNU General Public License as published by\n"
00201     "the Free Software Foundation; either version 2 of the License, or\n"
00202     "(at your option) any later version.\n\n"
00203     "This program is distributed in the hope that it will be useful,\n"
00204     "but WITHOUT ANY WARRANTY; without even the implied warranty of\n"
00205     "MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the\n"
00206     "GNU General Public License for more details.\n\n"
00207     "You should have received a copy of the GNU General Public License\n"
00208     "along with this program; if not, write to the Free Software Foundation,\n"
00209     "Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA\n",
00210                     vimos_science_create,
00211                     vimos_science_exec,
00212                     vimos_science_destroy);
00213 
00214     cpl_pluginlist_append(list, plugin);
00215     
00216     return 0;
00217 }
00218 
00219 
00230 static int vimos_science_create(cpl_plugin *plugin)
00231 {
00232     cpl_recipe    *recipe;
00233     cpl_parameter *p;
00234 
00235 
00236     /* 
00237      * Check that the plugin is part of a valid recipe 
00238      */
00239 
00240     if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE) 
00241         recipe = (cpl_recipe *)plugin;
00242     else 
00243         return -1;
00244 
00245     /* 
00246      * Create the parameters list in the cpl_recipe object 
00247      */
00248 
00249     recipe->parameters = cpl_parameterlist_new(); 
00250 
00251 
00252     /*
00253      * Dispersion
00254      */
00255 
00256     p = cpl_parameter_new_value("fors.vimos_science.dispersion",
00257                                 CPL_TYPE_DOUBLE,
00258                                 "Resampling step (Angstrom/pixel)",
00259                                 "fors.vimos_science",
00260                                 0.0);
00261     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "dispersion");
00262     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00263     cpl_parameterlist_append(recipe->parameters, p);
00264 
00265     /*
00266      * Sky lines alignment
00267      */
00268 
00269     p = cpl_parameter_new_value("fors.vimos_science.skyalign",
00270                                 CPL_TYPE_INT,
00271                                 "Polynomial order for sky lines alignment, "
00272                                 "or -1 to avoid alignment",
00273                                 "fors.vimos_science",
00274                                 0);
00275     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "skyalign");
00276     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00277     cpl_parameterlist_append(recipe->parameters, p);
00278 
00279     /*
00280      * Line catalog table column containing the sky reference wavelengths
00281      */
00282 
00283     p = cpl_parameter_new_value("fors.vimos_science.wcolumn",
00284                                 CPL_TYPE_STRING,
00285                                 "Name of sky line catalog table column "
00286                                 "with wavelengths",
00287                                 "fors.vimos_science",
00288                                 "WLEN");
00289     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "wcolumn");
00290     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00291     cpl_parameterlist_append(recipe->parameters, p);
00292 
00293     /*
00294      * Start wavelength for spectral extraction
00295      */
00296 
00297     p = cpl_parameter_new_value("fors.vimos_science.startwavelength",
00298                                 CPL_TYPE_DOUBLE,
00299                                 "Start wavelength in spectral extraction",
00300                                 "fors.vimos_science",
00301                                 0.0);
00302     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "startwavelength");
00303     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00304     cpl_parameterlist_append(recipe->parameters, p);
00305 
00306     /*
00307      * End wavelength for spectral extraction
00308      */
00309 
00310     p = cpl_parameter_new_value("fors.vimos_science.endwavelength",
00311                                 CPL_TYPE_DOUBLE,
00312                                 "End wavelength in spectral extraction",
00313                                 "fors.vimos_science",
00314                                 0.0);
00315     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "endwavelength");
00316     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00317     cpl_parameterlist_append(recipe->parameters, p);
00318 
00319     /*
00320      * Reference wavelength for wavelength calibration
00321      */
00322 
00323     p = cpl_parameter_new_value("fors.vimos_science.reference",
00324                                 CPL_TYPE_DOUBLE,
00325                                 "Reference wavelength for calibration",
00326                                 "fors.vimos_science",
00327                                 0.0);
00328     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "reference");
00329     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00330     cpl_parameterlist_append(recipe->parameters, p);
00331 
00332     /*
00333      * Flux conservation
00334      */
00335 
00336     p = cpl_parameter_new_value("fors.vimos_science.flux",
00337                                 CPL_TYPE_BOOL,
00338                                 "Apply flux conservation",
00339                                 "fors.vimos_science",
00340                                 TRUE);
00341     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "flux");
00342     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00343     cpl_parameterlist_append(recipe->parameters, p);
00344 
00345     /*
00346      * Apply flat field
00347      */
00348 
00349     p = cpl_parameter_new_value("fors.vimos_science.flatfield",
00350                                 CPL_TYPE_BOOL,
00351                                 "Apply flat field",
00352                                 "fors.vimos_science",
00353                                 TRUE);
00354     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "flatfield");
00355     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00356     cpl_parameterlist_append(recipe->parameters, p);
00357 
00358     /*
00359      * Global sky subtraction
00360      */
00361 
00362     p = cpl_parameter_new_value("fors.vimos_science.skyglobal",
00363                                 CPL_TYPE_BOOL,
00364                                 "Subtract global sky spectrum from CCD",
00365                                 "fors.vimos_science",
00366                                 FALSE);
00367     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "skyglobal");
00368     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00369     cpl_parameterlist_append(recipe->parameters, p);
00370 
00371     /*
00372      * Local sky subtraction on extracted spectra
00373      */
00374 
00375     p = cpl_parameter_new_value("fors.vimos_science.skymedian",
00376                                 CPL_TYPE_BOOL,
00377                                 "Sky subtraction from extracted slit spectra",
00378                                 "fors.vimos_science",
00379                                 FALSE);
00380     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "skymedian");
00381     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00382     cpl_parameterlist_append(recipe->parameters, p);
00383 
00384     /*
00385      * Local sky subtraction on CCD spectra
00386      */
00387 
00388     p = cpl_parameter_new_value("fors.vimos_science.skylocal",
00389                                 CPL_TYPE_BOOL,
00390                                 "Sky subtraction from CCD slit spectra",
00391                                 "fors.vimos_science",
00392                                 TRUE);
00393     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "skylocal");
00394     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00395     cpl_parameterlist_append(recipe->parameters, p);
00396 
00397     /*
00398      * Cosmic rays removal
00399      */
00400 
00401     p = cpl_parameter_new_value("fors.vimos_science.cosmics",
00402                                 CPL_TYPE_BOOL,
00403                                 "Eliminate cosmic rays hits (only if global "
00404                                 "sky subtraction is also requested)",
00405                                 "fors.vimos_science",
00406                                 FALSE);
00407     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "cosmics");
00408     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00409     cpl_parameterlist_append(recipe->parameters, p);
00410 
00411     /*
00412      * Slit margin
00413      */
00414 
00415     p = cpl_parameter_new_value("fors.vimos_science.slit_margin",
00416                                 CPL_TYPE_INT,
00417                                 "Number of pixels to exclude at each slit "
00418                                 "in object detection and extraction",
00419                                 "fors.vimos_science",
00420                                 3);
00421     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "slit_margin");
00422     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00423     cpl_parameterlist_append(recipe->parameters, p);
00424 
00425     /*
00426      * Extraction radius
00427      */
00428 
00429     p = cpl_parameter_new_value("fors.vimos_science.ext_radius",
00430                                 CPL_TYPE_INT,
00431                                 "Maximum extraction radius for detected "
00432                                 "objects (pixel)",
00433                                 "fors.vimos_science",
00434                                 6);
00435     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "ext_radius");
00436     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00437     cpl_parameterlist_append(recipe->parameters, p);
00438 
00439     /*
00440      * Contamination radius
00441      */
00442 
00443     p = cpl_parameter_new_value("fors.vimos_science.cont_radius",
00444                                 CPL_TYPE_INT,
00445                                 "Minimum distance at which two objects "
00446                                 "of equal luminosity do not contaminate "
00447                                 "each other (pixel)",
00448                                 "fors.vimos_science",
00449                                 0);
00450     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "cont_radius");
00451     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00452     cpl_parameterlist_append(recipe->parameters, p);
00453 
00454     /*
00455      * Object extraction method
00456      */
00457 
00458     p = cpl_parameter_new_value("fors.vimos_science.ext_mode",
00459                                 CPL_TYPE_INT,
00460                                 "Object extraction method: 0 = aperture, "
00461                                 "1 = Horne optimal extraction",
00462                                 "fors.vimos_science",
00463                                 1);
00464     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "ext_mode");
00465     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00466     cpl_parameterlist_append(recipe->parameters, p);
00467 
00468     /*
00469      * Normalise output by exposure time
00470      */
00471 
00472     p = cpl_parameter_new_value("fors.vimos_science.time_normalise",
00473                                 CPL_TYPE_BOOL,
00474                                 "Normalise output spectra by the exposure time",
00475                                 "fors.vimos_science",
00476                                 TRUE);
00477     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "time_normalise");
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 vimos_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 vimos_science(recipe->parameters, recipe->frames);
00503 }
00504 
00505 
00514 static int vimos_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 vimos_science(cpl_parameterlist *parlist, cpl_frameset *frameset)
00539 {
00540 
00541     const char *recipe = "vimos_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     double      reference;
00554     int         flux;
00555     int         flatfield;
00556     int         skyglobal;
00557     int         skylocal;
00558     int         skymedian;
00559     int         cosmics;
00560     int         slit_margin;
00561     int         ext_radius;
00562     int         cont_radius;
00563     int         ext_mode;
00564     int         time_normalise;
00565 
00566     /*
00567      * CPL objects
00568      */
00569 
00570     cpl_imagelist    *all_science;
00571     cpl_image       **images;
00572 
00573     cpl_image        *bias           = NULL;
00574     cpl_image        *norm_flat      = NULL;
00575     cpl_image        *spectra        = NULL;
00576     cpl_image        *rectified      = NULL;
00577     cpl_image        *coordinate     = NULL;
00578     cpl_image        *rainbow        = NULL;
00579     cpl_image        *mapped         = NULL;
00580     cpl_image        *mapped_sky     = NULL;
00581     cpl_image        *mapped_cleaned = NULL;
00582     cpl_image        *smapped        = NULL;
00583     cpl_image        *wavemap        = NULL;
00584     cpl_image        *skymap         = NULL;
00585     cpl_image        *skylocalmap    = NULL;
00586     cpl_image        *dummy          = NULL;
00587 
00588     cpl_table        *grism_table    = NULL;
00589     cpl_table        *overscans      = NULL;
00590     cpl_table        *wavelengths    = NULL;
00591     cpl_table        *idscoeff       = NULL;
00592     cpl_table        *slits          = NULL;
00593     cpl_table        *maskslits      = NULL;
00594     cpl_table        *polytraces     = NULL;
00595     cpl_table        *offsets        = NULL;
00596     cpl_table        *sky            = NULL;
00597 
00598     cpl_vector       *lines          = NULL;
00599 
00600     cpl_propertylist *header         = NULL;
00601     cpl_propertylist *save_header    = NULL;
00602 
00603     /*
00604      * Auxiliary variables
00605      */
00606 
00607     char        version[80];
00608     const char *science_tag;
00609     const char *master_norm_flat_tag;
00610     const char *disp_coeff_tag;
00611     const char *disp_coeff_sky_tag;
00612     const char *wavelength_map_tag;
00613     const char *wavelength_map_sky_tag;
00614     const char *curv_coeff_tag;
00615     const char *slit_location_tag;
00616     const char *reduced_science_tag;
00617     const char *reduced_sky_tag;
00618     const char *reduced_error_tag;
00619     const char *mapped_science_tag;
00620     const char *unmapped_science_tag;
00621     const char *mapped_science_sky_tag;
00622     const char *mapped_sky_tag;
00623     const char *unmapped_sky_tag;
00624     const char *global_sky_spectrum_tag;
00625     const char *object_table_tag;
00626     const char *skylines_offsets_tag;
00627     const char *key_gris_name;
00628     const char *key_gris_id;
00629     const char *key_filt_name;
00630     const char *key_filt_id;
00631     int         mos;
00632     int         treat_as_lss = 0;
00633     int         nslits;
00634     int         nscience;
00635     int         quadrant;
00636     double     *xpos;
00637     double     *exptime = NULL;
00638     double      alltime;
00639     double      mxpos;
00640     double      mean_rms;
00641     int         nlines;
00642     double     *line;
00643     int         nx, ny;
00644     double      gain;
00645     double      ron;
00646     int         standard;
00647     int         highres;
00648     int         rotate = 1;
00649     int         rotate_back = -1;
00650     int         i;
00651 
00652 
00653     snprintf(version, 80, "%s-%s", PACKAGE, PACKAGE_VERSION);
00654 
00655     cpl_msg_set_indentation(2);
00656 
00657     if (dfs_files_dont_exist(frameset))
00658         vimos_science_exit(NULL);
00659 
00660 
00661     /* 
00662      * Get configuration parameters
00663      */
00664 
00665     cpl_msg_info(recipe, "Recipe %s configuration parameters:", recipe);
00666     cpl_msg_indent_more();
00667 
00668     if (cpl_frameset_count_tags(frameset, "GRISM_TABLE") > 1)
00669         vimos_science_exit("Too many in input: GRISM_TABLE");
00670 
00671     grism_table = dfs_load_table(frameset, "GRISM_TABLE", 1);
00672 
00673     dispersion = dfs_get_parameter_double(parlist, 
00674                     "fors.vimos_science.dispersion", grism_table);
00675 
00676     if (dispersion <= 0.0)
00677         vimos_science_exit("Invalid resampling step");
00678 
00679     skyalign = dfs_get_parameter_int(parlist, 
00680                     "fors.vimos_science.skyalign", NULL);
00681 
00682     if (skyalign > 2)
00683         vimos_science_exit("Max polynomial degree for sky alignment is 2");
00684 
00685     wcolumn = dfs_get_parameter_string(parlist, 
00686                     "fors.vimos_science.wcolumn", NULL);
00687 
00688     startwavelength = dfs_get_parameter_double(parlist, 
00689                     "fors.vimos_science.startwavelength", grism_table);
00690     if (startwavelength < 3000.0 || startwavelength > 13000.0)
00691         vimos_science_exit("Invalid wavelength");
00692 
00693     endwavelength = dfs_get_parameter_double(parlist, 
00694                     "fors.vimos_science.endwavelength", grism_table);
00695     if (endwavelength < 3000.0 || endwavelength > 13000.0)
00696         vimos_science_exit("Invalid wavelength");
00697 
00698     if (endwavelength - startwavelength <= 0.0)
00699         vimos_science_exit("Invalid wavelength interval");
00700 
00701     reference = dfs_get_parameter_double(parlist,
00702                 "fors.vimos_science.reference", grism_table);
00703 
00704     if (reference < startwavelength || reference > endwavelength)
00705         vimos_science_exit("Invalid reference wavelength");
00706 
00707     flux = dfs_get_parameter_bool(parlist, "fors.vimos_science.flux", NULL);
00708 
00709     flatfield = dfs_get_parameter_bool(parlist, "fors.vimos_science.flatfield", 
00710                                        NULL);
00711 
00712     skyglobal = dfs_get_parameter_bool(parlist, "fors.vimos_science.skyglobal", 
00713                                        NULL);
00714     skylocal  = dfs_get_parameter_bool(parlist, "fors.vimos_science.skylocal", 
00715                                        NULL);
00716     skymedian = dfs_get_parameter_bool(parlist, "fors.vimos_science.skymedian", 
00717                                        NULL);
00718 
00719     if (skylocal && skyglobal)
00720         vimos_science_exit("Cannot do both local and global sky subtraction");
00721 
00722     if (skylocal && skymedian)
00723         vimos_science_exit("Cannot do sky subtraction both on extracted "
00724                            "and non-extracted spectra");
00725 
00726     cosmics = dfs_get_parameter_bool(parlist, 
00727                                      "fors.vimos_science.cosmics", NULL);
00728 
00729     if (cosmics)
00730         if (!(skyglobal || skylocal))
00731             vimos_science_exit("Cosmic rays correction requires "
00732                               "either skylocal=true or skyglobal=true");
00733 
00734     slit_margin = dfs_get_parameter_int(parlist, 
00735                                         "fors.vimos_science.slit_margin",
00736                                         NULL);
00737     if (slit_margin < 0)
00738         vimos_science_exit("Value must be zero or positive");
00739 
00740     ext_radius = dfs_get_parameter_int(parlist, 
00741                                        "fors.vimos_science.ext_radius",
00742                                        NULL);
00743     if (ext_radius < 0)
00744         vimos_science_exit("Value must be zero or positive");
00745 
00746     cont_radius = dfs_get_parameter_int(parlist, 
00747                                         "fors.vimos_science.cont_radius",
00748                                         NULL);
00749     if (cont_radius < 0)
00750         vimos_science_exit("Value must be zero or positive");
00751 
00752     ext_mode = dfs_get_parameter_int(parlist, "fors.vimos_science.ext_mode",
00753                                        NULL);
00754     if (ext_mode < 0 || ext_mode > 1)
00755         vimos_science_exit("Invalid object extraction mode");
00756 
00757     time_normalise = dfs_get_parameter_bool(parlist, 
00758                              "fors.vimos_science.time_normalise", NULL);
00759 
00760     cpl_table_delete(grism_table); grism_table = NULL;
00761 
00762     if (cpl_error_get_code())
00763         vimos_science_exit("Failure getting the configuration parameters");
00764 
00765     
00766     /* 
00767      * Check input set-of-frames
00768      */
00769 
00770     cpl_msg_indent_less();
00771     cpl_msg_info(recipe, "Check input set-of-frames:");
00772     cpl_msg_indent_more();
00773 
00774     if (!dfs_equal_keyword(frameset, "ESO OCS CON QUAD"))
00775         vimos_science_exit("Input frames are not from the same quadrant");
00776 
00777     mos = cpl_frameset_count_tags(frameset, "MOS_SCIENCE");
00778     standard = 0;
00779 
00780     if (mos == 0) {
00781         mos = cpl_frameset_count_tags(frameset, "MOS_STANDARD");
00782         standard = 1;
00783     }
00784 
00785     if (mos == 0)
00786         vimos_science_exit("Missing input scientific frame");
00787 
00788     nscience = mos;
00789 
00790     if (standard) {
00791         science_tag            = "MOS_STANDARD";
00792         reduced_science_tag    = "MOS_STANDARD_REDUCED";
00793         unmapped_science_tag   = "MOS_UNMAPPED_STANDARD";
00794         mapped_science_tag     = "MOS_STANDARD_EXTRACTED";
00795         mapped_science_sky_tag = "MOS_STANDARD_SKY_EXTRACTED";
00796     }
00797     else {
00798         science_tag            = "MOS_SCIENCE";
00799         reduced_science_tag    = "MOS_SCIENCE_REDUCED";
00800         unmapped_science_tag   = "MOS_UNMAPPED_SCIENCE";
00801         mapped_science_tag     = "MOS_SCIENCE_EXTRACTED";
00802         mapped_science_sky_tag = "MOS_SCIENCE_SKY_EXTRACTED";
00803     }
00804 
00805     reduced_sky_tag         = "MOS_SKY_REDUCED";
00806     reduced_error_tag       = "MOS_ERROR_REDUCED";
00807     master_norm_flat_tag    = "MOS_MASTER_SCREEN_FLAT";
00808     disp_coeff_sky_tag      = "MOS_DISP_COEFF_SKY";
00809     wavelength_map_sky_tag  = "MOS_WAVELENGTH_MAP_SKY";
00810     disp_coeff_tag          = "MOS_DISP_COEFF";
00811     wavelength_map_tag      = "WAVELENGTH_MAP_MXU";
00812     curv_coeff_tag          = "MOS_CURV_COEFF";
00813     slit_location_tag       = "MOS_SLIT_LOCATION";
00814     skylines_offsets_tag    = "MOS_SKYLINES_OFFSETS_SLIT";
00815     mapped_sky_tag          = "MOS_SCIENCE_SKY";
00816     unmapped_sky_tag        = "MOS_UNMAPPED_SKY";
00817     global_sky_spectrum_tag = "MOS_GLOBAL_SKY_SPECTRUM";
00818     object_table_tag        = "OBJECT_TABLE";
00819 
00820     if (cpl_frameset_count_tags(frameset, "MASTER_BIAS") == 0)
00821         vimos_science_exit("Missing required input: MASTER_BIAS");
00822 
00823     if (cpl_frameset_count_tags(frameset, "MASTER_BIAS") > 1)
00824         vimos_science_exit("Too many in input: MASTER_BIAS");
00825 
00826     if (skyalign >= 0)
00827         if (cpl_frameset_count_tags(frameset, "SKY_LINE_CATALOG") > 1)
00828             vimos_science_exit("Too many in input: SKY_LINE_CATALOG");
00829 
00830     if (cpl_frameset_count_tags(frameset, disp_coeff_tag) == 0) {
00831         cpl_msg_error(recipe, "Missing required input: %s", disp_coeff_tag);
00832         vimos_science_exit(NULL);
00833     }
00834 
00835     if (cpl_frameset_count_tags(frameset, disp_coeff_tag) > 1) {
00836         cpl_msg_error(recipe, "Too many in input: %s", disp_coeff_tag);
00837         vimos_science_exit(NULL);
00838     }
00839 
00840     if (cpl_frameset_count_tags(frameset, master_norm_flat_tag) > 1) {
00841         if (flatfield) {
00842             cpl_msg_error(recipe, "Too many in input: %s", 
00843                           master_norm_flat_tag);
00844             vimos_science_exit(NULL);
00845         }
00846         else {
00847             cpl_msg_warning(recipe, "%s in input are ignored, "
00848                             "since flat field correction was not requested", 
00849                             master_norm_flat_tag);
00850         }
00851     }
00852 
00853     if (cpl_frameset_count_tags(frameset, master_norm_flat_tag) == 1) {
00854         if (!flatfield) {
00855             cpl_msg_warning(recipe, "%s in input is ignored, "
00856                             "since flat field correction was not requested", 
00857                             master_norm_flat_tag);
00858         }
00859     }
00860 
00861     if (cpl_frameset_count_tags(frameset, master_norm_flat_tag) == 0) {
00862         if (flatfield) {
00863             cpl_msg_error(recipe, "Flat field correction was requested, "
00864                           "but no %s are found in input",
00865                           master_norm_flat_tag);
00866             vimos_science_exit(NULL);
00867         }
00868     }
00869 
00870     cpl_msg_indent_less();
00871 
00872 
00873     /*
00874      * Loading input data
00875      */
00876 
00877     exptime = cpl_calloc(nscience, sizeof(double));
00878 
00879     if (nscience > 1) {
00880 
00881         cpl_msg_info(recipe, "Load %d scientific frames and median them...",
00882                      nscience);
00883         cpl_msg_indent_more();
00884 
00885         all_science = cpl_imagelist_new();
00886 
00887         header = dfs_load_header(frameset, science_tag, 0);
00888 
00889         if (header == NULL)
00890             vimos_science_exit("Cannot load scientific frame header");
00891 
00892         alltime = exptime[0] = cpl_propertylist_get_double(header, "EXPTIME");
00893 
00894         if (cpl_error_get_code() != CPL_ERROR_NONE)
00895             vimos_science_exit("Missing keyword EXPTIME in scientific "
00896                               "frame header");
00897 
00898         cpl_propertylist_delete(header); header = NULL;
00899 
00900         cpl_msg_info(recipe, "Scientific frame 1 exposure time: %.2f s", 
00901                      exptime[0]);
00902 
00903         for (i = 1; i < nscience; i++) {
00904 
00905             header = dfs_load_header(frameset, NULL, 0);
00906 
00907             if (header == NULL)
00908                 vimos_science_exit("Cannot load scientific frame header");
00909     
00910             exptime[i] = cpl_propertylist_get_double(header, "EXPTIME");
00911 
00912             alltime += exptime[i];
00913     
00914             if (cpl_error_get_code() != CPL_ERROR_NONE)
00915                 vimos_science_exit("Missing keyword EXPTIME in scientific "
00916                                   "frame header");
00917     
00918             cpl_propertylist_delete(header); header = NULL;
00919 
00920             cpl_msg_info(recipe, "Scientific frame %d exposure time: %.2f s", 
00921                          i+1, exptime[i]);
00922         }
00923 
00924         spectra = dfs_load_image(frameset, science_tag, CPL_TYPE_FLOAT, 0, 0);
00925 
00926         if (spectra == NULL)
00927             vimos_science_exit("Cannot load scientific frame");
00928 
00929         cpl_image_divide_scalar(spectra, exptime[0]);
00930         cpl_imagelist_set(all_science, spectra, 0);
00931 
00932         for (i = 1; i < nscience; i++) {
00933 
00934             spectra = dfs_load_image(frameset, NULL, CPL_TYPE_FLOAT, 0, 0);
00935 
00936             if (spectra) {
00937                 cpl_image_divide_scalar(spectra, exptime[i]);
00938                 cpl_imagelist_set(all_science, spectra, i);
00939             }
00940             else
00941                 vimos_science_exit("Cannot load scientific frame");
00942 
00943         }
00944 
00945         spectra = cpl_imagelist_collapse_median_create(all_science);
00946         cpl_image_multiply_scalar(spectra, alltime);
00947 
00948         cpl_imagelist_delete(all_science);
00949     }
00950     else {
00951         cpl_msg_info(recipe, "Load scientific exposure...");
00952         cpl_msg_indent_more();
00953 
00954         header = dfs_load_header(frameset, science_tag, 0);
00955 
00956         if (header == NULL)
00957             vimos_science_exit("Cannot load scientific frame header");
00958 
00959         alltime = exptime[0] = cpl_propertylist_get_double(header, "EXPTIME");
00960 
00961         if (cpl_error_get_code() != CPL_ERROR_NONE)
00962             vimos_science_exit("Missing keyword EXPTIME in scientific "
00963                               "frame header");
00964 
00965         cpl_propertylist_delete(header); header = NULL;
00966 
00967         cpl_msg_info(recipe, "Scientific frame exposure time: %.2f s", 
00968                      exptime[0]);
00969 
00970         spectra = dfs_load_image(frameset, science_tag, CPL_TYPE_FLOAT, 0, 0);
00971     }
00972 
00973     if (spectra == NULL)
00974         vimos_science_exit("Cannot load scientific frame");
00975 
00976     cpl_free(exptime); exptime = NULL;
00977 
00978     cpl_msg_indent_less();
00979 
00980 
00981     /*
00982      * Get some info from header
00983      */
00984 
00985     header = dfs_load_header(frameset, science_tag, 0);
00986 
00987     if (header == NULL)
00988         vimos_science_exit("Cannot load scientific frame header");
00989 
00990     quadrant = cpl_propertylist_get_int(header, "ESO OCS CON QUAD");
00991 
00992     switch (quadrant) {
00993     case 1:
00994         key_gris_name = "ESO INS GRIS1 NAME";
00995         key_gris_id = "ESO INS GRIS1 ID";
00996         key_filt_name = "ESO INS FILT1 NAME";
00997         key_filt_id = "ESO INS FILT1 ID";
00998         break;
00999     case 2: 
01000         key_gris_name = "ESO INS GRIS2 NAME";
01001         key_gris_id = "ESO INS GRIS2 ID";
01002         key_filt_name = "ESO INS FILT2 NAME";
01003         key_filt_id = "ESO INS FILT2 ID";
01004         break;
01005     case 3: 
01006         key_gris_name = "ESO INS GRIS3 NAME";
01007         key_gris_id = "ESO INS GRIS3 ID";
01008         key_filt_name = "ESO INS FILT3 NAME";
01009         key_filt_id = "ESO INS FILT3 ID";
01010         break;
01011     case 4:
01012         key_gris_name = "ESO INS GRIS4 NAME";
01013         key_gris_id = "ESO INS GRIS4 ID";
01014         key_filt_name = "ESO INS FILT4 NAME";
01015         key_filt_id = "ESO INS FILT4 ID";
01016         break;
01017     }
01018 
01019 /*
01020     if (!dfs_equal_keyword(frameset, key_gris_id))
01021         vimos_science_exit("Input frames are not from the same grism");
01022 
01023     if (!dfs_equal_keyword(frameset, key_filt_id))
01024         vimos_science_exit("Input frames are not from the same filter");
01025 */
01026 
01027     gain = cpl_propertylist_get_double(header, "ESO DET OUT1 CONAD");
01028 
01029     if (cpl_error_get_code() != CPL_ERROR_NONE)
01030         vimos_science_exit("Missing keyword ESO DET OUT1 CONAD in scientific "
01031                           "frame header");
01032 
01033     cpl_msg_info(recipe, "The gain factor is: %.2f e-/ADU", gain);
01034 
01035     ron = cpl_propertylist_get_double(header, "ESO DET OUT1 RON");
01036 
01037     if (cpl_error_get_code() != CPL_ERROR_NONE)
01038         vimos_science_exit("Missing keyword ESO DET OUT1 RON in scientific "
01039                           "frame header");
01040 
01041     ron /= gain;     /* Convert from electrons to ADU */
01042 
01043     cpl_msg_info(recipe, "The read-out-noise is: %.2f ADU", ron);
01044 
01045     maskslits = mos_load_slits_vimos(header);
01046 
01047     /*
01048      * Check if all slits have the same X offset: in such case,
01049      * treat the observation as a long-slit one!
01050      */
01051 
01052     mxpos = cpl_table_get_column_median(maskslits, "xtop");
01053     xpos = cpl_table_get_data_double(maskslits, "xtop");
01054     nslits = cpl_table_get_nrow(maskslits);
01055      
01056     treat_as_lss = 1;
01057     for (i = 0; i < nslits; i++) { 
01058         if (fabs(mxpos-xpos[i]) > 0.01) {
01059             treat_as_lss = 0;
01060             break;
01061         }
01062     }
01063 
01064     cpl_table_delete(maskslits); maskslits = NULL;
01065 
01066     if (treat_as_lss) {
01067         cpl_msg_warning(recipe, "All MOS slits have the same offset: %.2f\n"
01068                         "The LSS data reduction strategy is applied!",
01069                         mxpos);
01070         skylines_offsets_tag   = "MOS_SKYLINES_OFFSETS_LONG";
01071     }
01072 
01073     if (treat_as_lss) {
01074         if (skylocal) {
01075             if (cosmics)
01076                 vimos_science_exit("Cosmic rays correction for long-slit-like "
01077                                   "data requires --skyglobal=true");
01078             skymedian = skylocal;
01079             skylocal = 0;
01080         }
01081     }
01082     else {
01083         if (cpl_frameset_count_tags(frameset, curv_coeff_tag) == 0) {
01084             cpl_msg_error(recipe, "Missing required input: %s", curv_coeff_tag);
01085             vimos_science_exit(NULL);
01086         }
01087 
01088         if (cpl_frameset_count_tags(frameset, curv_coeff_tag) > 1) {
01089             cpl_msg_error(recipe, "Too many in input: %s", curv_coeff_tag);
01090             vimos_science_exit(NULL);
01091         }
01092 
01093         if (cpl_frameset_count_tags(frameset, slit_location_tag) == 0) {
01094             cpl_msg_error(recipe, "Missing required input: %s", 
01095                           slit_location_tag);
01096             vimos_science_exit(NULL);
01097         }
01098         
01099         if (cpl_frameset_count_tags(frameset, slit_location_tag) > 1) {
01100             cpl_msg_error(recipe, "Too many in input: %s", slit_location_tag);
01101             vimos_science_exit(NULL);
01102         }
01103     }
01104 
01105     /* Leave the header on for the next step... */
01106 
01107 
01108     /*
01109      * Remove the master bias
01110      */
01111 
01112     cpl_msg_info(recipe, "Remove the master bias...");
01113 
01114     bias = dfs_load_image(frameset, "MASTER_BIAS", CPL_TYPE_FLOAT, 0, 1);
01115 
01116     if (bias == NULL)
01117         vimos_science_exit("Cannot load master bias");
01118 
01119     overscans = mos_load_overscans_vimos(header, 1);
01120     cpl_propertylist_delete(header); header = NULL;
01121     dummy = mos_remove_bias(spectra, bias, overscans);
01122     cpl_image_delete(spectra); spectra = dummy; dummy = NULL;
01123     cpl_image_delete(bias); bias = NULL;
01124     cpl_table_delete(overscans); overscans = NULL;
01125 
01126     if (spectra == NULL)
01127         vimos_science_exit("Cannot remove bias from scientific frame");
01128 
01129     /*
01130      * Rotate frames horizontally with red to the right
01131      */
01132 
01133     cpl_image_turn(spectra, rotate);
01134 
01135     nx = cpl_image_get_size_x(spectra);
01136     ny = cpl_image_get_size_y(spectra);
01137 
01138     cpl_msg_indent_less();
01139     cpl_msg_info(recipe, "Load normalised flat field (if present)...");
01140     cpl_msg_indent_more();
01141 
01142     if (flatfield) {
01143 
01144         norm_flat = dfs_load_image(frameset, master_norm_flat_tag, 
01145                                    CPL_TYPE_FLOAT, 0, 1);
01146 
01147         if (norm_flat) {
01148             cpl_image_turn(norm_flat, rotate);
01149             cpl_msg_info(recipe, "Apply flat field correction...");
01150             if (cpl_image_divide(spectra, norm_flat) != CPL_ERROR_NONE) {
01151                 cpl_msg_error(recipe, "Failure of flat field correction: %s",
01152                               cpl_error_get_message());
01153                 vimos_science_exit(NULL);
01154             }
01155             cpl_image_delete(norm_flat); norm_flat = NULL;
01156         }
01157         else {
01158             cpl_msg_error(recipe, "Cannot load input %s for flat field "
01159                           "correction", master_norm_flat_tag);
01160             vimos_science_exit(NULL);
01161         }
01162 
01163     }
01164 
01165 
01166     if (skyalign >= 0) {
01167         cpl_msg_indent_less();
01168         cpl_msg_info(recipe, "Load input sky line catalog...");
01169         cpl_msg_indent_more();
01170 
01171         wavelengths = dfs_load_table(frameset, "SKY_LINE_CATALOG", 1);
01172 
01173         if (wavelengths) {
01174 
01175             /*
01176              * Cast the wavelengths into a (double precision) CPL vector
01177              */
01178 
01179             nlines = cpl_table_get_nrow(wavelengths);
01180 
01181             if (nlines == 0)
01182                 vimos_science_exit("Empty input sky line catalog");
01183 
01184             if (cpl_table_has_column(wavelengths, wcolumn) != 1) {
01185                 cpl_msg_error(recipe, "Missing column %s in input line "
01186                               "catalog table", wcolumn);
01187                 vimos_science_exit(NULL);
01188             }
01189 
01190             line = cpl_malloc(nlines * sizeof(double));
01191     
01192             for (i = 0; i < nlines; i++)
01193                 line[i] = cpl_table_get(wavelengths, wcolumn, i, NULL);
01194 
01195             cpl_table_delete(wavelengths); wavelengths = NULL;
01196 
01197             lines = cpl_vector_wrap(nlines, line);
01198         }
01199         else {
01200             cpl_msg_info(recipe, "No sky line catalog found in input - fine!");
01201         }
01202     }
01203 
01204 
01205     /*
01206      * Load the spectral curvature table, or provide a dummy one
01207      * in case of LSS-like data (single slit with flat curvature)
01208      */
01209 
01210     if (!treat_as_lss) {
01211         polytraces = dfs_load_table(frameset, curv_coeff_tag, 1);
01212         if (polytraces == NULL)
01213             vimos_science_exit("Cannot load spectral curvature table");
01214     }
01215 
01216 
01217     /*
01218      * Load the slit location table, or provide a dummy one in case
01219      * of LSS-like data (single slit spanning whole image)
01220      */
01221 
01222     if (treat_as_lss) {
01223         slits = cpl_table_new(1);
01224         cpl_table_new_column(slits, "slit_id", CPL_TYPE_INT);
01225         cpl_table_set_int(slits, "slit_id", 0, 1);
01226         cpl_table_new_column(slits, "position", CPL_TYPE_INT);
01227         cpl_table_set_int(slits, "position", 0, 0);
01228         cpl_table_new_column(slits, "length", CPL_TYPE_INT);
01229         cpl_table_set_int(slits, "length", 0, ny);
01230     }
01231     else {
01232         slits = dfs_load_table(frameset, slit_location_tag, 1);
01233         if (slits == NULL)
01234             vimos_science_exit("Cannot load slits location table");
01235         mos_rotate_slits(slits, -rotate, nx, ny);
01236     }
01237 
01238 
01239     /*
01240      * Load the wavelength calibration table
01241      */
01242 
01243     idscoeff = dfs_load_table(frameset, disp_coeff_tag, 1);
01244     if (idscoeff == NULL)
01245         vimos_science_exit("Cannot load wavelength calibration table");
01246 
01247     cpl_msg_indent_less();
01248     cpl_msg_info(recipe, "Processing scientific spectra...");
01249     cpl_msg_indent_more();
01250 
01251     /*
01252      * This one will also generate the spatial map from the spectral 
01253      * curvature table (in the case of multislit data)
01254      */
01255 
01256     if (treat_as_lss) {
01257         smapped = cpl_image_duplicate(spectra);
01258     }
01259     else {
01260         coordinate = cpl_image_new(nx, ny, CPL_TYPE_FLOAT);
01261 
01262         smapped = mos_spatial_calibration(spectra, slits, polytraces, reference,
01263                                           startwavelength, endwavelength,
01264                                           dispersion, flux, coordinate);
01265     }
01266 
01267 
01268     /*
01269      * Generate a rectified wavelength map from the wavelength calibration 
01270      * table
01271      */
01272 
01273     rainbow = mos_map_idscoeff(idscoeff, nx, reference, startwavelength, 
01274                                endwavelength);
01275 
01276     if (dispersion > 1.0)
01277         highres = 0;
01278     else
01279         highres = 1;
01280 
01281     if (skyalign >= 0) {
01282         if (skyalign) {
01283             cpl_msg_info(recipe, "Align wavelength solution to reference "
01284             "skylines applying %d order residual fit...", skyalign);
01285         }
01286         else {
01287             cpl_msg_info(recipe, "Align wavelength solution to reference "
01288             "skylines applying median offset...");
01289         }
01290 
01291         if (treat_as_lss) {
01292             offsets = mos_wavelength_align_lss(smapped, reference, 
01293                                                startwavelength, endwavelength, 
01294                                                idscoeff, lines, highres, 
01295                                                skyalign, rainbow, 4);
01296         }
01297         else {
01298             offsets = mos_wavelength_align(smapped, slits, reference, 
01299                                            startwavelength, endwavelength, 
01300                                            idscoeff, lines, highres, skyalign, 
01301                                            rainbow, 4);
01302         }
01303 
01304         cpl_vector_delete(lines); lines = NULL;
01305 
01306         if (offsets) {
01307             if (standard)
01308                 cpl_msg_warning(recipe, "Alignment of the wavelength solution "
01309                                 "to reference sky lines may be unreliable in "
01310                                 "this case!");
01311 
01312             if (dfs_save_table(frameset, offsets, skylines_offsets_tag, NULL, 
01313                                parlist, recipe, version))
01314                 vimos_science_exit(NULL);
01315 
01316             cpl_table_delete(offsets); offsets = NULL;
01317         }
01318         else {
01319             cpl_msg_warning(recipe, "Alignment of the wavelength solution "
01320                             "to reference sky lines could not be done!");
01321             skyalign = -1;
01322         }
01323 
01324     }
01325 
01326     if (treat_as_lss) {
01327         wavemap = rainbow;
01328         rainbow = NULL;
01329     }
01330     else {
01331         wavemap = mos_map_wavelengths(coordinate, rainbow, slits, 
01332                                       polytraces, reference, 
01333                                       startwavelength, endwavelength,
01334                                       dispersion);
01335     }
01336 
01337     cpl_image_delete(rainbow); rainbow = NULL;
01338     cpl_image_delete(coordinate); coordinate = NULL;
01339 
01340     /*
01341      * Here the wavelength calibrated slit spectra are created. This frame
01342      * contains sky_science.
01343      */
01344 
01345     mapped_sky = mos_wavelength_calibration(smapped, reference,
01346                                             startwavelength, endwavelength,
01347                                             dispersion, idscoeff, flux);
01348 
01349     cpl_msg_indent_less();
01350     cpl_msg_info(recipe, "Check applied wavelength against skylines...");
01351     cpl_msg_indent_more();
01352 
01353     mean_rms = mos_distortions_rms(mapped_sky, NULL, startwavelength,
01354                                    dispersion, 6, highres);
01355 
01356     cpl_msg_info(recipe, "Mean residual: %f", mean_rms);
01357 
01358     mean_rms = cpl_table_get_column_mean(idscoeff, "error");
01359 
01360     cpl_msg_info(recipe, "Mean model accuracy: %f pixel (%f A)",
01361                  mean_rms, mean_rms * dispersion);
01362 
01363     header = cpl_propertylist_new();
01364     cpl_propertylist_update_double(header, "CRPIX1", 1.0);
01365     cpl_propertylist_update_double(header, "CRPIX2", 1.0);
01366     cpl_propertylist_update_double(header, "CRVAL1", 
01367                                    startwavelength + dispersion/2);
01368     cpl_propertylist_update_double(header, "CRVAL2", 1.0);
01369     /* cpl_propertylist_update_double(header, "CDELT1", dispersion);
01370     cpl_propertylist_update_double(header, "CDELT2", 1.0); */
01371     cpl_propertylist_update_double(header, "CD1_1", dispersion);
01372     cpl_propertylist_update_double(header, "CD1_2", 0.0);
01373     cpl_propertylist_update_double(header, "CD2_1", 0.0);
01374     cpl_propertylist_update_double(header, "CD2_2", 1.0);
01375     cpl_propertylist_update_string(header, "CTYPE1", "LINEAR");
01376     cpl_propertylist_update_string(header, "CTYPE2", "PIXEL");
01377 
01378     if (time_normalise) {
01379         dummy = cpl_image_divide_scalar_create(mapped_sky, alltime);
01380         if (dfs_save_image(frameset, dummy, mapped_science_sky_tag, header, 
01381                            parlist, recipe, version))
01382             vimos_science_exit(NULL);
01383         cpl_image_delete(dummy); dummy = NULL;
01384     }
01385     else {
01386         if (dfs_save_image(frameset, mapped_sky, mapped_science_sky_tag, 
01387                            header, parlist, recipe, version))
01388             vimos_science_exit(NULL);
01389     }
01390 
01391     if (skyglobal == 0 && skymedian == 0 && skylocal == 0) {
01392         cpl_image_delete(mapped_sky); mapped_sky = NULL;
01393     }
01394 
01395     if (skyglobal || skylocal) {
01396 
01397         cpl_msg_indent_less();
01398 
01399         if (skyglobal) {
01400             cpl_msg_info(recipe, "Global sky determination...");
01401             cpl_msg_indent_more();
01402             skymap = cpl_image_new(nx, ny, CPL_TYPE_FLOAT);
01403             sky = mos_sky_map_super(spectra, wavemap, dispersion, 
01404                                     2.0, 50, skymap);
01405             if (sky)
01406                 cpl_image_subtract(spectra, skymap);
01407             else
01408                 cpl_image_delete(skymap); skymap = NULL;
01409         }
01410         else {
01411             cpl_msg_info(recipe, "Local sky determination...");
01412             cpl_msg_indent_more();
01413             skymap = mos_subtract_sky(spectra, slits, polytraces, reference,
01414                            startwavelength, endwavelength, dispersion);
01415         }
01416 
01417         if (skymap) {
01418             if (skyglobal) {
01419                 if (time_normalise)
01420                     cpl_table_divide_scalar(sky, "sky", alltime);
01421                 if (dfs_save_table(frameset, sky, global_sky_spectrum_tag, 
01422                                    NULL, parlist, recipe, version))
01423                     vimos_science_exit(NULL);
01424     
01425                 cpl_table_delete(sky); sky = NULL;
01426             }
01427 
01428             save_header = dfs_load_header(frameset, science_tag, 0);
01429 
01430             if (time_normalise)
01431                 cpl_image_divide_scalar(skymap, alltime);
01432             cpl_image_turn(skymap, rotate_back);
01433             if (dfs_save_image(frameset, skymap, unmapped_sky_tag,
01434                                save_header, parlist, recipe, version))
01435                 vimos_science_exit(NULL);
01436 
01437             cpl_image_delete(skymap); skymap = NULL;
01438 
01439             cpl_image_turn(spectra, rotate_back);
01440             if (dfs_save_image(frameset, spectra, unmapped_science_tag,
01441                                save_header, parlist, recipe, version))
01442                 vimos_science_exit(NULL);
01443             cpl_image_turn(spectra, rotate);
01444 
01445             cpl_propertylist_delete(save_header); save_header = NULL;
01446 
01447             if (cosmics) {
01448                 cpl_msg_info(recipe, "Removing cosmic rays...");
01449                 mos_clean_cosmics(spectra, gain, -1., -1.);
01450             }
01451 
01452             /*
01453              * The spatially rectified image, that contained the sky,
01454              * is replaced by a sky-subtracted spatially rectified image:
01455              */
01456 
01457             cpl_image_delete(smapped); smapped = NULL;
01458 
01459             if (treat_as_lss) {
01460                 smapped = cpl_image_duplicate(spectra);
01461             }
01462             else {
01463                 smapped = mos_spatial_calibration(spectra, slits, polytraces, 
01464                                                   reference, startwavelength, 
01465                                                   endwavelength, dispersion, 
01466                                                   flux, NULL);
01467             }
01468         }
01469         else {
01470             cpl_msg_warning(recipe, "Sky subtraction failure");
01471             if (cosmics)
01472                 cpl_msg_warning(recipe, "Cosmic rays removal not performed!");
01473             cosmics = skylocal = skyglobal = 0;
01474         }
01475     }
01476 
01477     cpl_image_delete(spectra); spectra = NULL;
01478     cpl_table_delete(polytraces); polytraces = NULL;
01479 
01480     if (skyalign >= 0) {
01481         save_header = dfs_load_header(frameset, science_tag, 0);
01482         cpl_image_turn(wavemap, rotate_back);
01483         if (dfs_save_image(frameset, wavemap, wavelength_map_sky_tag,
01484                            save_header, parlist, recipe, version))
01485             vimos_science_exit(NULL);
01486         cpl_propertylist_delete(save_header); save_header = NULL;
01487     }
01488 
01489     cpl_image_delete(wavemap); wavemap = NULL;
01490 
01491     mapped = mos_wavelength_calibration(smapped, reference,
01492                                         startwavelength, endwavelength,
01493                                         dispersion, idscoeff, flux);
01494 
01495     cpl_image_delete(smapped); smapped = NULL;
01496 
01497     if (skymedian) {
01498         cpl_msg_indent_less();
01499         cpl_msg_info(recipe, "Local sky determination...");
01500         cpl_msg_indent_more();
01501     
01502         skylocalmap = mos_sky_local_old(mapped, slits);       
01503         cpl_image_subtract(mapped, skylocalmap);
01504         cpl_image_delete(skylocalmap); skylocalmap = NULL;
01505     }
01506 
01507     if (skyglobal || skymedian || skylocal) {
01508 
01509         skylocalmap = cpl_image_subtract_create(mapped_sky, mapped);
01510 
01511         cpl_image_delete(mapped_sky); mapped_sky = NULL;
01512 
01513         if (time_normalise) {
01514             dummy = cpl_image_divide_scalar_create(skylocalmap, alltime);
01515             if (dfs_save_image(frameset, dummy, mapped_sky_tag, header,
01516                                parlist, recipe, version))
01517                 vimos_science_exit(NULL);
01518             cpl_image_delete(dummy); dummy = NULL;
01519         }
01520         else {
01521             if (dfs_save_image(frameset, skylocalmap, mapped_sky_tag, header,
01522                                parlist, recipe, version))
01523                 vimos_science_exit(NULL);
01524         }
01525 
01526         cpl_msg_indent_less();
01527         cpl_msg_info(recipe, "Object detection...");
01528         cpl_msg_indent_more();
01529 
01530         if (cosmics || nscience > 1) {
01531             dummy = mos_detect_objects(mapped, slits, slit_margin, ext_radius, 
01532                                        cont_radius);
01533         }
01534         else {
01535             mapped_cleaned = cpl_image_duplicate(mapped);
01536             mos_clean_cosmics(mapped_cleaned, gain, -1., -1.);
01537             dummy = mos_detect_objects(mapped_cleaned, slits, slit_margin, 
01538                                        ext_radius, cont_radius);
01539 
01540             cpl_image_delete(mapped_cleaned); mapped_cleaned = NULL;
01541         }
01542 
01543         cpl_image_delete(dummy); dummy = NULL;
01544 
01545         mos_rotate_slits(slits, rotate, ny, nx);
01546         if (dfs_save_table(frameset, slits, object_table_tag, NULL, parlist, 
01547                            recipe, version))
01548             vimos_science_exit(NULL);
01549 
01550         cpl_msg_indent_less();
01551         cpl_msg_info(recipe, "Object extraction...");
01552         cpl_msg_indent_more();
01553 
01554         images = mos_extract_objects(mapped, skylocalmap, slits, 
01555                                      ext_mode, ron, gain, 1);
01556 
01557         cpl_image_delete(skylocalmap); skylocalmap = NULL;
01558 
01559         if (images) {
01560             if (time_normalise)
01561                 cpl_image_divide_scalar(images[0], alltime);
01562             if (dfs_save_image(frameset, images[0], reduced_science_tag, header,
01563                                parlist, recipe, version))
01564                 vimos_science_exit(NULL);
01565             cpl_image_delete(images[0]);
01566     
01567             if (time_normalise)
01568                 cpl_image_divide_scalar(images[1], alltime);
01569             if (dfs_save_image(frameset, images[1], reduced_sky_tag, header,
01570                                parlist, recipe, version))
01571                 vimos_science_exit(NULL);
01572             cpl_image_delete(images[1]);
01573     
01574             if (time_normalise)
01575                 cpl_image_divide_scalar(images[2], alltime);
01576             if (dfs_save_image(frameset, images[2], reduced_error_tag, header,
01577                                parlist, recipe, version))
01578                 vimos_science_exit(NULL);
01579             cpl_image_delete(images[2]);
01580 
01581             cpl_free(images);
01582         }
01583         else {
01584             cpl_msg_warning(recipe, "No objects found: the products "
01585                             "%s, %s, and %s are not created", 
01586                             reduced_science_tag, reduced_sky_tag, 
01587                             reduced_error_tag);
01588         }
01589 
01590     }
01591 
01592     cpl_table_delete(slits); slits = NULL;
01593 
01594     if (skyalign >= 0) {
01595         if (dfs_save_table(frameset, idscoeff, disp_coeff_sky_tag, NULL, 
01596                            parlist, recipe, version))
01597             vimos_science_exit(NULL);
01598     }
01599 
01600     cpl_table_delete(idscoeff); idscoeff = NULL;
01601 
01602     if (skyglobal || skymedian || skylocal) {
01603         if (time_normalise)
01604             cpl_image_divide_scalar(mapped, alltime);
01605         if (dfs_save_image(frameset, mapped, mapped_science_tag, header, 
01606                            parlist, recipe, version))
01607             vimos_science_exit(NULL);
01608     }
01609 
01610     cpl_image_delete(mapped); mapped = NULL;
01611     cpl_propertylist_delete(header); header = NULL;
01612 
01613     if (cpl_error_get_code()) {
01614         cpl_msg_error(cpl_error_get_where(), cpl_error_get_message());
01615         vimos_science_exit(NULL);
01616     }
01617     else 
01618         return 0;
01619 }

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