CR2RE Pipeline Reference Manual 1.7.0
cr2res_obs_pol.c
1/*
2 * This file is part of the CR2RES Pipeline
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 <string.h>
29#include <math.h>
30#include <cpl.h>
31
32#include "cr2res_utils.h"
33#include "cr2res_idp.h"
34#include "cr2res_pol.h"
35#include "cr2res_nodding.h"
36#include "cr2res_calib.h"
37#include "cr2res_pfits.h"
38#include "cr2res_dfs.h"
39#include "cr2res_bpm.h"
40#include "cr2res_trace.h"
41#include "cr2res_extract.h"
42#include "cr2res_io.h"
43#include "cr2res_qc.h"
44
45/*-----------------------------------------------------------------------------
46 Define
47 -----------------------------------------------------------------------------*/
48
49#define RECIPE_STRING "cr2res_obs_pol"
50
51/*-----------------------------------------------------------------------------
52 Plugin registration
53 -----------------------------------------------------------------------------*/
54
55int cpl_plugin_get_info(cpl_pluginlist * list);
56
57/*-----------------------------------------------------------------------------
58 Private function prototypes
59 -----------------------------------------------------------------------------*/
60
61static int * cr2res_obs_pol_get_order_numbers(
62 const cpl_table ** extracted,
63 int next,
64 int * norders) ;
65static cpl_table * cr2res_obs_pol_relvar(
66 const hdrl_image * model,
67 const cpl_image * var,
68 int relative,
69 const cpl_table * tw,
70 int extract_height) ;
71static void cr2res_obs_pol_corr_accumulate(
72 cpl_table ** corr,
73 cpl_vector ** intens,
74 cpl_table ** flat_relvar,
75 cpl_table ** bkg_relvar,
76 int order,
77 int n) ;
78static void cr2res_obs_pol_corr_add(
79 cpl_table * pol_spec,
80 const cpl_table * corr,
81 int ngroups) ;
82static int cr2res_obs_pol_check_inputs_validity(
83 const cpl_frameset * rawframes,
84 cpl_size * ngroups) ;
85static int cr2res_obs_pol_reduce(
86 const cpl_frameset * rawframes,
87 const cpl_frame * trace_wave_frame,
88 const cpl_frame * detlin_frame,
89 const cpl_frame * master_dark_frame,
90 const cpl_frame * master_flat_frame,
91 const cpl_frame * bpm_frame,
92 const cpl_frame * blaze_frame,
93 int subtract_nolight_rows,
94 int subtract_interorder_column,
95 int cosmics,
96 int error_method,
97 int extract_oversample,
98 double extract_pclip,
99 int extract_swath_width,
100 int extract_height,
101 double extract_smooth_slit,
102 double extract_smooth_spec,
103 int save_group,
104 int reduce_det,
105 hdrl_image ** in_calib_1_a,
106 hdrl_image ** in_calib_2_a,
107 hdrl_image ** in_calib_3_a,
108 hdrl_image ** in_calib_4_a,
109 cpl_table ** trace_wave_1u_a,
110 cpl_table ** trace_wave_1d_a,
111 cpl_table ** trace_wave_2u_a,
112 cpl_table ** trace_wave_2d_a,
113 cpl_table ** trace_wave_3u_a,
114 cpl_table ** trace_wave_3d_a,
115 cpl_table ** trace_wave_4u_a,
116 cpl_table ** trace_wave_4d_a,
117 cpl_table ** extract1D_1u_a,
118 cpl_table ** extract1D_1d_a,
119 cpl_table ** extract1D_2u_a,
120 cpl_table ** extract1D_2d_a,
121 cpl_table ** extract1D_3u_a,
122 cpl_table ** extract1D_3d_a,
123 cpl_table ** extract1D_4u_a,
124 cpl_table ** extract1D_4d_a,
125 cpl_table ** pol_speca,
126 hdrl_image ** in_calib_1_b,
127 hdrl_image ** in_calib_2_b,
128 hdrl_image ** in_calib_3_b,
129 hdrl_image ** in_calib_4_b,
130 cpl_table ** trace_wave_1u_b,
131 cpl_table ** trace_wave_1d_b,
132 cpl_table ** trace_wave_2u_b,
133 cpl_table ** trace_wave_2d_b,
134 cpl_table ** trace_wave_3u_b,
135 cpl_table ** trace_wave_3d_b,
136 cpl_table ** trace_wave_4u_b,
137 cpl_table ** trace_wave_4d_b,
138 cpl_table ** extract1D_1u_b,
139 cpl_table ** extract1D_1d_b,
140 cpl_table ** extract1D_2u_b,
141 cpl_table ** extract1D_2d_b,
142 cpl_table ** extract1D_3u_b,
143 cpl_table ** extract1D_3d_b,
144 cpl_table ** extract1D_4u_b,
145 cpl_table ** extract1D_4d_b,
146 cpl_table ** pol_specb,
147 cpl_propertylist ** ext_plista,
148 cpl_propertylist ** ext_plistb) ;
149static int cr2res_obs_pol_reduce_one(
150 const cpl_frameset * rawframes,
151 const cpl_frameset * raw_background_frames,
152 const cpl_frame * trace_wave_frame,
153 const cpl_frame * detlin_frame,
154 const cpl_frame * master_dark_frame,
155 const cpl_frame * master_flat_frame,
156 const cpl_frame * bpm_frame,
157 const cpl_frame * blaze_frame,
158 int subtract_nolight_rows,
159 int subtract_interorder_column,
160 int cosmics,
161 int error_method,
162 int extract_oversample,
163 double extract_pclip,
164 int extract_swath_width,
165 int extract_height,
166 double extract_smooth_slit,
167 double extract_smooth_spec,
168 int save_group,
169 int reduce_det,
170 hdrl_image *** in_calib,
171 cpl_table *** tw,
172 cpl_table *** extract1D,
173 cpl_table ** pol_spec,
174 cpl_propertylist ** ext_plist) ;
175static int cr2res_obs_pol_create(cpl_plugin *);
176static int cr2res_obs_pol_exec(cpl_plugin *);
177static int cr2res_obs_pol_destroy(cpl_plugin *);
178static int cr2res_obs_pol(cpl_frameset *, const cpl_parameterlist *);
179
180/*-----------------------------------------------------------------------------
181 Static variables
182 -----------------------------------------------------------------------------*/
183
184static char cr2res_obs_pol_description[] = "\
185Polarimetry Observation \n\
186 The input raw frames are separated in A and B nodding positions. A \n\
187 and B frames are reduced separately. The frames are grouped by blocks \n\
188 of 4 frames. Each block generates 8 extractions. For each order \n\
189 found, 8 spectra are passed to the demodulation functions. \n\
190 The results are stored in polarimetry tables (1 per group of 4 \n\
191 frames). The polarimetry tables are then merged together if there are \n\
192 several group of 4 frames available. \n\
193 \n\
194 Inputs \n\
195 raw.fits " CR2RES_OBS_POLARIMETRY_OTHER_RAW " [4 to 4n] \n\
196 trace.fits " CR2RES_CAL_FLAT_TW_PROCATG " [1] \n\
197 or " CR2RES_CAL_FLAT_TW_MERGED_PROCATG " \n\
198 or " CR2RES_UTIL_TRACE_TW_PROCATG " \n\
199 or " CR2RES_UTIL_WAVE_TW_PROCATG " \n\
200 or " CR2RES_CAL_WAVE_TW_PROCATG " \n\
201 or " CR2RES_UTIL_SLIT_CURV_TW_PROCATG " \n\
202 detlin.fits " CR2RES_CAL_DETLIN_COEFFS_PROCATG " [0 to 1] \n\
203 bpm.fits " CR2RES_CAL_DARK_BPM_PROCATG " [0 to 1] \n\
204 or " CR2RES_CAL_FLAT_BPM_PROCATG " \n\
205 or " CR2RES_CAL_DETLIN_BPM_PROCATG " \n\
206 or " CR2RES_UTIL_BPM_MERGE_PROCATG " \n\
207 or " CR2RES_UTIL_BPM_SPLIT_PROCATG " \n\
208 master_dark.fits " CR2RES_CAL_DARK_MASTER_PROCATG " [0 to 1] \n\
209 master_flat.fits " CR2RES_CAL_FLAT_MASTER_PROCATG " [0 to 1] \n\
210 blaze.fits " CR2RES_CAL_FLAT_EXTRACT_1D_PROCATG " [0 to 1] \n\
211 \n\
212 Outputs \n\
213 cr2res_obs_pol_specA.fits " CR2RES_OBS_POL_SPECA_PROCATG " \n\
214 cr2res_obs_pol_specB.fits " CR2RES_OBS_POL_SPECB_PROCATG " \n\
215 \n\
216 Algorithm \n\
217 loop on detectors d: \n\
218 cr2res_obs_pol_reduce() \n\
219 -> pol_speca(d) \n\
220 -> pol_specb(d) \n\
221 Save pol_speca \n\
222 Save pol_specb \n\
223 \n\
224 cr2res_obs_pol_reduce(): \n\
225 Read the nodding positions of the raw frames \n\
226 Split the rawframes in A and B positions \n\
227 Call cr2res_obs_pol_reduce_one() successively on A and B \n\
228 -> pol_speca \n\
229 -> pol_specb \n\
230 \n\
231 cr2res_obs_pol_reduce_one(): \n\
232 Read the DITs \n\
233 Read the decker positions \n\
234 Load the image list \n\
235 Apply the calibrations to the image list \n\
236 Load the trace_wave \n\
237 For each group of 4 images g: \n\
238 Compute 8 extracted tables (2 per image): \n\
239 1u, 1d, 2u, 2d, 3u, 3d, 4u, 4d \n\
240 where u/d are for up and down \n\
241 1->4 is derived with cr2res_pol_sort_frames() \n\
242 The decker info is used to derive the 8 slit fractions\n\
243 Count norders the number of different orders in those 8 tables \n\
244 [Note : 1 extracted table has 1 spectrum per order] \n\
245 loop on orders o: \n\
246 Get the 8 spectra/wl/error for this order from \n\
247 1u, 1d, 2u, 2d, 3u, 3d, 4u, 4d \n\
248 Call the cr2res_pol_demod_stokes() \n\
249 -> demod_stokes(o, g) \n\
250 Call cr2res_pol_demod_null() \n\
251 -> demod_null(o, g) \n\
252 Call cr2res_pol_demod_intens() \n\
253 -> demod_intens(o, g) \n\
254 Create pol_spec(g) from demod_stokes(g), \n\
255 demod_null(g), \n\
256 demod_intens(g) \n\
257 Merge pol_spec(g) into pol_spec \n\
258 \n\
259 Library functions used \n\
260 cr2res_io_find_TRACE_WAVE() \n\
261 cr2res_io_find_BPM() \n\
262 cr2res_obs_pol_reduce() \n\
263 cr2res_nodding_read_positions() \n\
264 cr2res_combine_nodding_split_frames() \n\
265 cr2res_obs_pol_reduce_one() \n\
266 cr2res_io_read_dits() \n\
267 cr2res_io_read_decker_positions() \n\
268 cr2res_io_load_image_list_from_set() \n\
269 cr2res_calib_imagelist() \n\
270 cr2res_io_load_TRACE_WAVE() \n\
271 cr2res_pol_sort_frames() \n\
272 cr2res_trace_slit_fraction_create() \n\
273 cr2res_trace_new_slit_fraction() \n\
274 cr2res_extract_traces() \n\
275 cr2res_obs_pol_get_order_numbers() \n\
276 cr2res_pol_demod_stokes() \n\
277 cr2res_pol_demod_null() \n\
278 cr2res_pol_demod_intens() \n\
279 cr2res_pol_POL_SPEC_create() \n\
280 cr2res_pol_spec_pol_merge() \n\
281 cr2res_io_save_POL_SPEC() \n\
282";
283
284/*-----------------------------------------------------------------------------
285 Function code
286 -----------------------------------------------------------------------------*/
287
288/*----------------------------------------------------------------------------*/
298/*----------------------------------------------------------------------------*/
299int cpl_plugin_get_info(cpl_pluginlist * list)
300{
301 cpl_recipe * recipe = cpl_calloc(1, sizeof *recipe );
302 cpl_plugin * plugin = &recipe->interface;
303
304 if (cpl_plugin_init(plugin,
305 CPL_PLUGIN_API,
306 CR2RES_BINARY_VERSION,
307 CPL_PLUGIN_TYPE_RECIPE,
308 RECIPE_STRING,
309 "Polarimetry Observation recipe",
310 cr2res_obs_pol_description,
311 CR2RES_PIPELINE_AUTHORS,
312 PACKAGE_BUGREPORT,
313 cr2res_get_license(),
314 cr2res_obs_pol_create,
315 cr2res_obs_pol_exec,
316 cr2res_obs_pol_destroy)) {
317 cpl_msg_error(cpl_func, "Plugin initialization failed");
318 (void)cpl_error_set_where(cpl_func);
319 return 1;
320 }
321
322 if (cpl_pluginlist_append(list, plugin)) {
323 cpl_msg_error(cpl_func, "Error adding plugin to list");
324 (void)cpl_error_set_where(cpl_func);
325 return 1;
326 }
327
328 return 0;
329}
330
331/*----------------------------------------------------------------------------*/
339/*----------------------------------------------------------------------------*/
340static int cr2res_obs_pol_create(cpl_plugin * plugin)
341{
342 cpl_recipe * recipe ;
343 cpl_parameter * p ;
344
345 /* Check that the plugin is part of a valid recipe */
346 if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE)
347 recipe = (cpl_recipe *)plugin;
348 else
349 return -1;
350
351 /* Create the parameters list in the cpl_recipe object */
352 recipe->parameters = cpl_parameterlist_new();
353
354 /* Fill the parameters list */
355 p = cpl_parameter_new_value("cr2res.cr2res_obs_pol.subtract_nolight_rows",
356 CPL_TYPE_BOOL,
357 "Subtract median row from baffled region at detector bottom",
358 "cr2res.cr2res_obs_pol", FALSE);
359 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "subtract_nolight_rows");
360 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
361 cpl_parameterlist_append(recipe->parameters, p);
362
363 p = cpl_parameter_new_value(
364 "cr2res.cr2res_obs_pol.subtract_interorder_column", CPL_TYPE_BOOL,
365 "Subtract column-by-column fit to the pixel values between orders",
366 "cr2res.cr2res_obs_pol", TRUE);
367 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI,
368 "subtract_interorder_column");
369 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
370 cpl_parameterlist_append(recipe->parameters, p);
371
372 p = cpl_parameter_new_value("cr2res.cr2res_obs_pol.cosmics",
373 CPL_TYPE_BOOL, "Find and mark cosmic rays hits as bad",
374 "cr2res.cr2res_obs_pol", FALSE);
375 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "cosmics");
376 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
377 cpl_parameterlist_append(recipe->parameters, p);
378
379 p = cpl_parameter_new_enum("cr2res.cr2res_obs_pol.error_method",
380 CPL_TYPE_STRING, "The 1d extraction error calculation method",
381 "cr2res.cr2res_obs_pol",
382 "Horne", 2, "Poisson", "Horne");
383 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "error_method");
384 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
385 cpl_parameterlist_append(recipe->parameters, p);
386
387 p = cpl_parameter_new_value("cr2res.cr2res_obs_pol.extract_oversample",
388 CPL_TYPE_INT, "factor by which to oversample the extraction",
389 "cr2res.cr2res_obs_pol", 5);
390 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "extract_oversample");
391 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
392 cpl_parameterlist_append(recipe->parameters, p);
393
394 p = cpl_parameter_new_value("cr2res.cr2res_obs_pol.extract_pclip",
395 CPL_TYPE_DOUBLE, "percentage of high and low pixels to clip on first iteration",
396 "cr2res.cr2res_obs_pol", 0.1);
397 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "extract_pclip");
398 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
399 cpl_parameterlist_append(recipe->parameters, p);
400
401 p = cpl_parameter_new_value("cr2res.cr2res_obs_pol.extract_swath_width",
402 CPL_TYPE_INT, "The swath width", "cr2res.cr2res_obs_pol", 800);
403 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "extract_swath_width");
404 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
405 cpl_parameterlist_append(recipe->parameters, p);
406
407 p = cpl_parameter_new_value("cr2res.cr2res_obs_pol.extract_height",
408 CPL_TYPE_INT, "Extraction height", "cr2res.cr2res_obs_pol", -1);
409 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "extract_height");
410 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
411 cpl_parameterlist_append(recipe->parameters, p);
412
413 p = cpl_parameter_new_value("cr2res.cr2res_obs_pol.extract_smooth_spec",
414 CPL_TYPE_DOUBLE, "Smoothing along the spectrum",
415 "cr2res.cr2res_obs_pol", 0.0);
416 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "extract_smooth_spec");
417 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
418 cpl_parameterlist_append(recipe->parameters, p);
419
420 p = cpl_parameter_new_value("cr2res.cr2res_obs_pol.extract_smooth_slit",
421 CPL_TYPE_DOUBLE, "Smoothing along the slit",
422 "cr2res.cr2res_obs_pol", 2.0);
423 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "extract_smooth_slit");
424 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
425 cpl_parameterlist_append(recipe->parameters, p);
426
427 p = cpl_parameter_new_value("cr2res.cr2res_obs_pol.detector", CPL_TYPE_INT,
428 "Only reduce the specified detector", "cr2res.cr2res_obs_pol", 0);
429 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "detector");
430 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
431 cpl_parameterlist_append(recipe->parameters, p);
432
433 p = cpl_parameter_new_value("cr2res.cr2res_obs_nodding.create_idp",
434 CPL_TYPE_BOOL, "Flag to produce IDP files",
435 "cr2res.cr2res_obs_nodding", FALSE);
436 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "idp");
437 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
438 cpl_parameterlist_append(recipe->parameters, p);
439
440 p = cpl_parameter_new_value("cr2res.cr2res_obs_pol.save_group",
441 CPL_TYPE_INT,
442 "Save extra files for the specified group number (0: no save)",
443 "cr2res.cr2res_obs_pol", 0);
444 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI, "save_group");
445 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
446 cpl_parameterlist_append(recipe->parameters, p);
447
448 return 0;
449}
450
451/*----------------------------------------------------------------------------*/
457/*----------------------------------------------------------------------------*/
458static int cr2res_obs_pol_exec(cpl_plugin * plugin)
459{
460 cpl_recipe *recipe;
461
462 /* Get the recipe out of the plugin */
463 if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE)
464 recipe = (cpl_recipe *)plugin;
465 else return -1;
466
467 return cr2res_obs_pol(recipe->frames, recipe->parameters);
468}
469
470/*----------------------------------------------------------------------------*/
476/*----------------------------------------------------------------------------*/
477static int cr2res_obs_pol_destroy(cpl_plugin * plugin)
478{
479 cpl_recipe *recipe;
480
481 /* Get the recipe out of the plugin */
482 if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE)
483 recipe = (cpl_recipe *)plugin;
484 else return -1 ;
485
486 cpl_parameterlist_delete(recipe->parameters);
487 return 0 ;
488}
489
490/*----------------------------------------------------------------------------*/
497/*----------------------------------------------------------------------------*/
498static int cr2res_obs_pol(
499 cpl_frameset * frameset,
500 const cpl_parameterlist * parlist)
501{
502 const cpl_parameter * param ;
503 int extract_oversample, extract_swath_width, cosmics,
504 extract_height, reduce_det, create_idp,
505 subtract_nolight_rows,
506 subtract_interorder_column, save_group,
507 error_method ;
508 double extract_smooth_slit, extract_smooth_spec,
509 extract_pclip ;
510 double barycorr;
511 cpl_frameset * rawframes ;
512 cpl_frameset * raw_flat_frames ;
513 const cpl_frame * trace_wave_frame ;
514 const cpl_frame * detlin_frame ;
515 const cpl_frame * master_dark_frame ;
516 const cpl_frame * master_flat_frame ;
517 const cpl_frame * bpm_frame ;
518 const cpl_frame * blaze_frame ;
519 cpl_propertylist * qc_main ;
520 hdrl_image * in_calib_1_a[CR2RES_NB_DETECTORS] ;
521 hdrl_image * in_calib_2_a[CR2RES_NB_DETECTORS] ;
522 hdrl_image * in_calib_3_a[CR2RES_NB_DETECTORS] ;
523 hdrl_image * in_calib_4_a[CR2RES_NB_DETECTORS] ;
524 hdrl_image * in_calib_1_b[CR2RES_NB_DETECTORS] ;
525 hdrl_image * in_calib_2_b[CR2RES_NB_DETECTORS] ;
526 hdrl_image * in_calib_3_b[CR2RES_NB_DETECTORS] ;
527 hdrl_image * in_calib_4_b[CR2RES_NB_DETECTORS] ;
528 cpl_table * tw_1u_a[CR2RES_NB_DETECTORS] ;
529 cpl_table * tw_1d_a[CR2RES_NB_DETECTORS] ;
530 cpl_table * tw_2u_a[CR2RES_NB_DETECTORS] ;
531 cpl_table * tw_2d_a[CR2RES_NB_DETECTORS] ;
532 cpl_table * tw_3u_a[CR2RES_NB_DETECTORS] ;
533 cpl_table * tw_3d_a[CR2RES_NB_DETECTORS] ;
534 cpl_table * tw_4u_a[CR2RES_NB_DETECTORS] ;
535 cpl_table * tw_4d_a[CR2RES_NB_DETECTORS] ;
536 cpl_table * tw_1u_b[CR2RES_NB_DETECTORS] ;
537 cpl_table * tw_1d_b[CR2RES_NB_DETECTORS] ;
538 cpl_table * tw_2u_b[CR2RES_NB_DETECTORS] ;
539 cpl_table * tw_2d_b[CR2RES_NB_DETECTORS] ;
540 cpl_table * tw_3u_b[CR2RES_NB_DETECTORS] ;
541 cpl_table * tw_3d_b[CR2RES_NB_DETECTORS] ;
542 cpl_table * tw_4u_b[CR2RES_NB_DETECTORS] ;
543 cpl_table * tw_4d_b[CR2RES_NB_DETECTORS] ;
544 cpl_table * extract1D_1u_a[CR2RES_NB_DETECTORS] ;
545 cpl_table * extract1D_1d_a[CR2RES_NB_DETECTORS] ;
546 cpl_table * extract1D_2u_a[CR2RES_NB_DETECTORS] ;
547 cpl_table * extract1D_2d_a[CR2RES_NB_DETECTORS] ;
548 cpl_table * extract1D_3u_a[CR2RES_NB_DETECTORS] ;
549 cpl_table * extract1D_3d_a[CR2RES_NB_DETECTORS] ;
550 cpl_table * extract1D_4u_a[CR2RES_NB_DETECTORS] ;
551 cpl_table * extract1D_4d_a[CR2RES_NB_DETECTORS] ;
552 cpl_table * extract1D_1u_b[CR2RES_NB_DETECTORS] ;
553 cpl_table * extract1D_1d_b[CR2RES_NB_DETECTORS] ;
554 cpl_table * extract1D_2u_b[CR2RES_NB_DETECTORS] ;
555 cpl_table * extract1D_2d_b[CR2RES_NB_DETECTORS] ;
556 cpl_table * extract1D_3u_b[CR2RES_NB_DETECTORS] ;
557 cpl_table * extract1D_3d_b[CR2RES_NB_DETECTORS] ;
558 cpl_table * extract1D_4u_b[CR2RES_NB_DETECTORS] ;
559 cpl_table * extract1D_4d_b[CR2RES_NB_DETECTORS] ;
560 cpl_table * pol_speca[CR2RES_NB_DETECTORS] ;
561 cpl_table * pol_specb[CR2RES_NB_DETECTORS] ;
562 cpl_propertylist * ext_plista[CR2RES_NB_DETECTORS] ;
563 cpl_propertylist * ext_plistb[CR2RES_NB_DETECTORS] ;
564 char * out_file;
565 cpl_table * eop_table ;
566 int det_nr, save_products ;
567
568 /* Initialise */
569
570 /* RETRIEVE INPUT PARAMETERS */
571 param = cpl_parameterlist_find_const(parlist,
572 "cr2res.cr2res_obs_pol.subtract_nolight_rows");
573 subtract_nolight_rows = cpl_parameter_get_bool(param);
574 param = cpl_parameterlist_find_const(parlist,
575 "cr2res.cr2res_obs_pol.subtract_interorder_column");
576 subtract_interorder_column = cpl_parameter_get_bool(param);
577 param = cpl_parameterlist_find_const(parlist,
578 "cr2res.cr2res_obs_pol.cosmics");
579 cosmics = cpl_parameter_get_bool(param);
580 param = cpl_parameterlist_find_const(parlist,
581 "cr2res.cr2res_obs_pol.error_method");
582 if(strcmp(cpl_parameter_get_string(param),"Horne") == 0)
583 error_method = CR2RES_EXTRACT_ERROR_HORNE;
584 else
585 error_method = CR2RES_EXTRACT_ERROR_POISSON;
586 param = cpl_parameterlist_find_const(parlist,
587 "cr2res.cr2res_obs_pol.extract_oversample");
588 extract_oversample = cpl_parameter_get_int(param);
589 param = cpl_parameterlist_find_const(parlist,
590 "cr2res.cr2res_obs_pol.extract_pclip");
591 extract_pclip = cpl_parameter_get_double(param);
592 param = cpl_parameterlist_find_const(parlist,
593 "cr2res.cr2res_obs_pol.extract_swath_width");
594 extract_swath_width = cpl_parameter_get_int(param);
595 param = cpl_parameterlist_find_const(parlist,
596 "cr2res.cr2res_obs_pol.extract_height");
597 extract_height = cpl_parameter_get_int(param);
598 param = cpl_parameterlist_find_const(parlist,
599 "cr2res.cr2res_obs_pol.extract_smooth_slit");
600 extract_smooth_slit = cpl_parameter_get_double(param);
601 param = cpl_parameterlist_find_const(parlist,
602 "cr2res.cr2res_obs_pol.extract_smooth_spec");
603 extract_smooth_spec = cpl_parameter_get_double(param);
604 param = cpl_parameterlist_find_const(parlist,
605 "cr2res.cr2res_obs_pol.detector");
606 reduce_det = cpl_parameter_get_int(param);
607 param = cpl_parameterlist_find_const(parlist,
608 "cr2res.cr2res_obs_nodding.create_idp");
609 create_idp = cpl_parameter_get_bool(param);
610 param = cpl_parameterlist_find_const(parlist,
611 "cr2res.cr2res_obs_pol.save_group");
612 save_group = cpl_parameter_get_int(param);
613
614 /* Identify the RAW and CALIB frames in the input frameset */
615 if (cr2res_dfs_set_groups(frameset)) {
616 cpl_msg_error(__func__, "Cannot identify RAW and CALIB frames") ;
617 cpl_error_set(__func__, CPL_ERROR_ILLEGAL_INPUT) ;
618 return -1 ;
619 }
620
621 /* Get Calibration frames */
622 trace_wave_frame = cr2res_io_find_TRACE_WAVE(frameset) ;
623 if (trace_wave_frame == NULL) {
624 cpl_msg_error(__func__, "Could not find TRACE_WAVE frame") ;
625 return -1 ;
626 }
627 detlin_frame = cpl_frameset_find_const(frameset,
628 CR2RES_CAL_DETLIN_COEFFS_PROCATG);
629 master_dark_frame = cpl_frameset_find_const(frameset,
630 CR2RES_CAL_DARK_MASTER_PROCATG) ;
631 master_flat_frame = cpl_frameset_find_const(frameset,
632 CR2RES_CAL_FLAT_MASTER_PROCATG) ;
633 bpm_frame = cr2res_io_find_BPM(frameset) ;
634 blaze_frame = cpl_frameset_find_const(frameset,
635 CR2RES_CAL_FLAT_EXTRACT_1D_PROCATG) ;
636
637 /* Get the RAW Frames */
638 rawframes = cr2res_extract_frameset(frameset,
639 CR2RES_OBS_POLARIMETRY_OTHER_RAW) ;
640 if (rawframes == NULL) {
641 cpl_msg_error(__func__, "Could not find RAW frames") ;
642 return -1 ;
643 }
644
645 /* Get the RAW flat frames */
646 raw_flat_frames = cr2res_extract_frameset(frameset, CR2RES_FLAT_RAW) ;
647
648 /* Loop on the detectors */
649 for (det_nr=1 ; det_nr<=CR2RES_NB_DETECTORS ; det_nr++) {
650 /* Initialise */
651 in_calib_1_a[det_nr-1] = NULL ;
652 in_calib_2_a[det_nr-1] = NULL ;
653 in_calib_3_a[det_nr-1] = NULL ;
654 in_calib_4_a[det_nr-1] = NULL ;
655 in_calib_1_b[det_nr-1] = NULL ;
656 in_calib_2_b[det_nr-1] = NULL ;
657 in_calib_3_b[det_nr-1] = NULL ;
658 in_calib_4_b[det_nr-1] = NULL ;
659 tw_1u_a[det_nr-1] = NULL ;
660 tw_1d_a[det_nr-1] = NULL ;
661 tw_2u_a[det_nr-1] = NULL ;
662 tw_2d_a[det_nr-1] = NULL ;
663 tw_3u_a[det_nr-1] = NULL ;
664 tw_3d_a[det_nr-1] = NULL ;
665 tw_4u_a[det_nr-1] = NULL ;
666 tw_4d_a[det_nr-1] = NULL ;
667 tw_1u_b[det_nr-1] = NULL ;
668 tw_1d_b[det_nr-1] = NULL ;
669 tw_2u_b[det_nr-1] = NULL ;
670 tw_2d_b[det_nr-1] = NULL ;
671 tw_3u_b[det_nr-1] = NULL ;
672 tw_3d_b[det_nr-1] = NULL ;
673 tw_4u_b[det_nr-1] = NULL ;
674 tw_4d_b[det_nr-1] = NULL ;
675 extract1D_1u_a[det_nr-1] = NULL ;
676 extract1D_1d_a[det_nr-1] = NULL ;
677 extract1D_2u_a[det_nr-1] = NULL ;
678 extract1D_2d_a[det_nr-1] = NULL ;
679 extract1D_3u_a[det_nr-1] = NULL ;
680 extract1D_3d_a[det_nr-1] = NULL ;
681 extract1D_4u_a[det_nr-1] = NULL ;
682 extract1D_4d_a[det_nr-1] = NULL ;
683 extract1D_1u_b[det_nr-1] = NULL ;
684 extract1D_1d_b[det_nr-1] = NULL ;
685 extract1D_2u_b[det_nr-1] = NULL ;
686 extract1D_2d_b[det_nr-1] = NULL ;
687 extract1D_3u_b[det_nr-1] = NULL ;
688 extract1D_3d_b[det_nr-1] = NULL ;
689 extract1D_4u_b[det_nr-1] = NULL ;
690 extract1D_4d_b[det_nr-1] = NULL ;
691 pol_speca[det_nr-1] = NULL ;
692 pol_specb[det_nr-1] = NULL ;
693 ext_plista[det_nr-1] = NULL ;
694 ext_plistb[det_nr-1] = NULL ;
695
696 /* Compute only one detector */
697 if (reduce_det != 0 && det_nr != reduce_det) continue ;
698
699 cpl_msg_info(__func__, "Process Detector %d", det_nr) ;
700 cpl_msg_indent_more() ;
701
702 /* Call the reduction function */
703 if (cr2res_obs_pol_reduce(rawframes,
704 trace_wave_frame, detlin_frame, master_dark_frame,
705 master_flat_frame, bpm_frame, blaze_frame,
706 subtract_nolight_rows, subtract_interorder_column,
707 cosmics, error_method, extract_oversample, extract_pclip,
708 extract_swath_width,
709 extract_height, extract_smooth_slit,
710 extract_smooth_spec, save_group, det_nr,
711 &(in_calib_1_a[det_nr-1]),
712 &(in_calib_2_a[det_nr-1]),
713 &(in_calib_3_a[det_nr-1]),
714 &(in_calib_4_a[det_nr-1]),
715 &(tw_1u_a[det_nr-1]),
716 &(tw_1d_a[det_nr-1]),
717 &(tw_2u_a[det_nr-1]),
718 &(tw_2d_a[det_nr-1]),
719 &(tw_3u_a[det_nr-1]),
720 &(tw_3d_a[det_nr-1]),
721 &(tw_4u_a[det_nr-1]),
722 &(tw_4d_a[det_nr-1]),
723 &(extract1D_1u_a[det_nr-1]),
724 &(extract1D_1d_a[det_nr-1]),
725 &(extract1D_2u_a[det_nr-1]),
726 &(extract1D_2d_a[det_nr-1]),
727 &(extract1D_3u_a[det_nr-1]),
728 &(extract1D_3d_a[det_nr-1]),
729 &(extract1D_4u_a[det_nr-1]),
730 &(extract1D_4d_a[det_nr-1]),
731 &(pol_speca[det_nr-1]),
732 &(in_calib_1_b[det_nr-1]),
733 &(in_calib_2_b[det_nr-1]),
734 &(in_calib_3_b[det_nr-1]),
735 &(in_calib_4_b[det_nr-1]),
736 &(tw_1u_b[det_nr-1]),
737 &(tw_1d_b[det_nr-1]),
738 &(tw_2u_b[det_nr-1]),
739 &(tw_2d_b[det_nr-1]),
740 &(tw_3u_b[det_nr-1]),
741 &(tw_3d_b[det_nr-1]),
742 &(tw_4u_b[det_nr-1]),
743 &(tw_4d_b[det_nr-1]),
744 &(extract1D_1u_b[det_nr-1]),
745 &(extract1D_1d_b[det_nr-1]),
746 &(extract1D_2u_b[det_nr-1]),
747 &(extract1D_2d_b[det_nr-1]),
748 &(extract1D_3u_b[det_nr-1]),
749 &(extract1D_3d_b[det_nr-1]),
750 &(extract1D_4u_b[det_nr-1]),
751 &(extract1D_4d_b[det_nr-1]),
752 &(pol_specb[det_nr-1]),
753 &(ext_plista[det_nr-1]),
754 &(ext_plistb[det_nr-1])) == -1) {
755 cpl_msg_warning(__func__, "Failed to reduce detector %d",
756 det_nr);
757 cpl_error_reset() ;
758 }
759 cpl_msg_indent_less() ;
760 }
761
762 /* Add ESO.DRS.TMID in the Main Header */
763 qc_main = cpl_propertylist_new();
764 cpl_propertylist_append_double(qc_main,
765 CR2RES_HEADER_DRS_TMID,
766 cr2res_utils_get_center_mjd(rawframes)) ;
767
768 /* Add barycentric correction */
769 eop_table = cr2res_io_get_eop_table() ;
770 if (eop_table != NULL) {
771 double ra, dec, mjd_obs, geolon, geolat, geoelev;
772 cpl_propertylist * plist ;
773 plist=cpl_propertylist_load(cpl_frame_get_filename(
774 cpl_frameset_get_position_const(rawframes, 0)), 0) ;
775
776 ra = cpl_propertylist_get_double(plist, "RA") ;
777 dec = cpl_propertylist_get_double(plist, "DEC") ;
778 mjd_obs = cpl_propertylist_get_double(plist, "MJD-OBS") ;
779 geolon = cpl_propertylist_get_double(plist, "ESO TEL GEOLON") ;
780 geolat = cpl_propertylist_get_double(plist, "ESO TEL GEOLAT") ;
781 geoelev = cpl_propertylist_get_double(plist, "ESO TEL GEOELEV") ;
782
783 cpl_propertylist_delete(plist) ;
784
785 barycorr = 0.0 ;
786 if (!cpl_error_get_code()) {
787 double mjd_cen;
788 mjd_cen = cr2res_utils_get_center_mjd(rawframes) ;
789 hdrl_barycorr_compute(ra, dec, eop_table, mjd_obs,
790 (mjd_cen-mjd_obs)*24*3600, geolon, geolat, geoelev,
791 0.0, 0.0, 0.0, 0.0, &barycorr);
792
793 cpl_msg_info(__func__, "Barycentric correction: %g m/s",
794 barycorr);
795 } else {
796 cpl_msg_info(__func__, "Cannot derive Barycentric correction");
797 cpl_error_reset() ;
798 }
799 cpl_table_delete(eop_table) ;
800 cpl_propertylist_append_double(qc_main, CR2RES_HEADER_DRS_BARYCORR,
801 barycorr);
802 }
803
804 if (error_method == CR2RES_EXTRACT_ERROR_HORNE)
805 cpl_propertylist_append_string(qc_main, CR2RES_HEADER_DRS_ERRMETHOD,
806 "Horne") ;
807 else
808 cpl_propertylist_append_string(qc_main, CR2RES_HEADER_DRS_ERRMETHOD,
809 "Poisson") ;
810
811 /* Add QC NUMSAT */
812 cpl_propertylist_append_int(qc_main,
813 CR2RES_HEADER_QC_NUMSAT,
814 cr2res_qc_numsat(rawframes)) ;
815
816 /* Save Extra Products if at least one is returned */
817 save_products = 0 ;
818 for (det_nr=1 ; det_nr<=CR2RES_NB_DETECTORS ; det_nr++) {
819 if (in_calib_1_a[det_nr-1] != NULL) save_products = 1 ;
820 }
821
822 /* Save only the used RAW - fill rawframes with CALIBS */
823 if (trace_wave_frame != NULL)
824 cpl_frameset_insert(rawframes,
825 cpl_frame_duplicate(trace_wave_frame)) ;
826 if (detlin_frame != NULL)
827 cpl_frameset_insert(rawframes,
828 cpl_frame_duplicate(detlin_frame)) ;
829 if (master_dark_frame != NULL)
830 cpl_frameset_insert(rawframes,
831 cpl_frame_duplicate(master_dark_frame)) ;
832 if (master_flat_frame!= NULL)
833 cpl_frameset_insert(rawframes,
834 cpl_frame_duplicate(master_flat_frame)) ;
835 if (bpm_frame!= NULL)
836 cpl_frameset_insert(rawframes,
837 cpl_frame_duplicate(bpm_frame)) ;
838 if (blaze_frame!= NULL)
839 cpl_frameset_insert(rawframes,
840 cpl_frame_duplicate(blaze_frame)) ;
841
842 if (save_products == 1) {
843 /* Calibrated images a */
844 out_file = cpl_sprintf("%s_in_calibA_1.fits", RECIPE_STRING) ;
845 cr2res_io_save_CALIBRATED(out_file, frameset, rawframes, parlist,
846 in_calib_1_a, NULL, ext_plista,
847 CR2RES_OBS_POL_CALIB_A_PROCATG, RECIPE_STRING);
848 cpl_free(out_file);
849 out_file = cpl_sprintf("%s_in_calibA_2.fits", RECIPE_STRING) ;
850 cr2res_io_save_CALIBRATED(out_file, frameset, rawframes, parlist,
851 in_calib_2_a, NULL, ext_plista,
852 CR2RES_OBS_POL_CALIB_A_PROCATG, RECIPE_STRING);
853 cpl_free(out_file);
854 out_file = cpl_sprintf("%s_in_calibA_3.fits", RECIPE_STRING) ;
855 cr2res_io_save_CALIBRATED(out_file, frameset, rawframes, parlist,
856 in_calib_3_a, NULL, ext_plista,
857 CR2RES_OBS_POL_CALIB_A_PROCATG, RECIPE_STRING);
858 cpl_free(out_file);
859 out_file = cpl_sprintf("%s_in_calibA_4.fits", RECIPE_STRING) ;
860 cr2res_io_save_CALIBRATED(out_file, frameset, rawframes, parlist,
861 in_calib_4_a, NULL, ext_plista,
862 CR2RES_OBS_POL_CALIB_A_PROCATG, RECIPE_STRING);
863 cpl_free(out_file);
864
865 /* TW Nodding a */
866 out_file = cpl_sprintf("%s_twA_1u.fits", RECIPE_STRING) ;
867 cr2res_io_save_EXTRACT_1D(out_file, frameset, rawframes, parlist,
868 tw_1u_a, NULL, ext_plista, CR2RES_OBS_POL_TWA_PROCATG,
869 RECIPE_STRING);
870 cpl_free(out_file);
871 out_file = cpl_sprintf("%s_twA_1d.fits", RECIPE_STRING) ;
872 cr2res_io_save_EXTRACT_1D(out_file, frameset, rawframes, parlist,
873 tw_1d_a, NULL, ext_plista, CR2RES_OBS_POL_TWA_PROCATG,
874 RECIPE_STRING);
875 cpl_free(out_file);
876 out_file = cpl_sprintf("%s_twA_2u.fits", RECIPE_STRING) ;
877 cr2res_io_save_EXTRACT_1D(out_file, frameset, rawframes, parlist,
878 tw_2u_a, NULL, ext_plista, CR2RES_OBS_POL_TWA_PROCATG,
879 RECIPE_STRING);
880 cpl_free(out_file);
881 out_file = cpl_sprintf("%s_twA_2d.fits", RECIPE_STRING) ;
882 cr2res_io_save_EXTRACT_1D(out_file, frameset, rawframes, parlist,
883 tw_2d_a, NULL, ext_plista, CR2RES_OBS_POL_TWA_PROCATG,
884 RECIPE_STRING);
885 cpl_free(out_file);
886 out_file = cpl_sprintf("%s_twA_3u.fits", RECIPE_STRING) ;
887 cr2res_io_save_EXTRACT_1D(out_file, frameset, rawframes, parlist,
888 tw_3u_a, NULL, ext_plista, CR2RES_OBS_POL_TWA_PROCATG,
889 RECIPE_STRING);
890 cpl_free(out_file);
891 out_file = cpl_sprintf("%s_twA_3d.fits", RECIPE_STRING) ;
892 cr2res_io_save_EXTRACT_1D(out_file, frameset, rawframes, parlist,
893 tw_3d_a, NULL, ext_plista, CR2RES_OBS_POL_TWA_PROCATG,
894 RECIPE_STRING);
895 cpl_free(out_file);
896 out_file = cpl_sprintf("%s_twA_4u.fits", RECIPE_STRING) ;
897 cr2res_io_save_EXTRACT_1D(out_file, frameset, rawframes, parlist,
898 tw_4u_a, NULL, ext_plista, CR2RES_OBS_POL_TWA_PROCATG,
899 RECIPE_STRING);
900 cpl_free(out_file);
901 out_file = cpl_sprintf("%s_twA_4d.fits", RECIPE_STRING) ;
902 cr2res_io_save_EXTRACT_1D(out_file, frameset, rawframes, parlist,
903 tw_4d_a, NULL, ext_plista, CR2RES_OBS_POL_TWA_PROCATG,
904 RECIPE_STRING);
905 cpl_free(out_file);
906
907 /* Extracted Nodding A */
908 out_file = cpl_sprintf("%s_extractedA_1u.fits", RECIPE_STRING) ;
909 cr2res_io_save_EXTRACT_1D(out_file, frameset, rawframes, parlist,
910 extract1D_1u_a, NULL, ext_plista,
911 CR2RES_OBS_POL_EXTRACTA_PROCATG, RECIPE_STRING);
912 cpl_free(out_file);
913 out_file = cpl_sprintf("%s_extractedA_1d.fits", RECIPE_STRING) ;
914 cr2res_io_save_EXTRACT_1D(out_file, frameset, rawframes, parlist,
915 extract1D_1d_a, NULL, ext_plista,
916 CR2RES_OBS_POL_EXTRACTA_PROCATG, RECIPE_STRING);
917 cpl_free(out_file);
918 out_file = cpl_sprintf("%s_extractedA_2u.fits", RECIPE_STRING) ;
919 cr2res_io_save_EXTRACT_1D(out_file, frameset, rawframes, parlist,
920 extract1D_2u_a, NULL, ext_plista,
921 CR2RES_OBS_POL_EXTRACTA_PROCATG, RECIPE_STRING);
922 cpl_free(out_file);
923 out_file = cpl_sprintf("%s_extractedA_2d.fits", RECIPE_STRING) ;
924 cr2res_io_save_EXTRACT_1D(out_file, frameset, rawframes, parlist,
925 extract1D_2d_a, NULL, ext_plista,
926 CR2RES_OBS_POL_EXTRACTA_PROCATG, RECIPE_STRING);
927 cpl_free(out_file);
928 out_file = cpl_sprintf("%s_extractedA_3u.fits", RECIPE_STRING) ;
929 cr2res_io_save_EXTRACT_1D(out_file, frameset, rawframes, parlist,
930 extract1D_3u_a, NULL, ext_plista,
931 CR2RES_OBS_POL_EXTRACTA_PROCATG, RECIPE_STRING);
932 cpl_free(out_file);
933 out_file = cpl_sprintf("%s_extractedA_3d.fits", RECIPE_STRING) ;
934 cr2res_io_save_EXTRACT_1D(out_file, frameset, rawframes, parlist,
935 extract1D_3d_a, NULL, ext_plista,
936 CR2RES_OBS_POL_EXTRACTA_PROCATG, RECIPE_STRING);
937 cpl_free(out_file);
938 out_file = cpl_sprintf("%s_extractedA_4u.fits", RECIPE_STRING) ;
939 cr2res_io_save_EXTRACT_1D(out_file, frameset, rawframes, parlist,
940 extract1D_4u_a, NULL, ext_plista,
941 CR2RES_OBS_POL_EXTRACTA_PROCATG, RECIPE_STRING);
942 cpl_free(out_file);
943 out_file = cpl_sprintf("%s_extractedA_4d.fits", RECIPE_STRING) ;
944 cr2res_io_save_EXTRACT_1D(out_file, frameset, rawframes, parlist,
945 extract1D_4d_a, NULL, ext_plista,
946 CR2RES_OBS_POL_EXTRACTA_PROCATG, RECIPE_STRING);
947 cpl_free(out_file);
948
949 /* Calibrated images b */
950 out_file = cpl_sprintf("%s_in_calibB_1.fits", RECIPE_STRING) ;
951 cr2res_io_save_CALIBRATED(out_file, frameset, rawframes, parlist,
952 in_calib_1_b, NULL, ext_plista,
953 CR2RES_OBS_POL_CALIB_B_PROCATG, RECIPE_STRING);
954 cpl_free(out_file);
955 out_file = cpl_sprintf("%s_in_calibB_2.fits", RECIPE_STRING) ;
956 cr2res_io_save_CALIBRATED(out_file, frameset, rawframes, parlist,
957 in_calib_2_b, NULL, ext_plista,
958 CR2RES_OBS_POL_CALIB_B_PROCATG, RECIPE_STRING);
959 cpl_free(out_file);
960 out_file = cpl_sprintf("%s_in_calibB_3.fits", RECIPE_STRING) ;
961 cr2res_io_save_CALIBRATED(out_file, frameset, rawframes, parlist,
962 in_calib_3_b, NULL, ext_plista,
963 CR2RES_OBS_POL_CALIB_B_PROCATG, RECIPE_STRING);
964 cpl_free(out_file);
965 out_file = cpl_sprintf("%s_in_calibB_4.fits", RECIPE_STRING) ;
966 cr2res_io_save_CALIBRATED(out_file, frameset, rawframes, parlist,
967 in_calib_4_b, NULL, ext_plista,
968 CR2RES_OBS_POL_CALIB_B_PROCATG, RECIPE_STRING);
969 cpl_free(out_file);
970
971 /* TW Nodding b */
972 out_file = cpl_sprintf("%s_twB_1u.fits", RECIPE_STRING) ;
973 cr2res_io_save_EXTRACT_1D(out_file, frameset, rawframes, parlist,
974 tw_1u_b, NULL, ext_plistb, CR2RES_OBS_POL_TWB_PROCATG,
975 RECIPE_STRING);
976 cpl_free(out_file);
977 out_file = cpl_sprintf("%s_twB_1d.fits", RECIPE_STRING) ;
978 cr2res_io_save_EXTRACT_1D(out_file, frameset, rawframes, parlist,
979 tw_1d_b, NULL, ext_plistb, CR2RES_OBS_POL_TWB_PROCATG,
980 RECIPE_STRING);
981 cpl_free(out_file);
982 out_file = cpl_sprintf("%s_twB_2u.fits", RECIPE_STRING) ;
983 cr2res_io_save_EXTRACT_1D(out_file, frameset, rawframes, parlist,
984 tw_2u_b, NULL, ext_plistb, CR2RES_OBS_POL_TWB_PROCATG,
985 RECIPE_STRING);
986 cpl_free(out_file);
987 out_file = cpl_sprintf("%s_twB_2d.fits", RECIPE_STRING) ;
988 cr2res_io_save_EXTRACT_1D(out_file, frameset, rawframes, parlist,
989 tw_2d_b, NULL, ext_plistb, CR2RES_OBS_POL_TWB_PROCATG,
990 RECIPE_STRING);
991 cpl_free(out_file);
992 out_file = cpl_sprintf("%s_twB_3u.fits", RECIPE_STRING) ;
993 cr2res_io_save_EXTRACT_1D(out_file, frameset, rawframes, parlist,
994 tw_3u_b, NULL, ext_plistb, CR2RES_OBS_POL_TWB_PROCATG,
995 RECIPE_STRING);
996 cpl_free(out_file);
997 out_file = cpl_sprintf("%s_twB_3d.fits", RECIPE_STRING) ;
998 cr2res_io_save_EXTRACT_1D(out_file, frameset, rawframes, parlist,
999 tw_3d_b, NULL, ext_plistb, CR2RES_OBS_POL_TWB_PROCATG,
1000 RECIPE_STRING);
1001 cpl_free(out_file);
1002 out_file = cpl_sprintf("%s_twB_4u.fits", RECIPE_STRING) ;
1003 cr2res_io_save_EXTRACT_1D(out_file, frameset, rawframes, parlist,
1004 tw_4u_b, NULL, ext_plistb, CR2RES_OBS_POL_TWB_PROCATG,
1005 RECIPE_STRING);
1006 cpl_free(out_file);
1007 out_file = cpl_sprintf("%s_twB_4d.fits", RECIPE_STRING) ;
1008 cr2res_io_save_EXTRACT_1D(out_file, frameset, rawframes, parlist,
1009 tw_4d_b, NULL, ext_plistb, CR2RES_OBS_POL_TWB_PROCATG,
1010 RECIPE_STRING);
1011 cpl_free(out_file);
1012
1013 /* Extracted Nodding B */
1014 out_file = cpl_sprintf("%s_extractedB_1u.fits", RECIPE_STRING) ;
1015 cr2res_io_save_EXTRACT_1D(out_file, frameset, rawframes, parlist,
1016 extract1D_1u_b, NULL, ext_plistb,
1017 CR2RES_OBS_POL_EXTRACTB_PROCATG, RECIPE_STRING);
1018 cpl_free(out_file);
1019 out_file = cpl_sprintf("%s_extractedB_1d.fits", RECIPE_STRING) ;
1020 cr2res_io_save_EXTRACT_1D(out_file, frameset, rawframes, parlist,
1021 extract1D_1d_b, NULL, ext_plistb,
1022 CR2RES_OBS_POL_EXTRACTB_PROCATG, RECIPE_STRING);
1023 cpl_free(out_file);
1024 out_file = cpl_sprintf("%s_extractedB_2u.fits", RECIPE_STRING) ;
1025 cr2res_io_save_EXTRACT_1D(out_file, frameset, rawframes, parlist,
1026 extract1D_2u_b, NULL, ext_plistb,
1027 CR2RES_OBS_POL_EXTRACTB_PROCATG, RECIPE_STRING);
1028 cpl_free(out_file);
1029 out_file = cpl_sprintf("%s_extractedB_2d.fits", RECIPE_STRING) ;
1030 cr2res_io_save_EXTRACT_1D(out_file, frameset, rawframes, parlist,
1031 extract1D_2d_b, NULL, ext_plistb,
1032 CR2RES_OBS_POL_EXTRACTB_PROCATG, RECIPE_STRING);
1033 cpl_free(out_file);
1034 out_file = cpl_sprintf("%s_extractedB_3u.fits", RECIPE_STRING) ;
1035 cr2res_io_save_EXTRACT_1D(out_file, frameset, rawframes, parlist,
1036 extract1D_3u_b, NULL, ext_plistb,
1037 CR2RES_OBS_POL_EXTRACTB_PROCATG, RECIPE_STRING);
1038 cpl_free(out_file);
1039 out_file = cpl_sprintf("%s_extractedB_3d.fits", RECIPE_STRING) ;
1040 cr2res_io_save_EXTRACT_1D(out_file, frameset, rawframes, parlist,
1041 extract1D_3d_b, NULL, ext_plistb,
1042 CR2RES_OBS_POL_EXTRACTB_PROCATG, RECIPE_STRING);
1043 cpl_free(out_file);
1044 out_file = cpl_sprintf("%s_extractedB_4u.fits", RECIPE_STRING) ;
1045 cr2res_io_save_EXTRACT_1D(out_file, frameset, rawframes, parlist,
1046 extract1D_4u_b, NULL, ext_plistb,
1047 CR2RES_OBS_POL_EXTRACTB_PROCATG, RECIPE_STRING);
1048 cpl_free(out_file);
1049 out_file = cpl_sprintf("%s_extractedB_4d.fits", RECIPE_STRING) ;
1050 cr2res_io_save_EXTRACT_1D(out_file, frameset, rawframes, parlist,
1051 extract1D_4d_b, NULL, ext_plistb,
1052 CR2RES_OBS_POL_EXTRACTB_PROCATG, RECIPE_STRING);
1053 cpl_free(out_file);
1054 }
1055
1056 /* Polarimetry Spectra */
1057 out_file = cpl_sprintf("%s_pol_specA.fits", RECIPE_STRING) ;
1058 cr2res_io_save_POL_SPEC(out_file, frameset, rawframes, parlist,
1059 pol_speca, qc_main, ext_plista, CR2RES_OBS_POL_SPECA_PROCATG,
1060 RECIPE_STRING) ;
1061 if (create_idp) {
1062 cr2res_idp_save(out_file, frameset, rawframes, parlist,
1063 pol_speca, qc_main, ext_plista,
1064 CR2RES_OBS_POL_EXTRACTA_IDP_PROCATG,
1065 RECIPE_STRING);
1066 }
1067 cpl_free(out_file);
1068
1069 out_file = cpl_sprintf("%s_pol_specB.fits", RECIPE_STRING) ;
1070 cr2res_io_save_POL_SPEC(out_file, frameset, rawframes, parlist,
1071 pol_specb, qc_main, ext_plistb, CR2RES_OBS_POL_SPECB_PROCATG,
1072 RECIPE_STRING) ;
1073 if (create_idp) {
1074 cr2res_idp_save(out_file, frameset, rawframes, parlist,
1075 pol_specb, qc_main, ext_plistb,
1076 CR2RES_OBS_POL_EXTRACTB_IDP_PROCATG,
1077 RECIPE_STRING);
1078 }
1079 //cpl_msg_info(__func__,"ERR: %d", cpl_error_get_code());
1080 cpl_free(out_file);
1081
1082 /* Free */
1083 cpl_frameset_delete(rawframes) ;
1084 cpl_propertylist_delete(qc_main) ;
1085 if (raw_flat_frames != NULL) cpl_frameset_delete(raw_flat_frames) ;
1086 for (det_nr=1 ; det_nr<=CR2RES_NB_DETECTORS ; det_nr++) {
1087 if (in_calib_1_a[det_nr-1] != NULL)
1088 hdrl_image_delete(in_calib_1_a[det_nr-1]) ;
1089 if (in_calib_2_a[det_nr-1] != NULL)
1090 hdrl_image_delete(in_calib_2_a[det_nr-1]) ;
1091 if (in_calib_3_a[det_nr-1] != NULL)
1092 hdrl_image_delete(in_calib_3_a[det_nr-1]) ;
1093 if (in_calib_4_a[det_nr-1] != NULL)
1094 hdrl_image_delete(in_calib_4_a[det_nr-1]) ;
1095 if (in_calib_1_b[det_nr-1] != NULL)
1096 hdrl_image_delete(in_calib_1_b[det_nr-1]) ;
1097 if (in_calib_2_b[det_nr-1] != NULL)
1098 hdrl_image_delete(in_calib_2_b[det_nr-1]) ;
1099 if (in_calib_3_b[det_nr-1] != NULL)
1100 hdrl_image_delete(in_calib_3_b[det_nr-1]) ;
1101 if (in_calib_4_b[det_nr-1] != NULL)
1102 hdrl_image_delete(in_calib_4_b[det_nr-1]) ;
1103 if (tw_1u_a[det_nr-1] != NULL) cpl_table_delete(tw_1u_a[det_nr-1]) ;
1104 if (tw_1d_a[det_nr-1] != NULL) cpl_table_delete(tw_1d_a[det_nr-1]) ;
1105 if (tw_2u_a[det_nr-1] != NULL) cpl_table_delete(tw_2u_a[det_nr-1]) ;
1106 if (tw_2d_a[det_nr-1] != NULL) cpl_table_delete(tw_2d_a[det_nr-1]) ;
1107 if (tw_3u_a[det_nr-1] != NULL) cpl_table_delete(tw_3u_a[det_nr-1]) ;
1108 if (tw_3d_a[det_nr-1] != NULL) cpl_table_delete(tw_3d_a[det_nr-1]) ;
1109 if (tw_4u_a[det_nr-1] != NULL) cpl_table_delete(tw_4u_a[det_nr-1]) ;
1110 if (tw_4d_a[det_nr-1] != NULL) cpl_table_delete(tw_4d_a[det_nr-1]) ;
1111 if (tw_1u_b[det_nr-1] != NULL) cpl_table_delete(tw_1u_b[det_nr-1]) ;
1112 if (tw_1d_b[det_nr-1] != NULL) cpl_table_delete(tw_1d_b[det_nr-1]) ;
1113 if (tw_2u_b[det_nr-1] != NULL) cpl_table_delete(tw_2u_b[det_nr-1]) ;
1114 if (tw_2d_b[det_nr-1] != NULL) cpl_table_delete(tw_2d_b[det_nr-1]) ;
1115 if (tw_3u_b[det_nr-1] != NULL) cpl_table_delete(tw_3u_b[det_nr-1]) ;
1116 if (tw_3d_b[det_nr-1] != NULL) cpl_table_delete(tw_3d_b[det_nr-1]) ;
1117 if (tw_4u_b[det_nr-1] != NULL) cpl_table_delete(tw_4u_b[det_nr-1]) ;
1118 if (tw_4d_b[det_nr-1] != NULL) cpl_table_delete(tw_4d_b[det_nr-1]) ;
1119
1120 if (extract1D_1u_a[det_nr-1] != NULL)
1121 cpl_table_delete(extract1D_1u_a[det_nr-1]) ;
1122 if (extract1D_1d_a[det_nr-1] != NULL)
1123 cpl_table_delete(extract1D_1d_a[det_nr-1]) ;
1124 if (extract1D_2u_a[det_nr-1] != NULL)
1125 cpl_table_delete(extract1D_2u_a[det_nr-1]) ;
1126 if (extract1D_2d_a[det_nr-1] != NULL)
1127 cpl_table_delete(extract1D_2d_a[det_nr-1]) ;
1128 if (extract1D_3u_a[det_nr-1] != NULL)
1129 cpl_table_delete(extract1D_3u_a[det_nr-1]) ;
1130 if (extract1D_3d_a[det_nr-1] != NULL)
1131 cpl_table_delete(extract1D_3d_a[det_nr-1]) ;
1132 if (extract1D_4u_a[det_nr-1] != NULL)
1133 cpl_table_delete(extract1D_4u_a[det_nr-1]) ;
1134 if (extract1D_4d_a[det_nr-1] != NULL)
1135 cpl_table_delete(extract1D_4d_a[det_nr-1]) ;
1136 if (extract1D_1u_b[det_nr-1] != NULL)
1137 cpl_table_delete(extract1D_1u_b[det_nr-1]) ;
1138 if (extract1D_1d_b[det_nr-1] != NULL)
1139 cpl_table_delete(extract1D_1d_b[det_nr-1]) ;
1140 if (extract1D_2u_b[det_nr-1] != NULL)
1141 cpl_table_delete(extract1D_2u_b[det_nr-1]) ;
1142 if (extract1D_2d_b[det_nr-1] != NULL)
1143 cpl_table_delete(extract1D_2d_b[det_nr-1]) ;
1144 if (extract1D_3u_b[det_nr-1] != NULL)
1145 cpl_table_delete(extract1D_3u_b[det_nr-1]) ;
1146 if (extract1D_3d_b[det_nr-1] != NULL)
1147 cpl_table_delete(extract1D_3d_b[det_nr-1]) ;
1148 if (extract1D_4u_b[det_nr-1] != NULL)
1149 cpl_table_delete(extract1D_4u_b[det_nr-1]) ;
1150 if (extract1D_4d_b[det_nr-1] != NULL)
1151 cpl_table_delete(extract1D_4d_b[det_nr-1]) ;
1152 if (pol_speca[det_nr-1] != NULL) cpl_table_delete(pol_speca[det_nr-1]) ;
1153 if (pol_specb[det_nr-1] != NULL) cpl_table_delete(pol_specb[det_nr-1]) ;
1154 if (ext_plista[det_nr-1] != NULL)
1155 cpl_propertylist_delete(ext_plista[det_nr-1]) ;
1156 if (ext_plistb[det_nr-1] != NULL)
1157 cpl_propertylist_delete(ext_plistb[det_nr-1]) ;
1158 }
1159
1160 return (int)cpl_error_get_code();
1161}
1162
1163/*----------------------------------------------------------------------------*/
1196/*----------------------------------------------------------------------------*/
1197static int cr2res_obs_pol_reduce(
1198 const cpl_frameset * rawframes,
1199 const cpl_frame * trace_wave_frame,
1200 const cpl_frame * detlin_frame,
1201 const cpl_frame * master_dark_frame,
1202 const cpl_frame * master_flat_frame,
1203 const cpl_frame * bpm_frame,
1204 const cpl_frame * blaze_frame,
1205 int subtract_nolight_rows,
1206 int subtract_interorder_column,
1207 int cosmics,
1208 int error_method,
1209 int extract_oversample,
1210 double extract_pclip,
1211 int extract_swath_width,
1212 int extract_height,
1213 double extract_smooth_slit,
1214 double extract_smooth_spec,
1215 int save_group,
1216 int reduce_det,
1217 hdrl_image ** in_calib_1_a,
1218 hdrl_image ** in_calib_2_a,
1219 hdrl_image ** in_calib_3_a,
1220 hdrl_image ** in_calib_4_a,
1221 cpl_table ** trace_wave_1u_a,
1222 cpl_table ** trace_wave_1d_a,
1223 cpl_table ** trace_wave_2u_a,
1224 cpl_table ** trace_wave_2d_a,
1225 cpl_table ** trace_wave_3u_a,
1226 cpl_table ** trace_wave_3d_a,
1227 cpl_table ** trace_wave_4u_a,
1228 cpl_table ** trace_wave_4d_a,
1229 cpl_table ** extract1D_1u_a,
1230 cpl_table ** extract1D_1d_a,
1231 cpl_table ** extract1D_2u_a,
1232 cpl_table ** extract1D_2d_a,
1233 cpl_table ** extract1D_3u_a,
1234 cpl_table ** extract1D_3d_a,
1235 cpl_table ** extract1D_4u_a,
1236 cpl_table ** extract1D_4d_a,
1237 cpl_table ** pol_speca,
1238 hdrl_image ** in_calib_1_b,
1239 hdrl_image ** in_calib_2_b,
1240 hdrl_image ** in_calib_3_b,
1241 hdrl_image ** in_calib_4_b,
1242 cpl_table ** trace_wave_1u_b,
1243 cpl_table ** trace_wave_1d_b,
1244 cpl_table ** trace_wave_2u_b,
1245 cpl_table ** trace_wave_2d_b,
1246 cpl_table ** trace_wave_3u_b,
1247 cpl_table ** trace_wave_3d_b,
1248 cpl_table ** trace_wave_4u_b,
1249 cpl_table ** trace_wave_4d_b,
1250 cpl_table ** extract1D_1u_b,
1251 cpl_table ** extract1D_1d_b,
1252 cpl_table ** extract1D_2u_b,
1253 cpl_table ** extract1D_2d_b,
1254 cpl_table ** extract1D_3u_b,
1255 cpl_table ** extract1D_3d_b,
1256 cpl_table ** extract1D_4u_b,
1257 cpl_table ** extract1D_4d_b,
1258 cpl_table ** pol_specb,
1259 cpl_propertylist ** ext_plista,
1260 cpl_propertylist ** ext_plistb)
1261{
1262 cpl_frameset * rawframes_a ;
1263 cpl_frameset * rawframes_b ;
1264 cr2res_nodding_pos * nod_positions ;
1265 hdrl_image ** in_calib_a_loc ;
1266 hdrl_image ** in_calib_b_loc ;
1267 cpl_table ** trace_wave_a_loc ;
1268 cpl_table ** trace_wave_b_loc ;
1269 cpl_table ** extract1D_a_loc ;
1270 cpl_table ** extract1D_b_loc ;
1271 cpl_table * pol_speca_loc ;
1272 cpl_table * pol_specb_loc ;
1273 cpl_propertylist * ext_plista_loc ;
1274 cpl_propertylist * ext_plistb_loc ;
1275 cpl_size ngroups ;
1276
1277 /* Check Inputs */
1278 if (pol_speca == NULL || pol_specb == NULL || ext_plista == NULL ||
1279 ext_plistb == NULL || rawframes == NULL ||
1280 trace_wave_frame == NULL) return -1 ;
1281
1282 if (cr2res_obs_pol_check_inputs_validity(rawframes, &ngroups) != 1) {
1283 cpl_msg_error(__func__, "Validity Checks failed") ;
1284 return -1 ;
1285 }
1286
1287 /* Initialize */
1288 *in_calib_1_a = *in_calib_2_a = *in_calib_3_a = *in_calib_4_a = NULL ;
1289 *trace_wave_1u_a = *trace_wave_1d_a = *trace_wave_2u_a =
1290 *trace_wave_2d_a = *trace_wave_3u_a = *trace_wave_3d_a =
1291 *trace_wave_4u_a = *trace_wave_4d_a = NULL ;
1292 *extract1D_1u_a = *extract1D_1d_a = *extract1D_2u_a =
1293 *extract1D_2d_a = *extract1D_3u_a = *extract1D_3d_a =
1294 *extract1D_4u_a = *extract1D_4d_a = NULL ;
1295 *in_calib_1_b = *in_calib_2_b = *in_calib_3_b = *in_calib_4_b = NULL ;
1296 *trace_wave_1u_b = *trace_wave_1d_b = *trace_wave_2u_b =
1297 *trace_wave_2d_b = *trace_wave_3u_b = *trace_wave_3d_b =
1298 *trace_wave_4u_b = *trace_wave_4d_b = NULL ;
1299 *extract1D_1u_b = *extract1D_1d_b = *extract1D_2u_b =
1300 *extract1D_2d_b = *extract1D_3u_b = *extract1D_3d_b =
1301 *extract1D_4u_b = *extract1D_4d_b = NULL ;
1302
1303 /* Get the Nodding positions */
1304 nod_positions = cr2res_nodding_read_positions(rawframes) ;
1305
1306 /* Split the frames */
1307 if (cr2res_combine_nodding_split_frames(rawframes, nod_positions,
1308 &rawframes_a, &rawframes_b)) {
1309 cpl_msg_error(__func__, "Failed to split the nodding positions") ;
1310 cpl_free(nod_positions) ;
1311 return -1 ;
1312 }
1313 cpl_free(nod_positions) ;
1314
1315 /* Reduce A position */
1316 cpl_msg_info(__func__, "Compute Polarimetry for nodding A position") ;
1317 cpl_msg_indent_more() ;
1318 if (cr2res_obs_pol_reduce_one(rawframes_a, rawframes_b,
1319 trace_wave_frame, detlin_frame, master_dark_frame,
1320 master_flat_frame, bpm_frame, blaze_frame,
1321 subtract_nolight_rows, subtract_interorder_column, cosmics,
1322 error_method, extract_oversample, extract_pclip,
1323 extract_swath_width, extract_height,
1324 extract_smooth_slit, extract_smooth_spec, save_group,
1325 reduce_det,
1326 &in_calib_a_loc, &trace_wave_a_loc, &extract1D_a_loc,
1327 &pol_speca_loc, &ext_plista_loc) == -1) {
1328 cpl_msg_warning(__func__, "Failed to Reduce A nodding frames") ;
1329 }
1330 cpl_msg_indent_less() ;
1331
1332 /* Reduce B position */
1333 cpl_msg_info(__func__, "Compute Polarimetry for nodding B position") ;
1334 cpl_msg_indent_more() ;
1335 if (cr2res_obs_pol_reduce_one(rawframes_b, rawframes_a,
1336 trace_wave_frame, detlin_frame, master_dark_frame,
1337 master_flat_frame, bpm_frame, blaze_frame,
1338 subtract_nolight_rows, subtract_interorder_column, cosmics,
1339 error_method, extract_oversample, extract_pclip,
1340 extract_swath_width, extract_height,
1341 extract_smooth_slit, extract_smooth_spec, save_group,
1342 reduce_det,
1343 &in_calib_b_loc, &trace_wave_b_loc, &extract1D_b_loc,
1344 &pol_specb_loc, &ext_plistb_loc) == -1) {
1345 cpl_msg_warning(__func__, "Failed to Reduce B nodding frames") ;
1346 }
1347 cpl_msg_indent_less() ;
1348 if (rawframes_a != NULL) cpl_frameset_delete(rawframes_a);
1349 if (rawframes_b != NULL) cpl_frameset_delete(rawframes_b);
1350
1351 /* Return */
1352 if (in_calib_a_loc != NULL) {
1353 *in_calib_1_a = in_calib_a_loc[0] ;
1354 *in_calib_2_a = in_calib_a_loc[1] ;
1355 *in_calib_3_a = in_calib_a_loc[2] ;
1356 *in_calib_4_a = in_calib_a_loc[3] ;
1357 cpl_free(in_calib_a_loc) ;
1358 }
1359 if (trace_wave_a_loc != NULL) {
1360 *trace_wave_1u_a = trace_wave_a_loc[0] ;
1361 *trace_wave_1d_a = trace_wave_a_loc[1] ;
1362 *trace_wave_2u_a = trace_wave_a_loc[2] ;
1363 *trace_wave_2d_a = trace_wave_a_loc[3] ;
1364 *trace_wave_3u_a = trace_wave_a_loc[4] ;
1365 *trace_wave_3d_a = trace_wave_a_loc[5] ;
1366 *trace_wave_4u_a = trace_wave_a_loc[6] ;
1367 *trace_wave_4d_a = trace_wave_a_loc[7] ;
1368 cpl_free(trace_wave_a_loc) ;
1369 }
1370 if (extract1D_a_loc != NULL) {
1371 *extract1D_1u_a = extract1D_a_loc[0] ;
1372 *extract1D_1d_a = extract1D_a_loc[1] ;
1373 *extract1D_2u_a = extract1D_a_loc[2] ;
1374 *extract1D_2d_a = extract1D_a_loc[3] ;
1375 *extract1D_3u_a = extract1D_a_loc[4] ;
1376 *extract1D_3d_a = extract1D_a_loc[5] ;
1377 *extract1D_4u_a = extract1D_a_loc[6] ;
1378 *extract1D_4d_a = extract1D_a_loc[7] ;
1379 cpl_free(extract1D_a_loc) ;
1380 }
1381 if (in_calib_b_loc != NULL) {
1382 *in_calib_1_b = in_calib_b_loc[0] ;
1383 *in_calib_2_b = in_calib_b_loc[1] ;
1384 *in_calib_3_b = in_calib_b_loc[2] ;
1385 *in_calib_4_b = in_calib_b_loc[3] ;
1386 cpl_free(in_calib_b_loc) ;
1387 }
1388 if (trace_wave_b_loc != NULL) {
1389 *trace_wave_1u_b = trace_wave_b_loc[0] ;
1390 *trace_wave_1d_b = trace_wave_b_loc[1] ;
1391 *trace_wave_2u_b = trace_wave_b_loc[2] ;
1392 *trace_wave_2d_b = trace_wave_b_loc[3] ;
1393 *trace_wave_3u_b = trace_wave_b_loc[4] ;
1394 *trace_wave_3d_b = trace_wave_b_loc[5] ;
1395 *trace_wave_4u_b = trace_wave_b_loc[6] ;
1396 *trace_wave_4d_b = trace_wave_b_loc[7] ;
1397 cpl_free(trace_wave_b_loc) ;
1398 }
1399 if (extract1D_b_loc != NULL) {
1400 *extract1D_1u_b = extract1D_b_loc[0] ;
1401 *extract1D_1d_b = extract1D_b_loc[1] ;
1402 *extract1D_2u_b = extract1D_b_loc[2] ;
1403 *extract1D_2d_b = extract1D_b_loc[3] ;
1404 *extract1D_3u_b = extract1D_b_loc[4] ;
1405 *extract1D_3d_b = extract1D_b_loc[5] ;
1406 *extract1D_4u_b = extract1D_b_loc[6] ;
1407 *extract1D_4d_b = extract1D_b_loc[7] ;
1408 cpl_free(extract1D_b_loc) ;
1409 }
1410 *pol_speca = pol_speca_loc ;
1411 *pol_specb = pol_specb_loc ;
1412 *ext_plista = ext_plista_loc ;
1413 *ext_plistb = ext_plistb_loc ;
1414 return 0 ;
1415}
1416
1417/*----------------------------------------------------------------------------*/
1445/*----------------------------------------------------------------------------*/
1446static int cr2res_obs_pol_reduce_one(
1447 const cpl_frameset * rawframes,
1448 const cpl_frameset * raw_background_frames,
1449 const cpl_frame * trace_wave_frame,
1450 const cpl_frame * detlin_frame,
1451 const cpl_frame * master_dark_frame,
1452 const cpl_frame * master_flat_frame,
1453 const cpl_frame * bpm_frame,
1454 const cpl_frame * blaze_frame,
1455 int subtract_nolight_rows,
1456 int subtract_interorder_column,
1457 int cosmics,
1458 int error_method,
1459 int extract_oversample,
1460 double extract_pclip,
1461 int extract_swath_width,
1462 int extract_height,
1463 double extract_smooth_slit,
1464 double extract_smooth_spec,
1465 int save_group,
1466 int reduce_det,
1467 hdrl_image *** in_calib,
1468 cpl_table *** tw,
1469 cpl_table *** extract1D,
1470 cpl_table ** pol_spec,
1471 cpl_propertylist ** ext_plist)
1472{
1473 cpl_vector * dits ;
1474 cpl_vector * ndits ;
1475 cpl_vector * ndsamples ;
1476 cr2res_decker * decker_positions ;
1477 hdrl_imagelist * in ;
1478 hdrl_imagelist * in_calib_loc ;
1479 hdrl_image * mflat ;
1480 cpl_image * flat_relerr2 ;
1481 cpl_image * bkg_var ;
1482 cpl_table ** flat_relvar ;
1483 cpl_table ** bkg_relvar ;
1484 cpl_table * corr ;
1485 int do_backgr ;
1486 cpl_image * contrib;
1487 hdrl_image * backgr;
1488 cpl_table * blaze_table ;
1489 cpl_table * trace_wave ;
1490 //cpl_table * trace_wave_corrected ;
1491 cpl_table ** trace_wave_extract ;
1492 const char * fname ;
1493 char * decker_name ;
1494 cpl_table * slit_func ;
1495 hdrl_image * model_master ;
1496 hdrl_image ** input_images ;
1497 cpl_table ** pol_spec_one_group ;
1498 cpl_table ** extract_1d ;
1499 char * colname ;
1500 const double * pcol_data ;
1501 double * pvec_data ;
1502 cpl_vector ** intens ;
1503 cpl_vector ** wl ;
1504 cpl_vector ** errors ;
1505 cpl_vector ** demod_wl ;
1506 cpl_bivector ** demod_stokes ;
1507 cpl_bivector ** demod_null ;
1508 cpl_bivector ** demod_intens ;
1509 int * orders ;
1510 cpl_table * pol_spec_merged ;
1511 cpl_propertylist * plist ;
1512 cpl_propertylist * ext_plist_loc ;
1513 const char * first_fname ;
1514 cpl_size nframes, nspec_group, spec_size ;
1515 int ngroups, i, j, k, l, o, norders, frame_idx ;
1516 double gain, error_factor, blaze_norm ;
1517 int * order_idx_values ;
1518 int nb_order_idx_values, order_zp;
1519
1520 /* TODO, make parameters */
1521 int extract_niter = 10;
1522 double extract_kappa = 10;
1523
1524 /* Check Inputs */
1525 if (in_calib == NULL || tw == NULL || extract1D == NULL ||
1526 pol_spec == NULL || ext_plist == NULL || rawframes == NULL ||
1527 trace_wave_frame == NULL)
1528 return -1 ;
1529
1530 /* Initialise */
1531 *in_calib = NULL ;
1532 *tw = NULL ;
1533 *extract1D = NULL ;
1534 *pol_spec = NULL ;
1535 *ext_plist = NULL ;
1536 first_fname = cpl_frame_get_filename(
1537 cpl_frameset_get_position_const(rawframes, 0)) ;
1538
1539 /* Get the Gain */
1540 if (reduce_det == 1) gain = CR2RES_GAIN_CHIP1 ;
1541 else if (reduce_det == 2) gain = CR2RES_GAIN_CHIP2 ;
1542 else if (reduce_det == 3) gain = CR2RES_GAIN_CHIP3 ;
1543 else {
1544 cpl_msg_error(__func__, "Failed to get the Gain value") ;
1545 return -1 ;
1546 }
1547
1548 /* Get the order zeropoint */
1549 if ((plist = cpl_propertylist_load(cpl_frame_get_filename(trace_wave_frame),
1550 0)) == NULL) {
1551 cpl_msg_error(__func__, "Cannot read the ORDER_ZP from the input TW") ;
1552 return -1 ;
1553 }
1554 order_zp = cr2res_pfits_get_order_zp(plist) ;
1555 cpl_propertylist_delete(plist) ;
1556 if (cpl_error_get_code()) {
1557 cpl_msg_error(__func__, "Missing ORDER_ZP in the header - Skip") ;
1558 cpl_error_reset() ;
1559 /* Negative Zerop to log the fact that it is missing */
1560 order_zp = -100 ;
1561 }
1562
1563 /* Check number of frames */
1564 nframes = cpl_frameset_get_size(rawframes) ;
1565 if (nframes == 0 || nframes % CR2RES_POLARIMETRY_GROUP_SIZE) {
1566 cpl_msg_error(__func__,
1567 "Input number of frames is %"CPL_SIZE_FORMAT" and should be multiple of %d",
1568 nframes, CR2RES_POLARIMETRY_GROUP_SIZE) ;
1569 return -1 ;
1570 }
1571
1572 /* Load the DITs if necessary */
1573 if (master_dark_frame != NULL) dits = cr2res_io_read_dits(rawframes) ;
1574 else dits = NULL ;
1575 if (cpl_msg_get_level() == CPL_MSG_DEBUG && dits != NULL)
1576 cpl_vector_dump(dits, stdout) ;
1577
1578 /* Load the NDITs */
1579 ndits = cr2res_io_read_ndits(rawframes) ;
1580
1581 /* Set the error factor. Note that we don't multiply by 4 because */
1582 /* each the beams get extracted from each frame, then errors propagated on.*/
1583 if (error_method == CR2RES_EXTRACT_ERROR_HORNE)
1584 error_factor = -1.0;
1585 else
1586 error_factor = gain * cpl_vector_get(ndits, 0);
1587
1588 for (i=0; i<cpl_vector_get_size(ndits); i++){
1589 if (cpl_vector_get(ndits,i) != cpl_vector_get(ndits, 0))
1590 cpl_msg_warning(__func__, "Raw frames have different NDIT! "
1591 "Error spectrum will likely be scaled incorrectly.");
1592 }
1593
1594 /* Get the decker positions */
1595 if ((decker_positions = cr2res_io_read_decker_positions(rawframes))==NULL) {
1596 cpl_msg_error(__func__, "Cannot read the Decker positions") ;
1597 if (dits != NULL) cpl_vector_delete(dits) ;
1598 return -1 ;
1599 }
1600
1601 /* Load image list */
1602 cpl_msg_info(__func__, "Load the images") ;
1603 if ((in = cr2res_io_load_image_list_from_set(rawframes,
1604 reduce_det)) == NULL) {
1605 cpl_msg_error(__func__, "Cannot load images") ;
1606 if (dits != NULL) cpl_vector_delete(dits) ;
1607 if (ndits != NULL) cpl_vector_delete(ndits) ;
1608 cpl_free(decker_positions) ;
1609 return -1 ;
1610 }
1611 if (hdrl_imagelist_get_size(in) != cpl_frameset_get_size(rawframes)) {
1612 cpl_msg_error(__func__, "Inconsistent number of loaded images") ;
1613 if (dits != NULL) cpl_vector_delete(dits) ;
1614 if (ndits != NULL) cpl_vector_delete(ndits) ;
1615 cpl_free(decker_positions) ;
1616 hdrl_imagelist_delete(in) ;
1617 return -1 ;
1618 }
1619
1620 /* Calibrate the images */
1621 cpl_msg_info(__func__, "Apply the calibrations") ;
1622 cpl_msg_indent_more() ;
1623 ndsamples = cr2res_io_read_ndsamples(rawframes) ;
1624 in_calib_loc = cr2res_calib_imagelist(in, reduce_det, 0,
1625 subtract_nolight_rows, subtract_interorder_column, cosmics,
1626 master_flat_frame, master_dark_frame, bpm_frame, detlin_frame,
1627 dits, ndits, ndsamples) ;
1628 cpl_vector_delete(ndsamples) ;
1629 if (in_calib_loc == NULL) {
1630 cpl_msg_error(__func__, "Failed to apply the calibrations") ;
1631 cpl_msg_indent_less() ;
1632 if (dits != NULL) cpl_vector_delete(dits) ;
1633 if (ndits != NULL) cpl_vector_delete(ndits) ;
1634 cpl_free(decker_positions) ;
1635 hdrl_imagelist_delete(in) ;
1636 return -1 ;
1637 }
1638 cpl_msg_indent_less() ;
1639 hdrl_imagelist_delete(in) ;
1640 if (ndits != NULL) cpl_vector_delete(ndits) ;
1641 ndits = NULL ;
1642
1643 /* Flat errors are common to all frames and cancel in the Stokes */
1644 /* ratios: remove them here and add them to the intensity only */
1645 do_backgr = cpl_frameset_get_size(raw_background_frames) > 1 ;
1646 mflat = NULL ;
1647 flat_relerr2 = NULL ;
1648 bkg_var = NULL ;
1649 corr = NULL ;
1650 if (master_flat_frame != NULL) {
1651 if ((mflat = cr2res_io_load_MASTER_FLAT(
1652 cpl_frame_get_filename(master_flat_frame),
1653 reduce_det)) == NULL) {
1654 cpl_msg_error(__func__, "Cannot load the master flat") ;
1655 if (dits != NULL) cpl_vector_delete(dits) ;
1656 cpl_free(decker_positions) ;
1657 hdrl_imagelist_delete(in_calib_loc) ;
1658 return -1 ;
1659 }
1660 flat_relerr2 = cr2res_calib_flat_relerr2(mflat) ;
1661 for (i=0 ; i<hdrl_imagelist_get_size(in_calib_loc) ; i++)
1662 cr2res_calib_flat_var(hdrl_imagelist_get(in_calib_loc, i),
1663 flat_relerr2, -1.0) ;
1664 }
1665 /* The dark cancels in the background subtraction */
1666 if (do_backgr && dits != NULL &&
1667 cr2res_calib_dark_var(in_calib_loc, master_dark_frame,
1668 reduce_det, mflat, dits, -1.0)) {
1669 cpl_msg_error(__func__, "Cannot remove the master dark errors") ;
1670 cpl_vector_delete(dits) ;
1671 cpl_free(decker_positions) ;
1672 hdrl_imagelist_delete(in_calib_loc) ;
1673 if (mflat != NULL) hdrl_image_delete(mflat) ;
1674 if (flat_relerr2 != NULL) cpl_image_delete(flat_relerr2) ;
1675 return -1 ;
1676 }
1677 if (dits != NULL) {
1678 cpl_vector_delete(dits) ;
1679 dits = NULL ;
1680 }
1681
1682 /* Apply Background Correction using the other nodding position */
1683 if (do_backgr) {
1684
1685 hdrl_imagelist * in_backgr ;
1686 cpl_msg_info(__func__, "Apply the background Correction") ;
1687 cpl_msg_indent_more() ;
1688
1689 /* Load image list for BACKGROUND frames */
1690 cpl_msg_info(__func__, "Load the background images") ;
1691 if ((in = cr2res_io_load_image_list_from_set(raw_background_frames,
1692 reduce_det)) == NULL) {
1693 cpl_msg_error(__func__, "Cannot load background images") ;
1694 cpl_msg_indent_less() ;
1695 cpl_free(decker_positions) ;
1696 hdrl_imagelist_delete(in_calib_loc) ;
1697 if (mflat != NULL) hdrl_image_delete(mflat) ;
1698 if (flat_relerr2 != NULL) cpl_image_delete(flat_relerr2) ;
1699 return -1 ;
1700 }
1701 /* Load the DITs & NDITs */
1702 if (master_dark_frame != NULL)
1703 dits = cr2res_io_read_dits(raw_background_frames) ;
1704 else
1705 dits = NULL ;
1706 ndits = cr2res_io_read_ndits(raw_background_frames) ;
1707
1708 /* Calibrate the background same as science images */
1709 cpl_msg_info(__func__, "Apply the calibrations to background") ;
1710 cpl_msg_indent_more() ;
1711 ndsamples = cr2res_io_read_ndsamples(raw_background_frames) ;
1712 in_backgr = cr2res_calib_imagelist(in, reduce_det, 0,
1713 subtract_nolight_rows, 1, cosmics, master_flat_frame,
1714 master_dark_frame, bpm_frame, detlin_frame,
1715 dits, ndits, ndsamples) ;
1716 cpl_vector_delete(ndsamples) ;
1717 if (in_backgr == NULL) {
1718 cpl_msg_error(__func__,
1719 "Failed to apply the calibrations to background") ;
1720 cpl_msg_indent_less() ;
1721 cpl_msg_indent_less() ;
1722 if (dits != NULL) cpl_vector_delete(dits) ;
1723 if (ndits != NULL) cpl_vector_delete(ndits) ;
1724 cpl_free(decker_positions) ;
1725 hdrl_imagelist_delete(in) ;
1726 hdrl_imagelist_delete(in_calib_loc) ;
1727 if (mflat != NULL) hdrl_image_delete(mflat) ;
1728 if (flat_relerr2 != NULL) cpl_image_delete(flat_relerr2) ;
1729 return -1 ;
1730 }
1731 cpl_msg_indent_less() ;
1732 hdrl_imagelist_delete(in) ;
1733
1734 /* Same common flat and dark error removal as the science frames */
1735 if (flat_relerr2 != NULL)
1736 for (i=0 ; i<hdrl_imagelist_get_size(in_backgr) ; i++)
1737 cr2res_calib_flat_var(hdrl_imagelist_get(in_backgr, i),
1738 flat_relerr2, -1.0) ;
1739 if (dits != NULL && cr2res_calib_dark_var(in_backgr,
1740 master_dark_frame, reduce_det, mflat, dits, -1.0)) {
1741 cpl_msg_error(__func__, "Cannot remove the master dark errors") ;
1742 cpl_msg_indent_less() ;
1743 cpl_vector_delete(dits) ;
1744 if (ndits != NULL) cpl_vector_delete(ndits) ;
1745 cpl_free(decker_positions) ;
1746 hdrl_imagelist_delete(in_calib_loc) ;
1747 hdrl_imagelist_delete(in_backgr) ;
1748 if (mflat != NULL) hdrl_image_delete(mflat) ;
1749 if (flat_relerr2 != NULL) cpl_image_delete(flat_relerr2) ;
1750 return -1 ;
1751 }
1752
1753 /* Collapse background images */
1754 if (hdrl_imagelist_collapse_mean(in_backgr, &backgr, &contrib) != \
1755 CPL_ERROR_NONE) {
1756 cpl_msg_error(__func__, "Failed to collapse background") ;
1757 cpl_msg_indent_less() ;
1758 if (dits != NULL) cpl_vector_delete(dits) ;
1759 if (ndits != NULL) cpl_vector_delete(ndits) ;
1760 cpl_free(decker_positions) ;
1761 hdrl_imagelist_delete(in_calib_loc) ;
1762 hdrl_imagelist_delete(in_backgr) ;
1763 if (mflat != NULL) hdrl_image_delete(mflat) ;
1764 if (flat_relerr2 != NULL) cpl_image_delete(flat_relerr2) ;
1765 return -1;
1766 }
1767 if (contrib != NULL) cpl_image_delete(contrib) ;
1768 if (in_backgr != NULL) hdrl_imagelist_delete(in_backgr) ;
1769
1770 /* Subtract background from calibrated images */
1771 if (hdrl_imagelist_sub_image(in_calib_loc, backgr) != CPL_ERROR_NONE) {
1772 cpl_msg_error(__func__, "Failed to subtract background") ;
1773 cpl_msg_indent_less() ;
1774 if (dits != NULL) cpl_vector_delete(dits) ;
1775 if (ndits != NULL) cpl_vector_delete(ndits) ;
1776 cpl_free(decker_positions) ;
1777 hdrl_imagelist_delete(in_calib_loc);
1778 hdrl_image_delete(backgr);
1779 if (mflat != NULL) hdrl_image_delete(mflat) ;
1780 if (flat_relerr2 != NULL) cpl_image_delete(flat_relerr2) ;
1781 return -1;
1782 }
1783 /* The background is common to all frames: keep its variance */
1784 if (error_method == CR2RES_EXTRACT_ERROR_HORNE)
1785 bkg_var = cpl_image_power_create(
1786 hdrl_image_get_error_const(backgr), 2.0) ;
1787 hdrl_image_delete(backgr);
1788 cpl_msg_indent_less() ;
1789 } else {
1790 cpl_msg_warning(__func__, "No background subtraction");
1791 }
1792 if (dits != NULL) cpl_vector_delete(dits) ;
1793 if (ndits != NULL) cpl_vector_delete(ndits) ;
1794 if (mflat != NULL) hdrl_image_delete(mflat) ;
1795
1796 /* Load the trace wave */
1797 cpl_msg_info(__func__, "Load the TRACE WAVE") ;
1798 if ((trace_wave = cr2res_io_load_TRACE_WAVE(cpl_frame_get_filename(
1799 trace_wave_frame), reduce_det)) == NULL) {
1800 cpl_msg_error(__func__, "Failed to Load the traces file") ;
1801 cpl_free(decker_positions) ;
1802 hdrl_imagelist_delete(in_calib_loc) ;
1803 if (flat_relerr2 != NULL) cpl_image_delete(flat_relerr2) ;
1804 if (bkg_var != NULL) cpl_image_delete(bkg_var) ;
1805 return -1 ;
1806 }
1807
1808 /* Correct trace_wave with some provided raw flats */
1809 cpl_msg_warning(__func__, "TRACE ADJUST NOT YET IMPLEMENTED") ;
1810 /*if (raw_flat_frames != NULL) {
1811 cpl_msg_info(__func__, "Try to correct the reproducibility error") ;
1812 cpl_msg_indent_more() ;
1813 trace_wave_corrected = cr2res_trace_adjust(trace_wave, raw_flat_frames,
1814 reduce_det) ;
1815 if (trace_wave_corrected != NULL) {
1816 cpl_table_delete(trace_wave) ;
1817 trace_wave = trace_wave_corrected ;
1818 trace_wave_corrected = NULL ;
1819 }
1820 cpl_msg_indent_less() ;
1821 }*/
1822
1823 /* Load Blaze */
1824 blaze_table = NULL ;
1825 blaze_norm = 0 ;
1826 if (blaze_frame != NULL) {
1827 cpl_msg_info(__func__, "Load the BLAZE") ;
1828 if ((blaze_table = cr2res_io_load_EXTRACT_1D(cpl_frame_get_filename(
1829 blaze_frame), reduce_det)) == NULL) {
1830 cpl_msg_error(__func__, "Failed to Load the Blaze file") ;
1831 cpl_free(decker_positions) ;
1832 hdrl_imagelist_delete(in_calib_loc) ;
1833 cpl_table_delete(trace_wave) ;
1834 if (flat_relerr2 != NULL) cpl_image_delete(flat_relerr2) ;
1835 if (bkg_var != NULL) cpl_image_delete(bkg_var) ;
1836 return -1 ;
1837 }
1838 cpl_propertylist * blaze_plist ;
1839 blaze_plist = cpl_propertylist_load_regexp(cpl_frame_get_filename(
1840 blaze_frame),0,CR2RES_HEADER_QC_BLAZE_NORM,0);
1841 if(cpl_propertylist_get_size(blaze_plist)>0){
1842 blaze_norm = cpl_propertylist_get_double(blaze_plist, CR2RES_HEADER_QC_BLAZE_NORM);
1843 }
1844 else {
1845 blaze_norm = -1;
1846 cpl_msg_warning(__func__, "QC BLAZE NORM value not found, reverting to per trace normalization") ;
1847 }
1848 if(blaze_plist!=NULL){
1849 cpl_propertylist_delete(blaze_plist);
1850 }
1851 }
1852
1853 /* Compute the number of groups */
1854 ngroups = nframes/CR2RES_POLARIMETRY_GROUP_SIZE ;
1855 nspec_group = 2*CR2RES_POLARIMETRY_GROUP_SIZE ;
1856
1857 /* Allocate pol_spec_group containers */
1858 pol_spec_one_group = cpl_malloc(ngroups * sizeof(cpl_table*)) ;
1859 for (i = 0; i < ngroups; i++) pol_spec_one_group[i] = NULL;
1860
1861 /* Loop on the groups */
1862 for (i = 0; i < ngroups; i++) {
1863 int *pol_sorting;
1864 cpl_msg_info(__func__, "Process %d-group number %d/%d",
1865 CR2RES_POLARIMETRY_GROUP_SIZE, i+1, ngroups) ;
1866 cpl_msg_indent_more() ;
1867
1868 /* Compute the proper order of the frames group */
1869 if ((pol_sorting = cr2res_pol_sort_frames(
1870 cpl_frameset_get_position_const(rawframes,
1871 i * CR2RES_POLARIMETRY_GROUP_SIZE),
1872 cpl_frameset_get_position_const(rawframes,
1873 i * CR2RES_POLARIMETRY_GROUP_SIZE + 1),
1874 cpl_frameset_get_position_const(rawframes,
1875 i * CR2RES_POLARIMETRY_GROUP_SIZE + 2),
1876 cpl_frameset_get_position_const(rawframes,
1877 i * CR2RES_POLARIMETRY_GROUP_SIZE + 3))) == NULL) {
1878 cpl_msg_warning(__func__, "Failed to sort the files") ;
1879 continue ;
1880 }
1881 if (cpl_msg_get_level() == CPL_MSG_DEBUG) {
1882 for (j=0 ; j<CR2RES_POLARIMETRY_GROUP_SIZE ; j++) {
1883 if (pol_sorting[j] != j)
1884 cpl_msg_warning(__func__,
1885 "Frame #%d moved to position #%d",
1886 j+1, pol_sorting[j]+1) ;
1887 }
1888 }
1889
1890 /* Container for the current group calibrated images 1 2 3 4 */
1891 input_images = cpl_malloc(CR2RES_POLARIMETRY_GROUP_SIZE *
1892 sizeof(hdrl_image*));
1893 for (j = 0; j < CR2RES_POLARIMETRY_GROUP_SIZE; j++)
1894 input_images[j] = NULL ;
1895
1896 /* Container for extracted tables: 1u, 1d, 2u, 2d, 3u, 3d, 4u, 4d */
1897 extract_1d = cpl_malloc(nspec_group * sizeof(cpl_table*));
1898 trace_wave_extract = cpl_malloc(nspec_group * sizeof(cpl_table*));
1899 flat_relvar = cpl_malloc(nspec_group * sizeof(cpl_table*));
1900 bkg_relvar = cpl_malloc(nspec_group * sizeof(cpl_table*));
1901 for (j = 0; j < nspec_group; j++) {
1902 extract_1d[j] = NULL;
1903 trace_wave_extract[j] = NULL;
1904 flat_relvar[j] = NULL;
1905 bkg_relvar[j] = NULL;
1906 }
1907
1908 /* Loop on the frames */
1909 for (j=0 ; j<CR2RES_POLARIMETRY_GROUP_SIZE ; j++) {
1910 frame_idx = i * CR2RES_POLARIMETRY_GROUP_SIZE + pol_sorting[j] ;
1911 fname = cpl_frame_get_filename(cpl_frameset_get_position_const(
1912 rawframes,frame_idx)) ;
1913
1914 /* Store the input image */
1915 input_images[j] = hdrl_image_duplicate(
1916 hdrl_imagelist_get_const(in_calib_loc, frame_idx)) ;
1917
1918 /* Extract up */
1919 decker_name = cr2res_decker_print_position(
1920 decker_positions[frame_idx]) ;
1921 cpl_msg_info(__func__,
1922 "Extract Up Spectrum from %s (Det %d / Decker %s)",
1923 fname, reduce_det, decker_name) ;
1924 cpl_free(decker_name) ;
1925 cpl_msg_indent_more();
1926
1927 trace_wave_extract[2 * j] = cr2res_pol_get_beam_trace(trace_wave, decker_positions[frame_idx], 1);
1928
1929 /* Execute the extraction */
1930 cpl_msg_info(__func__, "Spectra Extraction") ;
1931 if (cr2res_extract_traces(input_images[j], trace_wave_extract[2*j],
1932 NULL, blaze_table, blaze_norm, -1, -1, CR2RES_EXTR_OPT_CURV,
1933 extract_height, extract_swath_width, extract_oversample, extract_pclip,
1934 extract_smooth_slit, extract_smooth_spec,
1935 extract_niter, extract_kappa, error_factor, 0, 0, 0,
1936 &(extract_1d[2*j]), &slit_func, &model_master) == -1) {
1937 cpl_msg_error(__func__, "Failed Extraction") ;
1938 extract_1d[2*j] = NULL ;
1939 cpl_table_delete(trace_wave_extract[2*j]) ;
1940 trace_wave_extract[2*j] = NULL ;
1941 } else {
1942 if (flat_relerr2 != NULL)
1943 flat_relvar[2*j] = cr2res_obs_pol_relvar(
1944 model_master, flat_relerr2, 1,
1945 trace_wave_extract[2*j], extract_height) ;
1946 if (bkg_var != NULL)
1947 bkg_relvar[2*j] = cr2res_obs_pol_relvar(
1948 model_master, bkg_var, 0,
1949 trace_wave_extract[2*j], extract_height) ;
1950 cpl_table_delete(slit_func) ;
1951 hdrl_image_delete(model_master) ;
1952 }
1953 cpl_msg_indent_less() ;
1954
1955 /* Extract Down */
1956 decker_name = cr2res_decker_print_position(
1957 decker_positions[frame_idx]) ;
1958 cpl_msg_info(__func__,
1959 "Extract Down Spectrum from %s (Det %d / Decker %s)",
1960 fname, reduce_det, decker_name) ;
1961 cpl_free(decker_name) ;
1962 cpl_msg_indent_more() ;
1963
1964 trace_wave_extract[2*j+1] = cr2res_pol_get_beam_trace(trace_wave, decker_positions[frame_idx], 2);
1965
1966 /* Execute the extraction */
1967 cpl_msg_info(__func__, "Spectra Extraction") ;
1968 if (cr2res_extract_traces(input_images[j],
1969 trace_wave_extract[2*j+1], NULL, blaze_table, blaze_norm, -1, -1,
1970 CR2RES_EXTR_OPT_CURV, extract_height,
1971 extract_swath_width, extract_oversample, extract_pclip,
1972 extract_smooth_slit, extract_smooth_spec,
1973 extract_niter, extract_kappa, error_factor, 0, 0, 0,
1974 &(extract_1d[2*j+1]), &slit_func, &model_master)== -1) {
1975 cpl_msg_error(__func__, "Failed Extraction") ;
1976 extract_1d[2*j+1] = NULL ;
1977 cpl_table_delete(trace_wave_extract[2*j+1]) ;
1978 trace_wave_extract[2*j+1] = NULL ;
1979 } else {
1980 if (flat_relerr2 != NULL)
1981 flat_relvar[2*j+1] = cr2res_obs_pol_relvar(
1982 model_master, flat_relerr2, 1,
1983 trace_wave_extract[2*j+1], extract_height) ;
1984 if (bkg_var != NULL)
1985 bkg_relvar[2*j+1] = cr2res_obs_pol_relvar(
1986 model_master, bkg_var, 0,
1987 trace_wave_extract[2*j+1], extract_height) ;
1988 cpl_table_delete(slit_func) ;
1989 hdrl_image_delete(model_master) ;
1990 }
1991 cpl_msg_indent_less() ;
1992 }
1993 cpl_free(pol_sorting) ;
1994
1995 /* How many orders */
1996 if ((orders = cr2res_obs_pol_get_order_numbers(
1997 (const cpl_table **)extract_1d,
1998 nspec_group,
1999 &norders)) == NULL) {
2000 cpl_msg_warning(__func__, "No Order found in the extracted tables");
2001 for (j = 0; j < CR2RES_POLARIMETRY_GROUP_SIZE; j++) {
2002 hdrl_image_delete(input_images[j]) ;
2003 input_images[j] = NULL ;
2004 }
2005 cpl_free(input_images) ;
2006 for (j=0 ; j<nspec_group ; j++) {
2007 if (extract_1d[j] != NULL)
2008 cpl_table_delete(extract_1d[j]) ;
2009 if (trace_wave_extract[j] != NULL)
2010 cpl_table_delete(trace_wave_extract[j]) ;
2011 if (flat_relvar[j] != NULL)
2012 cpl_table_delete(flat_relvar[j]) ;
2013 if (bkg_relvar[j] != NULL)
2014 cpl_table_delete(bkg_relvar[j]) ;
2015 }
2016 cpl_free(extract_1d) ;
2017 cpl_free(trace_wave_extract) ;
2018 cpl_free(flat_relvar) ;
2019 cpl_free(bkg_relvar) ;
2020 continue ;
2021 }
2022 cpl_msg_debug(__func__, "%d different orders found", norders) ;
2023
2024 /* Allocate data containers */
2025 demod_wl = cpl_malloc(norders * sizeof(cpl_vector*)) ;
2026 demod_stokes = cpl_malloc(norders * sizeof(cpl_bivector*)) ;
2027 demod_null = cpl_malloc(norders * sizeof(cpl_bivector*)) ;
2028 demod_intens = cpl_malloc(norders * sizeof(cpl_bivector*)) ;
2029 for (o = 0; o < norders; o++){
2030 demod_wl[o] = NULL;
2031 demod_stokes[o] = NULL;
2032 demod_null[o] = NULL;
2033 demod_intens[o] = NULL;
2034 }
2035
2036 intens = cpl_malloc(nspec_group * sizeof(cpl_vector*));
2037 wl = cpl_malloc(nspec_group * sizeof(cpl_vector*));
2038 errors = cpl_malloc(nspec_group * sizeof(cpl_vector*));
2039 for (k = 0; k < nspec_group; k++){
2040 intens[k] = NULL;
2041 wl[k] = NULL;
2042 errors[k] = NULL;
2043 }
2044
2045 /* Loop on the orders */
2046 for (o=0 ; o<norders ; o++) {
2047 cpl_msg_info(__func__, "Compute Polarimetry for order %d",
2048 orders[o]) ;
2049 /* Get the inputs for the demod functions calls */
2050 for (k=0 ; k<nspec_group ; k++) {
2051 spec_size = cpl_table_get_nrow(extract_1d[k]) ;
2052 /* Get the SPEC for this order/trace 1 */
2053 colname = cr2res_dfs_SPEC_colname(orders[o], 1) ;
2054 if (cpl_table_has_valid(extract_1d[k], colname)){
2055 if ((pcol_data = cpl_table_get_data_double_const(
2056 extract_1d[k], colname)) == NULL) {
2057 cpl_error_reset() ;
2058 } else {
2059 intens[k] = cpl_vector_new(spec_size) ;
2060 pvec_data = cpl_vector_get_data(intens[k]) ;
2061 for (l=0 ; l<spec_size ; l++)
2062 pvec_data[l] = pcol_data[l] ;
2063 }
2064 }
2065 cpl_free(colname) ;
2066
2067 /* Get the WAVELENGTH for this order/trace 1 */
2068 colname = cr2res_dfs_WAVELENGTH_colname(orders[o], 1) ;
2069 if (cpl_table_has_valid(extract_1d[k], colname)){
2070 if ((pcol_data = cpl_table_get_data_double_const(
2071 extract_1d[k], colname)) == NULL) {
2072 cpl_error_reset() ;
2073 } else {
2074 wl[k] = cpl_vector_new(spec_size) ;
2075 pvec_data = cpl_vector_get_data(wl[k]) ;
2076 for (l=0 ; l<spec_size ; l++)
2077 pvec_data[l] = pcol_data[l] ;
2078 }
2079 }
2080 cpl_free(colname) ;
2081
2082 /* Get the ERROR for this order/trace 1 */
2083 colname = cr2res_dfs_SPEC_ERR_colname(orders[o], 1) ;
2084 if (cpl_table_has_valid(extract_1d[k], colname)){
2085 if ((pcol_data = cpl_table_get_data_double_const(
2086 extract_1d[k], colname)) == NULL) {
2087 cpl_error_reset() ;
2088 } else {
2089 errors[k] = cpl_vector_new(spec_size) ;
2090 pvec_data = cpl_vector_get_data(errors[k]) ;
2091 for (l=0 ; l<spec_size ; l++)
2092 pvec_data[l] = pcol_data[l] ;
2093 }
2094 }
2095 cpl_free(colname) ;
2096 }
2097
2098 /* Keep the 1st wl of the 8 as reference for the output file*/
2099 if (wl[0] != NULL) demod_wl[o] = cpl_vector_duplicate(wl[0]) ;
2100 else demod_wl[o] = NULL ;
2101
2102 /* Call Library Demodulation functions */
2103 demod_stokes[o] = cr2res_pol_demod_stokes(intens, wl, errors,
2104 nspec_group) ;
2105 demod_null[o] = cr2res_pol_demod_null(intens, wl, errors,
2106 nspec_group) ;
2107 demod_intens[o] = cr2res_pol_demod_intens(intens, wl, errors,
2108 nspec_group) ;
2109 if ((flat_relerr2 != NULL || bkg_var != NULL) &&
2110 demod_intens[o] != NULL)
2111 cr2res_obs_pol_corr_accumulate(&corr, intens, flat_relvar,
2112 bkg_relvar, orders[o], nspec_group) ;
2113
2114 /* Free */
2115 for (k=0 ; k<nspec_group ; k++) {
2116 if (intens[k] != NULL) {
2117 cpl_vector_delete(intens[k]) ;
2118 intens[k] = NULL;
2119 }
2120 if (wl[k] != NULL) {
2121 cpl_vector_delete(wl[k]) ;
2122 wl[k] = NULL;
2123 }
2124 if (errors[k] != NULL) {
2125 cpl_vector_delete(errors[k]) ;
2126 errors[k] = NULL;
2127 }
2128 }
2129 }
2130 cpl_free(intens) ;
2131 cpl_free(wl) ;
2132 cpl_free(errors) ;
2133
2134 /* Save the extra products for this group */
2135 if (i+1==save_group) {
2136 *in_calib = input_images ;
2137 *tw = trace_wave_extract ;
2138 *extract1D = extract_1d ;
2139 } else {
2140 /* Free the input images */
2141 for (j = 0; j < CR2RES_POLARIMETRY_GROUP_SIZE; j++)
2142 if (input_images[j] != NULL)
2143 hdrl_image_delete(input_images[j]) ;
2144 cpl_free(input_images) ;
2145 /* Free the extracted tables / the TW */
2146 for (j=0 ; j<nspec_group ; j++) {
2147 if (trace_wave_extract[j] != NULL)
2148 cpl_table_delete(trace_wave_extract[j]) ;
2149 if (extract_1d[j] != NULL) cpl_table_delete(extract_1d[j]) ;
2150 }
2151 cpl_free(trace_wave_extract) ;
2152 cpl_free(extract_1d) ;
2153 }
2154 for (j=0 ; j<nspec_group ; j++) {
2155 if (flat_relvar[j] != NULL) cpl_table_delete(flat_relvar[j]) ;
2156 if (bkg_relvar[j] != NULL) cpl_table_delete(bkg_relvar[j]) ;
2157 }
2158 cpl_free(flat_relvar) ;
2159 cpl_free(bkg_relvar) ;
2160
2161 /* Create the pol_spec table */
2162 cpl_msg_info(__func__, "Create the POL_SPEC table for this group");
2163 pol_spec_one_group[i] = cr2res_pol_POL_SPEC_create(orders, demod_wl,
2164 demod_stokes, demod_null, demod_intens, norders) ;
2165
2166 /* Deallocate */
2167 for (o=0 ; o<norders ; o++) {
2168 if (demod_wl[o] != NULL) cpl_vector_delete(demod_wl[o]) ;
2169 if (demod_stokes[o] != NULL) cpl_bivector_delete(demod_stokes[o]) ;
2170 if (demod_null[o] != NULL) cpl_bivector_delete(demod_null[o]) ;
2171 if (demod_intens[o] != NULL) cpl_bivector_delete(demod_intens[o]) ;
2172 }
2173 cpl_free(demod_wl) ;
2174 cpl_free(demod_stokes) ;
2175 cpl_free(demod_null) ;
2176 cpl_free(demod_intens) ;
2177 cpl_free(orders) ;
2178 cpl_msg_indent_less() ;
2179 }
2180 if (blaze_table != NULL) cpl_table_delete(blaze_table) ;
2181 if (flat_relerr2 != NULL) cpl_image_delete(flat_relerr2) ;
2182 if (bkg_var != NULL) cpl_image_delete(bkg_var) ;
2183 cpl_free(decker_positions) ;
2184 hdrl_imagelist_delete(in_calib_loc) ;
2185
2186 /* Merge the groups together */
2187 if (ngroups > 1) {
2188 cpl_msg_info(__func__, "Merge the %d groups POL_SPEC tables into one",
2189 ngroups);
2190 pol_spec_merged = cr2res_pol_spec_pol_merge(
2191 (const cpl_table **)pol_spec_one_group, ngroups) ;
2192 } else {
2193 if (pol_spec_one_group[0] == NULL) {
2194 pol_spec_merged = NULL ;
2195 } else {
2196 pol_spec_merged = cpl_table_duplicate(pol_spec_one_group[0]) ;
2197 }
2198 }
2199
2200 if (corr != NULL) {
2201 if (pol_spec_merged != NULL)
2202 cr2res_obs_pol_corr_add(pol_spec_merged, corr, ngroups) ;
2203 cpl_table_delete(corr) ;
2204 }
2205
2206 /* Deallocate */
2207 for (i=0 ; i<ngroups ; i++) cpl_table_delete(pol_spec_one_group[i]) ;
2208 cpl_free(pol_spec_one_group) ;
2209
2210 /* Check */
2211 if (pol_spec_merged == NULL) {
2212 cpl_msg_error(__func__, "Cannot create the POL_SPEC table");
2213 cpl_table_delete(trace_wave) ;
2214 for (j=0 ; j<nspec_group ; j++) {
2215 cpl_table **temp_extract1D = extract1D[j];
2216 if (temp_extract1D != NULL && *temp_extract1D != NULL)
2217 cpl_table_delete(*temp_extract1D);
2218 cpl_table **temp_tw = tw[j];
2219 if (temp_tw != NULL && *temp_tw != NULL)
2220 cpl_table_delete(*temp_tw) ;
2221 }
2222 cpl_free(*tw) ; *tw = NULL ;
2223 cpl_free(*extract1D) ; *extract1D = NULL ;
2224 return -1 ;
2225 }
2226
2227 /* Store the extension header for product saving */
2228 ext_plist_loc = cpl_propertylist_load(first_fname,
2229 cr2res_io_get_ext_idx(first_fname, reduce_det, 1)) ;
2230
2231 cpl_msg_info(__func__,"Order zero-point: %d", order_zp );
2232 /* Real Orders in QCs */
2233 if (order_zp > 0) {
2234 /* Get the order numbers from the TW rows */
2235 order_idx_values = cr2res_trace_get_order_idx_values(trace_wave,
2236 &nb_order_idx_values);
2237 /* Compute the Real Order numbers and store them in QCs */
2238 for (i = 0; i < nb_order_idx_values; i++) {
2239 char *key_name;
2240 char *spec_name;
2241 char *err_name;
2242 cpl_vector *spec_vec;
2243 cpl_vector *err_vec;
2244 double snr;
2245 int order_real, order_idx, order_idxp;
2246 order_idx = order_idx_values[i] ;
2247 order_idxp = cr2res_io_convert_order_idx_to_idxp(order_idx) ;
2248 order_real = cr2res_order_idx_to_real(order_idx, order_zp) ;
2249 key_name = cpl_sprintf(CR2RES_HEADER_QC_REAL_ORDER, order_idxp) ;
2250 cpl_propertylist_append_int(ext_plist_loc, key_name, order_real) ;
2251 cpl_free(key_name) ;
2252
2253 /* Compute SNR and store QC*/
2254 key_name = cpl_sprintf(CR2RES_HEADER_QC_SNR, order_idxp);
2255 spec_name = cr2res_dfs_POL_INTENS_colname(order_idxp);
2256 err_name = cr2res_dfs_POL_INTENS_ERROR_colname(order_idxp);
2257
2258 spec_vec = cpl_vector_wrap(cpl_table_get_nrow(pol_spec_merged),
2259 cpl_table_get_data_double(pol_spec_merged, spec_name));
2260 err_vec = cpl_vector_wrap(cpl_table_get_nrow(pol_spec_merged),
2261 cpl_table_get_data_double(pol_spec_merged, err_name));
2262 snr = cr2res_qc_compute_snr(spec_vec, err_vec);
2263 cpl_propertylist_append_double(ext_plist_loc, key_name, snr);
2264
2265 cpl_vector_unwrap(spec_vec);
2266 cpl_vector_unwrap(err_vec);
2267 cpl_free(key_name);
2268 cpl_free(spec_name) ;
2269 cpl_free(err_name) ;
2270 }
2271 cpl_free(order_idx_values) ;
2272 }
2273 cpl_table_delete(trace_wave) ;
2274
2275 /* Return */
2276 *pol_spec = pol_spec_merged ;
2277 *ext_plist = ext_plist_loc ;
2278 return 0 ;
2279}
2280
2281/*----------------------------------------------------------------------------*/
2292/*----------------------------------------------------------------------------*/
2293static int * cr2res_obs_pol_get_order_numbers(
2294 const cpl_table ** extracted,
2295 int next,
2296 int * norders)
2297{
2298 const char * col_name ;
2299 char * col_type ;
2300 int * tmp_orders_list ;
2301 int count_orders ;
2302 int * out_orders ;
2303 cpl_size j;
2304 int i, k, order, trace_nb, max_possible_spectra, new_order ;
2305
2306 /* Check entries */
2307 if (extracted == NULL || norders == NULL) return NULL ;
2308 for (i=0 ; i<next ; i++)
2309 if (extracted[i] == NULL) return NULL ;
2310
2311 /* Initialize */
2312 max_possible_spectra = 1000 ;
2313
2314 /* Allocate orders list */
2315 tmp_orders_list = cpl_malloc(max_possible_spectra * sizeof(int)) ;
2316
2317 /* Count the different orders */
2318 count_orders = 0 ;
2319 /* Loop over all columns */
2320 for (i = 0; i < next; i++) {
2321 cpl_size ncols;
2322 cpl_array *col_names;
2323 col_names = cpl_table_get_column_names(extracted[i]);
2324 ncols = cpl_table_get_ncol(extracted[i]) ;
2325 for (j=0 ; j<ncols ; j++) {
2326 col_name = cpl_array_get_string(col_names, j);
2327 col_type = cr2res_dfs_SPEC_colname_parse(col_name, &order,
2328 &trace_nb) ;
2329 if (col_type != NULL && !strcmp(col_type, CR2RES_COL_SPEC_SUFFIX)) {
2330 /* Is the current order a new one ? */
2331 new_order = 1 ;
2332 for (k=0 ; k<count_orders ; k++)
2333 if (tmp_orders_list[k] == order)
2334 new_order = 0 ;
2335
2336 /* Current order not yet stored */
2337 if (new_order) {
2338 tmp_orders_list[count_orders] = order ;
2339 count_orders ++ ;
2340 }
2341 }
2342 if (col_type != NULL) cpl_free(col_type) ;
2343 }
2344 cpl_array_delete(col_names) ;
2345 }
2346
2347 /* Allocate and fill output array */
2348 out_orders = cpl_malloc(count_orders * sizeof(int)) ;
2349 for (i=0 ; i<count_orders ; i++) out_orders[i] = tmp_orders_list[i] ;
2350
2351 /* Free and return */
2352 cpl_free(tmp_orders_list) ;
2353 *norders = count_orders ;
2354 return out_orders ;
2355}
2356
2357/*----------------------------------------------------------------------------*/
2372/*----------------------------------------------------------------------------*/
2373static cpl_table * cr2res_obs_pol_relvar(
2374 const hdrl_image * model,
2375 const cpl_image * var,
2376 int relative,
2377 const cpl_table * tw,
2378 int extract_height)
2379{
2380 const cpl_image * mod_ima ;
2381 const cpl_mask * bpm ;
2382 const cpl_binary * pbpm ;
2383 const double * pmod ;
2384 const double * pvar ;
2385 cpl_table * out ;
2386 cpl_size nx, ny, r, x, y ;
2387
2388 if (model == NULL || var == NULL || tw == NULL) return NULL ;
2389 mod_ima = hdrl_image_get_image_const(model) ;
2390 nx = cpl_image_get_size_x(mod_ima) ;
2391 ny = cpl_image_get_size_y(mod_ima) ;
2392 pmod = cpl_image_get_data_double_const(mod_ima) ;
2393 pvar = cpl_image_get_data_double_const(var) ;
2394 if (pmod == NULL || pvar == NULL ||
2395 cpl_image_get_size_x(var) != nx ||
2396 cpl_image_get_size_y(var) != ny) return NULL ;
2397 bpm = cpl_image_get_bpm_const(var) ;
2398 pbpm = bpm != NULL ? cpl_mask_get_data_const(bpm) : NULL ;
2399
2400 out = cpl_table_new(nx) ;
2401 for (r=0 ; r<cpl_table_get_nrow(tw) ; r++) {
2402 cpl_vector * ycen ;
2403 double * pcol ;
2404 char * colname ;
2405 int order, trace_nb, height ;
2406
2407 order = cpl_table_get(tw, CR2RES_COL_ORDER, r, NULL) ;
2408 trace_nb = cpl_table_get(tw, CR2RES_COL_TRACENB, r, NULL) ;
2409 height = extract_height > 0 ? extract_height :
2410 cr2res_trace_get_height(tw, order, trace_nb) ;
2411 colname = cr2res_dfs_SPEC_colname(order, trace_nb) ;
2412 if (height <= 0 || cpl_table_has_column(out, colname) ||
2413 (ycen = cr2res_trace_get_ycen(tw, order, trace_nb, nx))
2414 == NULL) {
2415 cpl_free(colname) ;
2416 cpl_error_reset() ;
2417 continue ;
2418 }
2419 cpl_table_new_column(out, colname, CPL_TYPE_DOUBLE) ;
2420 cpl_table_fill_column_window_double(out, colname, 0, nx, NAN) ;
2421 pcol = cpl_table_get_data_double(out, colname) ;
2422
2423 for (x=0 ; x<nx ; x++) {
2424 /* One extra row each side; the model is zero outside the trace */
2425 cpl_size ymin = (cpl_size)floor(cpl_vector_get(ycen, x)
2426 - height/2.0) - 1 ;
2427 cpl_size ymax = (cpl_size)ceil(cpl_vector_get(ycen, x)
2428 + height/2.0) + 1 ;
2429 double m2 = 0.0, num = 0.0 ;
2430 if (ymin < 1) ymin = 1 ;
2431 if (ymax > ny) ymax = ny ;
2432 for (y=ymin ; y<=ymax ; y++) {
2433 cpl_size idx = (y-1) * nx + x ;
2434 double m = pmod[idx] ;
2435 double v = pvar[idx] ;
2436 if (pbpm != NULL && pbpm[idx]) continue ;
2437 if (!isfinite(m) || !isfinite(v)) continue ;
2438 if (relative) v *= m * m ;
2439 m2 += m * m ;
2440 num += m * m * v ;
2441 }
2442 if (m2 > 0.0) pcol[x] = num / (m2 * m2) ;
2443 }
2444 cpl_vector_delete(ycen) ;
2445 cpl_free(colname) ;
2446 }
2447 return out ;
2448}
2449
2450/* Up/down flat sums, up/down background sums, background sum of squares */
2451static const char * cr2res_obs_pol_corr_sfx[5] = {"FU","FD","BU","BD","B2"} ;
2452
2453/* Relative variance at column x, or -1 if not available */
2454static double cr2res_obs_pol_relvar_at(
2455 const cpl_table * relvar,
2456 const char * colname,
2457 cpl_size x)
2458{
2459 double rv ;
2460 if (relvar == NULL || !cpl_table_has_column(relvar, colname) ||
2461 x >= cpl_table_get_nrow(relvar)) return -1.0 ;
2462 rv = cpl_table_get_double(relvar, colname, x, NULL) ;
2463 return isfinite(rv) ? rv : -1.0 ;
2464}
2465
2466/*----------------------------------------------------------------------------*/
2483/*----------------------------------------------------------------------------*/
2484static void cr2res_obs_pol_corr_accumulate(
2485 cpl_table ** corr,
2486 cpl_vector ** intens,
2487 cpl_table ** flat_relvar,
2488 cpl_table ** bkg_relvar,
2489 int order,
2490 int n)
2491{
2492 double * pc[5] ;
2493 char * base ;
2494 char * spec_col ;
2495 cpl_size size, x ;
2496 int c, k ;
2497
2498 if ((base = cr2res_dfs_POL_INTENS_ERROR_colname(order)) == NULL) return ;
2499 size = cpl_vector_get_size(intens[0]) ;
2500 if (*corr == NULL) *corr = cpl_table_new(size) ;
2501 if (cpl_table_get_nrow(*corr) < size) size = cpl_table_get_nrow(*corr) ;
2502 for (c=0 ; c<5 ; c++) {
2503 char * name = cpl_sprintf("%s_%s", base, cr2res_obs_pol_corr_sfx[c]) ;
2504 if (!cpl_table_has_column(*corr, name)) {
2505 cpl_table_new_column(*corr, name, CPL_TYPE_DOUBLE) ;
2506 cpl_table_fill_column_window_double(*corr, name, 0,
2507 cpl_table_get_nrow(*corr), 0.0) ;
2508 }
2509 pc[c] = cpl_table_get_data_double(*corr, name) ;
2510 cpl_free(name) ;
2511 }
2512 cpl_free(base) ;
2513
2514 spec_col = cr2res_dfs_SPEC_colname(order, 1) ;
2515 for (x=0 ; x<size ; x++) {
2516 for (k=0 ; k<n ; k++) {
2517 double rv, val ;
2518 if (x >= cpl_vector_get_size(intens[k])) continue ;
2519 val = cpl_vector_get(intens[k], x) ;
2520 if (!isfinite(val)) continue ;
2521 if (flat_relvar != NULL && (rv = cr2res_obs_pol_relvar_at(
2522 flat_relvar[k], spec_col, x)) >= 0.0)
2523 pc[k % 2][x] += val * sqrt(rv) ;
2524 if (bkg_relvar != NULL && (rv = cr2res_obs_pol_relvar_at(
2525 bkg_relvar[k], spec_col, x)) >= 0.0) {
2526 double s = fabs(val) * sqrt(rv) ;
2527 pc[2 + k % 2][x] += s ;
2528 pc[4][x] += s * s ;
2529 }
2530 }
2531 }
2532 cpl_free(spec_col) ;
2533}
2534
2535/*----------------------------------------------------------------------------*/
2544/*----------------------------------------------------------------------------*/
2545static void cr2res_obs_pol_corr_add(
2546 cpl_table * pol_spec,
2547 const cpl_table * corr,
2548 int ngroups)
2549{
2550 cpl_array * names ;
2551 cpl_size nrow, i, x ;
2552
2553 nrow = cpl_table_get_nrow(pol_spec) ;
2554 if (cpl_table_get_nrow(corr) < nrow) nrow = cpl_table_get_nrow(corr) ;
2555 names = cpl_table_get_column_names(pol_spec) ;
2556 for (i=0 ; i<cpl_array_get_size(names) ; i++) {
2557 const char * err_col = cpl_array_get_string(names, i) ;
2558 const double * pc[5] ;
2559 double * perr ;
2560 int c, found = 1 ;
2561
2562 for (c=0 ; c<5 ; c++) {
2563 char * name = cpl_sprintf("%s_%s", err_col,
2564 cr2res_obs_pol_corr_sfx[c]) ;
2565 if (cpl_table_has_column(corr, name))
2566 pc[c] = cpl_table_get_data_double_const(corr, name) ;
2567 else
2568 found = 0 ;
2569 cpl_free(name) ;
2570 }
2571 if (!found) continue ;
2572 perr = cpl_table_get_data_double(pol_spec, err_col) ;
2573 for (x=0 ; x<nrow ; x++) {
2574 double add = (pc[0][x] * pc[0][x] + pc[1][x] * pc[1][x] +
2575 pc[2][x] * pc[2][x] + pc[3][x] * pc[3][x] - pc[4][x]) /
2576 ((double)ngroups * ngroups) ;
2577 if (isfinite(perr[x]) && add > 0.0)
2578 perr[x] = sqrt(perr[x] * perr[x] + add) ;
2579 }
2580 }
2581 cpl_array_delete(names) ;
2582}
2583
2584/*----------------------------------------------------------------------------*/
2591/*----------------------------------------------------------------------------*/
2592static int cr2res_obs_pol_check_inputs_validity(
2593 const cpl_frameset * rawframes,
2594 cpl_size * ngroups)
2595{
2596 cr2res_nodding_pos * nod_positions ;
2597 cpl_frameset * raw_a ;
2598 cpl_frameset * raw_b ;
2599 cpl_size nframes, nframes_a, nframes_b, i ;
2600
2601 /* Check Inputs */
2602 if (rawframes == NULL || ngroups == NULL) return -1 ;
2603 nframes = cpl_frameset_get_size(rawframes) ;
2604
2605 /* Check that the MJD-OBS is increasing (only when several frames) */
2606 if (nframes > 1) {
2607 cpl_propertylist *plist;
2608 cpl_frame *cur_frame;
2609 const char *cur_fname;
2610 double mjd_obs_previous;
2611 /* Get the first MJDOBS */
2612 cur_frame = cpl_frameset_get_position((cpl_frameset*)rawframes, 0) ;
2613 cur_fname = cpl_frame_get_filename(cur_frame);
2614 plist = cpl_propertylist_load(cur_fname, 0) ;
2615 mjd_obs_previous = cr2res_pfits_get_mjd_obs(plist) ;
2616 cpl_propertylist_delete(plist) ;
2617 if (cpl_error_get_code() != CPL_ERROR_NONE) {
2618 cpl_msg_error(__func__, "Cannot access MJD-OBS") ;
2619 return -1 ;
2620 }
2621 /* Loop on the n-1 last frames */
2622 for (i=1 ; i<nframes ; i++) {
2623 double mjd_obs;
2624 /* Get the current MJDOBS */
2625 cur_frame = cpl_frameset_get_position((cpl_frameset*)rawframes, i) ;
2626 cur_fname = cpl_frame_get_filename(cur_frame);
2627 plist = cpl_propertylist_load(cur_fname, 0) ;
2628 mjd_obs = cr2res_pfits_get_mjd_obs(plist) ;
2629 cpl_propertylist_delete(plist) ;
2630 if (cpl_error_get_code() != CPL_ERROR_NONE) {
2631 cpl_msg_error(__func__, "Cannot access MJD-OBS") ;
2632 return -1 ;
2633 }
2634 if (mjd_obs <= mjd_obs_previous) {
2635 cpl_msg_error(__func__, "Rawframes MJD-OBS is not increasing") ;
2636 return 0 ;
2637 }
2638 mjd_obs_previous = mjd_obs ;
2639 }
2640 }
2641
2642 /* Get the nodding positions */
2643 nod_positions = cr2res_nodding_read_positions(rawframes) ;
2644 if (nod_positions == NULL) {
2645 cpl_msg_error(__func__, "Cannot read nodding positions") ;
2646 return -1 ;
2647 }
2648 if (cpl_msg_get_level() == CPL_MSG_DEBUG)
2649 for (i=0 ; i<nframes ; i++)
2650 cpl_msg_debug(__func__, "Frame %s - Nodding %c",
2651 cpl_frame_get_filename(
2652 cpl_frameset_get_position_const(rawframes,i)),
2653 cr2res_nodding_position_char(nod_positions[i])) ;
2654
2655 /* Split the frames */
2656 if (cr2res_combine_nodding_split_frames(rawframes, nod_positions,
2657 &raw_a, &raw_b)) {
2658 cpl_msg_error(__func__, "Cannot split the nodding positions") ;
2659 cpl_free(nod_positions) ;
2660 return -1 ;
2661 }
2662 cpl_free(nod_positions) ;
2663
2664 /* Check the sizes of A and B */
2665 nframes_a = cpl_frameset_get_size(raw_a) ;
2666 nframes_b = cpl_frameset_get_size(raw_b) ;
2667 if (nframes_a != nframes_b) {
2668 cpl_msg_error(__func__, "Number of A and B positions differ") ;
2669 cpl_frameset_delete(raw_a);
2670 cpl_frameset_delete(raw_b);
2671 return 0 ;
2672 }
2673 cpl_frameset_delete(raw_a);
2674 cpl_frameset_delete(raw_b);
2675
2676 /* number of frames must be modulo 4 */
2677 if (nframes_a % CR2RES_POLARIMETRY_GROUP_SIZE) {
2678 cpl_msg_error(__func__, "Require a multiple of 4 raw frames") ;
2679 return 0 ;
2680 }
2681
2682 /* Return */
2683 *ngroups = nframes_a/CR2RES_POLARIMETRY_GROUP_SIZE ;
2684 return 1 ;
2685}
2686
cr2res_calib_flat_relerr2
cpl_image * cr2res_calib_flat_relerr2(const hdrl_image *mflat)
Squared relative error of a master flat, (sigma_F/F)^2.
Definition cr2res_calib.c:596
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_calib_dark_var
int cr2res_calib_dark_var(hdrl_imagelist *list, const cpl_frame *master_dark_frame, int chip, const hdrl_image *mflat, const cpl_vector *dits, double sign)
Add or remove the master dark variance of cr2res_calib_image()
Definition cr2res_calib.c:645
cr2res_calib_flat_var
void cr2res_calib_flat_var(hdrl_image *ima, const cpl_image *relerr2, double sign)
Add or remove the flat-field variance term of an image.
Definition cr2res_calib.c:618
cr2res_dfs_SPEC_ERR_colname
char * cr2res_dfs_SPEC_ERR_colname(int order_idx, int trace)
Get the ERR column name for a given order/trace.
Definition cr2res_dfs.c:414
cr2res_dfs_POL_INTENS_ERROR_colname
char * cr2res_dfs_POL_INTENS_ERROR_colname(int order_idx)
Get the POL_INTENS_ERROR column name for a given order.
Definition cr2res_dfs.c:306
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_dfs_WAVELENGTH_colname
char * cr2res_dfs_WAVELENGTH_colname(int order_idx, int trace)
Get the WAVELENGTH column name for a given order/trace.
Definition cr2res_dfs.c:396
cr2res_dfs_SPEC_colname_parse
char * cr2res_dfs_SPEC_colname_parse(const char *colname, int *order_idx, int *trace)
Parse a column name ORDER_TRACE_TYPE format.
Definition cr2res_dfs.c:505
cr2res_dfs_POL_INTENS_colname
char * cr2res_dfs_POL_INTENS_colname(int order_idx)
Get the POL_INTENS column name for a given order.
Definition cr2res_dfs.c:290
cr2res_dfs_SPEC_colname
char * cr2res_dfs_SPEC_colname(int order_idx, int trace)
Get the SPEC column name for a given order/trace.
Definition cr2res_dfs.c:378
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_idp_save
int cr2res_idp_save(const char *filename, cpl_frameset *allframes, cpl_frameset *rawframes, const cpl_parameterlist *parlist, cpl_table **tables, cpl_propertylist *main_qc_plist, cpl_propertylist **ext_plist, const char *procatg, const char *recipe)
Save an IDP file.
Definition cr2res_idp.c:98
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_convert_order_idx_to_idxp
int cr2res_io_convert_order_idx_to_idxp(int order_idx)
Convert the order_idx to the order_idxp.
Definition cr2res_io.c:613
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_eop_table
cpl_table * cr2res_io_get_eop_table()
Get the eop_table.
Definition cr2res_io.c:120
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_save_CALIBRATED
int cr2res_io_save_CALIBRATED(const char *filename, cpl_frameset *allframes, cpl_frameset *inframes, const cpl_parameterlist *parlist, hdrl_image **calib_collapsed, const cpl_propertylist *qc_list, cpl_propertylist **ext_plist, const char *procatg, const char *recipe)
Save a CALIBRATED frame.
Definition cr2res_io.c:1638
cr2res_io_read_decker_positions
cr2res_decker * cr2res_io_read_decker_positions(const cpl_frameset *in)
Get the decker positions from a frame set.
Definition cr2res_io.c:521
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_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_io_save_POL_SPEC
int cr2res_io_save_POL_SPEC(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 Polarimetry spectrum.
Definition cr2res_io.c:2069
cr2res_io_load_MASTER_FLAT
hdrl_image * cr2res_io_load_MASTER_FLAT(const char *filename, int detector)
Load an hdrl image from a MASTER_FLAT.
Definition cr2res_io.c:1096
cr2res_pfits_get_mjd_obs
double cr2res_pfits_get_mjd_obs(const cpl_propertylist *plist)
find out the MJD-OBS value
Definition cr2res_pfits.c:217
cr2res_pol_demod_intens
cpl_bivector * cr2res_pol_demod_intens(cpl_vector **intens, cpl_vector **wl, cpl_vector **errors, int n)
Combine extracted spectra into Intensity spectrum.
Definition cr2res_pol.c:533
cr2res_pol_POL_SPEC_create
cpl_table * cr2res_pol_POL_SPEC_create(int *orders, cpl_vector **wl, cpl_bivector **stokes, cpl_bivector **null, cpl_bivector **intens, int norders)
Create the POL_SPEC table to be saved.
Definition cr2res_pol.c:678
cr2res_pol_spec_pol_merge
cpl_table * cr2res_pol_spec_pol_merge(const cpl_table **pol_spec_list, int pol_spec_nb)
Merge several POL_SPEC tables together, by averaging.
Definition cr2res_pol.c:899
cr2res_pol_sort_frames
int * cr2res_pol_sort_frames(const cpl_frame *frame1, const cpl_frame *frame2, const cpl_frame *frame3, const cpl_frame *frame4)
Compute the positions of the passed frames.
Definition cr2res_pol.c:815
cr2res_pol_get_beam_trace
cpl_table * cr2res_pol_get_beam_trace(const cpl_table *tw_in, cr2res_decker decker_position, int up_or_down)
Compute the traces for the polarimetric beams.
Definition cr2res_pol.c:958
cr2res_pol_demod_stokes
cpl_bivector * cr2res_pol_demod_stokes(cpl_vector **intens, cpl_vector **wl, cpl_vector **errors, int n)
Demodulate extracted spectra into Stokes parameter.
Definition cr2res_pol.c:106
cr2res_pol_demod_null
cpl_bivector * cr2res_pol_demod_null(cpl_vector **intens, cpl_vector **wl, cpl_vector **errors, int n)
Demodulate extracted spectra into Null spectrum.
Definition cr2res_pol.c:415
cr2res_qc_numsat
int cr2res_qc_numsat(const cpl_frameset *frameset)
Calculate the number of saturated pixels.
Definition cr2res_qc.c:1299
cr2res_qc_compute_snr
double cr2res_qc_compute_snr(cpl_vector *spec, cpl_vector *err)
Computes the SNR of one spectrum.
Definition cr2res_qc.c:1065
cr2res_trace_get_ycen
cpl_vector * cr2res_trace_get_ycen(const cpl_table *trace, int order_idx, int trace_nb, int size)
Retrieves the middle (All) polynomial from trace table and evaluates.
Definition cr2res_trace.c:769
cr2res_trace_get_height
int cr2res_trace_get_height(const cpl_table *trace, cpl_size order_idx, cpl_size trace_nb)
Computes the average height (pix) of an order, from trace polys.
Definition cr2res_trace.c:805
cr2res_trace_get_order_idx_values
int * cr2res_trace_get_order_idx_values(const cpl_table *trace, int *nb_order_idx_values)
Count and return the different order_idx values from a TW table.
Definition cr2res_trace.c:431
cr2res_decker_print_position
char * cr2res_decker_print_position(cr2res_decker dpos)
Get the decker position string for display.
Definition cr2res_utils.c:1276
cr2res_order_idx_to_real
int cr2res_order_idx_to_real(int order_idx, int order_zp)
Convert the order_idx into order_real.
Definition cr2res_utils.c:89
cr2res_combine_nodding_split_frames
int cr2res_combine_nodding_split_frames(const cpl_frameset *in, const cr2res_nodding_pos *positions, cpl_frameset **pos_a, cpl_frameset **pos_b)
Split A/B positions in 2 framesets.
Definition cr2res_nodding.c:100
cr2res_nodding_position_char
char cr2res_nodding_position_char(cr2res_nodding_pos pos)
Get the nodding position character for display.
Definition cr2res_nodding.c:83
cr2res_utils_get_center_mjd
double cr2res_utils_get_center_mjd(const cpl_frameset *frameset)
Calculate the middle of the exposures in the frameset.
Definition cr2res_utils.c:2289
cr2res_nodding_read_positions
cr2res_nodding_pos * cr2res_nodding_read_positions(const cpl_frameset *in)
Get the nodding positions from a frame set.
Definition cr2res_nodding.c:50
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
hdrl_barycorr_compute
cpl_error_code hdrl_barycorr_compute(double ra, double dec, const cpl_table *eop_table, double mjdobs, double time_to_mid_exposure, double longitude, double latitude, double elevation, double pressure, double temperature, double humidity, double wavelength, double *barycorr)
Derives the barycentric correction using the erfa function eraApco13(). The latter For a terrestrial ...
Definition hdrl_barycorr.c:198
hdrl_image_duplicate
hdrl_image * hdrl_image_duplicate(const hdrl_image *himg)
copy hdrl_image
Definition hdrl_image.c:391
hdrl_image_get_error_const
const cpl_image * hdrl_image_get_error_const(const hdrl_image *himg)
get error as cpl image
Definition hdrl_image.c:144
hdrl_image_get_image_const
const cpl_image * hdrl_image_get_image_const(const hdrl_image *himg)
get data as cpl image
Definition hdrl_image.c:118
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_const
const hdrl_image * hdrl_imagelist_get_const(const hdrl_imagelist *himlist, cpl_size inum)
Get an image from a list of images.
Definition hdrl_imagelist_io.c:229
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_sub_image
cpl_error_code hdrl_imagelist_sub_image(hdrl_imagelist *himlist, const hdrl_image *himg)
Subtract an image from an image list.
Definition hdrl_imagelist_basic.c:214
hdrl_imagelist_get
hdrl_image * hdrl_imagelist_get(const hdrl_imagelist *himlist, cpl_size inum)
Get an image from a list of images.
Definition hdrl_imagelist_io.c:210
Generated on Tue Oct 13 2026 12:27:26 for CR2RE Pipeline Reference Manual by doxygen 1.13.2