CR2RE Pipeline Reference Manual 1.7.0
cr2res_cal_wave.c
1/*
2
3 * Copyright (C) 2002,2003 European Southern Observatory
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 2 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02111-1307 USA
18 */
19
20#ifdef HAVE_CONFIG_H
21#include <config.h>
22#endif
23
24/*-----------------------------------------------------------------------------
25 Includes
26 -----------------------------------------------------------------------------*/
27
28#include <locale.h>
29#include <string.h>
30
31#include <cpl.h>
32#include <math.h>
33#include "hdrl.h"
34
35#include "cr2res_utils.h"
36#include "cr2res_calib.h"
37#include "cr2res_pfits.h"
38#include "cr2res_dfs.h"
39#include "cr2res_io.h"
40#include "cr2res_qc.h"
41#include "cr2res_slit_curv.h"
42#include "cr2res_extract.h"
43#include "cr2res_trace.h"
44#include "cr2res_wave.h"
45
46/*-----------------------------------------------------------------------------
47 Define
48 -----------------------------------------------------------------------------*/
49
50#define RECIPE_STRING "cr2res_cal_wave"
51#define CRIRES_PI 3.1415926535897932384626433832795029L
52
53/*-----------------------------------------------------------------------------
54 Plugin registration
55 -----------------------------------------------------------------------------*/
56
57int cpl_plugin_get_info(cpl_pluginlist * list);
58
59/*-----------------------------------------------------------------------------
60 Private function prototypes
61 -----------------------------------------------------------------------------*/
62
63static int cr2res_cal_wave_qc_une_flux(
64 const cpl_table ** extracted,
65 const cpl_propertylist ** ext_plist,
66 double cwlen,
67 cpl_propertylist * plist) ;
68static int cr2res_cal_wave_qc_fpi(
69 const cpl_table ** extracted,
70 const cpl_table ** tw,
71 const cpl_propertylist ** ext_plist,
72 double cwlen,
73 cpl_propertylist * plist) ;
74static int cr2res_cal_wave_qc_tilt(
75 const cpl_table ** tws,
76 double cwlen,
77 cpl_propertylist * plist) ;
78static int cr2res_cal_wave_reduce(
79 const cpl_frameset * rawframes_une,
80 const cpl_frameset * rawframes_fpet,
81 const cpl_frame * detlin_frame,
82 const cpl_frame * master_dark_frame,
83 const cpl_frame * master_flat_frame,
84 const cpl_frame * bpm_frame,
85 const cpl_frame * trace_wave_frame,
86 const cpl_frame * lines_frame,
87 int reduce_det,
88 int reduce_order,
89 int reduce_trace,
90 int slit_degree,
91 int subtract_nolight_rows,
92 int cosmics,
93 cr2res_collapse collapse,
94 int ext_height,
95 int ext_swath_width,
96 int ext_oversample,
97 double ext_smooth_slit,
98 double ext_smooth_spec,
99 cr2res_wavecal_type wavecal_type,
100 int wl_degree,
101 double wl_start,
102 double wl_end,
103 double wl_err,
104 double wl_shift,
105 int log_flag,
106 int fallback_input_wavecal_flag,
107 int keep_higher_degrees_flag,
108 int clean_spectrum,
109 int display,
110 double display_wmin,
111 double display_wmax,
112 cpl_table ** out_trace_wave_une,
113 cpl_table ** lines_diagnostics_une,
114 cpl_table ** out_extracted_une,
115 hdrl_image ** out_wave_map_une,
116 cpl_propertylist ** ext_plist_une,
117 cpl_table ** out_trace_wave_fpet,
118 cpl_table ** lines_diagnostics_fpet,
119 cpl_table ** out_extracted_fpet,
120 hdrl_image ** out_wave_map_fpet,
121 cpl_propertylist ** ext_plist_fpet) ;
122static int cr2res_cal_wave_create(cpl_plugin *);
123static int cr2res_cal_wave_exec(cpl_plugin *);
124static int cr2res_cal_wave_destroy(cpl_plugin *);
125static int cr2res_cal_wave(cpl_frameset *, const cpl_parameterlist *);
126
127/*-----------------------------------------------------------------------------
128 Static variables
129 -----------------------------------------------------------------------------*/
130
131static char cr2res_cal_wave_description[] = "\
132Spectrum Extraction and Wavelength Calibration \n\
133 This recipe performs the extraction of the various orders along \n\
134 the provided traces, and the wavelength calibration of these \n\
135 extracted spectra. \n\
136 It can support different methods (--wl_method parameter): \n\
137 XCORR: Cross Correlation with a emission lines catalog (default) \n\
138 LINE1D: Line identification and fitting for each 1D spectra \n\
139 LINE2D: Line identification and fitting for all 1D spectra at once \n\
140 \n\
141 Inputs \n\
142 raw.fits " CR2RES_WAVE_UNE_RAW " [1 to n] \n\
143 or " CR2RES_WAVE_FPET_RAW " [0 to n] \n\
144 trace.fits " CR2RES_CAL_FLAT_TW_PROCATG " [1] \n\
145 or " CR2RES_CAL_FLAT_TW_MERGED_PROCATG " \n\
146 or " CR2RES_UTIL_TRACE_TW_PROCATG " \n\
147 or " CR2RES_UTIL_WAVE_TW_PROCATG " \n\
148 or " CR2RES_CAL_WAVE_TW_PROCATG " \n\
149 or " CR2RES_UTIL_SLIT_CURV_TW_PROCATG " \n\
150 detlin.fits " CR2RES_CAL_DETLIN_COEFFS_PROCATG " [0 to 1] \n\
151 bpm.fits " CR2RES_CAL_DARK_BPM_PROCATG " [0 to 1] \n\
152 or " CR2RES_CAL_FLAT_BPM_PROCATG " \n\
153 or " CR2RES_CAL_DETLIN_BPM_PROCATG " \n\
154 or " CR2RES_UTIL_BPM_MERGE_PROCATG " \n\
155 or " CR2RES_UTIL_BPM_SPLIT_PROCATG " \n\
156 master_dark.fits " CR2RES_CAL_DARK_MASTER_PROCATG " [0 to 1] \n\
157 master_flat.fits " CR2RES_CAL_FLAT_MASTER_PROCATG " [0 to 1] \n\
158 lines.fits " CR2RES_EMISSION_LINES_PROCATG " [0 to 1] \n\
159 \n\
160 Outputs \n\
161 cr2res_cal_wave_tw.fits " CR2RES_CAL_WAVE_TW_PROCATG" \n\
162 cr2res_cal_wave_wave_map.fits " CR2RES_CAL_WAVE_MAP_PROCATG" \n\
163 cr2res_cal_wave_lines_diagnostics.fits "
164 CR2RES_CAL_WAVE_LINES_DIAGNOSTICS_PROCATG"\n\
165 cr2res_cal_wave_extracted.fits " CR2RES_CAL_WAVE_EXTRACT_1D_PROCATG"\n\
166 \n\
167 Algorithm \n\
168 loop on detectors d: \n\
169 Call cr2res_cal_wave_reduce() \n\
170 -> out_trace_wave[une|fpet](d) \n\
171 -> lines_diagnostics[une|fpet](d) \n\
172 -> out_extracted[une|fpet](d) \n\
173 -> out_wave_map[une|fpet](d) \n\
174 Save out_trace_wave[une|fpet] \n\
175 Save lines_diagnostics[une|fpet] \n\
176 Save out_extracted[une|fpet] \n\
177 Save out_wave_map[une|fpet] \n\
178 \n\
179 cr2res_cal_wave_reduce(): \n\
180 If FPET RAW frame is present, compute the slit_curvature from it, \n\
181 use it in the following, and store it in the out_trace_wave_fpet. \n\
182 \n\
183 Successively for UNE and FPET RAW frames: \n\
184 Load the raw image list \n\
185 Apply the calibrations to the image list \n\
186 Collapse the image list \n\
187 Extract along the traces from the collapsed image \n\
188 Compute the Wavelength with cr2res_wave_apply() \n\
189 -> out_trace_wave \n\
190 -> lines_diagnostics \n\
191 -> out_extracted \n\
192 Compute the Wavelength map \n\
193 -> out_wave_map \n\
194 \n\
195 cr2res_wave_apply() \n\
196 loop on the traces t: \n\
197 Get the spectrum \n\
198 Get the Initial guess \n\
199 Switch on the required method: \n\
200 CR2RES_LINE2D: cr2res_wave_2d() \n\
201 CR2RES_LINE1D: cr2res_wave_line_fitting() \n\
202 CR2RES_ETALON: cr2res_wave_etalon() \n\
203 CR2RES_XCORR: cr2res_wave_xcorr() \n\
204 \n\
205 Library Functions used \n\
206 cr2res_io_find_TRACE_WAVE() \n\
207 cr2res_io_find_BPM() \n\
208 cr2res_io_read_dits() \n\
209 cr2res_io_load_image_list_from_set() \n\
210 cr2res_calib_imagelist() \n\
211 cr2res_io_load_TRACE_WAVE() \n\
212 cr2res_extract_traces() \n\
213 cr2res_wave_apply() \n\
214 cr2res_extract_EXTRACT1D_get_spectrum() \n\
215 cr2res_wave_estimate_compute() \n\
216 cr2res_wave_2d() \n\
217 cr2res_wave_line_fitting() \n\
218 cr2res_wave_etalon() \n\
219 cr2res_wave_xcorr() \n\
220 cr2res_wave_gen_wave_map() \n\
221 cr2res_io_save_TRACE_WAVE() \n\
222 cr2res_io_save_WAVE_MAP() \n\
223 cr2res_io_save_LINES_DIAGNOSTICS() \n\
224 cr2res_io_save_EXTRACT_1D() \n\
225" ;
226
227/*-----------------------------------------------------------------------------
228 Function code
229 -----------------------------------------------------------------------------*/
230
231/*----------------------------------------------------------------------------*/
241/*----------------------------------------------------------------------------*/
242int cpl_plugin_get_info(cpl_pluginlist * list)
243{
244 cpl_recipe * recipe = cpl_calloc(1, sizeof *recipe );
245 cpl_plugin * plugin = &recipe->interface;
246
247 if (cpl_plugin_init(plugin,
248 CPL_PLUGIN_API,
249 CR2RES_BINARY_VERSION,
250 CPL_PLUGIN_TYPE_RECIPE,
251 RECIPE_STRING,
252 "Wavelength Calibration",
253 cr2res_cal_wave_description,
254 CR2RES_PIPELINE_AUTHORS,
255 PACKAGE_BUGREPORT,
256 cr2res_get_license(),
257 cr2res_cal_wave_create,
258 cr2res_cal_wave_exec,
259 cr2res_cal_wave_destroy)) {
260 cpl_msg_error(cpl_func, "Plugin initialization failed");
261 (void)cpl_error_set_where(cpl_func);
262 return 1;
263 }
264
265 if (cpl_pluginlist_append(list, plugin)) {
266 cpl_msg_error(cpl_func, "Error adding plugin to list");
267 (void)cpl_error_set_where(cpl_func);
268 return 1;
269 }
270
271 return 0;
272}
273
274/*----------------------------------------------------------------------------*/
282/*----------------------------------------------------------------------------*/
283static int cr2res_cal_wave_create(cpl_plugin * plugin)
284{
285 cpl_recipe * recipe;
286 cpl_parameter * p;
287
288 /* Check that the plugin is part of a valid recipe */
289 if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE)
290 recipe = (cpl_recipe *)plugin;
291 else
292 return -1;
293
294 /* Create the parameters list in the cpl_recipe object */
295 recipe->parameters = cpl_parameterlist_new();
296
297 /* Fill the parameters list */
298 p = cpl_parameter_new_value("cr2res.cr2res_cal_wave.detector",
299 CPL_TYPE_INT, "Only reduce the specified detector",
300 "cr2res.cr2res_cal_wave", 0);
301 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "detector");
302 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
303 cpl_parameterlist_append(recipe->parameters, p);
304
305 p = cpl_parameter_new_value("cr2res.cr2res_cal_wave.order",
306 CPL_TYPE_INT, "Only reduce the specified order",
307 "cr2res.cr2res_cal_wave", -1);
308 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "order");
309 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
310 cpl_parameterlist_append(recipe->parameters, p);
311
312 p = cpl_parameter_new_value("cr2res.cr2res_cal_wave.trace_nb",
313 CPL_TYPE_INT, "Only reduce the specified trace number",
314 "cr2res.cr2res_cal_wave", -1);
315 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "trace_nb");
316 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
317 cpl_parameterlist_append(recipe->parameters, p);
318
319 p = cpl_parameter_new_value("cr2res.cr2res_cal_wave.slit_degree",
320 CPL_TYPE_INT, "Slit fitting Polynomial degree (1 or 2)",
321 "cr2res.cr2res_cal_wave", 2);
322 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "slit_degree");
323 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
324 cpl_parameterlist_append(recipe->parameters, p);
325
326 p = cpl_parameter_new_value("cr2res.cr2res_cal_wave.subtract_nolight_rows",
327 CPL_TYPE_BOOL, "Subtract the no-light rows.",
328 "cr2res.cr2res_cal_wave", FALSE);
329 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "subtract_nolight_rows");
330 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
331 cpl_parameterlist_append(recipe->parameters, p);
332
333 p = cpl_parameter_new_value("cr2res.cr2res_cal_wave.cosmics",
334 CPL_TYPE_BOOL, "Find and mark cosmic rays hits as bad",
335 "cr2res.cr2res_cal_wave", FALSE);
336 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "cosmics");
337 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
338 cpl_parameterlist_append(recipe->parameters, p);
339
340 p = cpl_parameter_new_value("cr2res.cr2res_cal_wave.collapse_method",
341 CPL_TYPE_STRING, "Collapse the input images (MEAN or MEDIAN)",
342 "cr2res.cr2res_cal_wave", "MEDIAN");
343 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "collapse_method");
344 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
345 cpl_parameterlist_append(recipe->parameters, p);
346
347 p = cpl_parameter_new_value("cr2res.cr2res_cal_wave.extract_oversample",
348 CPL_TYPE_INT, "factor by which to oversample the extraction",
349 "cr2res.cr2res_cal_wave", 7);
350 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "extract_oversample");
351 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
352 cpl_parameterlist_append(recipe->parameters, p);
353
354 p = cpl_parameter_new_value("cr2res.cr2res_cal_wave.extract_swath_width",
355 CPL_TYPE_INT, "The swath width", "cr2res.cr2res_cal_wave", 2048);
356 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "extract_swath_width");
357 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
358 cpl_parameterlist_append(recipe->parameters, p);
359
360 p = cpl_parameter_new_value("cr2res.cr2res_cal_wave.extract_height",
361 CPL_TYPE_INT, "Extraction height",
362 "cr2res.cr2res_cal_wave", 160);
363 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "extract_height");
364 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
365 cpl_parameterlist_append(recipe->parameters, p);
366
367 p = cpl_parameter_new_value("cr2res.cr2res_cal_wave.extract_smooth_slit",
368 CPL_TYPE_DOUBLE, "Smoothing along the slit",
369 "cr2res.cr2res_cal_wave", 15.0);
370 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "extract_smooth_slit");
371 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
372 cpl_parameterlist_append(recipe->parameters, p);
373
374 p = cpl_parameter_new_value("cr2res.cr2res_cal_wave.extract_smooth_spec",
375 CPL_TYPE_DOUBLE, "Smoothing along the spectrum",
376 "cr2res.cr2res_cal_wave", 2.0E-7);
377 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "extract_smooth_spec");
378 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
379 cpl_parameterlist_append(recipe->parameters, p);
380
381 p = cpl_parameter_new_value("cr2res.cr2res_cal_wave.wl_method",
382 CPL_TYPE_STRING,
383 "Wavelength Method (XCORR / LINE1D / LINE2D)",
384 "cr2res.cr2res_cal_wave", "XCORR");
385 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "wl_method");
386 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
387 cpl_parameterlist_append(recipe->parameters, p);
388
389 p = cpl_parameter_new_value("cr2res.cr2res_cal_wave.wl_shift",
390 CPL_TYPE_DOUBLE, "Wavelength shift (nm) to apply to the guess",
391 "cr2res.cr2res_cal_wave", 0.0);
392 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "wl_shift");
393 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
394 cpl_parameterlist_append(recipe->parameters, p);
395
396 p = cpl_parameter_new_value("cr2res.cr2res_cal_wave.wl_est",
397 CPL_TYPE_STRING, "Estimated wavelength start and end",
398 "cr2res.cr2res_cal_wave", "-1.0, -1.0");
399 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "wl_est");
400 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
401 cpl_parameterlist_append(recipe->parameters, p);
402
403 p = cpl_parameter_new_value("cr2res.cr2res_cal_wave.wl_err",
404 CPL_TYPE_DOUBLE, "Estimated wavelength error (in nm)",
405 "cr2res.cr2res_cal_wave", 0.04);
406 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "wl_err");
407 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
408 cpl_parameterlist_append(recipe->parameters, p);
409
410 p = cpl_parameter_new_value("cr2res.cr2res_cal_wave.wl_degree",
411 CPL_TYPE_INT, "Wavelength Polynomial degree",
412 "cr2res.cr2res_cal_wave", 0);
413 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "wl_degree");
414 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
415 cpl_parameterlist_append(recipe->parameters, p);
416
417 p = cpl_parameter_new_value("cr2res.cr2res_cal_wave.log",
418 CPL_TYPE_BOOL, "Flag for taking the Log() value of the lines",
419 "cr2res.cr2res_cal_wave", FALSE);
420 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "log");
421 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
422 cpl_parameterlist_append(recipe->parameters, p);
423
424 p = cpl_parameter_new_value("cr2res.cr2res_cal_wave.fallback_input_wavecal",
425 CPL_TYPE_BOOL, "Flag for using the input WL when no computation",
426 "cr2res.cr2res_cal_wave", FALSE);
427 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "fallback");
428 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
429 cpl_parameterlist_append(recipe->parameters, p);
430
431 p = cpl_parameter_new_value("cr2res.cr2res_cal_wave.keep_higher_degrees",
432 CPL_TYPE_BOOL,
433 "Flag for re-using higher degrees of first guess in Cross-Corr.",
434 "cr2res.cr2res_cal_wave", TRUE);
435 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "keep");
436 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
437 cpl_parameterlist_append(recipe->parameters, p);
438
439 p = cpl_parameter_new_value("cr2res.cr2res_cal_wave.save_intermediate",
440 CPL_TYPE_BOOL,
441 "Flag to save UNE results (if UNE and FPET are used as inputs)",
442 "cr2res.cr2res_cal_wave", FALSE);
443 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "save");
444 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
445 cpl_parameterlist_append(recipe->parameters, p);
446
447 p = cpl_parameter_new_value("cr2res.cr2res_cal_wave.clean_spectrum",
448 CPL_TYPE_BOOL, "Flag to automatically clean the missing lines",
449 "cr2res.cr2res_cal_wave", TRUE);
450 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "clean_spectrum");
451 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
452 cpl_parameterlist_append(recipe->parameters, p);
453
454 p = cpl_parameter_new_value("cr2res.cr2res_cal_wave.display",
455 CPL_TYPE_BOOL, "Flag for display",
456 "cr2res.cr2res_cal_wave", FALSE);
457 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "display");
458 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
459 cpl_parameterlist_append(recipe->parameters, p);
460
461 p = cpl_parameter_new_value("cr2res.cr2res_cal_wave.display_range",
462 CPL_TYPE_STRING,
463 "Wavelength range to display [start, end] (in nm)",
464 "cr2res.cr2res_cal_wave", "-1.0, -1.0");
465 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "display_range");
466 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
467 cpl_parameterlist_append(recipe->parameters, p);
468
469 return 0;
470}
471
472/*----------------------------------------------------------------------------*/
478/*----------------------------------------------------------------------------*/
479static int cr2res_cal_wave_exec(cpl_plugin * plugin)
480{
481 cpl_recipe *recipe;
482
483 /* Get the recipe out of the plugin */
484 if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE)
485 recipe = (cpl_recipe *)plugin;
486 else return -1;
487
488 return cr2res_cal_wave(recipe->frames, recipe->parameters);
489}
490
491/*----------------------------------------------------------------------------*/
497/*----------------------------------------------------------------------------*/
498static int cr2res_cal_wave_destroy(cpl_plugin * plugin)
499{
500 cpl_recipe *recipe;
501
502 /* Get the recipe out of the plugin */
503 if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE)
504 recipe = (cpl_recipe *)plugin;
505 else return -1 ;
506
507 cpl_parameterlist_delete(recipe->parameters);
508 return 0 ;
509}
510
511/*----------------------------------------------------------------------------*/
518/*----------------------------------------------------------------------------*/
519static int cr2res_cal_wave(
520 cpl_frameset * frameset,
521 const cpl_parameterlist * parlist)
522{
523 const cpl_parameter * param;
524 int reduce_det, reduce_order, reduce_trace,
525 slit_degree, ext_oversample, ext_swath_width,
526 ext_height, wl_degree, display, log_flag,
527 fallback_input_wavecal_flag,
528 save_intermediate_flag, keep_higher_degrees_flag,
529 clean_spectrum, subtract_nolight_rows,
530 cosmics;
531 double ext_smooth_slit, ext_smooth_spec, wl_start, wl_end,
532 wl_err, wl_shift, display_wmin,
533 display_wmax, central_wlen;
534 cr2res_collapse collapse ;
535 cr2res_wavecal_type wavecal_type ;
536 const char * sval ;
537 cpl_frameset * rawframes_une ;
538 cpl_frameset * rawframes_fpet ;
539 const cpl_frame * detlin_frame ;
540 const cpl_frame * master_dark_frame ;
541 const cpl_frame * master_flat_frame ;
542 const cpl_frame * bpm_frame ;
543 const cpl_frame * trace_wave_frame ;
544 const cpl_frame * lines_frame ;
545 cpl_propertylist * plist ;
546 char * out_file;
547 cpl_table * out_trace_wave_une[CR2RES_NB_DETECTORS] ;
548 cpl_table * lines_diagnostics_une[CR2RES_NB_DETECTORS] ;
549 cpl_table * out_extracted_une[CR2RES_NB_DETECTORS] ;
550 hdrl_image * out_wave_map_une[CR2RES_NB_DETECTORS] ;
551 cpl_propertylist * ext_plist_une[CR2RES_NB_DETECTORS] ;
552 cpl_table * out_trace_wave_fpet[CR2RES_NB_DETECTORS] ;
553 cpl_table * lines_diagnostics_fpet[CR2RES_NB_DETECTORS] ;
554 cpl_table * out_extracted_fpet[CR2RES_NB_DETECTORS] ;
555 hdrl_image * out_wave_map_fpet[CR2RES_NB_DETECTORS] ;
556 cpl_propertylist * ext_plist_fpet[CR2RES_NB_DETECTORS] ;
557 cpl_propertylist * qc_main ;
558 int det_nr, i;
559
560 /* Needed for sscanf() */
561 setlocale(LC_NUMERIC, "C");
562
563 /* Initialise */
564 wl_start = wl_end = -1.0 ;
565
566 collapse = CR2RES_COLLAPSE_UNSPECIFIED ;
567 display_wmin = display_wmax = -1.0 ;
568
569 int qc_sizes[] = {CR2RES_NB_DETECTORS,CR2RES_NB_DETECTORS};
570 cpl_propertylist** qc_plists[] = {ext_plist_une, ext_plist_fpet} ;
571 cpl_propertylist** qc_plists_une[] = {ext_plist_une} ;
572
573 /* RETRIEVE INPUT PARAMETERS */
574 param = cpl_parameterlist_find_const(parlist,
575 "cr2res.cr2res_cal_wave.detector");
576 reduce_det = cpl_parameter_get_int(param);
577 param = cpl_parameterlist_find_const(parlist,
578 "cr2res.cr2res_cal_wave.order");
579 reduce_order = cpl_parameter_get_int(param);
580 param = cpl_parameterlist_find_const(parlist,
581 "cr2res.cr2res_cal_wave.trace_nb");
582 reduce_trace = cpl_parameter_get_int(param);
583 param = cpl_parameterlist_find_const(parlist,
584 "cr2res.cr2res_cal_wave.slit_degree");
585 slit_degree = cpl_parameter_get_int(param);
586 param = cpl_parameterlist_find_const(parlist,
587 "cr2res.cr2res_cal_wave.subtract_nolight_rows");
588 subtract_nolight_rows = cpl_parameter_get_bool(param) ;
589 param = cpl_parameterlist_find_const(parlist,
590 "cr2res.cr2res_cal_wave.cosmics");
591 cosmics = cpl_parameter_get_bool(param);
592 param = cpl_parameterlist_find_const(parlist,
593 "cr2res.cr2res_cal_wave.collapse_method");
594 sval = cpl_parameter_get_string(param);
595 if (!strcmp(sval, "MEAN")) collapse = CR2RES_COLLAPSE_MEAN ;
596 else if (!strcmp(sval, "MEDIAN")) collapse = CR2RES_COLLAPSE_MEDIAN ;
597 if (collapse!=CR2RES_COLLAPSE_MEAN && collapse!=CR2RES_COLLAPSE_MEDIAN) {
598 cpl_msg_error(__func__, "Unsupported collapse method") ;
599 cpl_error_set(__func__, CPL_ERROR_ILLEGAL_INPUT) ;
600 return -1 ;
601 }
602 param = cpl_parameterlist_find_const(parlist,
603 "cr2res.cr2res_cal_wave.extract_oversample");
604 ext_oversample = cpl_parameter_get_int(param);
605 param = cpl_parameterlist_find_const(parlist,
606 "cr2res.cr2res_cal_wave.extract_swath_width");
607 ext_swath_width = cpl_parameter_get_int(param);
608 param = cpl_parameterlist_find_const(parlist,
609 "cr2res.cr2res_cal_wave.extract_height");
610 ext_height = cpl_parameter_get_int(param);
611 param = cpl_parameterlist_find_const(parlist,
612 "cr2res.cr2res_cal_wave.extract_smooth_slit");
613 ext_smooth_slit = cpl_parameter_get_double(param);
614 param = cpl_parameterlist_find_const(parlist,
615 "cr2res.cr2res_cal_wave.extract_smooth_spec");
616 ext_smooth_spec = cpl_parameter_get_double(param);
617 param = cpl_parameterlist_find_const(parlist,
618 "cr2res.cr2res_cal_wave.wl_method");
619 sval = cpl_parameter_get_string(param) ;
620 if (!strcmp(sval, "XCORR")) wavecal_type = CR2RES_XCORR ;
621 else if (!strcmp(sval, "LINE1D")) wavecal_type = CR2RES_LINE1D ;
622 else if (!strcmp(sval, "LINE2D")) wavecal_type = CR2RES_LINE2D ;
623 else {
624 cpl_msg_error(__func__, "Invalid Method specified");
625 cpl_error_set(__func__, CPL_ERROR_ILLEGAL_INPUT) ;
626 return -1;
627 }
628 param = cpl_parameterlist_find_const(parlist,
629 "cr2res.cr2res_cal_wave.wl_shift");
630 wl_shift = cpl_parameter_get_double(param) ;
631 param = cpl_parameterlist_find_const(parlist,
632 "cr2res.cr2res_cal_wave.wl_est");
633 sval = cpl_parameter_get_string(param) ;
634 if (sscanf(sval, "%lg,%lg", &wl_start, &wl_end) != 2) {
635 return -1 ;
636 }
637 param = cpl_parameterlist_find_const(parlist,
638 "cr2res.cr2res_cal_wave.wl_err");
639 wl_err = cpl_parameter_get_double(param) ;
640 param = cpl_parameterlist_find_const(parlist,
641 "cr2res.cr2res_cal_wave.wl_degree");
642 wl_degree = cpl_parameter_get_int(param);
643 param = cpl_parameterlist_find_const(parlist,
644 "cr2res.cr2res_cal_wave.log");
645 log_flag = cpl_parameter_get_bool(param) ;
646 param = cpl_parameterlist_find_const(parlist,
647 "cr2res.cr2res_cal_wave.fallback_input_wavecal");
648 fallback_input_wavecal_flag = cpl_parameter_get_bool(param) ;
649 param = cpl_parameterlist_find_const(parlist,
650 "cr2res.cr2res_cal_wave.keep_higher_degrees");
651 keep_higher_degrees_flag = cpl_parameter_get_bool(param) ;
652 param = cpl_parameterlist_find_const(parlist,
653 "cr2res.cr2res_cal_wave.save_intermediate");
654 save_intermediate_flag = cpl_parameter_get_bool(param) ;
655 param = cpl_parameterlist_find_const(parlist,
656 "cr2res.cr2res_cal_wave.clean_spectrum");
657 clean_spectrum = cpl_parameter_get_bool(param) ;
658 param = cpl_parameterlist_find_const(parlist,
659 "cr2res.cr2res_cal_wave.display");
660 display = cpl_parameter_get_bool(param) ;
661 param = cpl_parameterlist_find_const(parlist,
662 "cr2res.cr2res_cal_wave.display_range");
663 sval = cpl_parameter_get_string(param) ;
664 if (sscanf(sval, "%lg,%lg", &display_wmin, &display_wmax) != 2) {
665 return -1 ;
666 }
667
668 /* Check Parameters */
669 if (display==TRUE && (reduce_order==-1 || reduce_det==0)){
670 cpl_msg_warning(__func__,
671 "Option --display can only be used with --order"
672 " and --detector");
673 display=FALSE;
674 }
675 if (wl_degree < 0) {
676 cpl_msg_error(__func__, "The degree needs to be >= 0");
677 cpl_error_set(__func__, CPL_ERROR_ILLEGAL_INPUT) ;
678 return -1 ;
679 }
680 if (wl_degree == 0 && !keep_higher_degrees_flag) {
681 cpl_msg_error(__func__, "The degree 0 can only be used with --keep");
682 cpl_error_set(__func__, CPL_ERROR_ILLEGAL_INPUT) ;
683 return -1 ;
684 }
685 if (slit_degree != 1 && slit_degree != 2) {
686 cpl_msg_error(__func__, "The slit fit degree must be 1 or 2");
687 cpl_error_set(__func__, CPL_ERROR_ILLEGAL_INPUT) ;
688 return -1 ;
689 }
690
691 /* Identify the RAW and CALIB frames in the input frameset */
692 if (cr2res_dfs_set_groups(frameset) != CPL_ERROR_NONE) {
693 cpl_msg_error(__func__, "Cannot identify RAW and CALIB frames") ;
694 cpl_error_set(__func__, CPL_ERROR_ILLEGAL_INPUT) ;
695 return -1 ;
696 }
697
698 /* Retrieve calibration data */
699 trace_wave_frame = cr2res_io_find_TRACE_WAVE(frameset) ;
700 if (trace_wave_frame == NULL) {
701 cpl_msg_error(__func__, "Could not find TRACE_WAVE frame") ;
702 return -1 ;
703 }
704 detlin_frame = cpl_frameset_find_const(frameset,
705 CR2RES_CAL_DETLIN_COEFFS_PROCATG);
706 master_dark_frame = cpl_frameset_find_const(frameset,
707 CR2RES_CAL_DARK_MASTER_PROCATG) ;
708 master_flat_frame = cpl_frameset_find_const(frameset,
709 CR2RES_CAL_FLAT_MASTER_PROCATG) ;
710 bpm_frame = cr2res_io_find_BPM(frameset) ;
711 lines_frame = cpl_frameset_find_const(frameset,
712 CR2RES_EMISSION_LINES_PROCATG) ;
713
714 /* Get the RAW Frames */
715 rawframes_une = cr2res_extract_frameset(frameset, CR2RES_WAVE_UNE_RAW) ;
716 rawframes_fpet = cr2res_extract_frameset(frameset, CR2RES_WAVE_FPET_RAW) ;
717 if (rawframes_une == NULL && rawframes_fpet == NULL) {
718 cpl_msg_error(__func__, "Could not find RAW frames") ;
719 return -1 ;
720 }
721 if (lines_frame == NULL && rawframes_une != NULL) {
722 cpl_frameset_delete(rawframes_une) ;
723 if (rawframes_fpet!=NULL) cpl_frameset_delete(rawframes_fpet) ;
724 cpl_msg_error(__func__, "The emission lines catalog is needed");
725 cpl_error_set(__func__, CPL_ERROR_ILLEGAL_INPUT) ;
726 return -1 ;
727 }
728 if (reduce_order > -1 && wavecal_type == CR2RES_LINE2D) {
729 if (rawframes_une !=NULL) cpl_frameset_delete(rawframes_une) ;
730 if (rawframes_fpet!=NULL) cpl_frameset_delete(rawframes_fpet) ;
731 cpl_msg_error(__func__, "Limiting to one order with LINE2D impossible");
732 cpl_error_set(__func__, CPL_ERROR_ILLEGAL_INPUT) ;
733 return -1 ;
734 }
735 if (reduce_order > -1 && wavecal_type == CR2RES_ETALON) {
736 if (rawframes_une !=NULL) cpl_frameset_delete(rawframes_une) ;
737 if (rawframes_fpet!=NULL) cpl_frameset_delete(rawframes_fpet) ;
738 cpl_msg_error(__func__, "Limiting to one order with ETALON impossible");
739 cpl_error_set(__func__, CPL_ERROR_ILLEGAL_INPUT) ;
740 return -1 ;
741 }
742
743 /* Loop over the detectors */
744 for (det_nr=1 ; det_nr<=CR2RES_NB_DETECTORS ; det_nr++) {
745
746 /* Initialise */
747 out_trace_wave_une[det_nr-1] = NULL ;
748 lines_diagnostics_une[det_nr-1] = NULL ;
749 out_extracted_une[det_nr-1] = NULL ;
750 out_wave_map_une[det_nr-1] = NULL ;
751 ext_plist_une[det_nr-1] = NULL ;
752
753 out_trace_wave_fpet[det_nr-1] = NULL ;
754 lines_diagnostics_fpet[det_nr-1] = NULL ;
755 out_extracted_fpet[det_nr-1] = NULL ;
756 out_wave_map_fpet[det_nr-1] = NULL ;
757 ext_plist_fpet[det_nr-1] = NULL ;
758
759 /* Compute only one detector but not this one */
760 if (reduce_det != 0 && det_nr != reduce_det) {
761 /* This Detector will not be processed here */
762 /* The output trace wave contains the input one ... */
763 if (rawframes_une != NULL)
764 out_trace_wave_une[det_nr-1] = cr2res_io_load_TRACE_WAVE(
765 cpl_frame_get_filename(trace_wave_frame), det_nr) ;
766 if (rawframes_fpet != NULL)
767 out_trace_wave_fpet[det_nr-1] = cr2res_io_load_TRACE_WAVE(
768 cpl_frame_get_filename(trace_wave_frame), det_nr) ;
769 /* ... without the WL / WL_ERR */
770 /* ... unless fallback_input_wavecal_flag is set */
771
772 if (!fallback_input_wavecal_flag) {
773 /* Reset WL / WL_ERR */
774 if (rawframes_une != NULL) {
775 cpl_table_erase_column(out_trace_wave_une[det_nr-1],
776 CR2RES_COL_WAVELENGTH) ;
777 cpl_table_new_column_array(out_trace_wave_une[det_nr-1],
778 CR2RES_COL_WAVELENGTH, CPL_TYPE_DOUBLE, 2) ;
779 cpl_table_erase_column(out_trace_wave_une[det_nr-1],
780 CR2RES_COL_WAVELENGTH_ERROR) ;
781 cpl_table_new_column_array(out_trace_wave_une[det_nr-1],
782 CR2RES_COL_WAVELENGTH_ERROR, CPL_TYPE_DOUBLE, 2) ;
783 }
784 if (rawframes_fpet != NULL) {
785 cpl_table_erase_column(out_trace_wave_fpet[det_nr-1],
786 CR2RES_COL_WAVELENGTH) ;
787 cpl_table_new_column_array(out_trace_wave_fpet[det_nr-1],
788 CR2RES_COL_WAVELENGTH, CPL_TYPE_DOUBLE, 2) ;
789 cpl_table_erase_column(out_trace_wave_fpet[det_nr-1],
790 CR2RES_COL_WAVELENGTH_ERROR) ;
791 cpl_table_new_column_array(out_trace_wave_fpet[det_nr-1],
792 CR2RES_COL_WAVELENGTH_ERROR, CPL_TYPE_DOUBLE, 2) ;
793 }
794 }
795 continue ;
796 }
797
798 cpl_msg_info(__func__, "Process detector number %d", det_nr) ;
799 cpl_msg_indent_more() ;
800
801 /* For Testing the QCs without re-running all the time */
802 if (0) {
803 out_trace_wave_une[det_nr-1] =
804 cr2res_io_load_TRACE_WAVE(
805"/home/yjung/P_cr2res/test_cr2res_cal_wave/res/cr2res_cal_wave_tw_une.fits",
806 det_nr) ;
807 out_extracted_une[det_nr-1] = cr2res_io_load_EXTRACT_1D(
808"/home/yjung/P_cr2res/test_cr2res_cal_wave/res/cr2res_cal_wave_extracted_une.fits",
809 det_nr) ;
810 ext_plist_une[det_nr-1] = cpl_propertylist_load(
811"/home/yjung/P_cr2res/test_cr2res_cal_wave/res/cr2res_cal_wave_extracted_une.fits",
812 det_nr) ;
813 out_trace_wave_fpet[det_nr-1] =
814 cr2res_io_load_TRACE_WAVE(
815"/home/yjung/P_cr2res/test_cr2res_cal_wave/res/cr2res_cal_wave_tw_fpet.fits",
816 det_nr) ;
817 out_extracted_fpet[det_nr-1] =
818 cr2res_io_load_EXTRACT_1D(
819"/home/yjung/P_cr2res/test_cr2res_cal_wave/res/cr2res_cal_wave_extracted_fpet.fits",
820 det_nr) ;
821 ext_plist_fpet[det_nr-1] =
822 cpl_propertylist_load(
823"/home/yjung/P_cr2res/test_cr2res_cal_wave/res/cr2res_cal_wave_extracted_fpet.fits",
824 det_nr) ;
825 lines_diagnostics_une[det_nr-1] = NULL ;
826 out_wave_map_une[det_nr-1] = NULL ;
827 out_wave_map_fpet[det_nr-1] = NULL ;
828 lines_diagnostics_fpet[det_nr-1] = NULL ;
829 } else {
830
831 /* Call the reduction function */
832 if (cr2res_cal_wave_reduce(rawframes_une, rawframes_fpet, detlin_frame,
833 master_dark_frame, master_flat_frame, bpm_frame,
834 trace_wave_frame, lines_frame, det_nr, reduce_order,
835 reduce_trace, slit_degree, subtract_nolight_rows, cosmics,
836 collapse, ext_height, ext_swath_width, ext_oversample,
837 ext_smooth_slit, ext_smooth_spec, wavecal_type, wl_degree,
838 wl_start, wl_end, wl_err, wl_shift, log_flag,
839 fallback_input_wavecal_flag, keep_higher_degrees_flag,
840 clean_spectrum, display, display_wmin,
841 display_wmax,
842 &(out_trace_wave_une[det_nr-1]),
843 &(lines_diagnostics_une[det_nr-1]),
844 &(out_extracted_une[det_nr-1]),
845 &(out_wave_map_une[det_nr-1]),
846 &(ext_plist_une[det_nr-1]),
847 &(out_trace_wave_fpet[det_nr-1]),
848 &(lines_diagnostics_fpet[det_nr-1]),
849 &(out_extracted_fpet[det_nr-1]),
850 &(out_wave_map_fpet[det_nr-1]),
851 &(ext_plist_fpet[det_nr-1])) == -1) {
852 cpl_msg_warning(__func__, "Failed to reduce detector %d", det_nr);
853 cpl_error_reset() ;
854 }
855 }
856
857 cpl_msg_indent_less() ;
858 }
859
860 /* QC Computation */
861 qc_main = cpl_propertylist_new() ;
862
863 cr2res_qc_calculate_mean_and_rmsd(qc_plists, 2, qc_sizes,
864 CR2RES_HEADER_QC_WAVE_CENTWL, qc_main,
865 CR2RES_HEADER_QC_WAVE_CENTWL_AVG, CR2RES_HEADER_QC_WAVE_CENTWL_RMS);
866
867 cr2res_qc_calculate_mean_and_rmsd(qc_plists, 2, qc_sizes,
868 CR2RES_HEADER_QC_WAVE_DISPWL, qc_main,
869 CR2RES_HEADER_QC_WAVE_DISPWL_AVG, CR2RES_HEADER_QC_WAVE_DISPWL_RMS);
870
871 char* qc_une_keywords[] = {
872 CR2RES_HEADER_QC_WAVE_LAMP_EFFIC,
873 CR2RES_HEADER_QC_WAVE_RESOL,
874 CR2RES_HEADER_QC_WAVE_RESOL_FWHM,
875 CR2RES_HEADER_QC_WAVE_POS,
876 CR2RES_HEADER_QC_WAVE_BESTXCORR,
877 CR2RES_HEADER_QC_OVEREXPOSED
878};
879
880
881 int qc_une_keywords_size = sizeof(qc_une_keywords) / sizeof(qc_une_keywords[0]);
882
883 if (rawframes_une != NULL) {
884 plist = cpl_propertylist_load(cpl_frame_get_filename(
885 cpl_frameset_get_position_const(rawframes_une, 0)), 0) ;
886 central_wlen = cr2res_pfits_get_cwlen(plist) ;
887 cpl_propertylist_delete(plist);
888
889 /* QC.UNE.FLUX Computation */
890 if (cr2res_cal_wave_qc_une_flux(
891 (const cpl_table **)out_extracted_une,
892 (const cpl_propertylist **)ext_plist_une,
893 central_wlen, qc_main)) {
894 cpl_msg_warning(__func__, "QC.UNE.FLUX Computation failed") ;
895 }
896
897 /* Should get min and max order_id's from column of all detectors here*/
898 /*int nb_traces = cpl_table_get_nrow(out_trace_wave_une[0]);*/
899
900 int min_order = INT_MAX;
901 int max_order = INT_MIN;
902 int current_min, current_max;
903
904 for (int qc_i = 0; qc_i < CR2RES_NB_DETECTORS; qc_i++) {
905 if (out_trace_wave_une[qc_i] == NULL) {
906 continue;
907 }
908 current_min = cpl_table_get_column_min(out_trace_wave_une[qc_i], "Order");
909 if (cpl_error_get_code() != CPL_ERROR_NONE) {
910 continue;
911 }
912 if (current_min < min_order) {
913 min_order = current_min;
914 }
915
916 current_max = cpl_table_get_column_max(out_trace_wave_une[qc_i], "Order");
917 if (cpl_error_get_code() != CPL_ERROR_NONE) {
918 continue;
919 }
920 if (current_max > max_order) {
921 max_order = current_max;
922 }
923 }
924
925 int qc_j;
926 for (qc_j=0; qc_j<qc_une_keywords_size ; qc_j++) {
927 for (int order_id=min_order ; order_id<=max_order ; order_id++) {
928
929 int qc_sizes_une[] = {CR2RES_NB_DETECTORS};
930 char * qc_ref, * qc_res_avg, * qc_res_rmsd ;
931
932 if(strcmp(qc_une_keywords[qc_j], CR2RES_HEADER_QC_OVEREXPOSED)==0) {
933 qc_ref = cpl_sprintf("%s%02d",
934 qc_une_keywords[qc_j], order_id);
935 qc_res_avg = cpl_sprintf("%s%02d %s",
936 qc_une_keywords[qc_j], order_id, "AVG");
937 qc_res_rmsd = cpl_sprintf("%s%02d %s",
938 qc_une_keywords[qc_j], order_id, "RMS");;
939 }
940 else {
941 qc_ref = cpl_sprintf("%s-%02d-%02d",
942 qc_une_keywords[qc_j], order_id, 1);
943 qc_res_avg = cpl_sprintf("%s-%02d-%02d %s",
944 qc_une_keywords[qc_j], order_id, 1, "AVG");
945 qc_res_rmsd = cpl_sprintf("%s-%02d-%02d %s",
946 qc_une_keywords[qc_j], order_id, 1, "RMS");
947 }
948
949 cr2res_qc_calculate_mean_and_rmsd(qc_plists_une, 1, qc_sizes_une,
950 qc_ref, qc_main,
951 qc_res_avg, qc_res_rmsd);
952
953 cpl_free(qc_ref);
954 cpl_free(qc_res_avg);
955 cpl_free(qc_res_rmsd);
956 }
957 }
958
959 if (cpl_error_get_code() != CPL_ERROR_NONE) {
960 cpl_msg_error(__func__, "Failed to set QC params for UNE: %s",
961 cpl_error_get_message_default ( cpl_error_get_code())) ;
962 return -1 ;
963 }
964
965 }
966
967 if (rawframes_fpet != NULL) {
968 plist = cpl_propertylist_load(cpl_frame_get_filename(
969 cpl_frameset_get_position_const(rawframes_fpet, 0)), 0) ;
970 central_wlen = cr2res_pfits_get_cwlen(plist) ;
971 cpl_propertylist_delete(plist);
972
973 /* QC.TILT Computation */
974 if (cr2res_cal_wave_qc_tilt((const cpl_table **)out_trace_wave_fpet,
975 central_wlen, qc_main)) {
976 cpl_msg_warning(__func__, "QC.TILT Computation failed") ;
977 }
978 /* QC.FPI Computation */
979 if (cr2res_cal_wave_qc_fpi(
980 (const cpl_table **)out_extracted_fpet,
981 (const cpl_table **)out_trace_wave_fpet,
982 (const cpl_propertylist **)ext_plist_fpet,
983 central_wlen, qc_main)) {
984 cpl_msg_warning(__func__, "QC.FPI Computation failed") ;
985 }
986 }
987
988 /* Save only the used RAW ? : rawframes_xxx instead of 2nd frameset */
989 /* Beware that the calibration PRO RECi CAL will be missing */
990
991 /* Save Products UNE */
992 if (rawframes_une != NULL &&
993 (rawframes_fpet == NULL || // UNE is the main product
994 save_intermediate_flag)) {
995 out_file = cpl_sprintf("%s_tw_une.fits", RECIPE_STRING) ;
996 cr2res_io_save_TRACE_WAVE(out_file, frameset, frameset, parlist,
997 out_trace_wave_une, qc_main, ext_plist_une,
998 CR2RES_CAL_WAVE_TW_PROCATG, RECIPE_STRING) ;
999 cpl_free(out_file);
1000
1001 out_file = cpl_sprintf("%s_wave_map_une.fits", RECIPE_STRING) ;
1002 cr2res_io_save_WAVE_MAP(out_file, frameset, frameset, parlist,
1003 out_wave_map_une, qc_main, ext_plist_une,
1004 CR2RES_CAL_WAVE_MAP_PROCATG, RECIPE_STRING) ;
1005 cpl_free(out_file);
1006
1007 out_file = cpl_sprintf("%s_extracted_une.fits", RECIPE_STRING) ;
1008 cr2res_io_save_EXTRACT_1D(out_file, frameset, frameset, parlist,
1009 out_extracted_une, qc_main, ext_plist_une,
1010 CR2RES_CAL_WAVE_EXTRACT_1D_PROCATG, RECIPE_STRING) ;
1011 cpl_free(out_file);
1012
1013 out_file = cpl_sprintf("%s_lines_diagnostics_une.fits",
1014 RECIPE_STRING);
1015 cr2res_io_save_LINES_DIAGNOSTICS(out_file, frameset, frameset,
1016 parlist, lines_diagnostics_une, qc_main, ext_plist_une,
1017 CR2RES_CAL_WAVE_LINES_DIAGNOSTICS_PROCATG, RECIPE_STRING) ;
1018 cpl_free(out_file);
1019 }
1020
1021 if (rawframes_fpet != NULL) {
1022 /* Save Products FPET */
1023 out_file = cpl_sprintf("%s_tw_fpet.fits", RECIPE_STRING) ;
1024 cr2res_io_save_TRACE_WAVE(out_file, frameset, frameset, parlist,
1025 out_trace_wave_fpet, qc_main, ext_plist_fpet,
1026 CR2RES_CAL_WAVE_TW_PROCATG, RECIPE_STRING) ;
1027 cpl_free(out_file);
1028
1029 out_file = cpl_sprintf("%s_wave_map_fpet.fits", RECIPE_STRING) ;
1030 cr2res_io_save_WAVE_MAP(out_file, frameset, frameset, parlist,
1031 out_wave_map_fpet, qc_main, ext_plist_fpet,
1032 CR2RES_CAL_WAVE_MAP_PROCATG, RECIPE_STRING) ;
1033 cpl_free(out_file);
1034
1035 out_file = cpl_sprintf("%s_extracted_fpet.fits", RECIPE_STRING) ;
1036 cr2res_io_save_EXTRACT_1D(out_file, frameset, frameset, parlist,
1037 out_extracted_fpet, qc_main, ext_plist_fpet,
1038 CR2RES_CAL_WAVE_EXTRACT_1D_PROCATG, RECIPE_STRING) ;
1039 cpl_free(out_file);
1040
1041 out_file = cpl_sprintf("%s_lines_diagnostics_fpet.fits",
1042 RECIPE_STRING);
1043 cr2res_io_save_LINES_DIAGNOSTICS(out_file, frameset, frameset, parlist,
1044 lines_diagnostics_fpet, qc_main, ext_plist_fpet,
1045 CR2RES_CAL_WAVE_LINES_DIAGNOSTICS_PROCATG, RECIPE_STRING) ;
1046 cpl_free(out_file);
1047 }
1048
1049 /* Free and return */
1050 cpl_propertylist_delete(qc_main) ;
1051 cpl_frameset_delete(rawframes_une) ;
1052 if (rawframes_fpet!=NULL) cpl_frameset_delete(rawframes_fpet) ;
1053 for (i=0 ; i<CR2RES_NB_DETECTORS ; i++) {
1054 if (ext_plist_une[i] != NULL)
1055 cpl_propertylist_delete(ext_plist_une[i]) ;
1056 if (out_trace_wave_une[i] != NULL)
1057 cpl_table_delete(out_trace_wave_une[i]) ;
1058 if (lines_diagnostics_une[i] != NULL)
1059 cpl_table_delete(lines_diagnostics_une[i]) ;
1060 if (out_extracted_une[i] != NULL)
1061 cpl_table_delete(out_extracted_une[i]) ;
1062 if (out_wave_map_une[i] != NULL)
1063 hdrl_image_delete(out_wave_map_une[i]) ;
1064 if (ext_plist_fpet[i] != NULL)
1065 cpl_propertylist_delete(ext_plist_fpet[i]) ;
1066 if (out_trace_wave_fpet[i] != NULL)
1067 cpl_table_delete(out_trace_wave_fpet[i]) ;
1068 if (lines_diagnostics_fpet[i] != NULL)
1069 cpl_table_delete(lines_diagnostics_fpet[i]) ;
1070 if (out_extracted_fpet[i] != NULL)
1071 cpl_table_delete(out_extracted_fpet[i]) ;
1072 if (out_wave_map_fpet[i] != NULL)
1073 hdrl_image_delete(out_wave_map_fpet[i]) ;
1074 }
1075 return (int)cpl_error_get_code();
1076}
1077
1078/*----------------------------------------------------------------------------*/
1127/*----------------------------------------------------------------------------*/
1128static int cr2res_cal_wave_reduce(
1129 const cpl_frameset * rawframes_une,
1130 const cpl_frameset * rawframes_fpet,
1131 const cpl_frame * detlin_frame,
1132 const cpl_frame * master_dark_frame,
1133 const cpl_frame * master_flat_frame,
1134 const cpl_frame * bpm_frame,
1135 const cpl_frame * trace_wave_frame,
1136 const cpl_frame * lines_frame,
1137 int reduce_det,
1138 int reduce_order,
1139 int reduce_trace,
1140 int slit_degree,
1141 int subtract_nolight_rows,
1142 int cosmics,
1143 cr2res_collapse collapse,
1144 int ext_height,
1145 int ext_swath_width,
1146 int ext_oversample,
1147 double ext_smooth_slit,
1148 double ext_smooth_spec,
1149 cr2res_wavecal_type wavecal_type,
1150 int wl_degree,
1151 double wl_start,
1152 double wl_end,
1153 double wl_err,
1154 double wl_shift,
1155 int log_flag,
1156 int fallback_input_wavecal_flag,
1157 int keep_higher_degrees_flag,
1158 int clean_spectrum,
1159 int display,
1160 double display_wmin,
1161 double display_wmax,
1162 cpl_table ** out_trace_wave_une,
1163 cpl_table ** lines_diagnostics_une,
1164 cpl_table ** out_extracted_une,
1165 hdrl_image ** out_wave_map_une,
1166 cpl_propertylist ** ext_plist_une,
1167 cpl_table ** out_trace_wave_fpet,
1168 cpl_table ** lines_diagnostics_fpet,
1169 cpl_table ** out_extracted_fpet,
1170 hdrl_image ** out_wave_map_fpet,
1171 cpl_propertylist ** ext_plist_fpet)
1172{
1173 cpl_vector * ndits_fpet ;
1174 cpl_vector * ndsamples_fpet ;
1175 hdrl_imagelist * in_fpet ;
1176 hdrl_imagelist * in_fpet_calib ;
1177 hdrl_image * collapsed_une ;
1178 hdrl_image * collapsed_fpet ;
1179 cpl_image * contrib ;
1180 cpl_table * tw_in ;
1181 cpl_table * extracted_une ;
1182 cpl_table * slit_func_une ;
1183 hdrl_image * model_master_une ;
1184 cpl_table * extracted_fpet ;
1185 cpl_table * slit_func_fpet ;
1186 hdrl_image * model_master_fpet ;
1187 hdrl_image * wl_map_une_out ;
1188 cpl_table * tw_une_out ;
1189 cpl_table * lines_diagnostics_une_out ;
1190 cpl_table * extracted_une_out ;
1191 cpl_propertylist * qcs_une_out ;
1192 cpl_propertylist * plist_une_out ;
1193 hdrl_image * wl_map_fpet_out ;
1194 cpl_table * tw_fpet_out ;
1195 cpl_table * lines_diagnostics_fpet_out ;
1196 cpl_table * extracted_fpet_out ;
1197 cpl_propertylist * qcs_fpet_out ;
1198 cpl_propertylist * plist_fpet_out ;
1199 cpl_propertylist * plist ;
1200 const char * first_file ;
1201 cpl_polynomial * slit_polya ;
1202 cpl_polynomial * slit_polyb ;
1203 cpl_polynomial * slit_polyc ;
1204 double gain, error_factor ;
1205 int i, order, trace_id, ext_nr,
1206 zp_order_fpet, nb_traces,
1207 grat1_order_fpet ;
1208
1209 /* TODO, make parameters */
1210 int extract_niter = 10;
1211 double extract_kappa = 10;
1212
1213 /* Check Inputs */
1214 if (rawframes_fpet==NULL && rawframes_une==NULL) return -1 ;
1215 if (trace_wave_frame==NULL || out_trace_wave_une==NULL ||
1216 out_trace_wave_fpet==NULL || lines_diagnostics_une == NULL ||
1217 lines_diagnostics_fpet == NULL || out_extracted_une == NULL ||
1218 out_extracted_fpet == NULL || out_wave_map_une == NULL ||
1219 out_wave_map_fpet == NULL || ext_plist_une == NULL ||
1220 ext_plist_fpet == NULL) return -1 ;
1221 if (collapse!=CR2RES_COLLAPSE_MEAN && collapse!=CR2RES_COLLAPSE_MEDIAN)
1222 return -1 ;
1223
1224 /* Get the Gain */
1225 if (reduce_det == 1) gain = CR2RES_GAIN_CHIP1 ;
1226 else if (reduce_det == 2) gain = CR2RES_GAIN_CHIP2 ;
1227 else if (reduce_det == 3) gain = CR2RES_GAIN_CHIP3 ;
1228 else {
1229 cpl_msg_error(__func__, "Failed to get the Gain value") ;
1230 return -1 ;
1231 }
1232
1233
1234 /* Compute the Slit Curvature using the first passed FPET frame */
1235 tw_in = NULL ;
1236 if (rawframes_fpet != NULL) {
1237
1238 const cpl_frame * fpet_frame ;
1239 const char * fpet_fname ;
1240 hdrl_image * fpet_image ;
1241
1242 fpet_frame = cpl_frameset_get_position_const(rawframes_fpet, 0) ;
1243 fpet_fname = cpl_frame_get_filename(fpet_frame) ;
1244 cpl_msg_info(__func__,
1245 "Compute the Slit Curvature from the first FPET frame %s",
1246 fpet_fname) ;
1247 cpl_msg_indent_more() ;
1248
1249 /* Load the trace wave */
1250 cpl_msg_info(__func__, "Load the TRACE WAVE") ;
1251 if ((tw_in = cr2res_io_load_TRACE_WAVE(cpl_frame_get_filename(
1252 trace_wave_frame), reduce_det)) == NULL) {
1253 cpl_msg_error(__func__, "Failed to Load the traces file") ;
1254 cpl_msg_indent_less() ;
1255 return -1 ;
1256 }
1257
1258 /* Load the FPET image of this detector */
1259 cpl_msg_info(__func__, "Load the FPET image") ;
1260 if ((fpet_image = cr2res_io_load_image(fpet_fname, reduce_det))==NULL) {
1261 cpl_msg_error(__func__, "Failed to load image") ;
1262 cpl_table_delete(tw_in) ;
1263 cpl_msg_indent_less() ;
1264 return -1 ;
1265 }
1266
1267 /* Loop over the traces and get the slit curvature */
1268 nb_traces = cpl_table_get_nrow(tw_in) ;
1269 for (i=0 ; i<nb_traces ; i++) {
1270 /* Get Order and trace id */
1271 order = cpl_table_get(tw_in, CR2RES_COL_ORDER, i, NULL) ;
1272 trace_id = cpl_table_get(tw_in, CR2RES_COL_TRACENB, i, NULL) ;
1273
1274 /* Check if this order needs to be skipped */
1275 if (reduce_order > -1 && order != reduce_order) continue ;
1276
1277 /* Check if this trace needs to be skipped */
1278 if (reduce_trace > -1 && trace_id != reduce_trace) continue ;
1279
1280 cpl_msg_info(__func__, "Process Order %d/Trace %d",order,trace_id) ;
1281 cpl_msg_indent_more() ;
1282
1283 /* Call the Slit Curvature Computation */
1284 /* TODO : Should those become parameters ? */
1285 int height = 100 ;
1286 int window = 15 ;
1287 int change_degree = 1;
1288 if (cr2res_slit_curv_compute_order_trace(fpet_image, tw_in, order,
1289 trace_id, height, window, change_degree, slit_degree,
1290 &slit_polya, &slit_polyb, &slit_polyc)) {
1291 cpl_msg_error(__func__, "Failed to compute slit curvature") ;
1292 cpl_table_delete(tw_in) ;
1293 hdrl_image_delete(fpet_image) ;
1294 cpl_msg_indent_less() ;
1295 return -1 ;
1296 }
1297
1298 cpl_array * slit_array ;
1299 /* Fill the SLIT_CURVE_A/B/C for the current trace */
1300 slit_array = cr2res_convert_poly_to_array(slit_polya, 3) ;
1301 cpl_polynomial_delete(slit_polya) ;
1302 cpl_table_set_array(tw_in, CR2RES_COL_SLIT_CURV_A, i, slit_array) ;
1303 cpl_array_delete(slit_array) ;
1304 slit_array = cr2res_convert_poly_to_array(slit_polyb, 3) ;
1305 cpl_polynomial_delete(slit_polyb) ;
1306 cpl_table_set_array(tw_in, CR2RES_COL_SLIT_CURV_B, i, slit_array) ;
1307 cpl_array_delete(slit_array) ;
1308 slit_array = cr2res_convert_poly_to_array(slit_polyc, 3) ;
1309 cpl_polynomial_delete(slit_polyc) ;
1310 cpl_table_set_array(tw_in, CR2RES_COL_SLIT_CURV_C, i, slit_array) ;
1311 cpl_array_delete(slit_array) ;
1312
1313 cpl_msg_indent_less() ;
1314 }
1315 hdrl_image_delete(fpet_image) ;
1316 cpl_msg_indent_less() ;
1317 }
1318
1319 /*
1320 At this point : tw_in is NULL or it is the input TW with the
1321 newly computed slit curvature
1322 */
1323
1324 /* Reduce the UNE */
1325 if (rawframes_une != NULL) {
1326
1327 cpl_vector * dits_une ;
1328 cpl_vector * ndits_une ;
1329 cpl_vector * ndsamples_une ;
1330 hdrl_imagelist * in_une ;
1331 hdrl_imagelist * in_une_calib ;
1332 hdrl_image * first_image ;
1333 int zp_order_une, grat1_order_une ;
1334
1335 cpl_msg_info(__func__, "Reduce %"CPL_SIZE_FORMAT" UNE Frames",
1336 cpl_frameset_get_size(rawframes_une)) ;
1337 cpl_msg_indent_more() ;
1338
1339 /* Load the UNE DITs if necessary */
1340 if (master_dark_frame != NULL)
1341 dits_une = cr2res_io_read_dits(rawframes_une) ;
1342 else
1343 dits_une = NULL ;
1344 if (cpl_msg_get_level() == CPL_MSG_DEBUG && dits_une != NULL)
1345 cpl_vector_dump(dits_une, stdout) ;
1346
1347 /*Load the NDITs */
1348 ndits_une = cr2res_io_read_ndits(rawframes_une);
1349
1350 /* Set the error factor. */
1351 error_factor = gain * cpl_vector_get(ndits_une, 0) *
1352 cpl_frameset_get_size(rawframes_une) ;
1353
1354 for (i=0; i<cpl_vector_get_size(ndits_une); i++){
1355 if (cpl_vector_get(ndits_une,i) != cpl_vector_get(ndits_une, 0))
1356 cpl_msg_warning(__func__, "UNE raw frames have different NDIT! "
1357 "Error spectrum will likely be scaled incorrectly.");
1358 }
1359
1360 /* Load UNE image list */
1361 if ((in_une = cr2res_io_load_image_list_from_set(rawframes_une,
1362 reduce_det)) == NULL) {
1363 cpl_msg_error(__func__, "Cannot load images") ;
1364 if (dits_une != NULL) cpl_vector_delete(dits_une) ;
1365 if (tw_in != NULL) cpl_table_delete(tw_in);
1366 cpl_msg_indent_less() ;
1367 return -1 ;
1368 }
1369 if (hdrl_imagelist_get_size(in_une) !=
1370 cpl_frameset_get_size(rawframes_une)) {
1371 cpl_msg_error(__func__, "Inconsistent number of loaded images") ;
1372 if (dits_une != NULL) cpl_vector_delete(dits_une) ;
1373 hdrl_imagelist_delete(in_une) ;
1374 if (tw_in != NULL) cpl_table_delete(tw_in);
1375 cpl_msg_indent_less() ;
1376 return -1 ;
1377 }
1378 /* Calibrate the UNE images */
1379 ndsamples_une = cr2res_io_read_ndsamples(rawframes_une) ;
1380 in_une_calib = cr2res_calib_imagelist(in_une, reduce_det, 0,
1381 subtract_nolight_rows, 1, cosmics, master_flat_frame,
1382 master_dark_frame, bpm_frame, detlin_frame,
1383 dits_une, ndits_une, ndsamples_une) ;
1384 cpl_vector_delete(ndsamples_une) ;
1385 if (in_une_calib == NULL) {
1386 cpl_msg_error(__func__, "Failed to apply the calibrations") ;
1387 if (dits_une != NULL) cpl_vector_delete(dits_une) ;
1388 if (ndits_une != NULL) cpl_vector_delete(ndits_une) ;
1389 hdrl_imagelist_delete(in_une) ;
1390 if (tw_in != NULL) cpl_table_delete(tw_in);
1391 cpl_msg_indent_less() ;
1392 return -1 ;
1393 }
1394 hdrl_imagelist_delete(in_une) ;
1395 if (dits_une != NULL) cpl_vector_delete(dits_une) ;
1396 if (ndits_une != NULL) cpl_vector_delete(ndits_une) ;
1397
1398 /* Collapse UNE */
1399 contrib = NULL ;
1400 if (collapse == CR2RES_COLLAPSE_MEAN) {
1401 cpl_msg_info(__func__, "Collapse (Mean) the input UNE frames") ;
1402 cpl_msg_indent_more() ;
1403 hdrl_imagelist_collapse_mean(in_une_calib, &collapsed_une,
1404 &contrib) ;
1405 cpl_msg_indent_less() ;
1406 } else if (collapse == CR2RES_COLLAPSE_MEDIAN) {
1407 cpl_msg_info(__func__, "Collapse (Median) the input UNE frames") ;
1408 cpl_msg_indent_more() ;
1409 hdrl_imagelist_collapse_median(in_une_calib, &collapsed_une,
1410 &contrib) ;
1411 cpl_msg_indent_less() ;
1412 } else {
1413 /* Should never happen */
1414 collapsed_une = NULL ;
1415 contrib = NULL ;
1416 }
1417 hdrl_imagelist_delete(in_une_calib) ;
1418 if (contrib != NULL) cpl_image_delete(contrib) ;
1419 if (cpl_error_get_code() != CPL_ERROR_NONE) {
1420 cpl_msg_error(__func__, "Failed Collapse: %d",
1421 cpl_error_get_code());
1422 cpl_msg_indent_less() ;
1423 if (tw_in != NULL) cpl_table_delete(tw_in);
1424 return -1 ;
1425 }
1426
1427 /* Load the trace wave */
1428 cpl_msg_info(__func__, "Load the TRACE WAVE") ;
1429 if (tw_in == NULL)
1430 if ((tw_in = cr2res_io_load_TRACE_WAVE(cpl_frame_get_filename(
1431 trace_wave_frame), reduce_det)) == NULL) {
1432 cpl_msg_error(__func__, "Failed to Load the traces file") ;
1433 hdrl_image_delete(collapsed_une) ;
1434 cpl_msg_indent_less() ;
1435 return -1 ;
1436 }
1437
1438 /* Execute the extraction for UNE */
1439 cpl_msg_info(__func__, "Spectra Extraction UNE") ;
1440 cpl_msg_indent_more() ;
1441 if (cr2res_extract_traces(collapsed_une, tw_in, NULL, NULL, 0,
1442 reduce_order, reduce_trace, CR2RES_EXTR_OPT_CURV,
1443 ext_height, ext_swath_width, ext_oversample, 0.0,
1444 ext_smooth_slit, ext_smooth_spec,
1445 extract_niter, extract_kappa, error_factor, 0, 0, 0,
1446 &extracted_une, &slit_func_une, &model_master_une) == -1) {
1447 cpl_msg_error(__func__, "Failed to extract");
1448 hdrl_image_delete(collapsed_une) ;
1449 cpl_table_delete(tw_in) ;
1450 cpl_msg_indent_less() ;
1451 cpl_msg_indent_less() ;
1452 return -1 ;
1453 }
1454 cpl_msg_indent_less() ;
1455 cpl_table_delete(slit_func_une) ;
1456 hdrl_image_delete(model_master_une) ;
1457 hdrl_image_delete(collapsed_une);
1458
1459 /* Compute the Wavelength Calibration for UNE */
1460 first_file = cpl_frame_get_filename(
1461 cpl_frameset_get_position_const(rawframes_une, 0)) ;
1462 plist = cpl_propertylist_load(first_file, 0) ;
1463 zp_order_une = cr2res_pfits_get_order_zp(plist) ;
1464 grat1_order_une = cr2res_pfits_get_order(plist) ;
1465 cpl_propertylist_delete(plist);
1466
1467 cpl_msg_info(__func__, "Compute the Wavelength for UNE") ;
1468 cpl_msg_indent_more() ;
1469 if (cr2res_wave_apply(tw_in, extracted_une, lines_frame, reduce_order,
1470 reduce_trace, wavecal_type, wl_degree, -1, // xdegree
1471 wl_start, wl_end,
1472 wl_err, wl_shift, log_flag, fallback_input_wavecal_flag,
1473 keep_higher_degrees_flag, clean_spectrum,
1474 display, display_wmin, display_wmax, zp_order_une,
1475 grat1_order_une,
1476 &qcs_une_out,
1477 &lines_diagnostics_une_out,
1478 &extracted_une_out,
1479 &tw_une_out) || cpl_error_get_code()) {
1480 cpl_msg_error(__func__, "Failed to calibrate");
1481 cpl_table_delete(tw_in) ;
1482 cpl_table_delete(extracted_une) ;
1483 cpl_msg_indent_less() ;
1484 cpl_msg_indent_less() ;
1485 return -1 ;
1486 }
1487 cpl_msg_indent_less() ;
1488 cpl_table_delete(extracted_une) ;
1489
1490 /* Generate the Wave Map */
1491 wl_map_une_out = cr2res_wave_gen_wave_map(tw_une_out) ;
1492
1493 /* Load the extension header for saving UNE*/
1494 first_file = cpl_frame_get_filename(
1495 cpl_frameset_get_position_const(rawframes_une, 0)) ;
1496 ext_nr = cr2res_io_get_ext_idx(first_file, reduce_det, 1) ;
1497 plist_une_out = cpl_propertylist_load(first_file, ext_nr) ;
1498 if (plist_une_out == NULL) {
1499 cpl_propertylist_delete(qcs_une_out) ;
1500 cpl_table_delete(tw_in) ;
1501 cpl_table_delete(tw_une_out) ;
1502 cpl_table_delete(lines_diagnostics_une_out) ;
1503 cpl_table_delete(extracted_une_out) ;
1504 hdrl_image_delete(wl_map_une_out) ;
1505 cpl_msg_error(__func__, "Failed to load the plist") ;
1506 cpl_msg_indent_less() ;
1507 return -1 ;
1508 }
1509
1510 /* Load the first RAW image and compute QC OVEREXPOSED */
1511 first_image = cr2res_io_load_image(first_file,reduce_det) ;
1512
1513 /* QC.OVEREXPOSED */
1514 /* Loop on the traces */
1515 nb_traces = cpl_table_get_nrow(tw_in) ;
1516 for (i=0 ; i<nb_traces ; i++) {
1517 /* Get Order and trace id */
1518 order = cpl_table_get(tw_in, CR2RES_COL_ORDER, i, NULL) ;
1519 trace_id = cpl_table_get(tw_in, CR2RES_COL_TRACENB, i, NULL) ;
1520
1521 /* Check if this order needs to be skipped */
1522 if (reduce_order > -1 && order != reduce_order) continue ;
1523
1524 /* Check if this trace needs to be skipped */
1525 if (reduce_trace > -1 && trace_id != reduce_trace) continue ;
1526
1527 double qc_overexposed;
1528 qc_overexposed = cr2res_qc_overexposed(
1529 hdrl_image_get_image(first_image), tw_in, order) ;
1530 char * qc_name = cpl_sprintf("%s%02d",
1531 CR2RES_HEADER_QC_OVEREXPOSED, order) ;
1532 if (qcs_une_out == NULL) qcs_une_out = cpl_propertylist_new() ;
1533 cpl_propertylist_append_double(qcs_une_out, qc_name,qc_overexposed);
1534 cpl_free(qc_name) ;
1535 }
1536 hdrl_image_delete(first_image) ;
1537
1538 if (qcs_une_out != NULL) {
1539 cpl_propertylist_append(plist_une_out, qcs_une_out) ;
1540 cpl_propertylist_delete(qcs_une_out) ;
1541 }
1542 cpl_msg_indent_less() ;
1543 } else {
1544 tw_une_out = NULL ;
1545 lines_diagnostics_une_out = NULL ;
1546 extracted_une_out = NULL ;
1547 wl_map_une_out = NULL ;
1548 plist_une_out = NULL ;
1549 }
1550
1551 /*
1552 At this point : tw_in is NULL or it is the input TW with the
1553 newly computed slit curvature
1554 */
1555
1556 /* Reduce the FPET */
1557 if (rawframes_fpet != NULL && reduce_order != -1){
1558 cpl_msg_warning(__func__,
1559 "Etalon wavecal requires all orders to be used, skipping it.");
1560 if (tw_in != NULL) cpl_table_delete(tw_in);
1561 tw_fpet_out = NULL ;
1562 lines_diagnostics_fpet_out = NULL ;
1563 extracted_fpet_out = NULL ;
1564 wl_map_fpet_out = NULL ;
1565 plist_fpet_out = NULL ;
1566 } else if (rawframes_fpet != NULL) {
1567 cpl_msg_info(__func__, "Reduce %"CPL_SIZE_FORMAT" FPET Frames",
1568 cpl_frameset_get_size(rawframes_fpet)) ;
1569 cpl_msg_indent_more() ;
1570
1571 /* Load FPET image list */
1572 if ((in_fpet = cr2res_io_load_image_list_from_set(rawframes_fpet,
1573 reduce_det)) == NULL) {
1574 cpl_msg_error(__func__, "Cannot load images") ;
1575 cpl_msg_indent_less() ;
1576 if (plist_une_out != NULL) cpl_propertylist_delete(plist_une_out) ;
1577 if (tw_une_out != NULL) cpl_table_delete(tw_une_out) ;
1578 if (lines_diagnostics_une_out != NULL)
1579 cpl_table_delete(lines_diagnostics_une_out) ;
1580 if (extracted_une_out != NULL) cpl_table_delete(extracted_une_out) ;
1581 if (wl_map_une_out != NULL) hdrl_image_delete(wl_map_une_out) ;
1582 return -1 ;
1583 }
1584 if (hdrl_imagelist_get_size(in_fpet) !=
1585 cpl_frameset_get_size(rawframes_fpet)) {
1586 cpl_msg_error(__func__, "Inconsistent number of loaded images") ;
1587 cpl_msg_indent_less() ;
1588 hdrl_imagelist_delete(in_fpet) ;
1589 if (plist_une_out != NULL) cpl_propertylist_delete(plist_une_out) ;
1590 if (tw_une_out != NULL) cpl_table_delete(tw_une_out) ;
1591 if (lines_diagnostics_une_out != NULL)
1592 cpl_table_delete(lines_diagnostics_une_out) ;
1593 if (extracted_une_out != NULL) cpl_table_delete(extracted_une_out) ;
1594 if (wl_map_une_out != NULL) hdrl_image_delete(wl_map_une_out) ;
1595 return -1 ;
1596 }
1597
1598 /* Load FPET NDITs*/
1599 ndits_fpet = cr2res_io_read_ndits(rawframes_fpet);
1600
1601 /* Set the error factor. */
1602 error_factor = gain * cpl_vector_get(ndits_fpet, 0) *
1603 cpl_frameset_get_size(rawframes_une) ;
1604
1605 for (i=0; i<cpl_vector_get_size(ndits_fpet); i++){
1606 if (cpl_vector_get(ndits_fpet,i) != cpl_vector_get(ndits_fpet, 0))
1607 cpl_msg_warning(__func__, "FPET raw frames have different NDIT! "
1608 "Error spectrum will likely be scaled incorrectly.");
1609 }
1610
1611 /* Calibrate the FPET images */
1612 ndsamples_fpet = cr2res_io_read_ndsamples(rawframes_fpet) ;
1613 in_fpet_calib = cr2res_calib_imagelist(in_fpet, reduce_det, 0,
1614 subtract_nolight_rows, 0, cosmics, master_flat_frame,
1615 NULL, bpm_frame, detlin_frame, NULL,
1616 ndits_fpet, ndsamples_fpet) ;
1617 cpl_vector_delete(ndsamples_fpet) ;
1618 if (in_fpet_calib == NULL) {
1619 cpl_msg_error(__func__, "Failed to apply the calibrations") ;
1620 cpl_msg_indent_less() ;
1621 if (in_fpet != NULL) hdrl_imagelist_delete(in_fpet) ;
1622 if (plist_une_out != NULL) cpl_propertylist_delete(plist_une_out) ;
1623 if (tw_une_out != NULL) cpl_table_delete(tw_une_out) ;
1624 if (lines_diagnostics_une_out != NULL)
1625 cpl_table_delete(lines_diagnostics_une_out) ;
1626 if (extracted_une_out != NULL) cpl_table_delete(extracted_une_out) ;
1627 if (wl_map_une_out != NULL) hdrl_image_delete(wl_map_une_out) ;
1628 if (ndits_fpet != NULL) cpl_vector_delete(ndits_fpet) ;
1629 return -1 ;
1630 }
1631 hdrl_imagelist_delete(in_fpet) ;
1632 if (ndits_fpet != NULL) cpl_vector_delete(ndits_fpet) ;
1633
1634 /* Collapse FPET */
1635 contrib = NULL ;
1636 if (collapse == CR2RES_COLLAPSE_MEAN) {
1637 cpl_msg_info(__func__, "Collapse (Mean) the input FPET frames") ;
1638 cpl_msg_indent_more() ;
1639 hdrl_imagelist_collapse_mean(in_fpet_calib, &collapsed_fpet,
1640 &contrib) ;
1641 cpl_msg_indent_less() ;
1642 } else if (collapse == CR2RES_COLLAPSE_MEDIAN) {
1643 cpl_msg_info(__func__, "Collapse (Median) the input FPET frames") ;
1644 cpl_msg_indent_more() ;
1645 hdrl_imagelist_collapse_median(in_fpet_calib, &collapsed_fpet,
1646 &contrib) ;
1647 cpl_msg_indent_less() ;
1648 } else {
1649 /* Should never happen */
1650 collapsed_fpet = NULL ;
1651 contrib = NULL ;
1652 }
1653 hdrl_imagelist_delete(in_fpet_calib) ;
1654 if (contrib != NULL) cpl_image_delete(contrib) ;
1655 if (cpl_error_get_code() != CPL_ERROR_NONE) {
1656 cpl_msg_error(__func__,
1657 "Failed Collapse: %d", cpl_error_get_code()) ;
1658 if (collapsed_fpet != NULL) hdrl_image_delete(collapsed_fpet) ;
1659 if (plist_une_out != NULL) cpl_propertylist_delete(plist_une_out) ;
1660 if (tw_une_out != NULL) cpl_table_delete(tw_une_out) ;
1661 if (lines_diagnostics_une_out != NULL)
1662 cpl_table_delete(lines_diagnostics_une_out) ;
1663 if (extracted_une_out != NULL) cpl_table_delete(extracted_une_out) ;
1664 if (wl_map_une_out != NULL) hdrl_image_delete(wl_map_une_out) ;
1665 cpl_msg_indent_less() ;
1666 return -1 ;
1667 }
1668
1669 /* Get the trace wave */
1670 if (tw_une_out != NULL) {
1671 cpl_msg_info(__func__, "Use the UNE output TRACE WAVE") ;
1672 if (tw_in != NULL) cpl_table_delete(tw_in);
1673 tw_in = cpl_table_duplicate(tw_une_out) ;
1674 } else {
1675 cpl_msg_info(__func__, "Load the TRACE WAVE") ;
1676 if (tw_in == NULL)
1677 if ((tw_in = cr2res_io_load_TRACE_WAVE(cpl_frame_get_filename(
1678 trace_wave_frame), reduce_det)) == NULL) {
1679 cpl_msg_error(__func__, "Failed to Load the traces file") ;
1680 if (collapsed_fpet != NULL)
1681 hdrl_image_delete(collapsed_fpet) ;
1682 if (plist_une_out != NULL)
1683 cpl_propertylist_delete(plist_une_out) ;
1684 if (tw_une_out != NULL) cpl_table_delete(tw_une_out) ;
1685 if (lines_diagnostics_une_out != NULL)
1686 cpl_table_delete(lines_diagnostics_une_out) ;
1687 if (extracted_une_out != NULL)
1688 cpl_table_delete(extracted_une_out) ;
1689 if (wl_map_une_out != NULL)
1690 hdrl_image_delete(wl_map_une_out) ;
1691 cpl_msg_indent_less() ;
1692 return -1 ;
1693 }
1694 }
1695
1696 /* Execute the extraction for FPET */
1697 cpl_msg_info(__func__, "Spectra Extraction FPET") ;
1698 cpl_msg_indent_more() ;
1699 if (cr2res_extract_traces(collapsed_fpet, tw_in, NULL, NULL, 0,
1700 reduce_order, reduce_trace, CR2RES_EXTR_OPT_CURV,
1701 ext_height, ext_swath_width, ext_oversample, 0.0,
1702 ext_smooth_slit, ext_smooth_spec,
1703 extract_niter, extract_kappa, error_factor, 0, 0, 0,
1704 &extracted_fpet, &slit_func_fpet, &model_master_fpet)==-1) {
1705 cpl_msg_error(__func__, "Failed to extract");
1706 cpl_table_delete(tw_in) ;
1707 if (collapsed_fpet != NULL) hdrl_image_delete(collapsed_fpet) ;
1708 if (plist_une_out != NULL) cpl_propertylist_delete(plist_une_out) ;
1709 if (tw_une_out != NULL) cpl_table_delete(tw_une_out) ;
1710 if (lines_diagnostics_une_out != NULL)
1711 cpl_table_delete(lines_diagnostics_une_out) ;
1712 if (extracted_une_out != NULL) cpl_table_delete(extracted_une_out) ;
1713 if (wl_map_une_out != NULL) hdrl_image_delete(wl_map_une_out) ;
1714 cpl_msg_indent_less() ;
1715 cpl_msg_indent_less() ;
1716 return -1 ;
1717 }
1718 cpl_msg_indent_less() ;
1719 cpl_table_delete(slit_func_fpet) ;
1720 hdrl_image_delete(model_master_fpet) ;
1721 hdrl_image_delete(collapsed_fpet);
1722
1723 /* Compute the Wavelength Calibration for FPET */
1724 first_file = cpl_frame_get_filename(
1725 cpl_frameset_get_position_const(rawframes_fpet, 0)) ;
1726 plist = cpl_propertylist_load(first_file, 0) ;
1727 zp_order_fpet = cr2res_pfits_get_order_zp(plist) ;
1728 grat1_order_fpet = cr2res_pfits_get_order(plist) ;
1729 cpl_propertylist_delete(plist);
1730
1731 cpl_msg_info(__func__, "Compute the Wavelength for FPET") ;
1732 cpl_msg_indent_more() ;
1733 cpl_msg_info(__func__, "Running FPET with degrees=4,6") ;
1734
1735 if (cr2res_wave_apply(tw_in, extracted_fpet, NULL, reduce_order,
1736 reduce_trace, CR2RES_ETALON, 4, 6, wl_start, wl_end,
1737 wl_err, wl_shift, log_flag, fallback_input_wavecal_flag,
1738 FALSE, // --keep=FALSE
1739 clean_spectrum,
1740 display, display_wmin, display_wmax, zp_order_fpet,
1741 grat1_order_fpet,
1742 &qcs_fpet_out,
1743 &lines_diagnostics_fpet_out,
1744 &extracted_fpet_out,
1745 &tw_fpet_out) || cpl_error_get_code()) {
1746 cpl_msg_error(__func__, "Failed to calibrate");
1747 cpl_table_delete(extracted_fpet) ;
1748 cpl_table_delete(tw_in) ;
1749 if (plist_une_out != NULL) cpl_propertylist_delete(plist_une_out) ;
1750 if (tw_une_out != NULL) cpl_table_delete(tw_une_out) ;
1751 if (lines_diagnostics_une_out != NULL)
1752 cpl_table_delete(lines_diagnostics_une_out) ;
1753 if (extracted_une_out != NULL) cpl_table_delete(extracted_une_out) ;
1754 if (wl_map_une_out != NULL) hdrl_image_delete(wl_map_une_out) ;
1755 cpl_msg_indent_less() ;
1756 cpl_msg_indent_less() ;
1757 return -1 ;
1758 }
1759 cpl_msg_indent_less() ;
1760 cpl_table_delete(tw_in) ;
1761 cpl_table_delete(extracted_fpet) ;
1762
1763 /* Generate the Wave Map */
1764 wl_map_fpet_out = cr2res_wave_gen_wave_map(tw_fpet_out) ;
1765
1766 /* Load the extension header for saving FPET */
1767 first_file = cpl_frame_get_filename(
1768 cpl_frameset_get_position_const(rawframes_fpet, 0)) ;
1769 ext_nr = cr2res_io_get_ext_idx(first_file, reduce_det, 1) ;
1770 plist_fpet_out = cpl_propertylist_load(first_file, ext_nr) ;
1771 if (plist_fpet_out == NULL) {
1772 cpl_propertylist_delete(qcs_fpet_out) ;
1773 cpl_table_delete(tw_fpet_out) ;
1774 cpl_table_delete(lines_diagnostics_fpet_out) ;
1775 cpl_table_delete(extracted_fpet_out) ;
1776 hdrl_image_delete(wl_map_fpet_out) ;
1777 if (plist_une_out != NULL) cpl_propertylist_delete(plist_une_out) ;
1778 if (tw_une_out != NULL) cpl_table_delete(tw_une_out) ;
1779 if (lines_diagnostics_une_out != NULL)
1780 cpl_table_delete(lines_diagnostics_une_out) ;
1781 if (extracted_une_out != NULL) cpl_table_delete(extracted_une_out) ;
1782 if (wl_map_une_out != NULL) hdrl_image_delete(wl_map_une_out) ;
1783 cpl_msg_indent_less() ;
1784 cpl_msg_error(__func__, "Failed to load the plist") ;
1785 return -1 ;
1786 }
1787 if (qcs_fpet_out != NULL) {
1788 cpl_propertylist_append(plist_fpet_out, qcs_fpet_out) ;
1789 cpl_propertylist_delete(qcs_fpet_out) ;
1790 }
1791 cpl_msg_indent_less() ;
1792 } else {
1793 if (tw_in != NULL) cpl_table_delete(tw_in);
1794 tw_fpet_out = NULL ;
1795 lines_diagnostics_fpet_out = NULL ;
1796 extracted_fpet_out = NULL ;
1797 wl_map_fpet_out = NULL ;
1798 plist_fpet_out = NULL ;
1799 }
1800
1801 /* Return */
1802 *out_trace_wave_une = tw_une_out ;
1803 *lines_diagnostics_une = lines_diagnostics_une_out ;
1804 *out_extracted_une = extracted_une_out ;
1805 *out_wave_map_une = wl_map_une_out ;
1806 *ext_plist_une = plist_une_out ;
1807
1808 *out_trace_wave_fpet = tw_fpet_out ;
1809 *lines_diagnostics_fpet = lines_diagnostics_fpet_out ;
1810 *out_extracted_fpet = extracted_fpet_out ;
1811 *out_wave_map_fpet = wl_map_fpet_out ;
1812 *ext_plist_fpet = plist_fpet_out ;
1813 return 0 ;
1814}
1815
1816/*----------------------------------------------------------------------------*/
1832/*----------------------------------------------------------------------------*/
1833static int cr2res_cal_wave_qc_une_flux(
1834 const cpl_table ** extracted,
1835 const cpl_propertylist ** ext_plist,
1836 double cwlen,
1837 cpl_propertylist * plist)
1838{
1839 const cpl_table * ref_extr ;
1840 const cpl_propertylist * ref_plist ;
1841 cpl_bivector * spec ;
1842 cpl_bivector * spec_err ;
1843 cpl_bivector * intensities ;
1844 double * pintens;
1845 double qc_flux, min_intens, max_intens;
1846 int ref_det, corder_idx ;
1847 cpl_size k;
1848
1849 /* Check Inputs */
1850 if (extracted==NULL || ext_plist == NULL) return -1 ;
1851 for (k=0 ; k< CR2RES_NB_DETECTORS ; k++) {
1852 if (extracted[k] == NULL) return -1 ;
1853 if (ext_plist[k] == NULL) return -1 ;
1854 }
1855
1856 /* Initialize */
1857 qc_flux = -1.0 ;
1858 min_intens = 1000.0 ;
1859 max_intens = 37000.0 ;
1860
1861 /* Detector used depends on the cwlen */
1862 if (cwlen <= 1100.0) ref_det = 0 ;
1863 else ref_det = 1 ;
1864 ref_extr = extracted[ref_det] ;
1865
1866 ref_plist = ext_plist[ref_det] ;
1867
1868 /* Get central order index */
1869 corder_idx = cr2res_pfits_get_order_idx(ref_plist,
1870 CR2RES_DETECTOR_SIZE/2.0) ;
1871
1872 cpl_msg_debug(__func__, "Central wl: %g / Cent. Order IDX: %d",
1873 cwlen, corder_idx) ;
1874
1875 /* Load the spectrum */
1876 cr2res_extract_EXTRACT1D_get_spectrum(ref_extr, corder_idx, 1, &spec,
1877 &spec_err) ;
1878 cpl_bivector_delete(spec_err) ;
1879
1880 /* Plot the spectrum */
1881 if (cpl_msg_get_level() == CPL_MSG_DEBUG) {
1882 cpl_errorstate status;
1883 status = cpl_errorstate_get();
1884 cpl_plot_bivector(
1885 "set grid;set xlabel 'Wavelength (nm)';set ylabel 'Intensity';",
1886 "t 'UNE Spectrum' w lines", "",
1887 spec);
1888 cpl_error_reset();
1889 cpl_errorstate_set(status);
1890 }
1891
1892 /* Compute the lines statistics */
1893 if ((intensities = cr2res_qc_lines_intens_bgd(spec)) == NULL) {
1894 cpl_msg_warning(__func__, "Cannot get the lines statistics") ;
1895 }
1896 else {
1897 cpl_size nintensities, nval;
1898 double * pbgs ;
1899 nintensities = cpl_bivector_get_size(intensities) ;
1900 pintens = cpl_bivector_get_x_data(intensities) ;
1901 pbgs = cpl_bivector_get_y_data(intensities) ;
1902 /* Compute qc_flux of lines in [min_intens,max_intens] */
1903 nval = 0 ;
1904 for (k=0 ; k<nintensities ; k++) {
1905 double val;
1906 val = pintens[k] - pbgs[k] ;
1907 if (val < max_intens && val > min_intens) {
1908 qc_flux += val ;
1909 nval++;
1910 }
1911 }
1912 if (nval > 0) qc_flux /= nval ;
1913 cpl_bivector_delete(intensities) ;
1914 }
1915 cpl_bivector_delete(spec) ;
1916
1917 /* Store the QCs */
1918 cpl_propertylist_append_double(plist, CR2RES_HEADER_QC_UNE_FLUX,
1919 qc_flux) ;
1920
1921 return 0 ;
1922}
1923/*----------------------------------------------------------------------------*/
1946/*----------------------------------------------------------------------------*/
1947static int cr2res_cal_wave_qc_fpi(
1948 const cpl_table ** extracted,
1949 const cpl_table ** tw,
1950 const cpl_propertylist ** ext_plist,
1951 double cwlen,
1952 cpl_propertylist * plist)
1953{
1954 const cpl_table * ref_extr ;
1955 const cpl_table * ref_tw ;
1956 const cpl_propertylist * ref_plist ;
1957 double qc_contrast, qc_flux, qc_separation;
1958 int corder_idx, ref_det, first_fringe_idx,
1959 last_fringe_idx, nb_fringes_between ;
1960 cpl_bivector * spec ;
1961 cpl_bivector * spec_err ;
1962 cpl_vector * xpeaks ;
1963 cpl_vector * lpeaks ;
1964 cpl_vector * fpeaks ;
1965 cpl_vector * xmins ;
1966 cpl_vector * fmins ;
1967 cpl_vector * lmins ;
1968 cpl_vector * tmp_vec ;
1969 cpl_vector * flux_vec ;
1970 cpl_polynomial * wlpoly ;
1971 cpl_size k, npeaks, end_val_search;
1972
1973
1974
1975 /* Check Inputs */
1976 if (extracted==NULL || ext_plist == NULL) return -1 ;
1977 for (k=0 ; k< CR2RES_NB_DETECTORS ; k++) {
1978 if (extracted[k] == NULL) return -1 ;
1979 if (ext_plist[k] == NULL) return -1 ;
1980 }
1981
1982 /* Detector used depends on the cwlen */
1983 if (cwlen <= 1100.0) ref_det = 0 ;
1984 else ref_det = 1 ;
1985 ref_extr = extracted[ref_det] ;
1986 ref_tw = tw[ref_det] ;
1987 ref_plist = ext_plist[ref_det] ;
1988
1989 /* Get central order index */
1990 corder_idx = cr2res_pfits_get_order_idx(ref_plist,
1991 CR2RES_DETECTOR_SIZE/2.0) ;
1992
1993 cpl_msg_debug(__func__, "Central wl: %g / Cent. Order IDX: %d",
1994 cwlen, corder_idx) ;
1995
1996 /* Load the spectrum */
1997 cr2res_extract_EXTRACT1D_get_spectrum(ref_extr, corder_idx, 1, &spec,
1998 &spec_err) ;
1999 cpl_bivector_delete(spec_err) ;
2000
2001 /* Plot the spectrum */
2002 if (cpl_msg_get_level() == CPL_MSG_DEBUG) {
2003 cpl_errorstate status;
2004
2005 status = cpl_errorstate_get();
2006 cpl_plot_bivector(
2007 "set grid;set xlabel 'Wavelength (nm)';set ylabel 'Intensity';",
2008 "t 'FPI Spectrum' w lines", "",
2009 spec);
2010 cpl_error_reset();
2011 cpl_errorstate_set(status);
2012 }
2013
2014 /* Identify the peaks */
2015 if ((xpeaks=cr2res_wave_etalon_measure_fringes(cpl_bivector_get_y(spec)))
2016 == NULL) {
2017 cpl_msg_warning(__func__, "Cannot identify the peaks") ;
2018 cpl_bivector_delete(spec) ;
2019 return -1 ;
2020 }
2021 npeaks = cpl_vector_get_size(xpeaks) ;
2022
2023 if (cpl_msg_get_level() == CPL_MSG_DEBUG) {
2024 cpl_vector_dump(xpeaks, stdout) ;
2025 }
2026
2027 /* Get the wavelength solution for this order */
2028 wlpoly = cr2res_get_trace_wave_poly(ref_tw, CR2RES_COL_WAVELENGTH,
2029 corder_idx, 1) ;
2030 if (wlpoly == NULL) {
2031 cpl_msg_warning(__func__, "Cannot get the wavelength") ;
2032 cpl_bivector_delete(spec) ;
2033 cpl_vector_delete(xpeaks) ;
2034 return -1 ;
2035 }
2036
2037 /* Get the wavelengths at the peak positions */
2038 lpeaks = cr2res_polynomial_eval_vector(wlpoly, xpeaks);
2039
2040 /* Get the Min after each peak */
2041 xmins = cpl_vector_new(npeaks) ;
2042 fpeaks = cpl_vector_new(npeaks);
2043 fmins = cpl_vector_new(npeaks);
2044 for (k=0 ; k<npeaks ; k++) {
2045 cpl_size start_val_search, xminpos;
2046
2047 /* Min of the peak */
2048 start_val_search = (int)cpl_vector_get(xpeaks, k) ;
2049 if (k == npeaks-1) {
2050 end_val_search = start_val_search ;// Not used
2051 } else {
2052 end_val_search = (int)cpl_vector_get(xpeaks, k+1) ;
2053 }
2054 tmp_vec = cpl_vector_extract(cpl_bivector_get_y(spec),
2055 start_val_search-1, end_val_search-1, 1) ;
2056 xminpos = cpl_vector_get_minpos(tmp_vec) ;
2057 cpl_vector_set(fpeaks,k,cpl_vector_get(tmp_vec,0));
2058 cpl_vector_set(fmins,k,cpl_vector_get(tmp_vec,xminpos));
2059 cpl_vector_delete(tmp_vec) ;
2060
2061 cpl_vector_set(xmins, k, start_val_search + xminpos) ;
2062 }
2063
2064 /* Get the wavelengths at the min positions */
2065 lmins = cr2res_polynomial_eval_vector(wlpoly, xmins);
2066 cpl_polynomial_delete(wlpoly) ;
2067
2068 for (k=0 ; k<npeaks ; k++) {
2069 cpl_msg_debug(__func__,
2070 "Peak #%"CPL_SIZE_FORMAT" Pos: %g / WL: %g / Val :%g",
2071 k+1,
2072 cpl_vector_get(xpeaks, k),
2073 cpl_vector_get(lpeaks, k),
2074 cpl_vector_get(cpl_bivector_get_y(spec),
2075 (int)cpl_vector_get(xpeaks, k))) ;
2076 cpl_msg_debug(__func__,
2077 " Min Pos: %g / WL: %g / Val :%g",
2078 cpl_vector_get(xmins, k),
2079 cpl_vector_get(lmins, k),
2080 cpl_vector_get(cpl_bivector_get_y(spec),
2081 (int)cpl_vector_get(xmins, k))) ;
2082 }
2083 cpl_bivector_delete(spec) ;
2084 cpl_vector_delete(xpeaks) ;
2085 cpl_vector_delete(xmins) ;
2086
2087 /* CONTRAST */
2088 /* Ignore first and last peaks */
2089 tmp_vec = cpl_vector_new(npeaks-2) ;
2090 flux_vec = cpl_vector_new(npeaks-2) ;
2091 for (k=1 ; k<npeaks-1 ; k++) {
2092 double fpeak, fmin, contrast_val ;
2093 fpeak = cpl_vector_get(fpeaks, k) ;
2094 fmin = cpl_vector_get(fmins, k) ;
2095 contrast_val = (fpeak-fmin) / (fmin+fpeak) ;
2096 cpl_vector_set(tmp_vec, k-1, contrast_val) ;
2097 cpl_vector_set(flux_vec, k-1, fpeak) ;
2098 }
2099 qc_contrast = cpl_vector_get_mean(tmp_vec) ;
2100 qc_flux = cpl_vector_get_mean(flux_vec) ;
2101 cpl_vector_delete(tmp_vec) ;
2102 cpl_vector_delete(flux_vec) ;
2103 cpl_vector_delete(lmins) ;
2104
2105 /* SEPARATION */
2106 first_fringe_idx = 1 ; // Skip first
2107 last_fringe_idx = npeaks-2 ; // Skip Last
2108 nb_fringes_between = last_fringe_idx - first_fringe_idx + 1 - 2 ;
2109
2110 qc_separation = (cpl_vector_get(lpeaks, last_fringe_idx) -
2111 cpl_vector_get(lpeaks, first_fringe_idx)) /
2112 nb_fringes_between ;
2113
2114 cpl_vector_delete(lpeaks) ;
2115 cpl_vector_delete(fpeaks) ;
2116 cpl_vector_delete(fmins) ;
2117
2118 /* Store the QCs */
2119 cpl_propertylist_append_double(plist, CR2RES_HEADER_QC_FPI_CONTRAST,
2120 qc_contrast) ;
2121 cpl_propertylist_append_double(plist, CR2RES_HEADER_QC_FPI_FLUX,
2122 qc_flux) ;
2123 cpl_propertylist_append_double(plist, CR2RES_HEADER_QC_FPI_SEPARATION,
2124 qc_separation) ;
2125
2126 return 0 ;
2127}
2128
2129/*----------------------------------------------------------------------------*/
2146/*----------------------------------------------------------------------------*/
2147static int cr2res_cal_wave_qc_tilt(
2148 const cpl_table ** tws,
2149 double cwlen,
2150 cpl_propertylist * plist)
2151{
2152 const cpl_table * ref_tw ;
2153 double tilt_avg;
2154 int ref_x, ntilt ;
2155 cpl_size nrows, k ;
2156
2157 /* Check Inputs */
2158 if (tws==NULL) return -1 ;
2159 for (k=0 ; k< CR2RES_NB_DETECTORS ; k++)
2160 if (tws[k] == NULL) return -1 ;
2161
2162 /* Initialize */
2163 ref_x = 1024 ;
2164 tilt_avg = 0.0 ;
2165 ntilt = 0 ;
2166
2167 /* Detector used depends on the cwlen */
2168 if (cwlen <= 1100.0) ref_tw = tws[0] ;
2169 else ref_tw = tws[1] ;
2170
2171 nrows = cpl_table_get_nrow(ref_tw) ;
2172
2173 /* Loop on rows */
2174 for (k=0 ; k<nrows ; k++) {
2175
2176 const cpl_array * tmp_array ;
2177 cpl_polynomial * slit_poly_a ;
2178 cpl_polynomial * slit_poly_b ;
2179 cpl_polynomial * slit_poly_c ;
2180 cpl_polynomial * upper_poly ;
2181 cpl_polynomial * lower_poly ;
2182 cpl_polynomial * slit_curv_poly ;
2183 char * qc_name;
2184 double top_x, top_y, bot_x, bot_y;
2185 double tilt;
2186 int cur_order ;
2187
2188 cur_order = cpl_table_get(ref_tw, CR2RES_COL_ORDER, k, NULL) ;
2189
2190 /* Get the Slit curvature for the trace */
2191 tmp_array = cpl_table_get_array(ref_tw, CR2RES_COL_SLIT_CURV_A, k) ;
2192 slit_poly_a = cr2res_convert_array_to_poly(tmp_array) ;
2193 tmp_array = cpl_table_get_array(ref_tw, CR2RES_COL_SLIT_CURV_B, k) ;
2194 slit_poly_b = cr2res_convert_array_to_poly(tmp_array) ;
2195 tmp_array = cpl_table_get_array(ref_tw, CR2RES_COL_SLIT_CURV_C, k) ;
2196 slit_poly_c = cr2res_convert_array_to_poly(tmp_array) ;
2197
2198 /* Get the Upper Polynomial */
2199 tmp_array = cpl_table_get_array(ref_tw, CR2RES_COL_UPPER, k) ;
2200 upper_poly = cr2res_convert_array_to_poly(tmp_array) ;
2201
2202 /* Get the Lower Polynomial */
2203 tmp_array = cpl_table_get_array(ref_tw, CR2RES_COL_LOWER, k) ;
2204 lower_poly = cr2res_convert_array_to_poly(tmp_array) ;
2205
2206 /* Create the slit curvature polynomial at position ref_x */
2207 slit_curv_poly = cr2res_slit_curv_build_poly(slit_poly_a,
2208 slit_poly_b, slit_poly_c, ref_x) ;
2209 cpl_polynomial_delete(slit_poly_a) ;
2210 cpl_polynomial_delete(slit_poly_b) ;
2211 cpl_polynomial_delete(slit_poly_c) ;
2212
2213 top_y = cpl_polynomial_eval_1d(upper_poly, ref_x, NULL) ;
2214 bot_y = cpl_polynomial_eval_1d(lower_poly, ref_x, NULL) ;
2215 cpl_polynomial_delete(upper_poly) ;
2216 cpl_polynomial_delete(lower_poly) ;
2217
2218 top_x = cpl_polynomial_eval_1d(slit_curv_poly, top_y, NULL) ;
2219 bot_x = cpl_polynomial_eval_1d(slit_curv_poly, bot_y, NULL) ;
2220 cpl_polynomial_delete(slit_curv_poly) ;
2221
2222 /* Compute the tilt */
2223 tilt = atan((top_x-bot_x)/(top_y-bot_y))*180.0/CRIRES_PI;
2224 cpl_msg_debug(__func__,
2225 "TILT : %g / Upper (x,y)=(%g,%g) Lower (x,y)=(%g,%g)",
2226 tilt, top_x, top_y, bot_x, bot_y) ;
2227
2228 /* Store the tilt */
2229 qc_name = cpl_sprintf(CR2RES_HEADER_QC_TILT, cur_order) ;
2230 cpl_propertylist_append_double(plist, qc_name, tilt) ;
2231 cpl_free(qc_name) ;
2232
2233 /* update the tilt average */
2234 tilt_avg += tilt ;
2235 ntilt++ ;
2236 }
2237
2238 /* Store global tilt */
2239 if (ntilt > 0) {
2240 cpl_msg_debug(__func__, "TILT GLOBAL: %g", tilt_avg/ntilt) ;
2241 cpl_propertylist_append_double(plist, CR2RES_HEADER_QC_TILT_GLOBAL,
2242 tilt_avg/ntilt) ;
2243 }
2244
2245 return 0 ;
2246}
cr2res_calib_imagelist
hdrl_imagelist * cr2res_calib_imagelist(const hdrl_imagelist *in, int chip, int clean_bad, int subtract_nolight_rows, int subtract_interorder_column, int cosmics_corr, const cpl_frame *flat, const cpl_frame *dark, const cpl_frame *bpm, const cpl_frame *detlin, const cpl_vector *dits, const cpl_vector *ndits, const cpl_vector *ndsamples)
The images calibration routine for a given chip on a list.
Definition cr2res_calib.c:75
cr2res_dfs_set_groups
cpl_error_code cr2res_dfs_set_groups(cpl_frameset *set)
Set the group as RAW or CALIB in a frameset.
Definition cr2res_dfs.c:53
cr2res_extract_EXTRACT1D_get_spectrum
int cr2res_extract_EXTRACT1D_get_spectrum(const cpl_table *tab, int order, int trace_nb, cpl_bivector **spec, cpl_bivector **spec_err)
Get a Spectrum and its error from the EXTRACT_1D table.
Definition cr2res_extract.c:1289
cr2res_extract_traces
int cr2res_extract_traces(const hdrl_image *img, const cpl_table *traces, const cpl_table *slit_func_in, const cpl_table *blaze_table_in, float blaze_norm, int reduce_order, int reduce_trace, cr2res_extr_method extr_method, int extr_height, int swath_width, int oversample, double pclip, double smooth_slit, double smooth_spec, int niter, double kappa, double error_factor, int display, int disp_order_idx, int disp_trace, cpl_table **extracted, cpl_table **slit_func, hdrl_image **model_master)
Extracts all the passed traces at once.
Definition cr2res_extract.c:186
cr2res_io_load_EXTRACT_1D
cpl_table * cr2res_io_load_EXTRACT_1D(const char *filename, int detector)
Load a table from a EXTRACT_1D.
Definition cr2res_io.c:1326
cr2res_io_find_BPM
const cpl_frame * cr2res_io_find_BPM(const cpl_frameset *in)
Get the first CR2RES_BPM_DRSTYPE frame from a frameset.
Definition cr2res_io.c:369
cr2res_io_find_TRACE_WAVE
const cpl_frame * cr2res_io_find_TRACE_WAVE(const cpl_frameset *in)
Get the first CR2RES_TW_DRSTYPE frame from a frameset.
Definition cr2res_io.c:231
cr2res_io_read_ndsamples
cpl_vector * cr2res_io_read_ndsamples(const cpl_frameset *in)
Get the NDSAMPLES from a frame set.
Definition cr2res_io.c:492
cr2res_io_save_WAVE_MAP
int cr2res_io_save_WAVE_MAP(const char *filename, cpl_frameset *allframes, cpl_frameset *inframes, const cpl_parameterlist *parlist, hdrl_image **data, const cpl_propertylist *qc_list, cpl_propertylist **ext_plist, const char *procatg, const char *recipe)
Save a WAVE_MAP.
Definition cr2res_io.c:1946
cr2res_io_save_TRACE_WAVE
int cr2res_io_save_TRACE_WAVE(const char *filename, cpl_frameset *allframes, cpl_frameset *inframes, const cpl_parameterlist *parlist, cpl_table **tables, const cpl_propertylist *qc_list, cpl_propertylist **ext_plist, const char *procatg, const char *recipe)
Save a TRACE_WAVE.
Definition cr2res_io.c:1700
cr2res_io_read_dits
cpl_vector * cr2res_io_read_dits(const cpl_frameset *in)
Get the DITS from a frame set.
Definition cr2res_io.c:432
cr2res_io_save_EXTRACT_1D
int cr2res_io_save_EXTRACT_1D(const char *filename, cpl_frameset *allframes, cpl_frameset *inframes, const cpl_parameterlist *parlist, cpl_table **tables, const cpl_propertylist *qc_list, cpl_propertylist **ext_plist, const char *procatg, const char *recipe)
Save a 1D extracted spectrum.
Definition cr2res_io.c:1762
cr2res_io_load_image_list_from_set
hdrl_imagelist * cr2res_io_load_image_list_from_set(const cpl_frameset *in, int detector)
Load an hdrl image list from an images frameset.
Definition cr2res_io.c:838
cr2res_io_get_ext_idx
int cr2res_io_get_ext_idx(const char *filename, int detector, int data)
Get the wished extension number for a detector.
Definition cr2res_io.c:674
cr2res_io_load_image
hdrl_image * cr2res_io_load_image(const char *in, int detector)
Load an hdrl image from a image file.
Definition cr2res_io.c:734
cr2res_io_load_TRACE_WAVE
cpl_table * cr2res_io_load_TRACE_WAVE(const char *filename, int detector)
Load a table from a TRACE_WAVE.
Definition cr2res_io.c:1139
cr2res_io_save_LINES_DIAGNOSTICS
int cr2res_io_save_LINES_DIAGNOSTICS(const char *filename, cpl_frameset *allframes, cpl_frameset *inframes, const cpl_parameterlist *parlist, cpl_table **tables, const cpl_propertylist *qc_list, cpl_propertylist **ext_plist, const char *procatg, const char *recipe)
Save a LINES_DIAGNOSTICS.
Definition cr2res_io.c:1731
cr2res_io_read_ndits
cpl_vector * cr2res_io_read_ndits(const cpl_frameset *in)
Get the NDITs from a frame set.
Definition cr2res_io.c:462
cr2res_pfits_get_order_idx
int cr2res_pfits_get_order_idx(const cpl_propertylist *plist, double yposition)
find out the order_idx closest to the passed y position
Definition cr2res_pfits.c:426
cr2res_pfits_get_cwlen
double cr2res_pfits_get_cwlen(const cpl_propertylist *plist)
find out the CWLEN value
Definition cr2res_pfits.c:229
cr2res_qc_overexposed
double cr2res_qc_overexposed(const cpl_image *ima, const cpl_table *tw, int order_idx)
Computes the Overexposed fraction.
Definition cr2res_qc.c:796
cr2res_qc_lines_intens_bgd
cpl_bivector * cr2res_qc_lines_intens_bgd(const cpl_bivector *spec)
Computes the lines intensities / background.
Definition cr2res_qc.c:612
cr2res_slit_curv_compute_order_trace
int cr2res_slit_curv_compute_order_trace(const hdrl_image *img, const cpl_table *trace_wave, const int order, const int trace, const int height, const int window, const cpl_size change_degree, const int slit_degree, cpl_polynomial **slit_poly_a, cpl_polynomial **slit_poly_b, cpl_polynomial **slit_poly_c)
Get the slit curvature directly from the image.
Definition cr2res_slit_curv.c:133
cr2res_slit_curv_build_poly
cpl_polynomial * cr2res_slit_curv_build_poly(cpl_polynomial *slit_poly_a, cpl_polynomial *slit_poly_b, cpl_polynomial *slit_poly_c, cpl_size x)
Create the slit curvature polynomial for x position.
Definition cr2res_slit_curv.c:557
cr2res_get_trace_wave_poly
cpl_polynomial * cr2res_get_trace_wave_poly(const cpl_table *trace_wave, const char *poly_column, int order_idx, int trace_nb)
Get a polynomial from a TRACE_WAVE table.
Definition cr2res_trace.c:685
cr2res_convert_array_to_poly
cpl_polynomial * cr2res_convert_array_to_poly(const cpl_array *arr)
Convert an array to polynomial.
Definition cr2res_utils.c:1305
cr2res_polynomial_eval_vector
cpl_vector * cr2res_polynomial_eval_vector(const cpl_polynomial *poly, const cpl_vector *vec)
Evaluate a polynomial on a vector.
Definition cr2res_utils.c:1059
cr2res_convert_poly_to_array
cpl_array * cr2res_convert_poly_to_array(const cpl_polynomial *poly, int size)
Convert a polynomial to array.
Definition cr2res_utils.c:1339
cr2res_extract_frameset
cpl_frameset * cr2res_extract_frameset(const cpl_frameset *in, const char *tag)
Extract the frames with the given tag from a frameset.
Definition cr2res_utils.c:1179
cr2res_get_license
const char * cr2res_get_license(void)
Get the pipeline copyright and license.
Definition cr2res_utils.c:2362
cr2res_wave_gen_wave_map
hdrl_image * cr2res_wave_gen_wave_map(const cpl_table *trace_wave)
Compute the wavelength map from the trace_wave table.
Definition cr2res_wave.c:1608
cr2res_wave_apply
int cr2res_wave_apply(cpl_table *tw_in, cpl_table *spectra_tab, const cpl_frame *catalog_frame, int reduce_order, int reduce_trace, cr2res_wavecal_type wavecal_type, int degree, int xdegree, double wl_start, double wl_end, double wl_err, double wl_shift, int log_flag, int fallback_input_wavecal_flag, int keep_higher_degrees_flag, int clean_spectrum, int display, double display_wmin, double display_wmax, int zp_order, int grat1_order, cpl_propertylist **qcs, cpl_table **lines_diagnostics, cpl_table **extracted_out, cpl_table **trace_wave_out)
Apply the Wavelength Calibration.
Definition cr2res_wave.c:150
hdrl_image_get_image
cpl_image * hdrl_image_get_image(hdrl_image *himg)
get data as cpl image
Definition hdrl_image.c:105
hdrl_image_delete
void hdrl_image_delete(hdrl_image *himg)
delete hdrl_image
Definition hdrl_image.c:379
hdrl_imagelist_collapse_mean
cpl_error_code hdrl_imagelist_collapse_mean(const hdrl_imagelist *himlist, hdrl_image **out, cpl_image **contrib)
Mean collapsing of image list.
Definition hdrl_imagelist_basic.c:476
hdrl_imagelist_delete
void hdrl_imagelist_delete(hdrl_imagelist *himlist)
Free all memory used by a hdrl_imagelist object including the images.
Definition hdrl_imagelist_io.c:384
hdrl_imagelist_get_size
cpl_size hdrl_imagelist_get_size(const hdrl_imagelist *himlist)
Get the number of images in the imagelist.
Definition hdrl_imagelist_io.c:148
hdrl_imagelist_collapse_median
cpl_error_code hdrl_imagelist_collapse_median(const hdrl_imagelist *himlist, hdrl_image **out, cpl_image **contrib)
Median collapsing of image list.
Definition hdrl_imagelist_basic.c:524
Generated on Tue Oct 13 2026 12:27:26 for CR2RE Pipeline Reference Manual by doxygen 1.13.2