31#include "cr2res_utils.h"
32#include "cr2res_idp.h"
33#include "cr2res_calib.h"
34#include "cr2res_pfits.h"
35#include "cr2res_dfs.h"
36#include "cr2res_bpm.h"
37#include "cr2res_trace.h"
38#include "cr2res_extract.h"
46#define RECIPE_STRING "cr2res_obs_staring"
52int cpl_plugin_get_info(cpl_pluginlist * list);
58static cpl_frameset * cr2res_obs_staring_find_RAW(
59 const cpl_frameset * in) ;
60static int cr2res_obs_staring_check_inputs_validity(
61 const cpl_frameset * rawframes) ;
62static int cr2res_obs_staring_reduce(
63 const cpl_frameset * rawframes,
64 const cpl_frame * trace_wave_frame,
65 const cpl_frame * detlin_frame,
66 const cpl_frame * master_dark_frame,
67 const cpl_frame * master_flat_frame,
68 const cpl_frame * bpm_frame,
69 const cpl_frame * blaze_frame,
70 const cpl_array * slit_frac,
71 int subtract_nolight_rows,
72 int subtract_interorder_column,
75 int extract_oversample,
77 int extract_swath_width,
79 double extract_smooth_slit,
80 double extract_smooth_spec,
82 hdrl_image ** combined,
84 cpl_table ** slitfunc,
86 cpl_propertylist ** ext_plist) ;
87static int cr2res_obs_staring_create(cpl_plugin *);
88static int cr2res_obs_staring_exec(cpl_plugin *);
89static int cr2res_obs_staring_destroy(cpl_plugin *);
90static int cr2res_obs_staring(cpl_frameset *,
const cpl_parameterlist *);
96static char cr2res_obs_staring_description[] =
"\
97Staring Observation \n\
98 This recipe handles staring observations. \n\
101 raw.fits " CR2RES_OBS_STARING_OTHER_RAW
" [1 to n] \n\
102 or " CR2RES_OBS_STARING_JITTER_RAW
" [1 to n] \n\
103 or " CR2RES_OBS_STARING_WAVE_SKY_RAW
" [1 to n] \n\
104 trace.fits " CR2RES_CAL_FLAT_TW_PROCATG
" [1] \n\
105 or " CR2RES_CAL_FLAT_TW_MERGED_PROCATG
" \n\
106 or " CR2RES_UTIL_TRACE_TW_PROCATG
" \n\
107 or " CR2RES_UTIL_WAVE_TW_PROCATG
" \n\
108 or " CR2RES_CAL_WAVE_TW_PROCATG
" \n\
109 or " CR2RES_UTIL_SLIT_CURV_TW_PROCATG
" \n\
110 detlin.fits " CR2RES_CAL_DETLIN_COEFFS_PROCATG
" [0 to 1] \n\
111 bpm.fits " CR2RES_CAL_DARK_BPM_PROCATG
" [0 to 1] \n\
112 or " CR2RES_CAL_FLAT_BPM_PROCATG
" \n\
113 or " CR2RES_CAL_DETLIN_BPM_PROCATG
" \n\
114 or " CR2RES_UTIL_BPM_MERGE_PROCATG
" \n\
115 or " CR2RES_UTIL_BPM_SPLIT_PROCATG
" \n\
116 master_dark.fits " CR2RES_CAL_DARK_MASTER_PROCATG
" [0 to 1] \n\
117 master_flat.fits " CR2RES_CAL_FLAT_MASTER_PROCATG
" [0 to 1] \n\
118 blaze.fits " CR2RES_CAL_FLAT_EXTRACT_1D_PROCATG
" [0 to 1] \n\
121 cr2res_obs_staring_slitfunc.fits "
122 CR2RES_OBS_STARING_SLITFUNC_PROCATG
"\n\
123 cr2res_obs_staring_model.fits "
124 CR2RES_OBS_STARING_SLITMODEL_PROCATG
"\n\
125 cr2res_obs_staring_extracted.fits "
126 CR2RES_OBS_STARING_EXTRACT_PROCATG
"\n\
129 loop on detectors d: \n\
130 call cr2res_obs_staring_reduce() \n\
139 cr2res_obs_staring_reduce() \n\
140 Load the input raw frames in an image list \n\
141 Apply the calibrations to the image list \n\
142 Remove the master dark and flat errors from each image \n\
143 Collapse the image list \n\
144 Add the master dark and flat errors to the collapsed image \n\
145 Load the input trace wave \n\
146 Recompute a new trace wave with the specified slit fraction \n\
147 (--slit_frac) if needed \n\
148 Extract the spectra from the collapsed image \n\
153 Compute QC parameters \n\
155 Library functions used \n\
156 cr2res_io_find_TRACE_WAVE() \n\
157 cr2res_io_find_BPM() \n\
158 cr2res_obs_staring_reduce() \n\
159 cr2res_io_read_dits() \n\
160 cr2res_io_load_image_list_from_set() \n\
161 cr2res_calib_imagelist() \n\
162 cr2res_io_load_TRACE_WAVE() \n\
163 cr2res_extract_traces() \n\
164 cr2res_io_save_COMBINED() \n\
165 cr2res_io_save_EXTRACT_1D() \n\
166 cr2res_io_save_SLIT_FUNC() \n\
184int cpl_plugin_get_info(cpl_pluginlist * list)
186 cpl_recipe * recipe = cpl_calloc(1,
sizeof *recipe );
187 cpl_plugin * plugin = &recipe->interface;
189 if (cpl_plugin_init(plugin,
191 CR2RES_BINARY_VERSION,
192 CPL_PLUGIN_TYPE_RECIPE,
194 "Staring Observation recipe",
195 cr2res_obs_staring_description,
196 CR2RES_PIPELINE_AUTHORS,
199 cr2res_obs_staring_create,
200 cr2res_obs_staring_exec,
201 cr2res_obs_staring_destroy)) {
202 cpl_msg_error(cpl_func,
"Plugin initialization failed");
203 (void)cpl_error_set_where(cpl_func);
207 if (cpl_pluginlist_append(list, plugin)) {
208 cpl_msg_error(cpl_func,
"Error adding plugin to list");
209 (void)cpl_error_set_where(cpl_func);
225static int cr2res_obs_staring_create(cpl_plugin * plugin)
227 cpl_recipe * recipe ;
231 if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE)
232 recipe = (cpl_recipe *)plugin;
237 recipe->parameters = cpl_parameterlist_new();
241 p = cpl_parameter_new_value(
242 "cr2res.cr2res_obs_staring.subtract_nolight_rows",
244 "Subtract median row from baffled region at detector bottom",
245 "cr2res.cr2res_obs_staring", FALSE);
246 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI,
"subtract_nolight_rows");
247 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
248 cpl_parameterlist_append(recipe->parameters, p);
250 p = cpl_parameter_new_value(
251 "cr2res.cr2res_obs_staring.subtract_interorder_column",
253 "Subtract column-by-column fit to the pixel values between"
255 "cr2res.cr2res_obs_staring", TRUE);
256 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI,
257 "subtract_interorder_column");
258 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
259 cpl_parameterlist_append(recipe->parameters, p);
261 p = cpl_parameter_new_value(
"cr2res.cr2res_obs_staring.cosmics",
262 CPL_TYPE_BOOL,
"Find and mark cosmic rays hits as bad",
263 "cr2res.cr2res_obs_staring", FALSE);
264 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI,
"cosmics");
265 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
266 cpl_parameterlist_append(recipe->parameters, p);
268 p = cpl_parameter_new_enum(
"cr2res.cr2res_obs_staring.error_method",
269 CPL_TYPE_STRING,
"The 1d extraction error calculation method",
270 "cr2res.cr2res_obs_staring",
271 "Horne", 2,
"Poisson",
"Horne");
272 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI,
"error_method");
273 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
274 cpl_parameterlist_append(recipe->parameters, p);
276 p = cpl_parameter_new_value(
"cr2res.cr2res_obs_staring.slit_frac",
277 CPL_TYPE_STRING,
"Wished slit fraction",
278 "cr2res.cr2res_obs_staring",
"-1.0, -1.0");
279 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI,
"slit_frac");
280 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
281 cpl_parameterlist_append(recipe->parameters, p);
283 p = cpl_parameter_new_value(
"cr2res.cr2res_obs_staring.extract_oversample",
284 CPL_TYPE_INT,
"factor by which to oversample the extraction",
285 "cr2res.cr2res_obs_staring", 5);
286 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI,
"extract_oversample");
287 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
288 cpl_parameterlist_append(recipe->parameters, p);
290 p = cpl_parameter_new_value(
"cr2res.cr2res_obs_staring.extract_pclip",
291 CPL_TYPE_DOUBLE,
"percentage of high and low pixels to clip on first iteration",
292 "cr2res.cr2res_obs_staring", 0.1);
293 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI,
"extract_pclip");
294 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
295 cpl_parameterlist_append(recipe->parameters, p);
297 p = cpl_parameter_new_value(
"cr2res.cr2res_obs_staring.extract_swath_width",
298 CPL_TYPE_INT,
"The swath width",
"cr2res.cr2res_obs_staring", 800);
299 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI,
"extract_swath_width");
300 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
301 cpl_parameterlist_append(recipe->parameters, p);
303 p = cpl_parameter_new_value(
"cr2res.cr2res_obs_staring.extract_height",
304 CPL_TYPE_INT,
"Extraction height",
305 "cr2res.cr2res_obs_staring", -1);
306 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI,
"extract_height");
307 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
308 cpl_parameterlist_append(recipe->parameters, p);
310 p = cpl_parameter_new_value(
"cr2res.cr2res_obs_staring.extract_smooth_slit",
312 "Smoothing along the slit (1 for high S/N, 5 for low)",
313 "cr2res.cr2res_obs_staring", 2.0);
314 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI,
"extract_smooth_slit");
315 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
316 cpl_parameterlist_append(recipe->parameters, p);
318 p = cpl_parameter_new_value(
"cr2res.cr2res_obs_staring.extract_smooth_spec",
319 CPL_TYPE_DOUBLE,
"Smoothing along the spectrum",
320 "cr2res.cr2res_obs_staring", 0.0);
321 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI,
"extract_smooth_spec");
322 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
323 cpl_parameterlist_append(recipe->parameters, p);
325 p = cpl_parameter_new_value(
"cr2res.cr2res_obs_staring.detector",
326 CPL_TYPE_INT,
"Only reduce the specified detector",
327 "cr2res.cr2res_obs_staring", 0);
328 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI,
"detector");
329 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
330 cpl_parameterlist_append(recipe->parameters, p);
332 p = cpl_parameter_new_value(
"cr2res.cr2res_obs_staring.create_idp",
333 CPL_TYPE_BOOL,
"Flag to produce IDP files",
334 "cr2res.cr2res_obs_staring", FALSE);
335 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI,
"idp");
336 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
337 cpl_parameterlist_append(recipe->parameters, p);
339 p = cpl_parameter_new_value(
"cr2res.cr2res_obs_staring.display_order",
340 CPL_TYPE_INT,
"Apply the display for the specified order",
341 "cr2res.cr2res_obs_staring", 0);
342 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI,
"display_order");
343 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
344 cpl_parameterlist_append(recipe->parameters, p);
346 p = cpl_parameter_new_value(
"cr2res.cr2res_obs_staring.display_trace",
347 CPL_TYPE_INT,
"Apply the display for the specified trace",
348 "cr2res.cr2res_obs_staring", 0);
349 cpl_parameter_set_alias(p, CPL_PARAMETER_MODE_CLI,
"display_trace");
350 cpl_parameter_disable(p, CPL_PARAMETER_MODE_ENV);
351 cpl_parameterlist_append(recipe->parameters, p);
363static int cr2res_obs_staring_exec(cpl_plugin * plugin)
368 if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE)
369 recipe = (cpl_recipe *)plugin;
372 return cr2res_obs_staring(recipe->frames, recipe->parameters);
382static int cr2res_obs_staring_destroy(cpl_plugin * plugin)
387 if (cpl_plugin_get_type(plugin) == CPL_PLUGIN_TYPE_RECIPE)
388 recipe = (cpl_recipe *)plugin;
391 cpl_parameterlist_delete(recipe->parameters);
403static int cr2res_obs_staring(
404 cpl_frameset * frameset,
405 const cpl_parameterlist * parlist)
407 const cpl_parameter * param ;
408 int subtract_nolight_rows, subtract_interorder_column,
409 extract_oversample, create_idp, cosmics,
411 extract_swath_width, extract_height, reduce_det;
412 double extract_smooth_slit, extract_smooth_spec,
413 extract_pclip, slit_low, slit_up ;
415 cpl_array * slit_frac ;
416 cpl_frameset * rawframes ;
417 const cpl_frame * trace_wave_frame ;
418 const cpl_frame * detlin_frame ;
419 const cpl_frame * master_dark_frame ;
420 const cpl_frame * master_flat_frame ;
421 const cpl_frame * bpm_frame ;
422 const cpl_frame * blaze_frame ;
424 hdrl_image * combined[CR2RES_NB_DETECTORS] ;
425 cpl_table * extract[CR2RES_NB_DETECTORS] ;
426 cpl_table * slitfunc[CR2RES_NB_DETECTORS] ;
427 hdrl_image * model[CR2RES_NB_DETECTORS] ;
428 cpl_propertylist * qc_main ;
429 cpl_propertylist * ext_plist[CR2RES_NB_DETECTORS] ;
431 cpl_table * eop_table ;
435 param = cpl_parameterlist_find_const(parlist,
436 "cr2res.cr2res_obs_staring.subtract_nolight_rows");
437 subtract_nolight_rows = cpl_parameter_get_bool(param);
438 param = cpl_parameterlist_find_const(parlist,
439 "cr2res.cr2res_obs_staring.subtract_interorder_column");
440 subtract_interorder_column = cpl_parameter_get_bool(param);
441 param = cpl_parameterlist_find_const(parlist,
442 "cr2res.cr2res_obs_staring.cosmics");
443 cosmics = cpl_parameter_get_bool(param);
444 param = cpl_parameterlist_find_const(parlist,
445 "cr2res.cr2res_obs_staring.error_method");
446 if (strcmp(cpl_parameter_get_string(param),
"Horne") == 0)
447 error_method = CR2RES_EXTRACT_ERROR_HORNE;
449 error_method = CR2RES_EXTRACT_ERROR_POISSON;
450 param = cpl_parameterlist_find_const(parlist,
451 "cr2res.cr2res_obs_staring.extract_oversample");
452 extract_oversample = cpl_parameter_get_int(param);
453 param = cpl_parameterlist_find_const(parlist,
454 "cr2res.cr2res_obs_staring.extract_pclip");
455 extract_pclip = cpl_parameter_get_double(param);
456 param = cpl_parameterlist_find_const(parlist,
457 "cr2res.cr2res_obs_staring.extract_swath_width");
458 extract_swath_width = cpl_parameter_get_int(param);
459 param = cpl_parameterlist_find_const(parlist,
460 "cr2res.cr2res_obs_staring.extract_height");
461 extract_height = cpl_parameter_get_int(param);
462 param = cpl_parameterlist_find_const(parlist,
463 "cr2res.cr2res_obs_staring.extract_smooth_slit");
464 extract_smooth_slit = cpl_parameter_get_double(param);
465 param = cpl_parameterlist_find_const(parlist,
466 "cr2res.cr2res_obs_staring.extract_smooth_spec");
467 extract_smooth_spec = cpl_parameter_get_double(param);
468 param = cpl_parameterlist_find_const(parlist,
469 "cr2res.cr2res_obs_staring.detector");
470 reduce_det = cpl_parameter_get_int(param);
471 param = cpl_parameterlist_find_const(parlist,
472 "cr2res.cr2res_obs_staring.create_idp");
473 create_idp = cpl_parameter_get_bool(param);
480 param = cpl_parameterlist_find_const(parlist,
481 "cr2res.cr2res_obs_staring.slit_frac");
482 sval = cpl_parameter_get_string(param) ;
483 if (sscanf(sval,
"%lg,%lg", &slit_low, &slit_up) != 2) {
484 cpl_msg_error(__func__,
"Invalid Slit Fraction specified");
485 cpl_error_set(__func__, CPL_ERROR_ILLEGAL_INPUT) ;
490 if (slit_low >= 0.0 && slit_up >= 0.0 && slit_low <= 1.0 && slit_up <= 1.0
491 && slit_up > slit_low) {
492 slit_frac = cpl_array_new(3, CPL_TYPE_DOUBLE) ;
493 cpl_array_set(slit_frac, 0, slit_low) ;
494 cpl_array_set(slit_frac, 1, (slit_low+slit_up)/2.0) ;
495 cpl_array_set(slit_frac, 2, slit_up) ;
502 if (slit_frac != NULL) cpl_array_delete(slit_frac) ;
503 cpl_msg_error(__func__,
"Cannot identify RAW and CALIB frames") ;
504 cpl_error_set(__func__, CPL_ERROR_ILLEGAL_INPUT) ;
510 if (trace_wave_frame == NULL) {
511 if (slit_frac != NULL) cpl_array_delete(slit_frac) ;
512 cpl_msg_error(__func__,
"Could not find TRACE_WAVE frame") ;
515 detlin_frame = cpl_frameset_find_const(frameset,
516 CR2RES_CAL_DETLIN_COEFFS_PROCATG);
517 master_dark_frame = cpl_frameset_find_const(frameset,
518 CR2RES_CAL_DARK_MASTER_PROCATG) ;
519 master_flat_frame = cpl_frameset_find_const(frameset,
520 CR2RES_CAL_FLAT_MASTER_PROCATG) ;
522 blaze_frame = cpl_frameset_find_const(frameset,
523 CR2RES_CAL_FLAT_EXTRACT_1D_PROCATG) ;
526 rawframes = cr2res_obs_staring_find_RAW(frameset) ;
527 if (rawframes == NULL) {
528 if (slit_frac != NULL) cpl_array_delete(slit_frac) ;
529 cpl_msg_error(__func__,
"Could not find RAW frames") ;
534 for (det_nr=1 ; det_nr<=CR2RES_NB_DETECTORS ; det_nr++) {
536 combined[det_nr-1] = NULL ;
537 extract[det_nr-1] = NULL ;
538 slitfunc[det_nr-1] = NULL ;
539 model[det_nr-1] = NULL ;
540 ext_plist[det_nr-1] = NULL ;
543 if (reduce_det != 0 && det_nr != reduce_det) continue ;
545 cpl_msg_info(__func__,
"Process Detector %d", det_nr) ;
546 cpl_msg_indent_more() ;
549 if (cr2res_obs_staring_reduce(rawframes,
550 trace_wave_frame, detlin_frame, master_dark_frame,
551 master_flat_frame, bpm_frame, blaze_frame, slit_frac,
552 subtract_nolight_rows, subtract_interorder_column,
553 cosmics, error_method, extract_oversample, extract_pclip,
554 extract_swath_width, extract_height, extract_smooth_slit,
555 extract_smooth_spec, det_nr,
556 &(combined[det_nr-1]),
557 &(extract[det_nr-1]),
558 &(slitfunc[det_nr-1]),
560 &(ext_plist[det_nr-1])) == -1) {
561 cpl_msg_warning(__func__,
"Failed to reduce detector %d", det_nr);
563 cpl_msg_indent_less() ;
565 cpl_msg_indent_less() ;
567 if (slit_frac != NULL) cpl_array_delete(slit_frac) ;
572 qc_main = cpl_propertylist_new();
573 cpl_propertylist_append_double(qc_main,
574 CR2RES_HEADER_DRS_TMID,
579 if (eop_table != NULL) {
580 double ra, dec, mjd_obs, geolon, geolat, geoelev;
581 cpl_propertylist *plist;
582 plist=cpl_propertylist_load(cpl_frame_get_filename(
583 cpl_frameset_get_position_const(rawframes, 0)), 0) ;
585 ra = cpl_propertylist_get_double(plist,
"RA") ;
586 dec = cpl_propertylist_get_double(plist,
"DEC") ;
587 mjd_obs = cpl_propertylist_get_double(plist,
"MJD-OBS") ;
588 geolon = cpl_propertylist_get_double(plist,
"ESO TEL GEOLON") ;
589 geolat = cpl_propertylist_get_double(plist,
"ESO TEL GEOLAT") ;
590 geoelev = cpl_propertylist_get_double(plist,
"ESO TEL GEOELEV") ;
592 cpl_propertylist_delete(plist) ;
595 if (!cpl_error_get_code()) {
599 (mjd_cen-mjd_obs)*24*3600, geolon, geolat, geoelev,
600 0.0, 0.0, 0.0, 0.0, &barycorr);
602 cpl_msg_info(__func__,
"Barycentric correction: %g m/s",
606 cpl_msg_info(__func__,
"Cannot derive Barycentric correction");
609 cpl_table_delete(eop_table) ;
610 cpl_propertylist_append_double(qc_main, CR2RES_HEADER_DRS_BARYCORR,
614 if (error_method == CR2RES_EXTRACT_ERROR_HORNE) {
615 cpl_propertylist_append_string(qc_main, CR2RES_HEADER_DRS_ERRMETHOD,
618 cpl_propertylist_append_string(qc_main, CR2RES_HEADER_DRS_ERRMETHOD,
623 cpl_propertylist_append_int(qc_main,
624 CR2RES_HEADER_QC_NUMSAT,
628 if (trace_wave_frame != NULL)
629 cpl_frameset_insert(rawframes,
630 cpl_frame_duplicate(trace_wave_frame)) ;
631 if (detlin_frame != NULL)
632 cpl_frameset_insert(rawframes,
633 cpl_frame_duplicate(detlin_frame)) ;
634 if (master_dark_frame != NULL)
635 cpl_frameset_insert(rawframes,
636 cpl_frame_duplicate(master_dark_frame)) ;
637 if (master_flat_frame!= NULL)
638 cpl_frameset_insert(rawframes,
639 cpl_frame_duplicate(master_flat_frame)) ;
640 if (bpm_frame!= NULL)
641 cpl_frameset_insert(rawframes,
642 cpl_frame_duplicate(bpm_frame)) ;
643 if (blaze_frame!= NULL)
644 cpl_frameset_insert(rawframes,
645 cpl_frame_duplicate(blaze_frame)) ;
647 out_file = cpl_sprintf(
"%s_combined.fits", RECIPE_STRING) ;
649 combined, qc_main, ext_plist,
650 CR2RES_OBS_STARING_COMBINED_PROCATG, RECIPE_STRING) ;
653 out_file = cpl_sprintf(
"%s_slitfunc.fits", RECIPE_STRING) ;
655 slitfunc, qc_main, ext_plist, CR2RES_OBS_STARING_SLITFUNC_PROCATG,
659 out_file = cpl_sprintf(
"%s_model.fits", RECIPE_STRING) ;
661 model, qc_main, ext_plist, CR2RES_OBS_STARING_SLITMODEL_PROCATG,
665 out_file = cpl_sprintf(
"%s_extracted.fits", RECIPE_STRING) ;
667 qc_main, ext_plist, CR2RES_OBS_STARING_EXTRACT_PROCATG,
671 extract, qc_main, ext_plist,
672 CR2RES_OBS_STARING_IDP_PROCATG,
676 cpl_frameset_delete(rawframes) ;
679 cpl_propertylist_delete(qc_main) ;
680 for (det_nr=1 ; det_nr<=CR2RES_NB_DETECTORS ; det_nr++) {
681 if (combined[det_nr-1] != NULL)
683 if (extract[det_nr-1] != NULL)
684 cpl_table_delete(extract[det_nr-1]) ;
685 if (slitfunc[det_nr-1] != NULL)
686 cpl_table_delete(slitfunc[det_nr-1]) ;
687 if (model[det_nr-1] != NULL)
689 if (ext_plist[det_nr-1] != NULL)
690 cpl_propertylist_delete(ext_plist[det_nr-1]) ;
693 return (
int)cpl_error_get_code();
723static int cr2res_obs_staring_reduce(
724 const cpl_frameset * rawframes,
725 const cpl_frame * trace_wave_frame,
726 const cpl_frame * detlin_frame,
727 const cpl_frame * master_dark_frame,
728 const cpl_frame * master_flat_frame,
729 const cpl_frame * bpm_frame,
730 const cpl_frame * blaze_frame,
731 const cpl_array * slit_frac,
732 int subtract_nolight_rows,
733 int subtract_interorder_column,
736 int extract_oversample,
737 double extract_pclip,
738 int extract_swath_width,
740 double extract_smooth_slit,
741 double extract_smooth_spec,
743 hdrl_image ** combined,
744 cpl_table ** extract,
745 cpl_table ** slitfunc,
747 cpl_propertylist ** ext_plist)
749 hdrl_imagelist * in ;
750 hdrl_imagelist * in_calib ;
753 cpl_vector * ndsamples ;
754 cpl_table * blaze_table ;
755 cpl_table * trace_wave ;
756 cpl_table * trace_wave_new ;
757 hdrl_image * collapsed ;
759 cpl_image * flat_relerr2 ;
760 cpl_image * contrib ;
761 cpl_propertylist * plist ;
763 hdrl_image * model_master ;
764 cpl_table * slit_func ;
765 cpl_table * extracted ;
766 int * order_idx_values ;
768 cpl_array * fwhm_array ;
770 const char * first_fname ;
771 double qc_signal, qc_fwhm, gain, error_factor, blaze_norm ;
772 int order_zp, nb_order_idx_values,
773 order_idx, order_idxp ;
776 if (combined == NULL || extract == NULL || ext_plist == NULL ||
777 rawframes == NULL || trace_wave_frame == NULL)
return -1 ;
780 if (reduce_det == 1) gain = CR2RES_GAIN_CHIP1 ;
781 else if (reduce_det == 2) gain = CR2RES_GAIN_CHIP2 ;
782 else if (reduce_det == 3) gain = CR2RES_GAIN_CHIP3 ;
784 cpl_msg_error(__func__,
"Failed to get the Gain value") ;
789 if (cr2res_obs_staring_check_inputs_validity(rawframes) != 1) {
790 cpl_msg_error(__func__,
"Invalid Inputs") ;
795 first_fname = cpl_frame_get_filename(
796 cpl_frameset_get_position_const(rawframes, 0)) ;
799 if ((plist = cpl_propertylist_load(cpl_frame_get_filename(trace_wave_frame),
801 cpl_msg_error(__func__,
"Cannot read the ORDER_ZP from the input TW") ;
804 order_zp = cr2res_pfits_get_order_zp(plist) ;
805 cpl_propertylist_delete(plist) ;
806 if (cpl_error_get_code()) {
807 cpl_msg_error(__func__,
"Missing ORDER_ZP in the header - Skip") ;
816 if (cpl_msg_get_level() == CPL_MSG_DEBUG && dits != NULL)
817 cpl_vector_dump(dits, stdout) ;
824 reduce_det)) == NULL) {
825 cpl_msg_error(__func__,
"Cannot load images") ;
826 if (dits != NULL) cpl_vector_delete(dits) ;
830 cpl_msg_error(__func__,
"Inconsistent number of loaded images") ;
831 if (dits != NULL) cpl_vector_delete(dits) ;
837 if (error_method == CR2RES_EXTRACT_ERROR_HORNE) {
838 error_factor = -1.0 ;
840 error_factor = gain * cpl_vector_get(ndits, 0) *
841 cpl_frameset_get_size(rawframes) ;
844 for (i=0; i<cpl_vector_get_size(ndits); i++){
845 if (cpl_vector_get(ndits,i) != cpl_vector_get(ndits, 0))
846 cpl_msg_warning(__func__,
"Raw frames have different NDIT! "
847 "Error spectrum will likely be scaled incorrectly.");
853 subtract_nolight_rows, subtract_interorder_column, cosmics,
854 master_flat_frame, master_dark_frame, bpm_frame, detlin_frame,
855 dits, ndits, ndsamples) ;
856 cpl_vector_delete(ndsamples) ;
857 if (in_calib == NULL) {
858 cpl_msg_error(__func__,
"Failed to apply the calibrations") ;
859 if (dits != NULL) cpl_vector_delete(dits) ;
860 if (ndits != NULL) cpl_vector_delete(ndits) ;
865 if (ndits != NULL) cpl_vector_delete(ndits) ;
870 flat_relerr2 = NULL ;
871 if (master_flat_frame != NULL) {
873 cpl_frame_get_filename(master_flat_frame),
874 reduce_det)) == NULL) {
875 cpl_msg_error(__func__,
"Cannot load the master flat") ;
876 if (dits != NULL) cpl_vector_delete(dits) ;
883 flat_relerr2, -1.0) ;
886 reduce_det, mflat, dits, -1.0)) {
887 cpl_msg_error(__func__,
"Cannot remove the master dark errors") ;
888 cpl_vector_delete(dits) ;
891 if (flat_relerr2 != NULL) cpl_image_delete(flat_relerr2) ;
896 cpl_msg_info(__func__,
"Collapse") ;
897 cpl_msg_indent_more() ;
900 cpl_msg_error(__func__,
"Failed to Collapse") ;
902 if (dits != NULL) cpl_vector_delete(dits) ;
904 if (flat_relerr2 != NULL) cpl_image_delete(flat_relerr2) ;
905 cpl_msg_indent_less() ;
908 cpl_image_delete(contrib) ;
910 cpl_msg_indent_less() ;
912 if (flat_relerr2 != NULL) {
914 cpl_image_delete(flat_relerr2) ;
918 cpl_vector * mean_dit = cpl_vector_new(1) ;
920 cpl_vector_set(mean_dit, 0, cpl_vector_get_mean(dits)) ;
923 mflat, mean_dit, 1.0) ;
926 cpl_vector_delete(mean_dit) ;
927 cpl_vector_delete(dits) ;
929 cpl_msg_error(__func__,
"Cannot add the master dark errors") ;
938 cpl_msg_info(__func__,
"Load the TRACE WAVE") ;
940 trace_wave_frame), reduce_det)) == NULL) {
941 cpl_msg_error(__func__,
"Failed to Load the traces file") ;
947 if (slit_frac != NULL) {
949 trace_wave, slit_frac)) == NULL) {
950 cpl_msg_warning(__func__,
951 "Failed to compute the traces for user specified slit fraction") ;
954 cpl_table_delete(trace_wave) ;
955 trace_wave = trace_wave_new ;
956 trace_wave_new = NULL ;
963 if (blaze_frame != NULL) {
964 cpl_msg_info(__func__,
"Load the BLAZE") ;
966 blaze_frame), reduce_det)) == NULL) {
967 cpl_msg_error(__func__,
"Failed to Load the Blaze file") ;
969 cpl_table_delete(trace_wave) ;
972 cpl_propertylist * blaze_plist ;
973 blaze_plist = cpl_propertylist_load_regexp(cpl_frame_get_filename(
974 blaze_frame),0,CR2RES_HEADER_QC_BLAZE_NORM,0);
975 if(cpl_propertylist_get_size(blaze_plist)>0){
976 blaze_norm = cpl_propertylist_get_double(blaze_plist, CR2RES_HEADER_QC_BLAZE_NORM);
980 cpl_msg_warning(__func__,
"QC BLAZE NORM value not found, reverting to per trace normalization") ;
982 if(blaze_plist!=NULL){
983 cpl_propertylist_delete(blaze_plist);
988 int extract_niter = 10;
989 double extract_kappa = 10;
992 cpl_msg_info(__func__,
"Spectra Extraction") ;
994 CR2RES_EXTR_OPT_CURV, extract_height, extract_swath_width,
995 extract_oversample, extract_pclip, extract_smooth_slit, extract_smooth_spec,
996 extract_niter, extract_kappa, error_factor,
997 0, 0, 0, &extracted, &slit_func, &model_master) == -1) {
998 cpl_msg_error(__func__,
"Failed to extract");
1000 cpl_table_delete(trace_wave) ;
1003 if (blaze_table != NULL) cpl_table_delete(blaze_table) ;
1006 plist = cpl_propertylist_load(first_fname,
1011 cpl_propertylist_append_double(plist, CR2RES_HEADER_QC_SIGNAL, qc_signal) ;
1014 qc_snrs =
cr2res_qc_snr(trace_wave, extracted, &order_idx_values,
1015 &nb_order_idx_values) ;
1016 for (i=0 ; i<nb_order_idx_values ; i++) {
1017 order_idx = order_idx_values[i] ;
1019 key_name = cpl_sprintf(CR2RES_HEADER_QC_SNR, order_idxp) ;
1020 cpl_propertylist_append_double(plist, key_name, qc_snrs[i]) ;
1021 cpl_free(key_name) ;
1023 cpl_free(order_idx_values) ;
1028 &nb_order_idx_values);
1031 fwhm_array = cpl_array_new(nb_order_idx_values, CPL_TYPE_DOUBLE);
1032 for (i=0 ; i<nb_order_idx_values ; i++) {
1033 order_idx = order_idx_values[i] ;
1036 extract_oversample);
1038 key_name = cpl_sprintf(CR2RES_HEADER_QC_SLITFWHM_ORDER, order_idxp) ;
1039 cpl_propertylist_append_double(plist, key_name, qc_fwhm) ;
1040 cpl_free(key_name) ;
1041 cpl_array_set(fwhm_array, i, qc_fwhm) ;
1043 cpl_free(order_idx_values) ;
1044 qc_fwhm = cpl_array_get_median(fwhm_array);
1045 cpl_propertylist_append_double(plist, CR2RES_HEADER_QC_SLITFWHM_MED,
1047 if (qc_fwhm < 3.5) {
1048 cpl_msg_warning(__func__,
"Median FWHM of the PSF along the slit "
1049 "is %gpix, i.e. below the slit width. This means the slit "
1050 "is likely not evenly filled with light "
1051 "in the spectral direction. This can result in a "
1052 "wavelength offset between different positions along the slit,"
1053 " and with respect to calibrations."
1057 cpl_array_delete(fwhm_array) ;
1063 &nb_order_idx_values);
1066 for (i = 0; i < nb_order_idx_values; i++) {
1068 order_idx = order_idx_values[i] ;
1071 key_name = cpl_sprintf(CR2RES_HEADER_QC_REAL_ORDER, order_idxp) ;
1072 cpl_propertylist_append_int(plist, key_name, order_real) ;
1073 cpl_free(key_name) ;
1075 cpl_free(order_idx_values) ;
1077 cpl_table_delete(trace_wave) ;
1080 *combined = collapsed ;
1081 *extract = extracted ;
1082 *slitfunc = slit_func ;
1083 *model = model_master ;
1084 *ext_plist = plist ;
1096static int cr2res_obs_staring_check_inputs_validity(
1097 const cpl_frameset * rawframes)
1100 if (rawframes == NULL)
return -1 ;
1114static cpl_frameset * cr2res_obs_staring_find_RAW(
1115 const cpl_frameset * in)
1117 cpl_frameset * out ;
1120 if (in == NULL)
return NULL ;
cpl_image * cr2res_calib_flat_relerr2(const hdrl_image *mflat)
Squared relative error of a master flat, (sigma_F/F)^2.
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.
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()
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.
cpl_error_code cr2res_dfs_set_groups(cpl_frameset *set)
Set the group as RAW or CALIB in a frameset.
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.
cpl_table * cr2res_io_load_EXTRACT_1D(const char *filename, int detector)
Load a table from a EXTRACT_1D.
const cpl_frame * cr2res_io_find_BPM(const cpl_frameset *in)
Get the first CR2RES_BPM_DRSTYPE frame from a frameset.
int cr2res_io_save_SLIT_MODEL(const char *filename, cpl_frameset *allframes, cpl_frameset *inframes, const cpl_parameterlist *parlist, hdrl_image **data, const cpl_propertylist *qc_list, cpl_propertylist **ext_plist, const char *procatg, const char *recipe)
Save a SLIT_MODEL.
const cpl_frame * cr2res_io_find_TRACE_WAVE(const cpl_frameset *in)
Get the first CR2RES_TW_DRSTYPE frame from a frameset.
cpl_vector * cr2res_io_read_ndsamples(const cpl_frameset *in)
Get the NDSAMPLES from a frame set.
int cr2res_io_convert_order_idx_to_idxp(int order_idx)
Convert the order_idx to the order_idxp.
cpl_vector * cr2res_io_read_dits(const cpl_frameset *in)
Get the DITS from a frame set.
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.
hdrl_imagelist * cr2res_io_load_image_list_from_set(const cpl_frameset *in, int detector)
Load an hdrl image list from an images frameset.
cpl_table * cr2res_io_get_eop_table()
Get the eop_table.
int cr2res_io_get_ext_idx(const char *filename, int detector, int data)
Get the wished extension number for a detector.
int cr2res_io_save_SLIT_FUNC(const char *filename, cpl_frameset *allframes, cpl_frameset *inframes, const cpl_parameterlist *parlist, cpl_table **slit_func, const cpl_propertylist *qc_list, cpl_propertylist **ext_plist, const char *procatg, const char *recipe)
Save a SLIT_FUNC.
int cr2res_io_save_COMBINED(const char *filename, cpl_frameset *allframes, cpl_frameset *inframes, const cpl_parameterlist *parlist, hdrl_image **data, const cpl_propertylist *qc_list, cpl_propertylist **ext_plist, const char *procatg, const char *recipe)
Save a COMBINED.
cpl_table * cr2res_io_load_TRACE_WAVE(const char *filename, int detector)
Load a table from a TRACE_WAVE.
cpl_vector * cr2res_io_read_ndits(const cpl_frameset *in)
Get the NDITs from a frame set.
hdrl_image * cr2res_io_load_MASTER_FLAT(const char *filename, int detector)
Load an hdrl image from a MASTER_FLAT.
int cr2res_qc_numsat(const cpl_frameset *frameset)
Calculate the number of saturated pixels.
double cr2res_qc_obs_nodding_signal(const cpl_table *extracted)
Computes the integrated flux over part of the spectrum.
double * cr2res_qc_snr(const cpl_table *tw, const cpl_table *extracted, int **out_order_idx_values, int *out_nb_order_idx_values)
Computes the SNR of several spectra.
double cr2res_qc_obs_slit_psf(const cpl_table *slitfu, int order_idxp, int oversample)
Computes the FWHM of the PSF along the slit for a given order.
cpl_table * cr2res_trace_new_slit_fraction(const cpl_table *traces, const cpl_array *new_slit_fraction)
Recompute the traces at a newly specified slit fraction.
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.
int cr2res_order_idx_to_real(int order_idx, int order_zp)
Convert the order_idx into order_real.
double cr2res_utils_get_center_mjd(const cpl_frameset *frameset)
Calculate the middle of the exposures in the frameset.
cpl_frameset * cr2res_extract_frameset(const cpl_frameset *in, const char *tag)
Extract the frames with the given tag from a frameset.
const char * cr2res_get_license(void)
Get the pipeline copyright and license.
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 ...
void hdrl_image_delete(hdrl_image *himg)
delete hdrl_image
cpl_error_code hdrl_imagelist_collapse_mean(const hdrl_imagelist *himlist, hdrl_image **out, cpl_image **contrib)
Mean collapsing of image list.
hdrl_image * hdrl_imagelist_unset(hdrl_imagelist *himlist, cpl_size pos)
Remove an image from an imagelist.
cpl_error_code hdrl_imagelist_set(hdrl_imagelist *himlist, hdrl_image *himg, cpl_size pos)
Insert an image into an imagelist.
void hdrl_imagelist_delete(hdrl_imagelist *himlist)
Free all memory used by a hdrl_imagelist object including the images.
cpl_size hdrl_imagelist_get_size(const hdrl_imagelist *himlist)
Get the number of images in the imagelist.
hdrl_imagelist * hdrl_imagelist_new(void)
Create an empty imagelist.
hdrl_image * hdrl_imagelist_get(const hdrl_imagelist *himlist, cpl_size inum)
Get an image from a list of images.