fors_pmos_calib.c

00001 /* $Id: fors_pmos_calib.c,v 1.17 2009/05/04 10:04:08 cizzo Exp $
00002  *
00003  * This file is part of the FORS Data Reduction Pipeline
00004  * Copyright (C) 2006 European Southern Observatory
00005  *
00006  * This program is free software; you can redistribute it and/or modify
00007  * it under the terms of the GNU General Public License as published by
00008  * the Free Software Foundation; either version 2 of the License, or
00009  * (at your option) any later version.
00010  *
00011  * This program is distributed in the hope that it will be useful,
00012  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014  * GNU General Public License for more details.
00015  *
00016  * You should have received a copy of the GNU General Public License
00017  * along with this program; if not, write to the Free Software
00018  * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
00019  */
00020 
00021 /*
00022  * $Author: cizzo $
00023  * $Date: 2009/05/04 10:04:08 $
00024  * $Revision: 1.17 $
00025  * $Name:  $
00026  */
00027 
00028 #ifdef HAVE_CONFIG_H
00029 #include <config.h>
00030 #endif
00031 
00032 #include <math.h>
00033 #include <string.h>
00034 #include <ctype.h>
00035 #include <cpl.h>
00036 #include <moses.h>
00037 #include <fors_dfs.h>
00038 #include <fors_qc.h>
00039 
00040 #define OFFSET    50
00041 #define TOLERANCE 10
00042 
00043 static int fors_pmos_calib_create(cpl_plugin *);
00044 static int fors_pmos_calib_exec(cpl_plugin *);
00045 static int fors_pmos_calib_destroy(cpl_plugin *);
00046 static int fors_pmos_calib(cpl_parameterlist *, cpl_frameset *);
00047 
00048 static char fors_pmos_calib_description[] =
00049 "This recipe is used to identify reference lines on PMOS arc lamp\n"
00050 "exposures, and trace the spectral edges on the corresponding flat field\n"
00051 "exposures. This information is used to determine the spectral extraction\n"
00052 "mask to be applied in the scientific data reduction, performed with the\n"
00053 "recipe fors_science.\n"
00054 "This recipe accepts both FORS1 and FORS2 frames. The input arc lamps and\n"
00055 "flat field exposures are assumed to be obtained quasi-simultaneously, so\n"
00056 "that they would be described by exactly the same instrument distortions.\n"
00057 "A line catalog must be specified, containing the wavelengths of the\n"
00058 "reference arc lamp lines used for the wavelength calibration. A grism\n"
00059 "table (typically depending on the instrument mode, and in particular on\n"
00060 "the grism used) may also be specified: this table contains a default\n"
00061 "recipe parameter setting to control the way spectra are extracted for\n"
00062 "a specific instrument mode, as it is used for automatic run of the\n"
00063 "pipeline on Paranal and in Garching. If this table is specified, it\n"
00064 "will modify the default recipe parameter setting, with the exception of\n"
00065 "those parameters which have been explicitly modifyed on the command line.\n"
00066 "If a grism table is not specified, the input recipe parameters values\n"
00067 "will always be read from the command line, or from an esorex configuration\n"
00068 "file if present, or from their generic default values (that are rarely\n"
00069 "meaningful). Finally a master bias frame must be input to this recipe.\n" 
00070 "The products SPECTRA_DETECTION_PMOS, SLIT_MAP_PMOS, and\n" 
00071 "DISP_RESIDUALS_PMOS, are just created if the --check parameter is set to\n" 
00072 "true.\n\n" 
00073 "Input files:\n\n"
00074 "  DO category:              Type:       Explanation:          Required:\n"
00075 "  SCREEN_FLAT_PMOS          Raw         Flat field exposures     Y\n"
00076 "  LAMP_PMOS                 Raw         Arc lamp exposure        Y\n"
00077 "  MASTER_BIAS or BIAS       Calib       Bias frame               Y\n"
00078 "  MASTER_LINECAT            Calib       Line catalog             Y\n"
00079 "  GRISM_TABLE               Calib       Grism table              .\n"
00080 "  MASTER_DISTORTION_TABLE   Calib       Master distortions table Y\n\n"
00081 "Output files:\n\n" 
00082 "  DO category:              Data type:  Explanation:\n"
00083 "  MASTER_SCREEN_FLAT_PMOS   FITS image  Combined (sum) flat field\n"
00084 "  MASTER_NORM_FLAT_PMOS     FITS image  Normalised flat field\n"
00085 "  MAPPED_SCREEN_FLAT_PMOS   FITS image  Wavelength calibrated flat field\n"
00086 "  MAPPED_NORM_FLAT_PMOS     FITS image  Wavelength calibrated normalised flat\n"
00087 "  REDUCED_LAMP_PMOS         FITS image  Wavelength calibrated arc spectrum\n"
00088 "  DISP_COEFF_PMOS           FITS table  Inverse dispersion coefficients\n"
00089 "  DISP_RESIDUALS_PMOS       FITS image  Residuals in wavelength calibration\n"
00090 "  DISP_RESIDUALS_TABLE_PMOS FITS table  Residuals in wavelength calibration\n"
00091 "  DELTA_IMAGE_PMOS          FITS image  Offset vs linear wavelength calib\n"
00092 "  WAVELENGTH_MAP_PMOS       FITS image  Wavelength for each pixel on CCD\n"
00093 "  SPECTRA_DETECTION_PMOS    FITS image  Check for preliminary detection\n"
00094 "  SLIT_MAP_PMOS             FITS image  Map of central wavelength on CCD\n"
00095 "  CURV_TRACES_PMOS          FITS table  Spectral curvature traces\n"
00096 "  CURV_COEFF_PMOS           FITS table  Spectral curvature coefficients\n"
00097 "  SPATIAL_MAP_PMOS          FITS image  Spatial position along slit on CCD\n"
00098 "  SPECTRAL_RESOLUTION_PMOS  FITS table  Resolution at reference arc lines\n"
00099 "  SLIT_LOCATION_PMOS        FITS table  Slits on product frames and CCD\n\n";
00100 
00101 #define fors_pmos_calib_exit(message)              \
00102 {                                             \
00103 if (message) cpl_msg_error(recipe, message);  \
00104 cpl_free(instrume);                           \
00105 cpl_free(pipefile);                           \
00106 cpl_free(fiterror);                           \
00107 cpl_free(fitlines);                           \
00108 cpl_image_delete(bias);                       \
00109 cpl_image_delete(master_bias);                \
00110 cpl_image_delete(coordinate);                 \
00111 cpl_image_delete(checkwave);                  \
00112 cpl_image_delete(flat);                       \
00113 cpl_image_delete(master_flat);                \
00114 cpl_image_delete(norm_flat);                  \
00115 cpl_image_delete(mapped_flat);                \
00116 cpl_image_delete(mapped_nflat);               \
00117 cpl_image_delete(rainbow);                    \
00118 cpl_image_delete(rectified);                  \
00119 cpl_image_delete(residual);                   \
00120 cpl_image_delete(smo_flat);                   \
00121 cpl_image_delete(spatial);                    \
00122 cpl_image_delete(spectra);                    \
00123 cpl_image_delete(wavemap);                    \
00124 cpl_image_delete(delta);                      \
00125 cpl_image_delete(rect_flat);                  \
00126 cpl_image_delete(rect_nflat);                 \
00127 cpl_image_delete(mapped_flat);                \
00128 cpl_image_delete(mapped_nflat);               \
00129 cpl_mask_delete(refmask);                     \
00130 cpl_propertylist_delete(header);              \
00131 cpl_propertylist_delete(save_header);         \
00132 cpl_propertylist_delete(qclist);              \
00133 cpl_table_delete(grism_table);                \
00134 cpl_table_delete(idscoeff);                   \
00135 cpl_table_delete(idscoeff_all);               \
00136 cpl_table_delete(restable);                   \
00137 cpl_table_delete(maskslits);                  \
00138 cpl_table_delete(overscans);                  \
00139 cpl_table_delete(traces);                     \
00140 cpl_table_delete(polytraces);                 \
00141 cpl_table_delete(slits);                      \
00142 cpl_table_delete(restab);                     \
00143 cpl_table_delete(global);                     \
00144 cpl_table_delete(wavelengths);                \
00145 cpl_vector_delete(lines);                     \
00146 cpl_msg_indent_less();                        \
00147 return -1;                                    \
00148 }
00149 
00150 #define fors_pmos_calib_exit_memcheck(message)       \
00151 {                                               \
00152 if (message) cpl_msg_info(recipe, message);     \
00153 printf("free instrume (%p)\n", instrume);       \
00154 cpl_free(instrume);                             \
00155 printf("free pipefile (%p)\n", pipefile);       \
00156 cpl_free(pipefile);                             \
00157 printf("free fiterror (%p)\n", fiterror);       \
00158 cpl_free(fiterror);                             \
00159 printf("free fitlines (%p)\n", fitlines);       \
00160 cpl_free(fitlines);                             \
00161 printf("free bias (%p)\n", bias);               \
00162 cpl_image_delete(bias);                         \
00163 printf("free master_bias (%p)\n", master_bias); \
00164 cpl_image_delete(master_bias);                  \
00165 printf("free coordinate (%p)\n", coordinate);   \
00166 cpl_image_delete(coordinate);                   \
00167 printf("free checkwave (%p)\n", checkwave);     \
00168 cpl_image_delete(checkwave);                    \
00169 printf("free flat (%p)\n", flat);               \
00170 cpl_image_delete(flat);                         \
00171 printf("free master_flat (%p)\n", master_flat); \
00172 cpl_image_delete(master_flat);                  \
00173 printf("free norm_flat (%p)\n", norm_flat);     \
00174 cpl_image_delete(norm_flat);                    \
00175 printf("free mapped_flat (%p)\n", mapped_flat); \
00176 cpl_image_delete(mapped_flat);                  \
00177 printf("free mapped_nflat (%p)\n", mapped_nflat); \
00178 cpl_image_delete(mapped_nflat);                 \
00179 printf("free rainbow (%p)\n", rainbow);         \
00180 cpl_image_delete(rainbow);                      \
00181 printf("free rectified (%p)\n", rectified);     \
00182 cpl_image_delete(rectified);                    \
00183 printf("free residual (%p)\n", residual);       \
00184 cpl_image_delete(residual);                     \
00185 printf("free smo_flat (%p)\n", smo_flat);       \
00186 cpl_image_delete(smo_flat);                     \
00187 printf("free spatial (%p)\n", spatial);         \
00188 cpl_image_delete(spatial);                      \
00189 printf("free spectra (%p)\n", spectra);         \
00190 cpl_image_delete(spectra);                      \
00191 printf("free wavemap (%p)\n", wavemap);         \
00192 cpl_image_delete(wavemap);                      \
00193 printf("free delta (%p)\n", delta);             \
00194 cpl_image_delete(delta);                        \
00195 printf("free rect_flat (%p)\n", rect_flat);     \
00196 cpl_image_delete(rect_flat);                    \
00197 printf("free rect_nflat (%p)\n", rect_nflat);   \
00198 cpl_image_delete(rect_nflat);                   \
00199 printf("free refmask (%p)\n", refmask);         \
00200 cpl_mask_delete(refmask);                       \
00201 printf("free header (%p)\n", header);           \
00202 cpl_propertylist_delete(header);                \
00203 printf("free save_header (%p)\n", save_header); \
00204 cpl_propertylist_delete(save_header);           \
00205 printf("free qclist (%p)\n", qclist);           \
00206 cpl_propertylist_delete(qclist);                \
00207 printf("free grism_table (%p)\n", grism_table); \
00208 cpl_table_delete(grism_table);                  \
00209 printf("free idscoeff (%p)\n", idscoeff);       \
00210 cpl_table_delete(idscoeff);                     \
00211 printf("free idscoeff_all (%p)\n", idscoeff_all);  \
00212 cpl_table_delete(idscoeff_all);                 \
00213 printf("free restable (%p)\n", restable);       \
00214 cpl_table_delete(restable);                     \
00215 printf("free maskslits (%p)\n", maskslits);     \
00216 cpl_table_delete(maskslits);                    \
00217 printf("free overscans (%p)\n", overscans);     \
00218 cpl_table_delete(overscans);                    \
00219 printf("free traces (%p)\n", traces);           \
00220 cpl_table_delete(traces);                       \
00221 printf("free polytraces (%p)\n", polytraces);   \
00222 cpl_table_delete(polytraces);                   \
00223 printf("free slits (%p)\n", slits);             \
00224 cpl_table_delete(slits);                        \
00225 printf("free restab (%p)\n", restab);           \
00226 cpl_table_delete(restab);                       \
00227 printf("free global (%p)\n", global);           \
00228 cpl_table_delete(global);                       \
00229 printf("free wavelengths (%p)\n", wavelengths); \
00230 cpl_table_delete(wavelengths);                  \
00231 printf("free lines (%p)\n", lines);             \
00232 cpl_vector_delete(lines);                       \
00233 cpl_msg_indent_less();                          \
00234 return 0;                                       \
00235 }
00236 
00237 
00249 int cpl_plugin_get_info(cpl_pluginlist *list)
00250 {
00251     cpl_recipe *recipe = cpl_calloc(1, sizeof *recipe );
00252     cpl_plugin *plugin = &recipe->interface;
00253 
00254     cpl_plugin_init(plugin,
00255                     CPL_PLUGIN_API,
00256                     FORS_BINARY_VERSION,
00257                     CPL_PLUGIN_TYPE_RECIPE,
00258                     "fors_pmos_calib",
00259                     "Determination of the extraction mask",
00260                     fors_pmos_calib_description,
00261                     "Carlo Izzo",
00262                     PACKAGE_BUGREPORT,
00263     "This file is currently part of the FORS Instrument Pipeline\n"
00264     "Copyright (C) 2002-2006 European Southern Observatory\n\n"
00265     "This program is free software; you can redistribute it and/or modify\n"
00266     "it under the terms of the GNU General Public License as published by\n"
00267     "the Free Software Foundation; either version 2 of the License, or\n"
00268     "(at your option) any later version.\n\n"
00269     "This program is distributed in the hope that it will be useful,\n"
00270     "but WITHOUT ANY WARRANTY; without even the implied warranty of\n"
00271     "MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the\n"
00272     "GNU General Public License for more details.\n\n"
00273     "You should have received a copy of the GNU General Public License\n"
00274     "along with this program; if not, write to the Free Software Foundation,\n"
00275     "Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA\n",
00276                     fors_pmos_calib_create,
00277                     fors_pmos_calib_exec,
00278                     fors_pmos_calib_destroy);
00279 
00280     cpl_pluginlist_append(list, plugin);
00281     
00282     return 0;
00283 }
00284 
00285 
00296 static int fors_pmos_calib_create(cpl_plugin *plugin)
00297 {
00298     cpl_recipe    *recipe;
00299     cpl_parameter *p;
00300 
00301 
00302     /* 
00303      * Check that the plugin is part of a valid recipe 
00304      */
00305 
00306     if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE) 
00307         recipe = (cpl_recipe *)plugin;
00308     else 
00309         return -1;
00310 
00311     /* 
00312      * Create the parameters list in the cpl_recipe object 
00313      */
00314 
00315     recipe->parameters = cpl_parameterlist_new(); 
00316 
00317 
00318     /*
00319      * Dispersion
00320      */
00321 
00322     p = cpl_parameter_new_value("fors.fors_pmos_calib.dispersion",
00323                                 CPL_TYPE_DOUBLE,
00324                                 "Expected spectral dispersion (Angstrom/pixel)",
00325                                 "fors.fors_pmos_calib",
00326                                 0.0);
00327     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "dispersion");
00328     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00329     cpl_parameterlist_append(recipe->parameters, p);
00330 
00331     /*
00332      * Peak detection level
00333      */
00334 
00335     p = cpl_parameter_new_value("fors.fors_pmos_calib.peakdetection",
00336                                 CPL_TYPE_DOUBLE,
00337                                 "Initial peak detection threshold (ADU)",
00338                                 "fors.fors_pmos_calib",
00339                                 0.0);
00340     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "peakdetection");
00341     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00342     cpl_parameterlist_append(recipe->parameters, p);
00343 
00344     /* 
00345      * Degree of wavelength calibration polynomial
00346      */
00347 
00348     p = cpl_parameter_new_value("fors.fors_pmos_calib.wdegree",
00349                                 CPL_TYPE_INT,
00350                                 "Degree of wavelength calibration polynomial",
00351                                 "fors.fors_pmos_calib",
00352                                 0);
00353     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "wdegree");
00354     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00355     cpl_parameterlist_append(recipe->parameters, p);
00356 
00357     /*
00358      * Reference lines search radius
00359      */
00360 
00361     p = cpl_parameter_new_value("fors.fors_pmos_calib.wradius",
00362                                 CPL_TYPE_INT,
00363                                 "Search radius if iterating pattern-matching "
00364                                 "with first-guess method",
00365                                 "fors.fors_pmos_calib",
00366                                 4);
00367     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "wradius");
00368     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00369     cpl_parameterlist_append(recipe->parameters, p);
00370 
00371     /*
00372      * Rejection threshold in dispersion relation polynomial fitting
00373      */
00374 
00375     p = cpl_parameter_new_value("fors.fors_pmos_calib.wreject",
00376                                 CPL_TYPE_DOUBLE,
00377                                 "Rejection threshold in dispersion "
00378                                 "relation fit (pixel)",
00379                                 "fors.fors_pmos_calib",
00380                                 0.7);
00381     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "wreject");
00382     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00383     cpl_parameterlist_append(recipe->parameters, p);
00384 
00385     /*
00386      * Line catalog table column containing the reference wavelengths
00387      */
00388 
00389     p = cpl_parameter_new_value("fors.fors_pmos_calib.wcolumn",
00390                                 CPL_TYPE_STRING,
00391                                 "Name of line catalog table column "
00392                                 "with wavelengths",
00393                                 "fors.fors_pmos_calib",
00394                                 "WLEN");
00395     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "wcolumn");
00396     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00397     cpl_parameterlist_append(recipe->parameters, p);
00398 
00399     /*
00400      * Degree of spectral curvature polynomial
00401      */
00402 
00403     p = cpl_parameter_new_value("fors.fors_pmos_calib.cdegree",
00404                                 CPL_TYPE_INT,
00405                                 "Degree of spectral curvature polynomial",
00406                                 "fors.fors_pmos_calib",
00407                                 0);
00408     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "cdegree");
00409     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00410     cpl_parameterlist_append(recipe->parameters, p);
00411 
00412     /*
00413      * Curvature solution interpolation
00414      */
00415  
00416     p = cpl_parameter_new_value("fors.fors_pmos_calib.cmode",
00417                                 CPL_TYPE_INT,
00418                                 "Interpolation mode of curvature solution "
00419                                 "(0 = no "
00420                                 "interpolation, 1 = fill gaps, 2 = global "
00421                                 "model)",
00422                                 "fors.fors_pmos_calib",
00423                                 1);
00424     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "cmode");
00425     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00426     cpl_parameterlist_append(recipe->parameters, p);
00427 
00428     /*
00429      * Start wavelength for spectral extraction
00430      */
00431 
00432     p = cpl_parameter_new_value("fors.fors_pmos_calib.startwavelength",
00433                                 CPL_TYPE_DOUBLE,
00434                                 "Start wavelength in spectral extraction",
00435                                 "fors.fors_pmos_calib",
00436                                 0.0);
00437     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "startwavelength");
00438     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00439     cpl_parameterlist_append(recipe->parameters, p);
00440 
00441     /*
00442      * End wavelength for spectral extraction
00443      */
00444 
00445     p = cpl_parameter_new_value("fors.fors_pmos_calib.endwavelength",
00446                                 CPL_TYPE_DOUBLE,
00447                                 "End wavelength in spectral extraction",
00448                                 "fors.fors_pmos_calib",
00449                                 0.0);
00450     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "endwavelength");
00451     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00452     cpl_parameterlist_append(recipe->parameters, p);
00453 
00454     /*
00455      * Degree of flat field fitting polynomial along dispersion direction
00456      */
00457 
00458     p = cpl_parameter_new_value("fors.fors_pmos_calib.ddegree",
00459                                 CPL_TYPE_INT,
00460                                 "Degree of flat field fitting polynomial "
00461                                 "along dispersion direction",
00462                                 "fors.fors_pmos_calib",
00463                                 -1);
00464     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "ddegree");
00465     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00466     cpl_parameterlist_append(recipe->parameters, p);
00467 
00468     /*
00469      * Smooth box radius for flat field along dispersion direction
00470      */
00471 
00472     p = cpl_parameter_new_value("fors.fors_pmos_calib.dradius",
00473                                 CPL_TYPE_INT,
00474                                 "Smooth box radius for flat field along "
00475                                 "dispersion direction",
00476                                 "fors.fors_pmos_calib",
00477                                 10);
00478     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "dradius");
00479     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00480     cpl_parameterlist_append(recipe->parameters, p);
00481 
00482     /*
00483      * Computation of QC1 parameters
00484      */
00485 
00486     p = cpl_parameter_new_value("fors.fors_pmos_calib.qc",
00487                                 CPL_TYPE_BOOL,
00488                                 "Compute QC1 parameters",
00489                                 "fors.fors_pmos_calib",
00490                                 TRUE);
00491     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "qc");
00492     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00493     cpl_parameterlist_append(recipe->parameters, p);
00494 
00495     /*
00496      * Create check products
00497      */
00498 
00499     p = cpl_parameter_new_value("fors.fors_pmos_calib.check",
00500                                 CPL_TYPE_BOOL,
00501                                 "Create intermediate products",
00502                                 "fors.fors_pmos_calib",
00503                                 FALSE);
00504     cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "check");
00505     cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
00506     cpl_parameterlist_append(recipe->parameters, p);
00507 
00508     return 0;
00509 }
00510 
00511 
00520 static int fors_pmos_calib_exec(cpl_plugin *plugin)
00521 {
00522     cpl_recipe *recipe;
00523     
00524     if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE) 
00525         recipe = (cpl_recipe *)plugin;
00526     else 
00527         return -1;
00528 
00529     return fors_pmos_calib(recipe->parameters, recipe->frames);
00530 }
00531 
00532 
00541 static int fors_pmos_calib_destroy(cpl_plugin *plugin)
00542 {
00543     cpl_recipe *recipe;
00544     
00545     if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE) 
00546         recipe = (cpl_recipe *)plugin;
00547     else 
00548         return -1;
00549 
00550     cpl_parameterlist_delete(recipe->parameters); 
00551 
00552     return 0;
00553 }
00554 
00555 
00565 static int fors_pmos_calib(cpl_parameterlist *parlist, cpl_frameset *frameset)
00566 {
00567 
00568     const char *recipe = "fors_pmos_calib";
00569 
00570 
00571     /*
00572      * Input parameters
00573      */
00574 
00575     double      dispersion;
00576     double      peakdetection;
00577     int         wdegree;
00578     int         wradius;
00579     double      wreject;
00580     const char *wcolumn;
00581     int         cdegree;
00582     int         cmode;
00583     double      startwavelength;
00584     double      endwavelength;
00585     int         ddegree;
00586     int         dradius;
00587     int         qc;
00588     int         check;
00589 
00590     /*
00591      * CPL objects
00592      */
00593 
00594     cpl_imagelist    *biases       = NULL;
00595     cpl_image        *bias         = NULL;
00596     cpl_image        *master_bias  = NULL;
00597     cpl_image        *multi_bias   = NULL;
00598     cpl_image        *flat         = NULL;
00599     cpl_image        *master_flat  = NULL;
00600     cpl_image        *smo_flat     = NULL;
00601     cpl_image        *norm_flat    = NULL;
00602     cpl_image        *spectra      = NULL;
00603     cpl_image        *wavemap      = NULL;
00604     cpl_image        *delta        = NULL;
00605     cpl_image        *residual     = NULL;
00606     cpl_image        *checkwave    = NULL;
00607     cpl_image        *rectified    = NULL;
00608     cpl_image        *dummy        = NULL;
00609     cpl_image        *refimage     = NULL;
00610     cpl_image        *coordinate   = NULL;
00611     cpl_image        *rainbow      = NULL;
00612     cpl_image        *spatial      = NULL;
00613     cpl_image        *rect_flat    = NULL;
00614     cpl_image        *rect_nflat   = NULL;
00615     cpl_image        *mapped_flat  = NULL;
00616     cpl_image        *mapped_nflat = NULL;
00617 
00618     cpl_mask         *refmask      = 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        *idscoeff_all = NULL;
00625     cpl_table        *restable     = NULL;
00626     cpl_table        *slits        = NULL;
00627     cpl_table        *positions    = NULL;
00628     cpl_table        *maskslits    = NULL;
00629     cpl_table        *traces       = NULL;
00630     cpl_table        *polytraces   = NULL;
00631     cpl_table        *restab       = NULL;
00632     cpl_table        *global       = NULL;
00633 
00634     cpl_vector       *lines        = NULL;
00635 
00636     cpl_propertylist *header       = NULL;
00637     cpl_propertylist *save_header  = NULL;
00638     cpl_propertylist *qclist       = NULL;
00639 
00640     /*
00641      * Auxiliary variables
00642      */
00643 
00644     char    version[80];
00645     const char   *arc_tag;
00646     const char   *flat_tag;
00647     const char   *master_screen_flat_tag;
00648     const char   *master_norm_flat_tag;
00649     const char   *reduced_lamp_tag;
00650     const char   *disp_residuals_tag;
00651     const char   *disp_coeff_tag;
00652     const char   *wavelength_map_tag;
00653     const char   *spectra_detection_tag;
00654     const char   *spectral_resolution_tag;
00655     const char   *slit_map_tag;
00656     const char   *curv_traces_tag;
00657     const char   *curv_coeff_tag;
00658     const char   *spatial_map_tag;
00659     const char   *slit_location_tag;
00660     const char   *master_distortion_tag = "MASTER_DISTORTION_TABLE";
00661     const char   *disp_residuals_table_tag;
00662     const char   *delta_image_tag;
00663     const char   *mapped_screen_flat_tag;
00664     const char   *mapped_norm_flat_tag;
00665     const char   *keyname;
00666     int     pmos;
00667     int     same_offset = 0;
00668     int     nslits;
00669     float  *data;
00670     double *xpos;
00671     double  mxpos;
00672     double  mean_rms;
00673     double  mean_rms_err;
00674     double  alltime;
00675     int     nflats;
00676     int     nbias;
00677     int     nlines;
00678     int     rebin;
00679     double *line;
00680     double *fiterror = NULL;
00681     int    *fitlines = NULL;
00682     int     nx, ny;
00683     double  reference;
00684     double  gain;
00685     int     ccd_xsize, ccd_ysize;
00686     int     i, j;
00687 
00688     char   *instrume = NULL;
00689     char   *pipefile = NULL;
00690 
00691 
00692     snprintf(version, 80, "%s-%s", PACKAGE, PACKAGE_VERSION);
00693 
00694     cpl_msg_set_indentation(2);
00695 
00696     if (dfs_files_dont_exist(frameset))
00697         fors_pmos_calib_exit(NULL);
00698 
00699     /* 
00700      * Get configuration parameters
00701      */
00702 
00703     cpl_msg_info(recipe, "Recipe %s configuration parameters:", recipe);
00704     cpl_msg_indent_more();
00705 
00706     if (cpl_frameset_count_tags(frameset, "GRISM_TABLE") > 1)
00707         fors_pmos_calib_exit("Too many in input: GRISM_TABLE");
00708 
00709     grism_table = dfs_load_table(frameset, "GRISM_TABLE", 1);
00710 
00711     dispersion = dfs_get_parameter_double(parlist, 
00712                     "fors.fors_pmos_calib.dispersion", grism_table);
00713 
00714     if (dispersion <= 0.0)
00715         fors_pmos_calib_exit("Invalid spectral dispersion value");
00716 
00717     peakdetection = dfs_get_parameter_double(parlist, 
00718                     "fors.fors_pmos_calib.peakdetection", grism_table);
00719     if (peakdetection <= 0.0)
00720         fors_pmos_calib_exit("Invalid peak detection level");
00721 
00722     wdegree = dfs_get_parameter_int(parlist, 
00723                     "fors.fors_pmos_calib.wdegree", grism_table);
00724 
00725     if (wdegree < 1)
00726         fors_pmos_calib_exit("Invalid polynomial degree");
00727 
00728     if (wdegree > 5)
00729         fors_pmos_calib_exit("Max allowed polynomial degree is 5");
00730 
00731     wradius = dfs_get_parameter_int(parlist, "fors.fors_pmos_calib.wradius", NULL);
00732 
00733     if (wradius < 0)
00734         fors_pmos_calib_exit("Invalid search radius");
00735 
00736     wreject = dfs_get_parameter_double(parlist, "fors.fors_pmos_calib.wreject", NULL);
00737 
00738     if (wreject <= 0.0)
00739         fors_pmos_calib_exit("Invalid rejection threshold");
00740 
00741     wcolumn = dfs_get_parameter_string(parlist, "fors.fors_pmos_calib.wcolumn", NULL);
00742 
00743     cdegree = dfs_get_parameter_int(parlist, 
00744                     "fors.fors_pmos_calib.cdegree", grism_table);
00745 
00746     if (cdegree < 1)
00747         fors_pmos_calib_exit("Invalid polynomial degree");
00748 
00749     if (cdegree > 5)
00750         fors_pmos_calib_exit("Max allowed polynomial degree is 5");
00751 
00752     cmode = dfs_get_parameter_int(parlist, "fors.fors_pmos_calib.cmode", NULL);
00753 
00754     if (cmode < 0 || cmode > 2)
00755         fors_pmos_calib_exit("Invalid curvature solution interpolation mode");
00756 
00757     startwavelength = dfs_get_parameter_double(parlist, 
00758                     "fors.fors_pmos_calib.startwavelength", grism_table);
00759     if (startwavelength > 1.0)
00760         if (startwavelength < 3000.0 || startwavelength > 13000.0)
00761             fors_pmos_calib_exit("Invalid wavelength");
00762 
00763     endwavelength = dfs_get_parameter_double(parlist, 
00764                     "fors.fors_pmos_calib.endwavelength", grism_table);
00765     if (endwavelength > 1.0) {
00766         if (endwavelength < 3000.0 || endwavelength > 13000.0)
00767             fors_pmos_calib_exit("Invalid wavelength");
00768         if (startwavelength < 1.0)
00769             fors_pmos_calib_exit("Invalid wavelength interval");
00770     }
00771 
00772     if (startwavelength > 1.0)
00773         if (endwavelength - startwavelength <= 0.0)
00774             fors_pmos_calib_exit("Invalid wavelength interval");
00775 
00776     ddegree = dfs_get_parameter_int(parlist, "fors.fors_pmos_calib.ddegree", NULL);
00777     dradius = dfs_get_parameter_int(parlist, "fors.fors_pmos_calib.dradius", NULL);
00778 
00779     if (dradius < 1)
00780         fors_pmos_calib_exit("Invalid smoothing box radius");
00781 
00782     qc = dfs_get_parameter_bool(parlist, "fors.fors_pmos_calib.qc", NULL);
00783 
00784     check = dfs_get_parameter_bool(parlist, "fors.fors_pmos_calib.check", NULL);
00785 
00786     cpl_table_delete(grism_table); grism_table = NULL;
00787 
00788     if (cpl_error_get_code())
00789         fors_pmos_calib_exit("Failure getting the configuration parameters");
00790 
00791 
00792     /* 
00793      * Check input set-of-frames
00794      */
00795 
00796     cpl_msg_indent_less();
00797     cpl_msg_info(recipe, "Check input set-of-frames:");
00798     cpl_msg_indent_more();
00799 
00800     if (!dfs_equal_keyword(frameset, "ESO INS GRIS1 ID")) 
00801         fors_pmos_calib_exit("Input frames are not from the same grism");
00802 
00803     if (!dfs_equal_keyword(frameset, "ESO INS FILT1 ID")) 
00804         fors_pmos_calib_exit("Input frames are not from the same filter");
00805 
00806     if (!dfs_equal_keyword(frameset, "ESO DET CHIP1 ID")) 
00807         fors_pmos_calib_exit("Input frames are not from the same chip");
00808 
00809     pmos = cpl_frameset_count_tags(frameset, "LAMP_PMOS");
00810 
00811     if (pmos == 0)
00812         fors_pmos_calib_exit("Missing input arc lamp frame");
00813 
00814     if (pmos) {
00815         cpl_msg_info(recipe, "PMOS data found");
00816         arc_tag                  = "LAMP_PMOS";
00817         flat_tag                 = "SCREEN_FLAT_PMOS";
00818         master_screen_flat_tag   = "MASTER_SCREEN_FLAT_PMOS";
00819         master_norm_flat_tag     = "MASTER_NORM_FLAT_PMOS";
00820         reduced_lamp_tag         = "REDUCED_LAMP_PMOS";
00821         disp_residuals_tag       = "DISP_RESIDUALS_PMOS";
00822         disp_coeff_tag           = "DISP_COEFF_PMOS";
00823         wavelength_map_tag       = "WAVELENGTH_MAP_PMOS";
00824         spectra_detection_tag    = "SPECTRA_DETECTION_PMOS";
00825         spectral_resolution_tag  = "SPECTRAL_RESOLUTION_PMOS";
00826         slit_map_tag             = "SLIT_MAP_PMOS";
00827         curv_traces_tag          = "CURV_TRACES_PMOS";
00828         curv_coeff_tag           = "CURV_COEFF_PMOS";
00829         spatial_map_tag          = "SPATIAL_MAP_PMOS";
00830         slit_location_tag        = "SLIT_LOCATION_PMOS";
00831         disp_residuals_table_tag = "DISP_RESIDUALS_TABLE_PMOS";
00832         delta_image_tag          = "DELTA_IMAGE_PMOS";
00833         mapped_screen_flat_tag   = "MAPPED_SCREEN_FLAT_PMOS";
00834         mapped_norm_flat_tag     = "MAPPED_NORM_FLAT_PMOS";
00835     }
00836 
00837     nbias = 0;
00838     if (cpl_frameset_count_tags(frameset, "MASTER_BIAS") == 0) {
00839         if (cpl_frameset_count_tags(frameset, "BIAS") == 0)
00840             fors_pmos_calib_exit("Missing required input: MASTER_BIAS or BIAS");
00841         nbias = cpl_frameset_count_tags(frameset, "BIAS");
00842     }
00843 
00844     if (cpl_frameset_count_tags(frameset, "MASTER_BIAS") > 1)
00845         fors_pmos_calib_exit("Too many in input: MASTER_BIAS");
00846 
00847     if (cpl_frameset_count_tags(frameset, "MASTER_LINECAT") == 0)
00848         fors_pmos_calib_exit("Missing required input: MASTER_LINECAT");
00849 
00850     if (cpl_frameset_count_tags(frameset, "MASTER_LINECAT") > 1)
00851         fors_pmos_calib_exit("Too many in input: MASTER_LINECAT");
00852 
00853     if (cpl_frameset_count_tags(frameset, "MASTER_LINECAT") == 0)
00854         fors_pmos_calib_exit("Missing required input: MASTER_LINECAT");
00855 
00856     if (cpl_frameset_count_tags(frameset, "MASTER_LINECAT") > 1)
00857         fors_pmos_calib_exit("Too many in input: MASTER_LINECAT");
00858 
00859     if (cpl_frameset_count_tags(frameset, "MASTER_DISTORTION_TABLE") == 0)
00860         fors_pmos_calib_exit("Missing required input: MASTER_DISTORTION_TABLE");
00861 
00862     if (cpl_frameset_count_tags(frameset, "MASTER_DISTORTION_TABLE") > 1)
00863         fors_pmos_calib_exit("Too many in input: MASTER_DISTORTION_TABLE");
00864 
00865     nflats = cpl_frameset_count_tags(frameset, flat_tag);
00866 
00867     if (nflats < 1) {
00868         cpl_msg_error(recipe, "Missing required input: %s", flat_tag);
00869         fors_pmos_calib_exit(NULL);
00870     }
00871 
00872     cpl_msg_indent_less();
00873 
00874     if (nflats > 1)
00875         cpl_msg_info(recipe, "Load %d flat field frames and sum them...",
00876                      nflats);
00877     else
00878         cpl_msg_info(recipe, "Load flat field exposure...");
00879 
00880     cpl_msg_indent_more();
00881 
00882     header = dfs_load_header(frameset, flat_tag, 0);
00883 
00884     if (header == NULL)
00885         fors_pmos_calib_exit("Cannot load flat field frame header");
00886 
00887     alltime = cpl_propertylist_get_double(header, "EXPTIME");
00888 
00889     if (cpl_error_get_code() != CPL_ERROR_NONE)
00890         fors_pmos_calib_exit("Missing keyword EXPTIME in flat field frame header");
00891 
00892     cpl_propertylist_delete(header);
00893 
00894     for (i = 1; i < nflats; i++) {
00895 
00896         header = dfs_load_header(frameset, NULL, 0);
00897 
00898         if (header == NULL)
00899             fors_pmos_calib_exit("Cannot load flat field frame header");
00900 
00901         alltime += cpl_propertylist_get_double(header, "EXPTIME");
00902 
00903         if (cpl_error_get_code() != CPL_ERROR_NONE)
00904             fors_pmos_calib_exit("Missing keyword EXPTIME in flat field "
00905                             "frame header");
00906 
00907         cpl_propertylist_delete(header);
00908 
00909     }
00910 
00911     master_flat = dfs_load_image(frameset, flat_tag, CPL_TYPE_FLOAT, 0, 0);
00912 
00913     if (master_flat == NULL)
00914         fors_pmos_calib_exit("Cannot load flat field");
00915 
00916     for (i = 1; i < nflats; i++) {
00917         flat = dfs_load_image(frameset, NULL, CPL_TYPE_FLOAT, 0, 0);
00918         if (flat) {
00919             cpl_image_add(master_flat, flat);
00920             cpl_image_delete(flat); flat = NULL;
00921         }
00922         else
00923             fors_pmos_calib_exit("Cannot load flat field");
00924     }
00925 
00926 /***
00927     if (nflats > 1)
00928         cpl_image_divide_scalar(master_flat, nflats);
00929 ***/
00930 
00931 
00932     /*
00933      * Get the reference wavelength and the rebin factor along the
00934      * dispersion direction from the arc lamp exposure
00935      */
00936 
00937     header = dfs_load_header(frameset, arc_tag, 0);
00938 
00939     if (header == NULL)
00940         fors_pmos_calib_exit("Cannot load arc lamp header");
00941 
00942     instrume = (char *)cpl_propertylist_get_string(header, "INSTRUME");
00943     if (instrume == NULL)
00944         fors_pmos_calib_exit("Missing keyword INSTRUME in arc lamp header");
00945 
00946     instrume = cpl_strdup(instrume);
00947 
00948     if (instrume[4] == '1')
00949         snprintf(version, 80, "%s/%s", "fors1", VERSION);
00950     if (instrume[4] == '2')
00951         snprintf(version, 80, "%s/%s", "fors2", VERSION);
00952 
00953     reference = cpl_propertylist_get_double(header, "ESO INS GRIS1 WLEN");
00954 
00955     if (cpl_error_get_code() != CPL_ERROR_NONE)
00956         fors_pmos_calib_exit("Missing keyword ESO INS GRIS1 WLEN in arc lamp "
00957                         "frame header");
00958 
00959     if (reference < 3000.0)   /* Perhaps in nanometers... */
00960         reference *= 10;
00961 
00962     if (reference < 3000.0 || reference > 13000.0) {
00963         cpl_msg_error(recipe, "Invalid central wavelength %.2f read from "
00964                       "keyword ESO INS GRIS1 WLEN in arc lamp frame header",
00965                       reference);
00966         fors_pmos_calib_exit(NULL);
00967     }
00968 
00969     cpl_msg_info(recipe, "The central wavelength is: %.2f", reference);
00970 
00971     rebin = cpl_propertylist_get_int(header, "ESO DET WIN1 BINX");
00972 
00973     if (cpl_error_get_code() != CPL_ERROR_NONE)
00974         fors_pmos_calib_exit("Missing keyword ESO DET WIN1 BINX in arc lamp "
00975                         "frame header");
00976 
00977     if (rebin != 1) {
00978         dispersion *= rebin;
00979         cpl_msg_warning(recipe, "The rebin factor is %d, and therefore the "
00980                         "working dispersion used is %f A/pixel", rebin, 
00981                         dispersion);
00982     }
00983 
00984     gain = cpl_propertylist_get_double(header, "ESO DET OUT1 CONAD");
00985 
00986     if (cpl_error_get_code() != CPL_ERROR_NONE)
00987         fors_pmos_calib_exit("Missing keyword ESO DET OUT1 CONAD in arc lamp "
00988                         "frame header");
00989 
00990     cpl_msg_info(recipe, "The gain factor is: %.2f e-/ADU", gain);
00991 
00992     if (pmos) {
00993         cpl_msg_info(recipe, "Produce mask slit position table...");
00994 
00995     maskslits = mos_load_slits_fors_mos(header);
00996 
00997         /*
00998          * Check if all slits have the same X offset: in such case, 
00999          * treat the observation as a long-slit one!
01000          */
01001 
01002         mxpos = cpl_table_get_column_median(maskslits, "xtop");
01003         xpos = cpl_table_get_data_double(maskslits, "xtop");
01004         nslits = cpl_table_get_nrow(maskslits);
01005 
01006         same_offset = 1;
01007         for (i = 0; i < nslits; i++) {
01008             if (fabs(mxpos-xpos[i]) > 0.01) {
01009                 same_offset = 0;
01010                 break;
01011             }
01012         }
01013 
01014         if (same_offset) {
01015             cpl_msg_info(recipe, "All slits have same offset: %.2f", mxpos);
01016         }
01017         else {
01018             cpl_msg_info(recipe, "All slits have same different offsets");
01019         }
01020     }
01021 
01022     /* Leave the header on for the next step... */
01023 
01024 
01025     /*
01026      * Remove the master bias
01027      */
01028 
01029     if (nbias) {
01030 
01031         /*
01032          * Set of raw BIASes in input, need to create master bias!
01033          */
01034 
01035         cpl_msg_info(recipe, "Generate the master from input raw biases...");
01036 
01037         if (nbias > 1) {
01038 
01039             biases = cpl_imagelist_new();
01040 
01041             bias = dfs_load_image(frameset, "BIAS", CPL_TYPE_FLOAT, 0, 0);
01042     
01043             if (bias == NULL)
01044                 fors_pmos_calib_exit("Cannot load bias frame");
01045 
01046             cpl_imagelist_set(biases, bias, 0);
01047     
01048             for (i = 1; i < nbias; i++) {
01049                 bias = dfs_load_image(frameset, NULL, CPL_TYPE_FLOAT, 0, 0);
01050                 if (bias)
01051                     cpl_imagelist_set(biases, bias, i);
01052                 else
01053                     fors_pmos_calib_exit("Cannot load bias frame");
01054             }
01055     
01056             master_bias = cpl_imagelist_collapse_median_create(biases);
01057 
01058             cpl_imagelist_delete(biases);
01059         }
01060         else {
01061             master_bias = dfs_load_image(frameset, "BIAS", 
01062                                          CPL_TYPE_FLOAT, 0, 1);
01063             if (master_bias == NULL)
01064                 fors_pmos_calib_exit("Cannot load bias");
01065         }
01066 
01067     }
01068     else {
01069         master_bias = dfs_load_image(frameset, "MASTER_BIAS", 
01070                                      CPL_TYPE_FLOAT, 0, 1);
01071         if (master_bias == NULL)
01072             fors_pmos_calib_exit("Cannot load master bias");
01073     }
01074 
01075     cpl_msg_info(recipe, "Remove the master bias...");
01076 
01077     overscans = mos_load_overscans_vimos(header, 1);
01078     cpl_propertylist_delete(header); header = NULL;
01079 
01080     if (nbias) {
01081         int xlow = cpl_table_get_int(overscans, "xlow", 0, NULL);
01082         int ylow = cpl_table_get_int(overscans, "ylow", 0, NULL);
01083         int xhig = cpl_table_get_int(overscans, "xhig", 0, NULL);
01084         int yhig = cpl_table_get_int(overscans, "yhig", 0, NULL);
01085         dummy = cpl_image_extract(master_bias, xlow+1, ylow+1, xhig, yhig);
01086         cpl_image_delete(master_bias); master_bias = dummy;
01087 
01088         if (dfs_save_image(frameset, master_bias, "MASTER_BIAS",
01089                            NULL, parlist, recipe, version))
01090             fors_pmos_calib_exit(NULL);
01091     }
01092 
01093     if (nflats > 1) {
01094         multi_bias = cpl_image_multiply_scalar_create(master_bias, nflats);
01095         dummy = mos_remove_bias(master_flat, multi_bias, overscans);
01096         cpl_image_delete(multi_bias);
01097     }
01098     else {
01099         dummy = mos_remove_bias(master_flat, master_bias, overscans);
01100     }
01101     cpl_image_delete(master_flat);
01102     master_flat = dummy;
01103 
01104     if (master_flat == NULL)
01105         fors_pmos_calib_exit("Cannot remove bias from flat field");
01106 
01107     wavelengths = dfs_load_table(frameset, "MASTER_LINECAT", 1);
01108 
01109     if (wavelengths == NULL)
01110     fors_pmos_calib_exit("Cannot load line catalog");
01111 
01112     /*
01113      * Cast the wavelengths into a (double precision) CPL vector
01114      */
01115 
01116     nlines = cpl_table_get_nrow(wavelengths);
01117 
01118     if (nlines == 0)
01119     fors_pmos_calib_exit("Empty input line catalog");
01120 
01121     if (cpl_table_has_column(wavelengths, wcolumn) != 1) {
01122     cpl_msg_error(recipe, "Missing column %s in input line catalog table",
01123               wcolumn);
01124     fors_pmos_calib_exit(NULL);
01125     }
01126 
01127     line = cpl_malloc(nlines * sizeof(double));
01128     
01129     for (i = 0; i < nlines; i++)
01130     line[i] = cpl_table_get(wavelengths, wcolumn, i, NULL);
01131 
01132     lines = cpl_vector_wrap(nlines, line);
01133 
01134     for (j = 0; j < pmos; j++) {
01135     int k;
01136 
01137     cpl_msg_indent_less();
01138     cpl_msg_info(recipe, "Processing arc lamp nb %d out of %d ...",
01139              j + 1, pmos);
01140     cpl_msg_indent_more();
01141 
01142     cpl_msg_info(recipe, "Load arc lamp exposure...");
01143     cpl_msg_indent_more();
01144 
01145     spectra = dfs_load_image(frameset, arc_tag, CPL_TYPE_FLOAT, 0, 0);
01146 
01147     /*
01148      * FIXME: Horrible workaround to avoid the problem because of the
01149      * multiple encapsulation of cpl_frameset_find() in different 
01150      * loading functions
01151      */
01152     for (k = 0; k < j; k ++) {
01153         cpl_image_delete(spectra);
01154         spectra = dfs_load_image(frameset, NULL, CPL_TYPE_FLOAT, 0, 0);
01155     }
01156 
01157     if (spectra == NULL)
01158         fors_pmos_calib_exit("Cannot load arc lamp exposure");
01159 
01160     cpl_msg_info(recipe, "Remove the master bias...");
01161 
01162     dummy = mos_remove_bias(spectra, master_bias, overscans);
01163     cpl_image_delete(spectra); spectra = dummy;
01164 
01165     if (spectra == NULL)
01166         fors_pmos_calib_exit("Cannot remove bias from arc lamp exposure");
01167 
01168     cpl_msg_indent_less();
01169     cpl_msg_info(recipe, "Load input line catalog...");
01170     cpl_msg_indent_more();
01171 
01172     /*
01173      * Here the PMOS calibration is carried out.
01174      */
01175 //
01176     if (mos_saturation_process(spectra))
01177         fors_pmos_calib_exit("Cannot process saturation");
01178 
01179     if (mos_subtract_background(spectra))
01180         fors_pmos_calib_exit("Cannot subtract the background");
01181 
01182     if (!j) {
01183         /*
01184          * Detecting spectra on the CCD
01185          */
01186 
01187         cpl_msg_indent_less();
01188         cpl_msg_info(recipe, "Detecting spectra on CCD...");
01189         cpl_msg_indent_more();
01190 
01191         ccd_xsize = nx = cpl_image_get_size_x(spectra);
01192         ccd_ysize = ny = cpl_image_get_size_y(spectra);
01193 
01194         refmask = cpl_mask_new(nx, ny);
01195 
01196         checkwave =
01197         mos_wavelength_calibration_raw(spectra, lines, dispersion, 
01198                            peakdetection, wradius, 
01199                            wdegree, wreject, reference,
01200                            &startwavelength, &endwavelength,
01201                            NULL, NULL, NULL, NULL, NULL, 
01202                            NULL, refmask);
01203 
01204         if (checkwave == NULL)
01205         fors_pmos_calib_exit("Wavelength calibration failure.");
01206 
01207         /*
01208          * Save check image to disk
01209          */
01210 
01211         header = cpl_propertylist_new();
01212         cpl_propertylist_update_double(header, "CRPIX1", 1.0);
01213         cpl_propertylist_update_double(header, "CRPIX2", 1.0);
01214         cpl_propertylist_update_double(header, "CRVAL1", 
01215                        startwavelength + dispersion/2);
01216         cpl_propertylist_update_double(header, "CRVAL2", 1.0);
01217         /* cpl_propertylist_update_double(header, "CDELT1", dispersion);
01218            cpl_propertylist_update_double(header, "CDELT2", 1.0); */
01219         cpl_propertylist_update_double(header, "CD1_1", dispersion);
01220         cpl_propertylist_update_double(header, "CD1_2", 0.0);
01221         cpl_propertylist_update_double(header, "CD2_1", 0.0);
01222         cpl_propertylist_update_double(header, "CD2_2", 1.0);
01223         cpl_propertylist_update_string(header, "CTYPE1", "LINEAR");
01224         cpl_propertylist_update_string(header, "CTYPE2", "PIXEL");
01225 
01226         if (check) {
01227         if (!j) {
01228             if(dfs_save_image_null(frameset, parlist,
01229                        spectra_detection_tag,
01230                        recipe, version)) {
01231             fors_pmos_calib_exit(NULL);
01232             }
01233         }
01234 
01235         if (dfs_save_image_ext(checkwave, spectra_detection_tag, header)) {
01236             fors_pmos_calib_exit(NULL);
01237         }
01238         }
01239 
01240         cpl_image_delete(checkwave); checkwave = NULL;
01241         cpl_propertylist_delete(header); header = NULL;
01242 
01243         if (mos_refmask_find_gaps(refmask, master_flat))
01244         fors_pmos_calib_exit("The gaps could not be found");
01245 
01246         cpl_msg_info(recipe,
01247              "Locate slits at reference wavelength on CCD...");
01248         slits = mos_locate_spectra(refmask);
01249 
01250         if (!slits) {
01251         cpl_msg_error(cpl_error_get_where(), cpl_error_get_message());
01252         fors_pmos_calib_exit("No slits could be detected!");
01253         }
01254 
01255             if (same_offset) {
01256             if (mos_check_slits(slits)) {
01257             fors_pmos_calib_exit("Some slits are missing. "
01258                                          "Cannot recover!");
01259             }
01260             }
01261 
01262         refimage = cpl_image_new_from_mask(refmask);
01263         cpl_mask_delete(refmask); refmask = NULL;
01264 
01265         if (check) {
01266         if (!j) {
01267             if(dfs_save_image_null(frameset, parlist,
01268                        slit_map_tag,
01269                        recipe, version)) {
01270             fors_pmos_calib_exit(NULL);
01271             }
01272         }
01273 
01274         save_header = dfs_load_header(frameset, arc_tag, 0);
01275 
01276         for (k = 0; k < j; k ++) {
01277             cpl_propertylist_delete(save_header);
01278             save_header = dfs_load_header(frameset, NULL, 0);
01279         }
01280 
01281         if (dfs_save_image_ext(refimage, slit_map_tag, save_header)) {
01282             fors_pmos_calib_exit(NULL);
01283         }
01284         cpl_propertylist_delete(save_header); save_header = NULL;
01285         }
01286 
01287         cpl_image_delete(refimage); refimage = NULL;
01288 
01289         if (same_offset == 0) { /* FIXLANDER: it was slit_ident */
01290 
01291         /*
01292          * Attempt slit identification: this recipe may continue even
01293          * in case of failed identification (i.e., the position table
01294          * is not produced, but an error is not set). In case of 
01295          * failure, the spectra would be still extracted, even if they
01296          * would not be associated to slits on the mask.
01297          * 
01298          * The reason for making the slit identification an user option
01299          * (via the parameter slit_ident) is to offer the possibility 
01300          * to avoid identifications that are only apparently successful
01301          * as it would happen in the case of an incorrect slit
01302          * description in the data header.
01303          */
01304 
01305         cpl_msg_indent_less();
01306         cpl_msg_info(recipe, 
01307                              "Attempt slit identification (optional)...");
01308         cpl_msg_indent_more();
01309 
01310         positions = mos_identify_slits(slits, maskslits, NULL);
01311 
01312         if (positions) {
01313             cpl_table_delete(slits);
01314             slits = positions;
01315 
01316             /*
01317              * Eliminate slits which are _entirely_ outside the CCD
01318              */
01319 
01320             cpl_table_and_selected_double(slits, "ybottom", 
01321                                                   CPL_GREATER_THAN, ny-1);
01322             cpl_table_or_selected_double(slits, "ytop", 
01323                                                   CPL_LESS_THAN, 0);
01324             cpl_table_erase_selected(slits);
01325 
01326             nslits = cpl_table_get_nrow(slits);
01327 
01328             if (nslits == 0)
01329             fors_pmos_calib_exit("No slits found on the CCD");
01330 
01331             cpl_msg_info(recipe, "%d slits are entirely or partially "
01332                  "contained in CCD", nslits);
01333 
01334         }
01335         else {
01336             same_offset = 1; /* FIXLANDER slit_ident = 0; */
01337             cpl_msg_info(recipe, "Global distortion model cannot be computed");
01338             if (cpl_error_get_code() != CPL_ERROR_NONE) {
01339             fors_pmos_calib_exit(NULL);
01340             }
01341         }
01342         }
01343 
01344 
01345         /*
01346          * Determination of spectral curvature
01347          */
01348 
01349         cpl_msg_indent_less();
01350         cpl_msg_info(recipe, "Determining spectral curvature...");
01351         cpl_msg_indent_more();
01352 
01353         cpl_msg_info(recipe, "Tracing master flat field spectra edges...");
01354         traces = mos_trace_flat(master_flat, slits, reference, 
01355                     startwavelength, endwavelength, dispersion);
01356 
01357         if (!traces)
01358         fors_pmos_calib_exit("Tracing failure");
01359 
01360         cpl_msg_info(recipe, "Fitting flat field spectra edges...");
01361         polytraces = mos_poly_trace(slits, traces, cdegree);
01362 
01363         if (!polytraces)
01364         fors_pmos_calib_exit("Trace fitting failure");
01365 
01366         if (cmode) {
01367         cpl_msg_info(recipe, 
01368                              "Computing global spectral curvature model...");
01369         mos_global_trace(slits, polytraces, cmode);
01370         }
01371 
01372 /* Old saving & appending:
01373 
01374         if (!j) {
01375         if (dfs_save_table(frameset, traces, curv_traces_tag, NULL,
01376                    parlist, recipe, version)) {
01377             fors_pmos_calib_exit(NULL);
01378         }
01379         } else {
01380         if (dfs_save_table_ext(traces, curv_traces_tag, NULL)) {
01381             fors_pmos_calib_exit(NULL);
01382         }
01383         }
01384 End of old saving & appending */
01385 
01386         if (!j) {
01387                 if(dfs_save_image_null(frameset, parlist, curv_traces_tag,
01388                                        recipe, version)) {
01389                     fors_pmos_calib_exit(NULL);
01390                 }
01391             }
01392 
01393             if (dfs_save_table_ext(traces, curv_traces_tag, NULL)) {
01394                 fors_pmos_calib_exit(NULL);
01395             }
01396 
01397         cpl_table_delete(traces); traces = NULL;
01398 
01399         coordinate = cpl_image_new(nx, ny, CPL_TYPE_FLOAT);
01400 
01401     }
01402 //
01403     spatial = mos_spatial_calibration(spectra, slits, polytraces,
01404                       reference, 
01405                       startwavelength, endwavelength, 
01406                       dispersion, 0, j ? NULL: coordinate);
01407 
01408     if (!j) {
01409 //
01410         if (same_offset) { /* FIXLANDER It was !slit_ident */
01411         cpl_image_delete(spectra); spectra = NULL;
01412         }
01413 
01414         /*
01415          * Flat field normalisation is done directly on the master flat
01416          * field (without spatial rectification first). The spectral
01417          * curvature model may be provided in input, in future releases.
01418          */
01419 
01420         cpl_msg_indent_less();
01421         cpl_msg_info(recipe, "Perform flat field normalisation...");
01422         cpl_msg_indent_more();
01423 
01424         norm_flat = cpl_image_duplicate(master_flat);
01425 
01426         smo_flat = mos_normalise_flat(norm_flat, coordinate, slits,
01427                       polytraces, reference,
01428                       startwavelength, endwavelength,
01429                       dispersion, dradius, ddegree);
01430 
01431         /* This may be a product */
01432         cpl_image_delete(smo_flat); smo_flat = NULL; 
01433 
01434  
01435         save_header = dfs_load_header(frameset, flat_tag, 0);
01436         cpl_propertylist_update_int(save_header, "ESO PRO DATANCOM",
01437                     nflats);
01438 
01439         rect_flat = mos_spatial_calibration(master_flat, slits, polytraces,
01440                         reference, startwavelength, 
01441                         endwavelength, dispersion, 0,
01442                         NULL);
01443         rect_nflat = mos_spatial_calibration(norm_flat, slits, polytraces, 
01444                          reference, startwavelength, 
01445                          endwavelength, dispersion, 0,
01446                          NULL);
01447 
01448 
01449         if (dfs_save_image(frameset, master_flat, master_screen_flat_tag,
01450                    save_header, parlist, recipe, version))
01451         fors_pmos_calib_exit(NULL);
01452 
01453 
01454         if (dfs_save_image(frameset, norm_flat, master_norm_flat_tag,
01455                    save_header, parlist, recipe, version))
01456         fors_pmos_calib_exit(NULL);
01457 
01458         cpl_image_delete(norm_flat); norm_flat = NULL;
01459         cpl_propertylist_delete(save_header); save_header = NULL;
01460 
01461     }
01462 //
01463 
01464     /*
01465      * Final wavelength calibration of spectra having their curvature
01466      * removed
01467      */
01468 
01469     cpl_msg_indent_less();
01470     cpl_msg_info(recipe, "Perform final wavelength calibration...");
01471     cpl_msg_indent_more();
01472 
01473     nx = cpl_image_get_size_x(spatial);
01474     ny = cpl_image_get_size_y(spatial);
01475 
01476     idscoeff = cpl_table_new(ny);
01477     restable = cpl_table_new(nlines);
01478     rainbow = cpl_image_new(nx, ny, CPL_TYPE_FLOAT);
01479     if (check)
01480         residual = cpl_image_new(nx, ny, CPL_TYPE_FLOAT);
01481     fiterror = cpl_calloc(ny, sizeof(double));
01482     fitlines = cpl_calloc(ny, sizeof(int));
01483 
01484     rectified = mos_wavelength_calibration_final(spatial, slits, lines, 
01485                              dispersion, peakdetection,
01486                              wradius, wdegree, wreject,
01487                              reference,
01488                              &startwavelength, 
01489                              &endwavelength, fitlines, 
01490                              fiterror, idscoeff,
01491                              rainbow, 
01492                              residual, restable);
01493 
01494 /*
01495   dfs_save_image(frameset, rainbow, "rainbow_calib", NULL, 
01496                  parlist, recipe, version);
01497 */
01498 
01499     if (rectified == NULL)
01500         fors_pmos_calib_exit("Wavelength calibration failure.");
01501 
01502 /* Old saving & appending:
01503 
01504     if (!j) {
01505         if (dfs_save_table(frameset, restable, disp_residuals_table_tag,
01506                    NULL, parlist, recipe, version)) {
01507         fors_pmos_calib_exit(NULL);
01508         }
01509     } else {
01510         if (dfs_save_table_ext(restable, disp_residuals_table_tag, NULL)) {
01511         fors_pmos_calib_exit(NULL);
01512         }
01513     }
01514 
01515 End of old saving & appending */
01516 
01517         if (!j) {
01518             if(dfs_save_image_null(frameset, parlist, disp_residuals_table_tag,
01519                                    recipe, version)) {
01520                 fors_pmos_calib_exit(NULL);
01521             }
01522         }
01523 
01524     header = dfs_load_header(frameset, arc_tag, 0);
01525 
01526     for (k = 0; k < j; k ++) {
01527         cpl_propertylist_delete(header);
01528         header = dfs_load_header(frameset, NULL, 0);
01529     }
01530 
01531         if (dfs_save_table_ext(restable, disp_residuals_table_tag, header)) {
01532             fors_pmos_calib_exit(NULL);
01533         }
01534 
01535     cpl_propertylist_delete(header);
01536 
01537     cpl_table_delete(restable); restable = NULL;
01538 
01539     cpl_table_wrap_double(idscoeff, fiterror, "error"); fiterror = NULL;
01540     cpl_table_set_column_unit(idscoeff, "error", "pixel");
01541     cpl_table_wrap_int(idscoeff, fitlines, "nlines"); fitlines = NULL;
01542 
01543     for (i = 0; i < ny; i++)
01544         if (!cpl_table_is_valid(idscoeff, "c0", i))
01545         cpl_table_set_invalid(idscoeff, "error", i);
01546 
01547     delta = mos_map_pixel(idscoeff, reference, startwavelength,
01548                   endwavelength, dispersion, 2);
01549 
01550     header = dfs_load_header(frameset, arc_tag, 0);
01551 
01552     for (k = 0; k < j; k ++) {
01553         cpl_propertylist_delete(header);
01554         header = dfs_load_header(frameset, NULL, 0);
01555     }
01556 
01557     cpl_propertylist_update_double(header, "CRPIX1", 1.0);
01558     cpl_propertylist_update_double(header, "CRPIX2", 1.0);
01559     cpl_propertylist_update_double(header, "CRVAL1",
01560                        startwavelength + dispersion/2);
01561     cpl_propertylist_update_double(header, "CRVAL2", 1.0);
01562     /* cpl_propertylist_update_double(header, "CDELT1", dispersion);
01563        cpl_propertylist_update_double(header, "CDELT2", 1.0); */
01564     cpl_propertylist_update_double(header, "CD1_1", dispersion);
01565     cpl_propertylist_update_double(header, "CD1_2", 0.0);
01566     cpl_propertylist_update_double(header, "CD2_1", 0.0);
01567     cpl_propertylist_update_double(header, "CD2_2", 1.0);
01568     cpl_propertylist_update_string(header, "CTYPE1", "LINEAR");
01569     cpl_propertylist_update_string(header, "CTYPE2", "PIXEL");
01570 
01571     if (!j) {
01572         if(dfs_save_image_null(frameset, parlist, delta_image_tag,
01573                    recipe, version)) {
01574         fors_pmos_calib_exit(NULL);
01575         }
01576     }
01577 
01578     if (dfs_save_image_ext(delta, delta_image_tag, header)) {
01579         fors_pmos_calib_exit(NULL);
01580     }
01581 
01582     cpl_image_delete(delta); delta = NULL;
01583     cpl_propertylist_delete(header); header = NULL;
01584 
01585     mean_rms = mos_distortions_rms(rectified, lines, startwavelength, 
01586                        dispersion, 6, 0);
01587 
01588     cpl_msg_info(recipe, "Mean residual: %f pixel", mean_rms);
01589 
01590     mean_rms = cpl_table_get_column_mean(idscoeff, "error");
01591     mean_rms_err = cpl_table_get_column_stdev(idscoeff, "error");
01592 
01593     cpl_msg_info(recipe, "Mean model accuracy: %f pixel (%f A)", 
01594              mean_rms, mean_rms * dispersion);
01595 
01596     restab = mos_resolution_table(rectified, startwavelength, dispersion, 
01597                       60000, lines);
01598 
01599     if (restab) {
01600         cpl_msg_info(recipe, "Mean spectral resolution: %.2f", 
01601              cpl_table_get_column_mean(restab, "resolution"));
01602         cpl_msg_info(recipe,
01603                   "Mean reference lines FWHM: %.2f +/- %.2f pixel",
01604                   cpl_table_get_column_mean(restab, "fwhm") / dispersion,
01605                   cpl_table_get_column_mean(restab, "fwhm_rms") / dispersion);
01606 
01607         if (qc) {
01608 
01609         header = dfs_load_header(frameset, arc_tag, 0)
01610 
01611 ;       for (k = 0; k < j; k ++) {
01612             cpl_propertylist_delete(header);
01613             header = dfs_load_header(frameset, NULL, 0);
01614         }
01615 
01616         if (header == NULL)
01617             fors_pmos_calib_exit("Cannot reload arc lamp header");
01618 
01619         qclist = cpl_propertylist_new();
01620 
01621         fors_qc_start_group(qclist, "2.0", instrume);
01622 
01623 
01624         /*
01625          * QC1 group header
01626          */
01627 
01628         if (fors_qc_write_string("PRO.CATG", spectral_resolution_tag,
01629                      "Product category", instrume))
01630             fors_pmos_calib_exit("Cannot write product category to "
01631                      "QC log file");
01632 
01633         if (fors_qc_keyword_to_paf(header, "ESO DPR TYPE", NULL,
01634                        "DPR type", instrume))
01635             fors_pmos_calib_exit("Missing keyword DPR TYPE in arc "
01636                      "lamp header");
01637 
01638         if (fors_qc_keyword_to_paf(header, "ESO TPL ID", NULL,
01639                        "Template", instrume))
01640             fors_pmos_calib_exit("Missing keyword TPL ID in arc "
01641                      "lamp header");
01642 
01643         if (fors_qc_keyword_to_paf(header, "ESO INS GRIS1 NAME", NULL,
01644                        "Grism name", instrume))
01645             fors_pmos_calib_exit("Missing keyword INS GRIS1 NAME in "
01646                      "arc lamp header");
01647 
01648         if (fors_qc_keyword_to_paf(header, "ESO INS GRIS1 ID", NULL,
01649                        "Grism identifier", instrume))
01650             fors_pmos_calib_exit("Missing keyword INS GRIS1 ID in arc "
01651                      "lamp header");
01652 
01653         if (cpl_propertylist_has(header, "ESO INS FILT1 NAME"))
01654             fors_qc_keyword_to_paf(header, "ESO INS FILT1 NAME", NULL,
01655                        "Filter name", instrume);
01656 
01657         if (fors_qc_keyword_to_paf(header, "ESO INS COLL NAME", NULL,
01658                        "Collimator name", instrume))
01659             fors_pmos_calib_exit("Missing keyword INS COLL NAME in arc "
01660                      "lamp header");
01661 
01662         if (fors_qc_keyword_to_paf(header, "ESO DET CHIP1 ID", NULL,
01663                        "Chip identifier", instrume))
01664             fors_pmos_calib_exit("Missing keyword DET CHIP1 ID in arc "
01665                      "lamp header");
01666 
01667         if (fors_qc_keyword_to_paf(header, "ARCFILE", NULL,
01668                        "Archive name of input data", 
01669                        instrume))
01670             fors_pmos_calib_exit("Missing keyword ARCFILE in arc "
01671                      "lamp header");
01672 
01673         cpl_propertylist_delete(header); header = NULL;
01674 
01675         pipefile = dfs_generate_filename_tfits(spectral_resolution_tag);
01676         if (fors_qc_write_string("PIPEFILE", pipefile,
01677                      "Pipeline product name", instrume))
01678             fors_pmos_calib_exit("Cannot write PIPEFILE to QC logfile");
01679         cpl_free(pipefile); pipefile = NULL;
01680 
01681 
01682         /*
01683          * QC1 parameters
01684          */
01685 
01686         keyname = "QC.PMOS.RESOLUTION";
01687 
01688         if (fors_qc_write_qc_double(qclist, 
01689                         cpl_table_get_column_mean(restab,
01690                                             "resolution"),
01691                         keyname,
01692                         "Angstrom",
01693                         "Mean spectral resolution",
01694                         instrume)) {
01695             fors_pmos_calib_exit("Cannot write mean spectral "
01696                                          "resolution to QC log file");
01697         }
01698 
01699         keyname = "QC.PMOS.RESOLUTION.RMS";
01700 
01701         if (fors_qc_write_qc_double(qclist, 
01702                         cpl_table_get_column_stdev(restab, 
01703                                             "resolution"),
01704                         keyname,
01705                         "Angstrom", 
01706                         "Scatter of spectral resolution",
01707                         instrume)) {
01708             fors_pmos_calib_exit("Cannot write spectral resolution "
01709                      "scatter to QC log file");
01710         }
01711 
01712         keyname = "QC.PMOS.RESOLUTION.NWAVE";
01713 
01714         if (fors_qc_write_qc_int(qclist, cpl_table_get_nrow(restab) -
01715                      cpl_table_count_invalid(restab, 
01716                                  "resolution"),
01717                      keyname,
01718                      NULL,
01719                      "Number of examined wavelengths "
01720                      "for resolution computation",
01721                      instrume)) {
01722             fors_pmos_calib_exit("Cannot write number of lines used "
01723                      "in spectral resolution computation "
01724                      "to QC log file");
01725         }
01726 
01727         keyname = "QC.PMOS.RESOLUTION.MEANRMS";
01728             
01729         if (fors_qc_write_qc_double(qclist,
01730                                             cpl_table_get_column_mean(restab,
01731                                             "resolution_rms"),
01732                         keyname, NULL,
01733                         "Mean error on spectral "
01734                         "resolution computation",
01735                         instrume)) {
01736             fors_pmos_calib_exit("Cannot write mean error in "
01737                      "spectral resolution computation "
01738                      "to QC log file");
01739         }
01740 
01741         keyname = "QC.PMOS.RESOLUTION.NLINES";
01742 
01743         if (fors_qc_write_qc_int(qclist,
01744                      cpl_table_get_column_mean(restab, 
01745                                                                    "nlines") *
01746                      cpl_table_get_nrow(restab),
01747                      keyname, NULL,
01748                      "Number of lines for spectral "
01749                      "resolution computation",
01750                      instrume)) {
01751             fors_pmos_calib_exit("Cannot write number of examined "
01752                      "wavelengths in spectral resolution "
01753                      "computation to QC log file");
01754         }
01755 
01756         fors_qc_end_group();
01757 
01758         }
01759 
01760 /* Old saving & appending:
01761 
01762         if (!j) {
01763         if (dfs_save_table(frameset, restab, spectral_resolution_tag, 
01764                                    qclist, parlist, recipe, version)) {
01765             fors_pmos_calib_exit(NULL);
01766         }
01767         } else {
01768         if (dfs_save_table_ext(restab, spectral_resolution_tag,
01769                        qclist)) {
01770             fors_pmos_calib_exit(NULL);
01771         }
01772         }
01773 
01774 End of old saving & appending */
01775 
01776             if (!j) {
01777                 if(dfs_save_image_null(frameset, parlist, 
01778                                        spectral_resolution_tag,
01779                                        recipe, version)) {
01780                     fors_pmos_calib_exit(NULL);
01781                 }
01782             }
01783 
01784         header = dfs_load_header(frameset, arc_tag, 0);
01785 
01786         for (k = 0; k < j; k ++) {
01787         cpl_propertylist_delete(header);
01788         header = dfs_load_header(frameset, NULL, 0);
01789         }
01790 
01791         cpl_propertylist_append(header, qclist);
01792 
01793             if (dfs_save_table_ext(restab, spectral_resolution_tag, header)) {
01794                 fors_pmos_calib_exit(NULL);
01795             }
01796 
01797         cpl_table_delete(restab); restab = NULL;
01798         cpl_propertylist_delete(qclist); qclist = NULL;
01799         cpl_propertylist_delete(header); header = NULL;
01800 
01801     }
01802     else
01803         fors_pmos_calib_exit("Cannot compute the spectral resolution table");
01804 
01805 /* Old saving & appending:
01806 
01807     if (!j) {
01808         if (dfs_save_table(frameset, idscoeff, disp_coeff_tag, NULL,
01809                    parlist, recipe, version)) {
01810         fors_pmos_calib_exit(NULL);
01811         }
01812     } else {
01813         if (dfs_save_table_ext(idscoeff, disp_coeff_tag, NULL)) {
01814         fors_pmos_calib_exit(NULL);
01815         }
01816     }
01817 
01818 End of old saving & appending */
01819 
01820         if (!j) {
01821             if(dfs_save_image_null(frameset, parlist, disp_coeff_tag,
01822                                    recipe, version)) {
01823                 fors_pmos_calib_exit(NULL);
01824             }
01825         }
01826 
01827     header = dfs_load_header(frameset, arc_tag, 0);
01828 
01829     for (k = 0; k < j; k ++) {
01830         cpl_propertylist_delete(header);
01831         header = dfs_load_header(frameset, NULL, 0);
01832     }
01833 
01834         if (dfs_save_table_ext(idscoeff, disp_coeff_tag, header)) {
01835             fors_pmos_calib_exit(NULL);
01836         }
01837 
01838     cpl_propertylist_delete(header);
01839 
01840     if (!j) {
01841         mapped_flat = mos_wavelength_calibration(rect_flat, reference,
01842                              startwavelength, 
01843                                                      endwavelength,
01844                              dispersion, idscoeff, 0);
01845 
01846         mapped_nflat = mos_wavelength_calibration(rect_nflat, reference,
01847                               startwavelength, 
01848                                                       endwavelength,
01849                               dispersion, idscoeff, 0);
01850 
01851         cpl_image_delete(rect_flat); rect_flat = NULL;
01852         cpl_image_delete(rect_nflat); rect_nflat = NULL;
01853     }
01854 
01855         /* Global removed */
01856 
01857     cpl_table_delete(idscoeff); idscoeff = NULL;
01858 
01859     header = dfs_load_header(frameset, arc_tag, 0);
01860 
01861     for (k = 0; k < j; k ++) {
01862         cpl_propertylist_delete(header);
01863         header = dfs_load_header(frameset, NULL, 0);
01864     }
01865 
01866     cpl_propertylist_update_double(header, "CRPIX1", 1.0);
01867     cpl_propertylist_update_double(header, "CRPIX2", 1.0);
01868     cpl_propertylist_update_double(header, "CRVAL1", 
01869                        startwavelength + dispersion/2);
01870     cpl_propertylist_update_double(header, "CRVAL2", 1.0);
01871     /* cpl_propertylist_update_double(header, "CDELT1", dispersion);
01872        cpl_propertylist_update_double(header, "CDELT2", 1.0); */
01873     cpl_propertylist_update_double(header, "CD1_1", dispersion);
01874     cpl_propertylist_update_double(header, "CD1_2", 0.0);
01875     cpl_propertylist_update_double(header, "CD2_1", 0.0);
01876     cpl_propertylist_update_double(header, "CD2_2", 1.0);
01877     cpl_propertylist_update_string(header, "CTYPE1", "LINEAR");
01878     cpl_propertylist_update_string(header, "CTYPE2", "PIXEL");
01879     cpl_propertylist_update_int(header, "ESO PRO DATANCOM", 1);
01880 
01881     if (!j) {
01882         if(dfs_save_image_null(frameset, parlist, reduced_lamp_tag,
01883                    recipe, version)) {
01884         fors_pmos_calib_exit(NULL);
01885         }
01886     }
01887 
01888     if (dfs_save_image_ext(rectified, reduced_lamp_tag, header)) {
01889         fors_pmos_calib_exit(NULL);
01890     }
01891 
01892     cpl_image_delete(rectified); rectified = NULL;
01893 
01894     cpl_propertylist_update_int(header, "ESO PRO DATANCOM", nflats);
01895 
01896     if (!j) {
01897         if (dfs_save_image(frameset, mapped_flat, mapped_screen_flat_tag,
01898                    header, parlist, recipe, version))
01899         fors_pmos_calib_exit(NULL);
01900         cpl_image_delete(mapped_flat); mapped_flat = NULL;
01901 
01902         if (dfs_save_image(frameset, mapped_nflat, mapped_norm_flat_tag,
01903                    header, parlist, recipe, version))
01904         fors_pmos_calib_exit(NULL);
01905         cpl_image_delete(mapped_nflat); mapped_nflat = NULL;
01906     }
01907 
01908     cpl_propertylist_delete(header); header = NULL;
01909 
01910     if (check) {
01911         save_header = dfs_load_header(frameset, arc_tag, 0);
01912         for (k = 0; k < j; k ++) {
01913         cpl_propertylist_delete(save_header);
01914         save_header = dfs_load_header(frameset, NULL, 0);
01915         }
01916 
01917         cpl_propertylist_update_double(save_header, "CRPIX2", 1.0);
01918         cpl_propertylist_update_double(save_header, "CRVAL2", 1.0);
01919         /* cpl_propertylist_update_double(save_header, "CDELT2", 1.0); */
01920         cpl_propertylist_update_double(save_header, "CD1_1", 1.0);
01921         cpl_propertylist_update_double(save_header, "CD1_2", 0.0);
01922         cpl_propertylist_update_double(save_header, "CD2_1", 0.0);
01923         cpl_propertylist_update_double(save_header, "CD2_2", 1.0);
01924         cpl_propertylist_update_string(save_header, "CTYPE1", "LINEAR");
01925         cpl_propertylist_update_string(save_header, "CTYPE2", "PIXEL");
01926 
01927         if (!j) {
01928         if(dfs_save_image_null(frameset, parlist, disp_residuals_tag,
01929                        recipe, version)) {
01930             fors_pmos_calib_exit(NULL);
01931         }
01932         }
01933 
01934         if (dfs_save_image_ext(residual, disp_residuals_tag, save_header)) {
01935         fors_pmos_calib_exit(NULL);
01936         }
01937 
01938         cpl_image_delete(residual); residual = NULL;
01939         cpl_propertylist_delete(save_header); save_header = NULL;
01940     }
01941 
01942     wavemap = mos_map_wavelengths(coordinate, rainbow, slits, polytraces, 
01943                       reference, startwavelength, endwavelength, 
01944                       dispersion);
01945 
01946     cpl_image_delete(rainbow); rainbow = NULL;
01947 
01948     save_header = dfs_load_header(frameset, arc_tag, 0);
01949 
01950     for (k = 0; k < j; k ++) {
01951         cpl_propertylist_delete(save_header);
01952         save_header = dfs_load_header(frameset, NULL, 0);
01953     }
01954 
01955     if (qc) {
01956         fors_qc_start_group(save_header, "2.0", instrume);
01957 
01958         /*
01959          * QC1 group header
01960          */
01961 
01962         if (fors_qc_write_string("PRO.CATG", wavelength_map_tag,
01963                      "Product category", instrume))
01964         fors_pmos_calib_exit("Cannot write product category to "
01965                      "QC log file");
01966 
01967         if (fors_qc_keyword_to_paf(save_header, "ESO DPR TYPE", NULL,
01968                        "DPR type", instrume))
01969         fors_pmos_calib_exit("Missing keyword DPR TYPE in arc "
01970                      "lamp header");
01971 
01972         if (fors_qc_keyword_to_paf(save_header, "ESO TPL ID", NULL,
01973                        "Template", instrume))
01974         fors_pmos_calib_exit("Missing keyword TPL ID in arc "
01975                      "lamp header");
01976 
01977         if (fors_qc_keyword_to_paf(save_header, "ESO INS GRIS1 NAME", NULL,
01978                        "Grism name", instrume))
01979         fors_pmos_calib_exit("Missing keyword INS GRIS1 NAME in arc "
01980                      "lamp header");
01981 
01982         if (fors_qc_keyword_to_paf(save_header, "ESO INS GRIS1 ID", NULL,
01983                        "Grism identifier", instrume))
01984         fors_pmos_calib_exit("Missing keyword INS GRIS1 ID in arc "
01985                      "lamp header");
01986 
01987         if (cpl_propertylist_has(save_header, "ESO INS FILT1 NAME"))
01988         fors_qc_keyword_to_paf(save_header, "ESO INS FILT1 NAME", NULL,
01989                        "Filter name", instrume);
01990 
01991         if (fors_qc_keyword_to_paf(save_header, "ESO INS COLL NAME", NULL,
01992                        "Collimator name", instrume))
01993         fors_pmos_calib_exit("Missing keyword INS COLL NAME in arc "
01994                      "lamp header");
01995 
01996         if (fors_qc_keyword_to_paf(save_header, "ESO DET CHIP1 ID", NULL,
01997                        "Chip identifier", instrume))
01998         fors_pmos_calib_exit("Missing keyword DET CHIP1 ID in arc "
01999                      "lamp header");
02000 
02001         if (fors_qc_keyword_to_paf(save_header, "ESO DET WIN1 BINX", NULL,
02002                        "Binning factor along X", instrume))
02003         fors_pmos_calib_exit("Missing keyword ESO DET WIN1 BINX "
02004                      "in frame header");
02005 
02006         if (fors_qc_keyword_to_paf(save_header, "ESO DET WIN1 BINY", NULL,
02007                        "Binning factor along Y", instrume))
02008         fors_pmos_calib_exit("Missing keyword ESO DET WIN1 BINY "
02009                      "in frame header");
02010 
02011         if (fors_qc_keyword_to_paf(save_header, "ARCFILE", NULL,
02012                        "Archive name of input data",
02013                        instrume))
02014         fors_pmos_calib_exit("Missing keyword ARCFILE in arc "
02015                      "lamp header");
02016 
02017         pipefile = dfs_generate_filename(wavelength_map_tag);
02018         if (fors_qc_write_string("PIPEFILE", pipefile,
02019                      "Pipeline product name", instrume))
02020         fors_pmos_calib_exit("Cannot write PIPEFILE to QC log file");
02021         cpl_free(pipefile); pipefile = NULL;
02022 
02023 
02024         /*
02025          * QC1 parameters
02026          */
02027 
02028         if (fors_qc_write_qc_double(save_header,
02029                     mean_rms,
02030                     "QC.WAVE.ACCURACY",
02031                     "pixel",
02032                     "Mean accuracy of wavecalib model",
02033                     instrume)) {
02034         fors_pmos_calib_exit("Cannot write mean wavelength calibration "
02035                      "accuracy to QC log file");
02036         }
02037 
02038 
02039         if (fors_qc_write_qc_double(save_header,
02040                     mean_rms_err,
02041                     "QC.WAVE.ACCURACY.ERROR",
02042                     "pixel",
02043                     "Error on accuracy of wavecalib model",
02044                     instrume)) {
02045         fors_pmos_calib_exit("Cannot write error on wavelength "
02046                      "calibration accuracy to QC log file");
02047         }
02048 
02049             if (same_offset && fabs(mxpos) < 0.05) { 
02050                                              /* Only if same offset is 0.0 */
02051 
02052                 data = cpl_image_get_data(wavemap);
02053 
02054                 if (fors_qc_write_qc_double(save_header,
02055                                            data[nx/2 + ccd_ysize*nx/2],
02056                                            "QC.PMOS.CENTRAL.WAVELENGTH",
02057                                            "Angstrom",
02058                                            "Wavelength at CCD center",
02059                                            instrume)) {
02060                     fors_pmos_calib_exit("Cannot write central wavelength "
02061                                          "to QC log file");
02062                 }
02063             }
02064 
02065         fors_qc_end_group();
02066 
02067     }
02068 
02069     if (!j) {
02070         if(dfs_save_image_null(frameset, parlist, wavelength_map_tag,
02071                    recipe, version)) {
02072         fors_pmos_calib_exit(NULL);
02073         }
02074     }
02075 
02076     if (dfs_save_image_ext(wavemap, wavelength_map_tag, save_header)) {
02077         fors_pmos_calib_exit(NULL);
02078     }
02079 
02080     cpl_image_delete(wavemap); wavemap = NULL;
02081 
02082     cpl_propertylist_erase_regexp(save_header, "^ESO QC ", 0);
02083 
02084     cpl_propertylist_delete(save_header); save_header = NULL;
02085 
02086     cpl_msg_indent_less();
02087 
02088     }
02089 
02090     if (dfs_save_image(frameset, coordinate, spatial_map_tag, save_header,
02091                parlist, recipe, version))
02092     fors_pmos_calib_exit(NULL);
02093 
02094     cpl_image_delete(coordinate); coordinate = NULL;
02095     cpl_propertylist_delete(save_header); save_header = NULL;
02096 
02097     header = NULL;
02098 
02099     if (qc) {
02100 
02101         double maxpos, maxneg, maxcurve, maxslope;
02102 
02103         header = dfs_load_header(frameset, arc_tag, 0);
02104 
02105         fors_qc_start_group(header, "2.0", instrume);
02106 
02107         /*
02108          * QC1 group header
02109          */
02110 
02111         if (fors_qc_write_string("PRO.CATG", curv_coeff_tag,
02112                                 "Product category", instrume))
02113             fors_pmos_calib_exit("Cannot write product category to "
02114                                  "QC log file");
02115 
02116         if (fors_qc_keyword_to_paf(header, "ESO DPR TYPE", NULL,
02117                                   "DPR type", instrume))
02118             fors_pmos_calib_exit("Missing keyword DPR TYPE in arc "
02119                                  "lamp header");
02120 
02121         if (fors_qc_keyword_to_paf(header, "ESO TPL ID", NULL,
02122                                   "Template", instrume))
02123             fors_pmos_calib_exit("Missing keyword TPL ID in arc "
02124                                  "lamp header");
02125         if (fors_qc_keyword_to_paf(header, "ARCFILE", NULL,
02126                                   "Archive name of input data",
02127                                   instrume))
02128             fors_pmos_calib_exit("Missing keyword ARCFILE in arc "
02129                                  "lamp header");
02130 
02131         pipefile = dfs_generate_filename(curv_coeff_tag);
02132         if (fors_qc_write_string("PIPEFILE", pipefile,
02133                                 "Pipeline product name", instrume))
02134             fors_pmos_calib_exit("Cannot write PIPEFILE to QC log file");
02135         cpl_free(pipefile); pipefile = NULL;
02136 
02137         /*
02138          * QC1 parameters
02139          */
02140 
02141         maxpos = fabs(cpl_table_get_column_max(polytraces, "c2"));
02142         maxneg = fabs(cpl_table_get_column_min(polytraces, "c2"));
02143         maxcurve = maxpos > maxneg ? maxpos : maxneg;
02144         if (fors_qc_write_qc_double(header,
02145                                    maxcurve,
02146                                    "QC.TRACE.MAX.CURVATURE",
02147                                    "Y pixel / X pixel ^2",
02148                                    "Max observed curvature in spectral tracing",
02149                                    instrume)) {
02150             fors_pmos_calib_exit("Cannot write max observed curvature in "
02151                                  "spectral tracing to QC log file");
02152         }
02153 
02154         maxpos = fabs(cpl_table_get_column_max(polytraces, "c1"));
02155         maxneg = fabs(cpl_table_get_column_min(polytraces, "c1"));
02156         maxslope = maxpos > maxneg ? maxpos : maxneg;
02157 
02158         if (fors_qc_write_qc_double(header,
02159                                    maxslope,
02160                                    "QC.TRACE.MAX.SLOPE",
02161                                    "Y pixel / X pixel",
02162                                    "Max observed slope in spectral tracing",
02163                                    instrume)) {
02164             fors_pmos_calib_exit("Cannot write max observed slope in spectral "
02165                                  "tracing to QC log file");
02166         }
02167 
02168         fors_qc_end_group();
02169     }
02170 
02171     if (dfs_save_table(frameset, polytraces, curv_coeff_tag, header,
02172                parlist, recipe, version)) {
02173     fors_pmos_calib_exit(NULL);
02174     }
02175 
02176     cpl_propertylist_delete(header); header = NULL;
02177     cpl_table_delete(polytraces); polytraces = NULL;
02178 
02179     /* FIXLANDER It was slit_ident == 0 and 
02180        it was in a different place above in the code */
02181 
02182     if (same_offset) {
02183     cpl_table * globaltbl = 
02184         dfs_load_table(frameset, master_distortion_tag, 1);
02185 
02186     cpl_table * slitpos   =
02187         mos_build_slit_location(globaltbl, maskslits, ccd_ysize);
02188 
02189     double * l_ytop = cpl_table_get_data_double(slitpos, "ytop");
02190     int    * l_id   = cpl_table_get_data_int(slitpos, "slit_id");
02191 
02192     int      npairs = cpl_table_get_nrow(slitpos);
02193 
02194     double * ytop   = cpl_table_get_data_double(slits, "ytop");
02195     double * ybot   = cpl_table_get_data_double(slits, "ybottom");
02196     int    * p_id;
02197 
02198     int      k;
02199 
02200     /* Just in case it has been modified */
02201     nslits = cpl_table_get_nrow(slits);
02202 
02203     cpl_table_new_column(slits, "pair_id", CPL_TYPE_INT);
02204     p_id   = cpl_table_get_data_int(slits, "pair_id");
02205 
02206     for (k = 0; k < npairs; k++) {
02207         int h;
02208 
02209         for (h = 0; h < nslits; h++) {
02210 
02211                 if (l_ytop[k] < ytop[h] && l_ytop[k] > ybot[h]) {
02212                     if (h + 1 < nslits) {
02213                 cpl_table_set_int(slits, "pair_id", h, l_id[k]);
02214                 cpl_table_set_int(slits, "pair_id", h + 1, l_id[k]);
02215                     }
02216                 }
02217 
02218 /***
02219         double difftop = ytop[h] - l_ytop[k];
02220 
02221         if (fabs(difftop + OFFSET) < TOLERANCE ||
02222             fabs(difftop - OFFSET) < TOLERANCE    )
02223             cpl_table_set_int(slits, "pair_id", h, l_id[k]);
02224 ***/
02225         }
02226     }
02227 
02228 /* %%% */
02229 
02230     cpl_table_fill_invalid_int(slits, "pair_id", -1);
02231 /*
02232     cpl_table_unselect_all(slits);
02233 
02234     for (k = 0; k < nslits; k++) {
02235         if (p_id[k] == -1) {
02236         cpl_table_select_row(slits, k);
02237         } else {
02238         if (k == 0) {
02239             if (p_id[k] != p_id[k + 1]) {
02240             cpl_table_select_row(slits, k);
02241             }
02242         } else if (k == nslits - 1) {
02243             if (p_id[k] != p_id[k - 1]) {
02244             cpl_table_select_row(slits, k);
02245             }
02246         } else {
02247             if (p_id[k] != p_id[k + 1] && p_id[k] != p_id[k - 1]) {
02248             cpl_table_select_row(slits, k);
02249             }
02250         }
02251         }
02252     }
02253 
02254     cpl_table_erase_selected(slits);
02255 */
02256     cpl_table_delete(slitpos);   slitpos   = NULL;
02257     cpl_table_delete(globaltbl); globaltbl = NULL;
02258 
02259         cpl_table_delete(maskslits); maskslits = NULL;
02260     }
02261 
02262     if (dfs_save_table(frameset, slits, slit_location_tag, NULL,
02263                parlist, recipe, version)) {
02264     fors_pmos_calib_exit(NULL);
02265     }
02266 
02267     cpl_table_delete(slits); slits = NULL;
02268 
02269     cpl_image_delete(spatial); spatial = NULL;
02270 
02271     cpl_free(instrume); instrume = NULL;
02272 
02273     cpl_table_delete(overscans); overscans = NULL;
02274     cpl_image_delete(master_bias); master_bias = NULL;
02275     cpl_image_delete(master_flat); master_flat = NULL;
02276 
02277     cpl_table_delete(wavelengths); wavelengths = NULL;
02278     cpl_vector_delete(lines); lines = NULL;
02279 
02280     if (cpl_error_get_code()) {
02281         cpl_msg_error(cpl_error_get_where(), cpl_error_get_message());
02282         fors_pmos_calib_exit(NULL);
02283     }
02284 
02285     return 0;
02286 }

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