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