recipes/fors_extract.c

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

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