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