GRAVI Pipeline Reference Manual 1.11.0
Loading...
Searching...
No Matches
gravi_wave.c
Go to the documentation of this file.
1/* $Id: gravi_calib.c,v 1.10 2012/03/23 15:10:40 nazouaoui Exp $
2 *
3 * This file is part of the GRAVI Pipeline
4 * Copyright (C) 2002,2003 European Southern Observatory
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program; if not, write to the Free Software
18 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 */
20
42
43/*-----------------------------------------------------------------------------
44 Includes
45 -----------------------------------------------------------------------------*/
46
47#ifdef HAVE_CONFIG_H
48#include <config.h>
49#endif
50
51#include <cpl.h>
52#include "cpl_wlcalib.h"
53#include "irplib_wavecal.h"
54#include <string.h>
55#include <stdio.h>
56#include <time.h>
57#include <math.h>
58#include <ctype.h>
59#include <complex.h>
60
61#include "gravi_data.h"
62#include "gravi_dfs.h"
63#include "gravi_pfits.h"
64#include "gravi_cpl.h"
65
66#include "gravi_utils.h"
67#include "gravi_ellipse.h"
68#include "gravi_signal.h"
69
70#include "gravi_preproc.h"
71#include "gravi_calib.h"
72#include "gravi_wave.h"
73#include "gravi_metrology.h"
74
75/*----------------------------------------------------------------------------
76 DEBUG
77 -----------------------------------------------------------------------------*/
78
79#define PLOT_WAVE_PHASE_VS_OPD 0
80#define WAVE_TO_PLOT 5
81
82/*-----------------------------------------------------------------------------
83 Private prototypes
84 -----------------------------------------------------------------------------*/
85
86cpl_table * gravi_opds_compute_guess (cpl_table * spectrumsc_table,
87 cpl_table * opdft_table,
88 cpl_table * met_table,
89 double dit_sc,
90 double lbd_met);
91
92cpl_table * gravi_opds_calibration (cpl_table * spectrum_table,
93 cpl_table * detector_table,
94 cpl_table * opdguess_table);
95
96cpl_error_code gravi_opds_correct_closures (cpl_table * opd_sc,
97 const char *name);
98
99cpl_vector * gravi_opds_fit_opdmet (cpl_table *opd_ft, double lbd_met);
100
101cpl_table * gravi_wave_fibre (cpl_table * spectrum_table,
102 cpl_table * detector_table,
103 cpl_table * opd_table);
104
105cpl_table * gravi_wave_fit_2d (cpl_table * wavefibre_table,
106 cpl_table * detector_table,
107 gravi_data * wave_param,
108 cpl_size fullstartx,
109 int spatial_order,
110 int spectral_order,
111 double * rms_residuals);
112
113cpl_array* my_ft_discrete(const cpl_vector *x, const cpl_vector *y, const cpl_vector *freqs);
114
115cpl_table * gravi_wave_fit_bandpass(cpl_table * spectrum_table,
116 cpl_table * detector_table,
117 cpl_table * opd_table);
118
119cpl_table * gravi_wave_fit_individual (cpl_table * wave_individual_table,
120 cpl_table * weight_individual_table,
121 cpl_table * wave_fitted_table,
122 cpl_table * opd_table,
123 cpl_table * spectrum_table,
124 cpl_table * detector_table,
125 double n0, double n1, double n2);
126
127cpl_error_code gravi_wave_correct_dispersion (cpl_table * wave_fibre,
128 double n0, double n1, double n2);
129
130cpl_imagelist * gravi_wave_test_image (cpl_table * wavedata_table,
131 cpl_table * wavefibre_table,
132 cpl_table * profile_table,
133 cpl_table * detector_table);
134
135
136
137
138/*-----------------------------------------------------------------------------
139 Functions code
140 -----------------------------------------------------------------------------*/
141
142/*----------------------------------------------------------------------------*/
153/*----------------------------------------------------------------------------*/
154
155cpl_table * gravi_compute_argon_wave (cpl_table * spectrum_table)
156{
157 gravi_msg_function_start(1);
158 cpl_ensure (spectrum_table, CPL_ERROR_NULL_INPUT, NULL);
159
160 /* Get the SC spectrums */
161 cpl_size nwave = cpl_table_get_column_depth (spectrum_table,"DATA1");
162 cpl_size nregion = gravi_spectrum_get_nregion (spectrum_table);
163 CPLCHECK_NUL ("Cannot get data");
164
165 cpl_msg_info (cpl_func, "spectrum_table = %lld",cpl_table_get_column_depth (spectrum_table, "DATA1"));
166
167 /* Prepare the outputs */
168 cpl_table * wave_data_sc = cpl_table_new (1);
169 CPLCHECK_NUL ("Cannot create output");
170
171 /* Init the list of lines*/
172 cpl_size nlines = 14;
173 double plines_str[] = {60,
174 160,
175 90,
176 1800,
177 50,
178 30,
179 30,
180 1500,
181 130,
182 510,
183 240,
184 169,
185 64,
186 73
187 };
188
189 double plines_lbd[] = {1.982291,
190 1.997118,
191 2.032256,
192 2.062186,
193 2.065277,
194 2.073922,
195 2.081672,
196 2.099184,
197 2.133871,
198 2.154009,
199 2.208321,
200 2.313952,
201 2.385154,
202 2.397306
203 };
204
205 cpl_vector * lines_lbd = cpl_vector_wrap (nlines,plines_lbd);
206 cpl_vector * lines_str = cpl_vector_wrap (nlines,plines_str);
207 cpl_bivector * lines = cpl_bivector_wrap_vectors (lines_lbd, lines_str);
208
209 /* Init the slit and lamp model */
210 cpl_wlcalib_slitmodel * model = cpl_wlcalib_slitmodel_new ();
211 cpl_wlcalib_slitmodel_set_wslit (model, 0.1);
212 cpl_wlcalib_slitmodel_set_wfwhm (model, 1.25);
213 cpl_wlcalib_slitmodel_set_threshold (model, 5.0);
214 cpl_wlcalib_slitmodel_set_catalog (model, lines);
215
216 /* Starting point of dispersion */
217 double values_hr[] = {1.97048,0.2228e-3,2.7728e-08,-4.72e-12};
218 double values_mr[] = {1.9697,0.2179e-2,2.5e-07,-5e-10};
219 double * values = nwave > 1000 ? values_hr : values_mr;
220
221 /* Fill the dispersion polynomial */
222 cpl_polynomial * dispersion0 = cpl_polynomial_new (1);
223 for (cpl_size order = 0 ; order < 4 ; order ++) {
224 cpl_polynomial_set_coeff (dispersion0, &order, values[order]);
225 }
226
227 /* QC parameters */
228 double minwave = -1e10;
229 double maxwave = +1e10;
230
231 cpl_vector * spectrum = cpl_vector_new (nwave);
232 cpl_vector * input_spectrum = cpl_vector_new (nwave);
233 cpl_vector * model_spectrum = cpl_vector_new (nwave);
234
235 /* Prepare output table */
236 cpl_table * fit_table = cpl_table_new (nregion);
237 gravi_table_new_column_array (fit_table, "DATA", "adu", CPL_TYPE_DOUBLE, nwave);
238 gravi_table_new_column_array (fit_table, "DATA_MODEL", "adu", CPL_TYPE_DOUBLE, nwave);
239 gravi_table_new_column_array (fit_table, "DATA_WAVE", "adu", CPL_TYPE_DOUBLE, nwave);
240
241 /* The starting point will be the fit of previous region */
242 cpl_polynomial * dispersion = cpl_polynomial_duplicate (dispersion0);
243
244 /* Loop on regions */
245 for (cpl_size reg = 0; reg < nregion ; reg ++ ) {
246 cpl_msg_info (cpl_func,"Fit region %lld over %lld", reg+1, nregion);
247
248 /* Copy the spectrum of this region into a vector */
249 const cpl_array * data = cpl_table_get_array (spectrum_table, GRAVI_DATA[reg], 0);
250 for (cpl_size wave=0; wave < nwave; wave ++) {
251 cpl_vector_set (spectrum, wave, log (CPL_MAX (cpl_array_get (data, wave, NULL), 1.0)));
252 cpl_vector_set (input_spectrum, wave, cpl_array_get (data, wave, NULL));
253 }
254
255 CPLCHECK_NUL ("Cannot prepare data");
256
257 /* Parameters: Explore a large hsize [pixel]
258 * only for the first region, since we keep
259 * the starting point */
260 int maxdeg = 3, nmaxima = 0, linelim = 99;
261 int maxite = 100 * maxdeg, maxfail = 10, maxcont = 10;
262 int hsize = (reg == 0 ? 40 : 5);
263 double pixtol = 1e-6, pixstep = 0.5, pxc = 0.0;
264
265 /* Use the log of spectrum to put
266 * less weight on the intensity */
267 irplib_polynomial_find_1d_from_correlation_all (dispersion, maxdeg, spectrum,
268 nmaxima, linelim,
269 (irplib_base_spectrum_model*)model,
270 &irplib_vector_fill_logline_spectrum,
271 pixtol, pixstep,
272 hsize, maxite, maxfail, maxcont,
273 CPL_FALSE,&pxc);
274 CPLCHECK_NUL ("Cannot fit polynomial");
275
276 /* Fill in the final fit */
277 irplib_vector_fill_line_spectrum (model_spectrum, dispersion, (void *)model);
278 cpl_polynomial_dump (dispersion, stdout);
279
280 /* Evaluate the polynomial for the output */
281 cpl_array * wave_data = cpl_array_new (nwave, CPL_TYPE_DOUBLE);
282 for (cpl_size pix = 0; pix < nwave ; pix++)
283 cpl_array_set (wave_data, pix, cpl_polynomial_eval_1d (dispersion, (double)pix, NULL) * 1e-6);
284
285 /* Check the limits for the QC */
286 minwave = CPL_MAX (cpl_array_get_min (wave_data), minwave);
287 maxwave = CPL_MIN (cpl_array_get_max (wave_data), maxwave);
288
289 /* Fill the WAVE_DATA_SC */
290 cpl_table_new_column_array (wave_data_sc, GRAVI_DATA[reg], CPL_TYPE_DOUBLE, nwave);
291 cpl_table_set_array (wave_data_sc, GRAVI_DATA[reg], 0, wave_data);
292 FREE (cpl_array_delete, wave_data);
293
294 /* Save the test tables */
295 cpl_array * tmp_array;
296 tmp_array = cpl_array_wrap_double (cpl_vector_get_data (input_spectrum), nwave);
297 cpl_table_set_array (fit_table, "DATA", reg, tmp_array);
298 cpl_array_unwrap (tmp_array);
299 tmp_array = cpl_array_wrap_double (cpl_vector_get_data (model_spectrum), nwave);
300 cpl_table_set_array (fit_table, "DATA_MODEL", reg, tmp_array);
301 cpl_array_unwrap (tmp_array);
302 tmp_array = cpl_array_new (nwave, CPL_TYPE_DOUBLE);
303 for (cpl_size wave = 0; wave < nwave; wave ++)
304 cpl_array_set (tmp_array, wave, cpl_polynomial_eval_1d (dispersion, (double)wave, NULL) * 1e-6);
305 cpl_table_set_array (fit_table, "DATA_WAVE", reg, tmp_array);
306 cpl_array_delete (tmp_array);
307
308 CPLCHECK_NUL ("Cannot set data");
309 } /* End loop on regions */
310
311 // cpl_table_save (fit_table, NULL, NULL, "fit_wave.fits", CPL_IO_CREATE);
312 FREE (cpl_table_delete, fit_table);
313
314 /* Delete */
315 FREE (cpl_vector_delete, spectrum);
316 FREE (cpl_vector_delete, input_spectrum);
317 FREE (cpl_vector_delete, model_spectrum);
318 FREE (cpl_polynomial_delete, dispersion0);
319 FREE (cpl_polynomial_delete, dispersion);
320
321 /* This memory desalocation does not work: */
322 // FREE (cpl_wlcalib_slitmodel_delete, model);
323
324 gravi_msg_function_exit(1);
325 return wave_data_sc;
326}
327
328/*----------------------------------------------------------------------------*/
358/*----------------------------------------------------------------------------*/
359
360cpl_error_code gravi_opds_correct_closures (cpl_table * phase_table,
361 const char *name)
362{
363 gravi_msg_function_start(1);
364 cpl_ensure_code (phase_table, CPL_ERROR_NULL_INPUT);
365 cpl_ensure_code (name, CPL_ERROR_NULL_INPUT);
366
367 int nclo = 4;
368 cpl_size nrow = cpl_table_get_nrow (phase_table) / GRAVI_NBASE;
369
370 /* Get the data */
371 double * phase = cpl_table_get_data_double (phase_table, name);
372 CPLCHECK_MSG ("Cannot get data");
373
374 /* Model and right-hand-side for the lineary system
375 * (unfilled matrix are 0.0) */
376 cpl_matrix * rhs_matrix = cpl_matrix_new (nrow * nclo, 1);
377 cpl_matrix * model_matrix = cpl_matrix_new (nrow * nclo, GRAVI_NBASE-1 + nclo);
378
379 for (cpl_size row = 0; row < nrow; row++) {
380 for (int clo = 0; clo < nclo; clo++) {
381 int b0 = GRAVI_CLO_BASE[clo][0];
382 int b1 = GRAVI_CLO_BASE[clo][1];
383 int b2 = GRAVI_CLO_BASE[clo][2];
384
385 /* Fill the rhs with measured closure */
386 double closure = phase[row*GRAVI_NBASE+b0] +
387 phase[row*GRAVI_NBASE+b1] -
388 phase[row*GRAVI_NBASE+b2];
389 cpl_matrix_set (rhs_matrix, row*nclo+clo, 0, closure);
390 CPLCHECK_MSG ("Cannot fill rhs_matrix");
391
392 /* Fill the k (unfilled are 0.0) */
393 cpl_matrix_set (model_matrix, row*nclo+clo, clo, 1.0);
394 CPLCHECK_MSG ("Cannot fill k in model_matrix");
395
396 /* Fill the f in of the model_matrix */
397 if (b0 < GRAVI_NBASE-1) cpl_matrix_set (model_matrix, row*nclo+clo, nclo+b0,
398 +phase[row*GRAVI_NBASE+b0]);
399 if (b1 < GRAVI_NBASE-1) cpl_matrix_set (model_matrix, row*nclo+clo, nclo+b1,
400 +phase[row*GRAVI_NBASE+b1]);
401 if (b2 < GRAVI_NBASE-1) cpl_matrix_set (model_matrix, row*nclo+clo, nclo+b2,
402 -phase[row*GRAVI_NBASE+b2]);
403 CPLCHECK_MSG ("Cannot fill phase in model_matrix");
404 }
405 } /* End loop on clo and rows */
406
407 /* Solve the system */
408 cpl_matrix * res_matrix = cpl_matrix_solve_normal (model_matrix, rhs_matrix);
409 CPLCHECK_MSG ("Cannot solve system");
410
411 /* Compute residuals */
412 cpl_matrix * residual_matrix = cpl_matrix_product_create (model_matrix, res_matrix);
413 cpl_matrix_subtract (residual_matrix, rhs_matrix);
414 double rms_fit = cpl_matrix_get_stdev (residual_matrix);
415
416 /* Correct baseline with (1-f). Last baseline is kept unchanged */
417 for (int base = 0; base < GRAVI_NBASE - 1; base ++) {
418 double f = cpl_matrix_get (res_matrix, nclo + base, 0);
419 cpl_msg_info (cpl_func,"correction factor f = 1 %+.20f", -1*f);
420 for (cpl_size row = 0; row < nrow; row++) phase[row*GRAVI_NBASE+base] *= 1 - f;
421 }
422
423 /* Dump the residual with their units */
424 const char * unit = cpl_table_get_column_unit (phase_table, name);
425 cpl_msg_info (cpl_func, "residual stdev = %.20g [%s]", rms_fit, unit);
426
427 /* Delete matrix */
428 FREE (cpl_matrix_delete, residual_matrix);
429 FREE (cpl_matrix_delete, res_matrix);
430 FREE (cpl_matrix_delete, model_matrix);
431 FREE (cpl_matrix_delete, rhs_matrix);
432
433 gravi_msg_function_exit(1);
434 return CPL_ERROR_NONE;
435}
436
437/*----------------------------------------------------------------------------*/
463/*----------------------------------------------------------------------------*/
464
465cpl_vector * gravi_opds_fit_opdmet (cpl_table * ft_table, double lbd_met)
466{
467 gravi_msg_function_start(1);
468 cpl_ensure (ft_table, CPL_ERROR_NULL_INPUT, NULL);
469
470
471 /* Get the number of acquisitions */
472 cpl_size nrow = cpl_table_get_nrow (ft_table) / GRAVI_NBASE;
473
474 /* Get the pointer to data */
475 double * opd_sc = cpl_table_get_data_double (ft_table, "OPD_SC");
476 double * opd_ft = cpl_table_get_data_double (ft_table, "OPD");
477 double * phase_met = cpl_table_get_data_double (ft_table, "PHASE_MET_FC");
478 CPLCHECK_NUL ("Cannot get data");
479
480 /* Number of valid rows */
481 cpl_size nrow_valid = 0;
482 for (cpl_size row = 0; row < nrow; row++) if (opd_sc[row*GRAVI_NBASE] != 0) nrow_valid++;
483 cpl_msg_info (cpl_func,"nrow_valid = %lld", nrow_valid);
484
485 cpl_ensure (nrow_valid > 100, CPL_ERROR_ILLEGAL_INPUT, NULL);
486
487 /* Model and right-hand-side for the lineary system
488 * (unfilled matrix are 0.0) */
489 cpl_matrix * rhs_matrix = cpl_matrix_new (nrow_valid*GRAVI_NBASE, 1);
490 cpl_matrix * model_matrix = cpl_matrix_new (nrow_valid*GRAVI_NBASE, GRAVI_NBASE + 2);
491
492 for (int base = 0; base < GRAVI_NBASE; base++) {
493 cpl_size row_valid = 0;
494
495 for (cpl_size row=0; row<nrow; row++) {
496 if (opd_sc[row*GRAVI_NBASE] == 0) continue;
497
498 int idv = row_valid * GRAVI_NBASE + base;
499 int id = row * GRAVI_NBASE + base;
500
501 /* Fill the OPD metrology */
502 cpl_matrix_set (rhs_matrix, idv, 0, phase_met[id] * lbd_met / CPL_MATH_2PI);
503 CPLCHECK_NUL ("Cannot set OPD");
504
505 /* Fill the model b and c */
506 cpl_matrix_set (model_matrix, idv, 0, 1*opd_sc[id]);
507 cpl_matrix_set (model_matrix, idv, 1, -1*opd_ft[id]);
508 CPLCHECK_NUL ("Cannot set SC or FT");
509
510 /* Fill the model Aij (unfilled matrix are 0.0) */
511 cpl_matrix_set (model_matrix, idv, 2 + base, 1.0);
512 CPLCHECK_NUL ("Cannot set the zero points");
513
514 row_valid++;
515 } /* End loop rows */
516 } /* End loop on base*/
517
518 /* Solve the system */
519 cpl_matrix * res_matrix = cpl_matrix_solve_normal (model_matrix, rhs_matrix);
520
521 /* Compute residuals */
522 cpl_matrix * residual_matrix = cpl_matrix_product_create (model_matrix, res_matrix);
523 cpl_matrix_subtract (residual_matrix, rhs_matrix);
524 double rms_fit = cpl_matrix_get_stdev (residual_matrix);
525 //cpl_plot_vector(NULL, NULL, NULL, cpl_vector_wrap(nrow_valid*GRAVI_NBASE,cpl_matrix_get_data(residual_matrix)));
526
527 /* Verbose */
528 cpl_msg_info (cpl_func, "coeff SC = %.20g ", cpl_matrix_get (res_matrix, 0, 0));
529 cpl_msg_info (cpl_func, "coeff FT = %.20g ", cpl_matrix_get (res_matrix, 1, 0));
530 cpl_msg_info (cpl_func, "residual stdev = %.20g [m]", rms_fit);
531
532 /* Set the Results */
533 cpl_vector * opd_coeff = cpl_vector_new(3);
534 cpl_vector_set (opd_coeff, GRAVI_SC, cpl_matrix_get (res_matrix, 0, 0));
535 cpl_vector_set (opd_coeff, GRAVI_FT, cpl_matrix_get (res_matrix, 1, 0));
536 cpl_vector_set (opd_coeff, 2, rms_fit);
537
538 /* Delete matrix */
539 FREE (cpl_matrix_delete, residual_matrix);
540 FREE (cpl_matrix_delete, res_matrix);
541 FREE (cpl_matrix_delete, model_matrix);
542 FREE (cpl_matrix_delete, rhs_matrix);
543
544 gravi_msg_function_exit(1);
545 return opd_coeff;
546}
547
548/*----------------------------------------------------------------------------*/
570/*----------------------------------------------------------------------------*/
571
572cpl_table * gravi_opds_compute_guess (cpl_table * spectrumsc_table,
573 cpl_table * ft_table,
574 cpl_table * vismet_table,
575 double dit_sc,
576 double lbd_met)
577{
578 gravi_msg_function_start(1);
579 cpl_ensure (spectrumsc_table, CPL_ERROR_NULL_INPUT, NULL);
580 cpl_ensure (ft_table, CPL_ERROR_NULL_INPUT, NULL);
581 cpl_ensure (vismet_table, CPL_ERROR_NULL_INPUT, NULL);
582 cpl_ensure (dit_sc>0, CPL_ERROR_ILLEGAL_INPUT, NULL);
583
584 int ntel = 4;
585
586 cpl_size nrow = cpl_table_get_nrow (spectrumsc_table);
587 cpl_size nrow_met = cpl_table_get_nrow (vismet_table) / ntel;
588 cpl_size nrow_ft = cpl_table_get_nrow (ft_table) / GRAVI_NBASE;
589 int * time_SC = cpl_table_get_data_int (spectrumsc_table, "TIME");
590 CPLCHECK_NUL ("Cannot get data");
591
592 /* Pointer on MET data */
593 int * time_MET = cpl_table_get_data_int (vismet_table, "TIME");
594 double * phase_MET = cpl_table_get_data_double (vismet_table, "PHASE_FC_DRS");
595 CPLCHECK_NUL ("Cannot get MET data");
596
597 /* Pointer on OPD data */
598 int * time_FT = cpl_table_get_data_int (ft_table, "TIME");
599 double * opd_FT = cpl_table_get_data_double (ft_table, "OPD");
600 CPLCHECK_NUL ("Cannot get FT data");
601
602 /* Create table */
603 cpl_table * guess_table = cpl_table_new (nrow * GRAVI_NBASE);
604 gravi_table_new_column (guess_table, "OPD", "m", CPL_TYPE_DOUBLE);
605
606 /* Loop on base */
607 for (int base = 0; base < GRAVI_NBASE; base++) {
608 int tel0 = GRAVI_BASE_TEL[base][0];
609 int tel1 = GRAVI_BASE_TEL[base][1];
610
611 /* Loop on the SC frames */
612 cpl_size row_met = 0, row_ft = 0;
613 for (cpl_size row = 0 ; row < nrow ; row++) {
614
615 /*
616 * Average the MET
617 */
618 int counter_met = 0;
619 double opd_met = 0.0;
620 while (time_MET[row_met*ntel] < (time_SC[row] + dit_sc/2.)) {
621 if ((time_MET[row_met*ntel] > (time_SC[row] - dit_sc/2.)) && (row_met < nrow_met)) {
622 opd_met += phase_MET[row_met*ntel+tel1] - phase_MET[row_met*ntel+tel0];
623 counter_met ++;
624 }
625
626 /* If not enough data to synchronize */
627 if (row_met > nrow_met - 2) {
628 cpl_msg_warning (cpl_func,"Not enough data to synchronise the MET with SC");
629 break;
630 }
631 row_met ++;
632 }
633 /* End sum the MET over the current SC DIT */
634
635 opd_met = opd_met * lbd_met / CPL_MATH_2PI / counter_met; // [m]
636
637 /*
638 * Average the FT
639 */
640 int counter_ft = 0;
641 double opd_ft = 0.0;
642 while (time_FT[row_ft*GRAVI_NBASE+base] < (time_SC[row] + dit_sc/2.)) {
643 if ((time_FT[row_ft*GRAVI_NBASE+base] > (time_SC[row] - dit_sc/2.)) && (row_ft < nrow_ft)) {
644 opd_ft += opd_FT[row_ft*GRAVI_NBASE+base];
645 counter_ft ++;
646 }
647
648 /* If not enough data to synchronize */
649 if (row_ft > nrow_ft - 2) {
650 cpl_msg_warning (cpl_func,"Not enough data to synchronise the FT with SC");
651 break;
652 }
653 row_ft ++;
654 }
655 /* End sum the FT over the current SC DIT */
656
657 opd_ft = opd_ft / counter_ft;
658
659 /* Set the total guess OPD as OPD_FT - OPD_MET */
660 cpl_table_set (guess_table, "OPD", row*GRAVI_NBASE+base, opd_ft - opd_met);
661 CPLCHECK_NUL ("Cannot compute opd guess");
662 } /* End loop on row */
663 } /* End loop on base */
664
665 /*
666 * Remove the mean of each base
667 */
668 cpl_msg_info (cpl_func, "Remove the mean from opdguess table");
669
670 for (int base = 0; base < GRAVI_NBASE; base++) {
671 double mean = gravi_table_get_column_mean (guess_table, "OPD", base, GRAVI_NBASE);
672 gravi_table_add_scalar (guess_table, "OPD", base, GRAVI_NBASE, -1*mean);
673 }
674
675 gravi_msg_function_exit(1);
676 return guess_table;
677}
678
679/*----------------------------------------------------------------------------*/
702/*----------------------------------------------------------------------------*/
703
704cpl_table * gravi_opds_calibration (cpl_table * spectrum_table,
705 cpl_table * detector_table,
706 cpl_table * guess_table)
707{
708 /* Verbose */
709 gravi_msg_function_start(1);
710 cpl_ensure (spectrum_table, CPL_ERROR_NULL_INPUT, NULL);
711 cpl_ensure (detector_table, CPL_ERROR_NULL_INPUT, NULL);
712
713 /* Verbose the phase algorithm selected */
714 if (USE_LINEAR_ELLIPSE) {
715 cpl_msg_info (cpl_func, "Solve linearly the ellipse X,Y");
716 }
717 else {
718 cpl_msg_info (cpl_func, "Solve non-linearly the ellipse X,Y");
719 }
720
721 /* Get the size of the vectors */
722 cpl_size nrow = cpl_table_get_nrow (spectrum_table);
723 cpl_size nwave = gravi_spectrum_get_nwave (spectrum_table);
724 cpl_size npol = gravi_spectrum_get_npol (spectrum_table);
725 CPLCHECK_NUL ("Cannot get size");
726
727 /* Create a fake OI_WAVE table, to average the different
728 * spectral channels [m] */
729
730 cpl_table ** oiwave_tables = cpl_calloc (npol, sizeof (cpl_table *));
731 for (int pol = 0; pol < npol; pol++) {
732 oiwave_tables[pol] = cpl_table_new (nwave);
733 cpl_table_new_column (oiwave_tables[pol], "EFF_WAVE", CPL_TYPE_DOUBLE);
734 for (cpl_size wave = 0; wave < nwave; wave++) {
735 double value = 1.97e-6 + (2.48e-6 - 1.97e-6) / nwave * wave;
736 cpl_table_set (oiwave_tables[pol], "EFF_WAVE", wave, value);
737 }
738 }
739
740 /* Create column in phase table */
741 cpl_table * output_table = cpl_table_new (GRAVI_NBASE * nrow);
742 gravi_table_new_column (output_table,"OPD", "m", CPL_TYPE_DOUBLE);
743 gravi_table_new_column (output_table,"TIME", "us", CPL_TYPE_INT);
744
745 /* Set the time */
746 for (cpl_size row = 0; row < nrow; row++) {
747 double value = cpl_table_get (spectrum_table, "TIME", row, NULL);
748 for (int base = 0; base < GRAVI_NBASE; base++)
749 cpl_table_set (output_table, "TIME", row*GRAVI_NBASE+base, value);
750 }
751
752 /* Loop on base */
753 for (int base = 0; base < GRAVI_NBASE; base ++) {
754
755 /* Check if a guess exists */
756 cpl_vector * opd_guess = NULL;
757 if (guess_table) {
758 opd_guess = cpl_vector_new (nrow);
759 for (cpl_size row = 0; row < nrow; row++) {
760 double value = cpl_table_get (guess_table, "OPD", row*GRAVI_NBASE+base, NULL);
761 cpl_vector_set (opd_guess, row, value);
762 }
763 }
764
765 /* Recover the mean OPD modulation (averaved over pol
766 * and channels) using the ellipses. In [m]. However since oiwave
767 * is not known... this is to a scaling factor correction. */
768 cpl_vector * mean_opd;
769 mean_opd = gravi_ellipse_meanopd_create (spectrum_table, detector_table,
770 oiwave_tables, opd_guess, base);
771 CPLCHECK_NUL ("Cannot compute opd");
772
773 /* Save the mean opd [m] */
774 for (cpl_size row = 0; row < nrow; row++) {
775 double value = cpl_vector_get (mean_opd, row);
776 cpl_table_set (output_table, "OPD", row*GRAVI_NBASE+base, value);
777 }
778 CPLCHECK_NUL ("Cannot set phase");
779
780 FREE (cpl_vector_delete, mean_opd);
781 FREE (cpl_vector_delete, opd_guess);
782
783 } /* End loop on base */
784
785 FREELOOP (cpl_table_delete, oiwave_tables, npol);
786
787 /* Verbose */
788 gravi_msg_function_exit(1);
789 return output_table;
790}
791
792/*----------------------------------------------------------------------------*/
824/*----------------------------------------------------------------------------*/
825
826cpl_error_code gravi_wave_compute_opds (gravi_data * spectrum_data,
827 cpl_table * met_table,
828 const cpl_parameterlist * parlist)
829{
830 gravi_msg_function_start(1);
831 cpl_ensure_code (spectrum_data, CPL_ERROR_NULL_INPUT);
832 cpl_ensure_code (met_table, CPL_ERROR_NULL_INPUT);
833
834 /* Get DIT */
835 cpl_propertylist * spectrum_header = gravi_data_get_header (spectrum_data);
836 double dit_sc = gravi_pfits_get_dit_sc (spectrum_header)*1e6; // [us]
837 CPLCHECK_MSG ("Cannot get DIT");
838
839 /* Get the input */
840 cpl_table * spectrumft_table = gravi_data_get_spectrum_data (spectrum_data, GRAVI_FT);
841 cpl_table * detectorft_table = gravi_data_get_imaging_detector (spectrum_data, GRAVI_FT);
842 cpl_table * spectrumsc_table = gravi_data_get_spectrum_data (spectrum_data, GRAVI_SC);
843 cpl_table * detectorsc_table = gravi_data_get_imaging_detector (spectrum_data, GRAVI_SC);
844 CPLCHECK_MSG ("Cannot get data");
845
846
847 /*
848 * Reduce the raw METROLOGY into VIS_MET
849 */
850 cpl_msg_info (cpl_func, "Compute the phase of MET_FC");
851
852 cpl_table * vismet_table = gravi_metrology_create (met_table, spectrum_header);
853 gravi_metrology_drs (met_table, vismet_table, spectrum_header, parlist);
854 /* also apply the TAC reduction */
855 gravi_metrology_tac(met_table, vismet_table, spectrum_header);
856
857 /* Compute the mean LBD_MET for this file */
858 double lbd_met = gravi_pfits_get_met_wavelength_mean (spectrum_header, met_table);
859
860 /*
861 * Compute the phase of SC and FT. For the SC, we use a guess
862 * of the opd modulation, computed from FT and MET, to unwrap.
863 */
864 cpl_msg_info (cpl_func, "Compute OPD of FT from ellipse");
865
866 cpl_table * ft_table;
867 ft_table = gravi_opds_calibration (spectrumft_table, detectorft_table, NULL);
868 gravi_opds_correct_closures (ft_table, "OPD");
869
870 cpl_msg_info (cpl_func, "Compute OPD of SC from ellipse");
871
872 cpl_table * guess_table;
873 guess_table = gravi_opds_compute_guess (spectrumsc_table, ft_table, vismet_table, dit_sc, lbd_met);
874
875 cpl_table * sc_table;
876 sc_table = gravi_opds_calibration (spectrumsc_table, detectorsc_table, guess_table);
877 gravi_opds_correct_closures (sc_table, "OPD");
878 FREE (cpl_table_delete, guess_table);
879
880 CPLCHECK_MSG ("Cannot calibrate phase");
881
882 /*
883 * Interpolate MET and SC at the sampling frequency of FT
884 * Results are saved in the ft_table table
885 */
886 cpl_msg_info (cpl_func, "Fit MET = a.SC - b.FT + c to get absolute modulation");
887
888 gravi_vis_create_met_ft (ft_table, vismet_table);
889 gravi_vis_create_opdsc_ft (ft_table, sc_table, dit_sc);
890
891 CPLCHECK_MSG ("Cannot resample SC or MET at the FT frequency");
892
893 /*
894 * Compute the scaling coefficients of OPDs by fitting:
895 * OPD_MET_ijt = a.OPD_SC_ijt - b.OPD_FT_ijt + c_ij
896 */
897 cpl_vector * coeff_vector = gravi_opds_fit_opdmet (ft_table, lbd_met);
898
899 CPLCHECK_MSG ("Cannot fit opdmet");
900
901 /* Set the CHI2 of the fit in the QC parameters */
902 cpl_propertylist_update_double (spectrum_header, QC_PHASECHI2,
903 cpl_vector_get (coeff_vector, 2));
904 cpl_propertylist_set_comment (spectrum_header, QC_PHASECHI2,
905 "chi2 of a.SC-b.FT+c=MET");
906
907 /* Add the OPD COEFF in QC parameters */
908 for (int type_data = 0; type_data < 2; type_data++) {
909 double tmp = cpl_vector_get (coeff_vector, type_data);
910 cpl_propertylist_update_float (spectrum_header, OPD_COEFF_SIGN(type_data), tmp);
911 cpl_propertylist_set_comment (spectrum_header, OPD_COEFF_SIGN(type_data), "wavelength correction");
912 }
913
914 /* Set a warning */
915 if (cpl_vector_get (coeff_vector, 2) > 1e-7) {
916 gravi_pfits_add_check (spectrum_header,"Residuals of fit MET=a.SC-b.FT+c are high");
917 }
918
919 if (cpl_vector_get (coeff_vector, 2) > 1.2e-7) {
920 cpl_msg_warning (cpl_func, "*************************************************");
921 cpl_msg_warning (cpl_func, "**** !!! residuals of the fit too high !!! ****");
922 cpl_msg_warning (cpl_func, "**** Residuals are:%7.0f nm ****",cpl_vector_get (coeff_vector, 2)*1e9);
923 cpl_msg_warning (cpl_func, "**** SC and RMN may be desynchronized ****");
924 cpl_msg_warning (cpl_func, "**** (or out of the envelope in LOW) ****");
925 cpl_msg_warning (cpl_func, "*************************************************");
926 }
927
928 CPLCHECK_MSG ("Cannot set QC parameter");
929
930 /*
931 * Correct opd from overall scaling by fitting MET
932 */
933 double coeff_sc = cpl_vector_get (coeff_vector, GRAVI_SC);
934 cpl_table_multiply_scalar (sc_table, "OPD", coeff_sc);
935
936 double coeff_ft = cpl_vector_get (coeff_vector, GRAVI_FT);
937 cpl_table_multiply_scalar (ft_table, "OPD", coeff_ft);
938
939 FREE (cpl_vector_delete, coeff_vector);
940 CPLCHECK_MSG ("Cannot correct OPDs from scaling coefficients");
941
942 /*
943 * Fill the output
944 */
945// if ((det_type == GRAVI_DET_SC || det_type == GRAVI_DET_ALL))
946// {
947 gravi_data_add_table (spectrum_data, NULL, "OPD_SC", sc_table);
948 CPLCHECK_MSG ("Cannot set OPD_SC table");
949// }
950// else
951// FREE (cpl_table_delete, sc_table);
952//
953//
954// if ((det_type == GRAVI_DET_FT || det_type == GRAVI_DET_ALL))
955// {
956 gravi_data_add_table (spectrum_data, NULL, "OPD_FT", ft_table);
957 CPLCHECK_MSG ("Cannot set OPD_FT table");
958// }
959// else
960// FREE (cpl_table_delete, ft_table);
961
962 gravi_data_add_table (spectrum_data, NULL, GRAVI_OI_VIS_MET_EXT, vismet_table);
963 CPLCHECK_MSG ("Cannot set OI_VIS_MET table");
964
965 gravi_msg_function_exit(1);
966 return CPL_ERROR_NONE;
967}
968
969/*----------------------------------------------------------------------------*/
993/*----------------------------------------------------------------------------*/
994
995cpl_table * gravi_wave_fibre (cpl_table * spectrum_table,
996 cpl_table * detector_table,
997 cpl_table * opd_table)
998{
999 gravi_msg_function_start(1);
1000 cpl_ensure (spectrum_table, CPL_ERROR_NULL_INPUT, NULL);
1001 cpl_ensure (detector_table, CPL_ERROR_NULL_INPUT, NULL);
1002 cpl_ensure (opd_table, CPL_ERROR_NULL_INPUT, NULL);
1003
1004 char name[100];
1005
1006 /* Create the output table */
1007 cpl_table * wave_fibre = cpl_table_new (1);
1008
1009 /* Get the number of wavelength, region, polarisation... */
1010 cpl_size nwave = cpl_table_get_column_depth (spectrum_table, "DATA1");
1011 cpl_size n_region = cpl_table_get_nrow (detector_table);
1012 cpl_size nrow = cpl_table_get_nrow (spectrum_table);
1013
1014 int npol = (n_region > 24 ? 2 : 1);
1015
1016 /*
1017 * Calibration of each polarization and base
1018 */
1019
1020 for (int pol = 0; pol < npol; pol++) {
1021 for (int base = 0; base < GRAVI_NBASE; base ++) {
1022 cpl_msg_info (cpl_func, "Compute wave fibre for pol %i over %i, base %i over %i",
1023 pol+1, npol, base+1, GRAVI_NBASE);
1024
1025 /* Get the index of the ABCD. */
1026 int iA = gravi_get_region (detector_table, base, 'A', pol);
1027 int iB = gravi_get_region (detector_table, base, 'B', pol);
1028 int iC = gravi_get_region (detector_table, base, 'C', pol);
1029 int iD = gravi_get_region (detector_table, base, 'D', pol);
1030 if (iA<0 || iB<0 || iC<0 || iD<0){
1031 cpl_msg_warning (cpl_func, "Don't found the A, B, C or D !!!");
1032 }
1033
1034 /* Sign of this baseline */
1035 double phi_sign = gravi_region_get_base_sign (detector_table, base);
1036
1037 /* Get the OPD into various flarous */
1038 cpl_matrix * opd_matrix = cpl_matrix_new (1, nrow);
1039 cpl_vector * opd_vector = cpl_vector_new (nrow);
1040 for (cpl_size row = 0; row < nrow; row ++ ) {
1041 double value = cpl_table_get (opd_table, "OPD", row*GRAVI_NBASE+base, NULL);
1042 cpl_matrix_set (opd_matrix, 0, row, value);
1043 cpl_vector_set (opd_vector, row, value);
1044 }
1045 CPLCHECK_NUL ("Cannot extract the OPD");
1046
1047 /* Create array to fill the wavelenghts */
1048 cpl_array * wavelength = cpl_array_new (nwave, CPL_TYPE_DOUBLE);
1049 cpl_array * wavechi2 = cpl_array_new (nwave, CPL_TYPE_DOUBLE);
1050 cpl_array_fill_window (wavelength, 0, nwave, 0.0);
1051 cpl_array_fill_window (wavechi2, 0, nwave, 1e10);
1052
1053 /*
1054 * Loop on spectral channel
1055 */
1056 for (cpl_size wave = 0; wave < nwave; wave++) {
1057
1058 cpl_vector * vector_T = NULL, * vector_X, * vector_Y;
1059
1060 /* Define the vector_X = C - A */
1061 vector_X = gravi_table_get_vector (spectrum_table, wave, GRAVI_DATA[iC]);
1062 vector_T = gravi_table_get_vector (spectrum_table, wave, GRAVI_DATA[iA]);
1063 cpl_vector_subtract (vector_X, vector_T);
1064 FREE (cpl_vector_delete, vector_T);
1065
1066 /* Define the vector_Y = D - B */
1067 vector_Y = gravi_table_get_vector (spectrum_table, wave, GRAVI_DATA[iD]);
1068 vector_T = gravi_table_get_vector (spectrum_table, wave, GRAVI_DATA[iB]);
1069 cpl_vector_subtract (vector_Y, vector_T);
1070 FREE (cpl_vector_delete, vector_T);
1071
1072 CPLCHECK_NUL ("Cannot compute vector_X or Y");
1073
1074 /* Compute envelope from OPD for this channel */
1075 cpl_vector * envelope_vector = gravi_compute_envelope (opd_vector, wave, nwave);
1076
1077 /* Compute the phase from the ellipse */
1078 cpl_vector * phase;
1079 phase = gravi_ellipse_phase_create (vector_X, vector_Y,
1080 envelope_vector);
1081 FREE (cpl_vector_delete, vector_X);
1082 FREE (cpl_vector_delete, vector_Y);
1083 FREE (cpl_vector_delete, envelope_vector);
1084
1085 /* If the computation of the ellipse fails, we continue with next wave */
1086 if (phase == NULL) {
1087 cpl_msg_warning (cpl_func, "Cannot compute wave for channel %lld and base %d", wave, base);
1088 continue;
1089 }
1090
1091 /* Multiply by the sign */
1092 cpl_vector_multiply_scalar (phase, phi_sign);
1093
1094 /* Unwrap phase from the OPD */
1095 double lbd_channel = 1.95e-6 + (2.46e-6 - 1.95e-6) / nwave * wave;
1096 gravi_vector_unwrap_with_guess (phase, opd_vector, CPL_MATH_2PI / lbd_channel);
1097
1098 /* Fit the slope of the phase versus OPD gives the
1099 * wavelength of the spectral element */
1100 const cpl_size mindeg = 0, maxdeg = 1;
1101 cpl_polynomial * fit_slope = cpl_polynomial_new (1);
1102 cpl_polynomial_fit (fit_slope, opd_matrix, NULL, phase, NULL, CPL_FALSE, &mindeg, &maxdeg);
1103
1104 /* Compute residuals */
1105 cpl_vector * residuals = cpl_vector_new (nwave);
1106 double rechisq;
1107 cpl_vector_fill_polynomial_fit_residual (residuals, phase, NULL, fit_slope, opd_matrix, &rechisq);
1108
1109 if(PLOT_WAVE_PHASE_VS_OPD && (wave == WAVE_TO_PLOT))
1110 {
1111 char gnuplot_str[200];
1112 sprintf (gnuplot_str, "set title 'Wavelength (base %d)'; set xlabel 'Phase [rad]'; set ylabel 'OPD (m)';", base);
1113 cpl_plot_vector (gnuplot_str, NULL, NULL, opd_vector);
1114 sprintf (gnuplot_str, "set title 'Wavelength residuals (base %d)'; set xlabel 'Phase [rad]'; set ylabel 'OPD (m)';", base);
1115 cpl_plot_vector (gnuplot_str, NULL, NULL, residuals);
1116 CPLCHECK_NUL ("Cannot plot OPD versus phase");
1117 }
1118
1119 CPLCHECK_NUL ("Cannot fit the OPD and phase");
1120
1121 /* Get the slope */
1122 const cpl_size pow_slope = 1;
1123 double slope = cpl_polynomial_get_coeff (fit_slope, &pow_slope);
1124
1125 /* Check slope sign. Should be positive */
1126 if (slope < 0.0 && wave == 0) {
1127 cpl_msg_warning (cpl_func, "Negative wavelength!! "
1128 "Report to DRS team");
1129 }
1130
1131 /* Set the wavelength */
1132 cpl_array_set (wavechi2, wave, sqrt(rechisq));
1133 cpl_array_set (wavelength, wave, CPL_MATH_2PI / fabs(slope));
1134
1135 /* Free memory */
1136 cpl_vector_delete (phase);
1137 cpl_vector_delete (residuals);
1138 cpl_polynomial_delete (fit_slope);
1139
1140 } /* End loop on wave */
1141
1142 cpl_matrix_delete (opd_matrix);
1143 cpl_vector_delete (opd_vector);
1144
1145 /* Set the wavelength array in the output table */
1146 sprintf (name, "BASE_%s_%s", GRAVI_BASE_NAME[base], GRAVI_POLAR(pol,npol));
1147 cpl_table_new_column_array (wave_fibre, name, CPL_TYPE_DOUBLE, nwave);
1148 cpl_table_set_column_unit (wave_fibre, name, "m");
1149 cpl_table_set_array (wave_fibre, name, 0, wavelength);
1150 cpl_array_delete (wavelength);
1151
1152 /* Set the chi2 array in the output table */
1153 sprintf (name, "BASE_%s_%s_CHI2", GRAVI_BASE_NAME[base], GRAVI_POLAR(pol,npol));
1154 cpl_table_new_column_array (wave_fibre, name, CPL_TYPE_DOUBLE, nwave);
1155 cpl_table_set_array (wave_fibre, name, 0, wavechi2);
1156 cpl_array_delete (wavechi2);
1157
1158 } /* End loop on base*/
1159 } /* End loop on polar */
1160
1161 gravi_msg_function_exit(1);
1162 return wave_fibre;
1163}
1164
1165/*----------------------------------------------------------------------------*/
1183/*----------------------------------------------------------------------------*/
1184
1185cpl_table * gravi_wave_fit_2d (cpl_table * wavefibre_table,
1186 cpl_table * detector_table,
1187 gravi_data * wave_param,
1188 cpl_size fullstartx,
1189 int spatial_order,
1190 int spectral_order,
1191 double * rms_residuals)
1192{
1193 gravi_msg_function_start(1);
1194 cpl_ensure (wavefibre_table, CPL_ERROR_NULL_INPUT, NULL);
1195 cpl_ensure (detector_table, CPL_ERROR_NULL_INPUT, NULL);
1196
1197 char name[100];
1198 *rms_residuals = 0;
1199
1200 /* Get numbers */
1201 cpl_size n_region = cpl_table_get_nrow (detector_table);
1202 int npol = n_region > 24 ? 2 : 1;
1203 cpl_size nwave = cpl_table_get_column_depth (wavefibre_table, npol > 1 ? "BASE_12_S" : "BASE_12_C");
1204 CPLCHECK_NUL ("Cannot get data");
1205
1206 /* get the wave param */
1207 cpl_propertylist * wave_param_plist = gravi_data_get_plist(wave_param,GRAVI_PRIMARY_HDR_EXT);
1208
1209 /* Odd index, for SC only */
1210 cpl_vector * odd_index = cpl_vector_new (nwave);
1211 for (int i = fullstartx; i < fullstartx + nwave; i++) {
1212 if (nwave > GRAVI_LBD_FTSC) cpl_vector_set (odd_index, i - fullstartx, ((i/64)%2 == 0) ? 0 : 1);
1213 else cpl_vector_set (odd_index, i - fullstartx, 0);
1214 }
1215
1216 CPLCHECK_NUL ("Cannot buid odd_index");
1217
1218 /* Save the 2D coeficient of each polar */
1219 cpl_polynomial ** coef_poly = cpl_calloc (npol, sizeof (cpl_polynomial*));
1220
1221 /* Loop on polarisation */
1222 for (int pol = 0; pol < npol; pol++) {
1223
1224 /* Prepare the 2D coordinates
1225 * and values to fit */
1226 cpl_vector * coord_X = cpl_vector_new (GRAVI_NBASE * nwave);
1227 cpl_vector * coord_Y = cpl_vector_new (GRAVI_NBASE * nwave);
1228
1229 cpl_vector * all_wavelength = cpl_vector_new (GRAVI_NBASE * nwave);
1230 cpl_vector * all_wavechi2 = cpl_vector_new (GRAVI_NBASE * nwave);
1231 cpl_vector * all_valid = cpl_vector_new (GRAVI_NBASE * nwave);
1232
1233 /*
1234 * Loop on base and wave to get all
1235 * wavelenght and coordinates
1236 */
1237 for (int base = 0; base < GRAVI_NBASE; base ++) {
1238
1239 /* Mean position of this baseline */
1240 int iA = gravi_get_region (detector_table, base, 'A', pol);
1241 int iB = gravi_get_region (detector_table, base, 'B', pol);
1242 int iC = gravi_get_region (detector_table, base, 'C', pol);
1243 int iD = gravi_get_region (detector_table, base, 'D', pol);
1244
1245 /* WAVE_FIBRE data */
1246 sprintf (name, "BASE_%s_%s", GRAVI_BASE_NAME[base], GRAVI_POLAR(pol,npol));
1247 cpl_array * wavelength = cpl_table_get_data_array (wavefibre_table, name)[0];
1248
1249 sprintf (name, "BASE_%s_%s_CHI2", GRAVI_BASE_NAME[base], GRAVI_POLAR(pol,npol));
1250 cpl_array * wavechi2 = cpl_table_get_data_array (wavefibre_table, name)[0];
1251
1252 /* Loop on wave */
1253 for (cpl_size wave = 0; wave < nwave; wave++) {
1254 cpl_size pos = base * nwave + wave;
1255
1256 /* Get the values */
1257 int nv = 0;
1258 double wave_value = cpl_array_get (wavelength, wave, &nv);
1259 double chi2_value = cpl_array_get (wavechi2, wave, &nv);
1260
1261 /* The FT accept all channel for the 2D fit */
1262 cpl_vector_set (all_valid, pos, 1);
1263
1264 if (nwave > 1000) {
1265 /* SC HIGH */
1266 if ((chi2_value > M_PI_4 ||
1267 wave_value < gravi_pfits_get_double_default(wave_param_plist, "ESO OIWAVE HIGH LBD MIN", 2.01e-6) ||
1268 wave_value > gravi_pfits_get_double_default(wave_param_plist, "ESO OIWAVE HIGH LBD MAX", 2.43e-6)))
1269 cpl_vector_set (all_valid, pos, 0);
1270 }
1271 else if (nwave > 100) {
1272 /* SC MEDIUM */
1273 if ((chi2_value > M_PI_4 ||
1274 wave_value < gravi_pfits_get_double_default(wave_param_plist, "ESO OIWAVE MED LBD MIN", 2.01e-6) ||
1275 wave_value > gravi_pfits_get_double_default(wave_param_plist, "ESO OIWAVE MED LBD MAX", 2.43e-6)))
1276 cpl_vector_set (all_valid, pos, 0);
1277 }
1278 else if (nwave > GRAVI_LBD_FTSC) {
1279 /* SC LOW */
1280 if ((chi2_value > M_PI_4 ||
1281 //wave_value < 2.01e-6 ||
1282 wave_value < gravi_pfits_get_double_default(wave_param_plist, "ESO OIWAVE LOW LBD MIN", 1.99e-6) ||
1283 wave_value > gravi_pfits_get_double_default(wave_param_plist, "ESO OIWAVE LOW LBD MAX", 2.5e-6) ||
1284 wave == 0 || wave == nwave-1))
1285 cpl_vector_set (all_valid, pos, 0);
1286 }
1287 else if (nwave == GRAVI_LBD_FTSC) {
1288 /* FT */
1289 if ((chi2_value > M_PI_4 ||
1290 //wave_value < 2.01e-6 ||
1291 wave_value < 1.99e-6 ||
1292 wave_value > 2.5e-6))
1293 cpl_vector_set (all_valid, pos, 0);
1294 }
1295
1296 /* Set the wavelength */
1297 cpl_vector_set (all_wavelength, pos, wave_value);
1298 cpl_vector_set (all_wavechi2, pos, chi2_value);
1299
1300 /* Set the X position as the mean of the 4 regions */
1301 cpl_vector_set (coord_X, pos, (double)(iA + iB + iC + iD) / 4.);
1302
1303 /* Set the Y position. Add a shift of 0.15 pixels
1304 * on odd output for SC detector */
1305 cpl_vector_set (coord_Y, pos, wave + cpl_vector_get (odd_index, wave)*0.15);
1306
1307 } /* End loop on wave */
1308 } /* End loop on base */
1309
1310 CPLCHECK_NUL ("Error get all wavelength");
1311
1312 /*
1313 * Reformat the valid point in vector and matrix
1314 */
1315 cpl_size nvalid = cpl_vector_get_sum (all_valid);
1316
1317 cpl_msg_info (cpl_func, "Remove %lld/%lld bad wave",
1318 GRAVI_NBASE*nwave - nvalid, GRAVI_NBASE * nwave);
1319
1320 cpl_vector * vector = cpl_vector_new (nvalid);
1321 cpl_matrix * matrix = cpl_matrix_new (2, nvalid);
1322
1323 for (cpl_size c = 0, i = 0 ; i < nwave * GRAVI_NBASE; i ++) {
1324 if (!cpl_vector_get (all_valid, i)) continue;
1325 cpl_vector_set (vector, c, cpl_vector_get (all_wavelength, i));
1326 cpl_matrix_set (matrix, 0, c, cpl_vector_get (coord_X, i));
1327 cpl_matrix_set (matrix, 1, c, cpl_vector_get (coord_Y, i));
1328 c++;
1329 }
1330
1331 CPLCHECK_NUL ("Error cropping all wavelength");
1332
1333 /*
1334 * Perform a 2D fit with a polynomial model
1335 * between the position and the wavelength = F(y, x)
1336 */
1337 cpl_size deg2d[2] = {2, 3};
1338 if ( (nwave < 20) && (nwave > 8) ) {deg2d[0] = 2; deg2d[1] = 7;} // FIXME Useless line ?
1339 deg2d[0] = spatial_order;
1340 deg2d[1] = spectral_order;
1341
1342 cpl_msg_info (cpl_func, "Fit a 2d polynomial {%lli..%lli} to the wavelengths map", deg2d[0], deg2d[1]);
1343
1344 cpl_polynomial * fit2d = cpl_polynomial_new (2);
1345 cpl_polynomial_fit (fit2d, matrix, NULL, vector, NULL, CPL_TRUE, NULL, deg2d);
1346 coef_poly[pol] = fit2d;
1347
1348 CPLCHECK_NUL ("Cannot fit 2D");
1349
1350 /*
1351 * Compute residuals
1352 */
1353 double rechisq = 0.0;
1354 cpl_vector * residuals = cpl_vector_new (nvalid);
1355 cpl_vector_fill_polynomial_fit_residual (residuals, vector, NULL, fit2d, matrix, &rechisq);
1356 *rms_residuals += cpl_vector_get_stdev(residuals)/npol;
1357 FREE (cpl_vector_delete, residuals);
1358 CPLCHECK_NUL ("Cannot compute residuals");
1359
1360
1361 FREE (cpl_matrix_delete, matrix);
1362 FREE (cpl_vector_delete, vector);
1363 FREE (cpl_vector_delete, all_wavelength);
1364 FREE (cpl_vector_delete, all_wavechi2);
1365 FREE (cpl_vector_delete, all_valid);
1366 FREE (cpl_vector_delete, coord_X);
1367 FREE (cpl_vector_delete, coord_Y);
1368
1369 } /* End loop on polarisation */
1370
1371
1372 /*
1373 * Create and fill the interpolated WAVE_DATA table
1374 */
1375
1376 cpl_table * wavedata_table = cpl_table_new (1);
1377 cpl_vector * pos = cpl_vector_new (2);
1378 cpl_array * value = cpl_array_new (nwave, CPL_TYPE_DOUBLE);
1379
1380 for (cpl_size region = 0 ; region < n_region; region ++) {
1381
1382 int pol = gravi_region_get_pol (detector_table, region);
1383 /* Loop on wave to evaluate the 2D polynome */
1384 for (cpl_size wave = 0; wave < nwave; wave ++) {
1385
1386 /* Evaluate */
1387 cpl_vector_set (pos, 0, region);
1388 cpl_vector_set (pos, 1, wave + cpl_vector_get (odd_index, wave)*0.15);
1389
1390 double result = cpl_polynomial_eval (coef_poly[pol], pos);
1391 cpl_array_set (value, wave, result);
1392
1393 }
1394 /* ensure cresent wavelength */
1395 double previous_wave = cpl_array_get(value, nwave/2, NULL);
1396 for (cpl_size wave_loop = nwave/2 ; wave_loop >= 0 ; wave_loop --){
1397 if (previous_wave < cpl_array_get(value, wave_loop, NULL))
1398 cpl_array_set(value, wave_loop, previous_wave);
1399 else previous_wave = cpl_array_get(value, wave_loop, NULL);
1400 }
1401
1402 previous_wave = cpl_array_get(value, nwave/2, NULL);
1403 for (cpl_size wave_loop = nwave/2 ; wave_loop < nwave ; wave_loop ++){
1404 if (previous_wave > cpl_array_get(value, wave_loop, NULL))
1405 cpl_array_set(value, wave_loop, previous_wave);
1406 else previous_wave = cpl_array_get(value, wave_loop, NULL);
1407 }
1408
1409 /* Add column */
1410 char * data_x = GRAVI_DATA[region];
1411 cpl_table_new_column_array (wavedata_table, data_x, CPL_TYPE_DOUBLE, nwave);
1412 cpl_table_set_array (wavedata_table, data_x, 0, value);
1413
1414 } /* End loop on regions */
1415
1416 /* Delete allocations */
1417 FREE (cpl_vector_delete, pos);
1418 FREE (cpl_array_delete, value);
1419 FREELOOP (cpl_polynomial_delete, coef_poly, npol);
1420 FREE (cpl_vector_delete, odd_index);
1421
1422 gravi_msg_function_exit(1);
1423 return wavedata_table;
1424}
1425
1426/*----------------------------------------------------------------------------*/
1427
1435/*----------------------------------------------------------------------------*/
1436cpl_array * my_ft_discrete(
1437 const cpl_vector *x, const cpl_vector *y, const cpl_vector *freqs
1438) {
1439 const cpl_size nx = cpl_vector_get_size(x);
1440 const cpl_size ny = cpl_vector_get_size(y);
1441 const cpl_size nf = cpl_vector_get_size(freqs);
1442
1443 cpl_ensure(nx == ny, CPL_ERROR_INCOMPATIBLE_INPUT, NULL);
1444
1445 cpl_matrix *kern_re = cpl_matrix_new(nx, nf);
1446 cpl_matrix *kern_im = cpl_matrix_new(nx, nf);
1447
1448 for (int j = 0; j < nx; j++) {
1449 double xval = cpl_vector_get(x, j);
1450 double yval = cpl_vector_get(y, j);
1451 for (int i = 0; i < nf; i++) {
1452 double fval = cpl_vector_get(freqs, i);
1453 double complex Cval = yval * cexp(-2.0 * M_PI * I * fval * xval);
1454 cpl_matrix_set(kern_re, j, i, creal(Cval));
1455 cpl_matrix_set(kern_im, j, i, cimag(Cval));
1456 }
1457 }
1458 CPLCHECK_NUL("Cannot compute fourier matrix");
1459
1460 cpl_array *result_re = cpl_array_new(nf, CPL_TYPE_DOUBLE);
1461 cpl_array *result_im = cpl_array_new(nf, CPL_TYPE_DOUBLE);
1462 for (int i = 0; i < nf; i++) {
1463 double re = 0;
1464 double im = 0;
1465 for (int j = 0; j < nx; j++) {
1466 re += cpl_matrix_get(kern_re, j, i);
1467 im += cpl_matrix_get(kern_im, j, i);
1468 }
1469 cpl_array_set(result_re, i, re);
1470 cpl_array_set(result_im, i, im);
1471 }
1472 CPLCHECK_NUL("Cannot collapse fourier matrix");
1473 cpl_array *result = cpl_array_new_complex_from_arrays(result_re, result_im);
1474
1475 cpl_matrix_delete(kern_re);
1476 cpl_matrix_delete(kern_im);
1477 cpl_array_delete(result_re);
1478 cpl_array_delete(result_im);
1479 return result;
1480} CPL_ATTR_ALLOC;
1481
1482/*----------------------------------------------------------------------------*/
1493/*----------------------------------------------------------------------------*/
1494cpl_table * gravi_wave_fit_bandpass(
1495 cpl_table * spectrum_table, cpl_table * detector_table, cpl_table * opd_table
1496) {
1497 gravi_msg_function_start(1);
1498 cpl_ensure (spectrum_table, CPL_ERROR_NULL_INPUT, NULL);
1499 cpl_ensure (detector_table, CPL_ERROR_NULL_INPUT, NULL);
1500 cpl_ensure (opd_table, CPL_ERROR_NULL_INPUT, NULL);
1501
1502 /* Get the number of wavelength, region, polarisation... */
1503 cpl_size nwave = cpl_table_get_column_depth (spectrum_table, "DATA1");
1504 cpl_size nrow = cpl_table_get_nrow (spectrum_table);
1505 cpl_size n_region = cpl_table_get_nrow (detector_table);
1506 const int npol = (n_region > 24 ? 2 : 1);
1507
1508 /* Create the output table */
1509 cpl_table * wave_bandpass = cpl_table_new (1);
1510
1511 /* Create grid containing the full range of valid wavenumbers */
1512 const double wave_min = 1.8 * 1e-6;
1513 const double wave_max = 2.8 * 1e-6;
1514 // const double delta_wavenumber = 0.001;
1515 // const cpl_size ngrid = ((1.0 / wave_min) - (1.0 / wave_max)) / delta_wavenumber;
1516 const cpl_size ngrid = 1000;
1517 const double delta_wavenumber = ((1.0 / wave_max) - (1.0 / wave_min)) / ngrid;
1518
1519 cpl_vector *wavenumber = cpl_vector_new (ngrid);
1520 for (int i = 0; i < ngrid; i++) {
1521 double wn = (1.0 / wave_min) + i * delta_wavenumber;
1522 cpl_vector_set (wavenumber, i, wn);
1523 }
1524
1525 /* Create all table columns first to put them in order */
1526 for (int reg = 0; reg < n_region; reg++) {
1527 cpl_table_new_column_array (wave_bandpass, GRAVI_DATA[reg], CPL_TYPE_DOUBLE, nwave);
1528 }
1529
1530 /*
1531 * Calibration of each polarization and base
1532 */
1533 for (int pol = 0; pol < npol; pol++) {
1534 for (int base = 0; base < GRAVI_NBASE; base ++) {
1535 cpl_msg_info (cpl_func, "Compute bandpass for pol %i over %i, base %i over %i",
1536 pol+1, npol, base+1, GRAVI_NBASE);
1537
1538 /* Get the OPD */
1539 cpl_vector * opd_vector = cpl_vector_new (nrow);
1540 for (cpl_size row = 0; row < nrow; row ++ ) {
1541 double value = cpl_table_get (opd_table, "OPD", row*GRAVI_NBASE+base, NULL);
1542 cpl_vector_set (opd_vector, row, value);
1543 }
1544 CPLCHECK_NUL ("Cannot extract the OPD");
1545
1546 for (char pha = 'A'; pha <= 'D'; pha++) {
1547 int reg = gravi_get_region (detector_table, base, pha, pol);
1548 if (reg < 0)
1549 cpl_msg_warning (cpl_func, "Don't found the A, B, C or D !!!");
1550 cpl_msg_debug(cpl_func, "REGION=%d", reg);
1551
1552 /* Compute FT for this REGION */
1553 cpl_array *reg_array = cpl_array_new(nwave, CPL_TYPE_DOUBLE);
1554 for (int wave = 0; wave < nwave; wave++) {
1555 cpl_vector *counts = gravi_table_get_vector (spectrum_table, wave, GRAVI_DATA[reg]);
1556 cpl_array *ftrans = my_ft_discrete (opd_vector, counts, wavenumber);
1557 if (!ftrans) {
1558 cpl_error_set(cpl_func, CPL_ERROR_ILLEGAL_OUTPUT);
1559 CPLCHECK_NUL("Cannot compute fourier transform");
1560 }
1561
1562 /* Take absolute value and find argmax */
1563 cpl_size argmax;
1564 cpl_array_abs (ftrans);
1565 cpl_array_get_maxpos (ftrans, &argmax);
1566
1567 if (argmax != 0 && argmax != ngrid - 1) {
1568 // /* Walk left to find turning point */
1569 // int im;
1570 // double x = cpl_array_get(ftrans, argmax, NULL);
1571 // for (im = argmax; im > 0; im--) {
1572 // double xm = cpl_array_get(ftrans, im, NULL);
1573 // double dx = xm - x;
1574 // if (dx > 0)
1575 // break;
1576 // x = xm;
1577 // }
1578
1579 // /* Walk right to find turning point */
1580 // int ip;
1581 // x = cpl_array_get(ftrans, argmax, NULL);
1582 // for (ip = argmax; ip < ngrid; ip++) {
1583 // double xp = cpl_array_get(ftrans, ip, NULL);
1584 // double dx = xp - x;
1585 // if (dx > 0)
1586 // break;
1587 // x = xp;
1588 // }
1589 // if (!(im < ip)) {
1590 // cpl_error_set(cpl_func, CPL_ERROR_ILLEGAL_OUTPUT);
1591 // CPLCHECK_NUL("Cannot find minima of central peak in bandpass");
1592 // }
1593
1594 // double avg = 0.0;
1595 // double wgt = 0.0;
1596 // for (int i = im; i < ip; i++) {
1597 // double wn = cpl_vector_get(wavenumber, i);
1598 // double w = cpl_array_get(ftrans, i, NULL);
1599 // avg += wn * w;
1600 // wgt += w;
1601 // }
1602 // avg /= wgt;
1603
1604 cpl_array_set (reg_array, wave, 1.0 / cpl_vector_get(wavenumber, argmax));
1605 // cpl_array_set (reg_array, wave, 1.0 / avg);
1606 } else {
1607 cpl_array_set (reg_array, wave, 0.0);
1608 }
1609
1610 cpl_vector_delete (counts);
1611 cpl_array_delete (ftrans);
1612 CPLCHECK_NUL ("Cannnot compute bandpass");
1613 }
1614 cpl_table_set_array(wave_bandpass, GRAVI_DATA[reg], 0, reg_array);
1615 cpl_array_delete(reg_array);
1616 }
1617
1618 cpl_vector_delete(opd_vector);
1619 }
1620 }
1621
1622 cpl_vector_delete(wavenumber);
1623
1624 /* Fix regions with no counts */
1625 for (int reg = 0; reg < n_region; reg++) {
1626 cpl_array *reg_arr = cpl_table_get_data_array(wave_bandpass, GRAVI_DATA[reg])[0];
1627 for (int wave = 0; wave < nwave; wave++) {
1628 cpl_vector *counts = gravi_table_get_vector (spectrum_table, wave, GRAVI_DATA[reg]);
1629 if (cpl_vector_get_sum(counts) == 0.0) {
1630 cpl_msg_warning(cpl_func, "region %s, pixel %d has no counts", GRAVI_DATA[reg], wave);
1631
1632 if (wave == 0)
1633 cpl_array_set(reg_arr, wave, cpl_array_get(reg_arr, 1, NULL));
1634 else if (wave == nwave - 1)
1635 cpl_array_set(reg_arr, wave, cpl_array_get(reg_arr, nwave - 2, NULL));
1636 else
1637 cpl_msg_error(cpl_func, "Invalid wave on intermediate pixel");
1638 }
1639 cpl_vector_delete (counts);
1640 }
1641 }
1642
1643 gravi_msg_function_exit(1);
1644 return wave_bandpass;
1645};
1646
1647/*----------------------------------------------------------------------------*/
1648/*
1649 * @brief Compute wavelength calibration using independent polynomial fitting.
1650 * @param TBD
1651 */
1652/*----------------------------------------------------------------------------*/
1653cpl_table * gravi_wave_fit_individual (cpl_table * wave_individual_table,
1654 cpl_table * weight_individual_table,
1655 cpl_table * wave_fitted_table,
1656 cpl_table * opd_table,
1657 cpl_table * spectrum_table,
1658 cpl_table * detector_table,
1659 double n0, double n1, double n2)
1660{
1661
1662 gravi_msg_function_start(1);
1663
1664 cpl_ensure (wave_individual_table, CPL_ERROR_NULL_INPUT, NULL);
1665 cpl_ensure (weight_individual_table, CPL_ERROR_NULL_INPUT, NULL);
1666 cpl_ensure (wave_fitted_table, CPL_ERROR_NULL_INPUT, NULL);
1667 cpl_ensure (opd_table, CPL_ERROR_NULL_INPUT, NULL);
1668 cpl_ensure (spectrum_table, CPL_ERROR_NULL_INPUT, NULL);
1669 cpl_ensure (detector_table, CPL_ERROR_NULL_INPUT, NULL);
1670
1671 /* Get the number of wavelength, region, polarisation... */
1672 cpl_size nwave = cpl_table_get_column_depth (spectrum_table, "DATA1");
1673 cpl_size n_region = cpl_table_get_nrow (detector_table);
1674 cpl_size nrow = cpl_table_get_nrow (spectrum_table);
1675 int npol = (n_region > 24 ? 2 : 1);
1676 cpl_size nwave_ref=3000;
1677 if (nwave<10) nwave_ref=600;
1678
1679 CPLCHECK_NUL ("Cannot buid odd_index");
1680
1681 /* Get OPD Table */
1682
1683
1684
1685 /* create temporary variables */
1686
1687 cpl_array * wave_individual_array = cpl_array_new (nwave, CPL_TYPE_DOUBLE);
1688 cpl_array * weight_individual_array = cpl_array_new (nwave, CPL_TYPE_DOUBLE);
1689 cpl_matrix * data_flux_matrix = cpl_matrix_new (nrow, nwave);
1690 cpl_matrix * vis_to_flux_matrix = cpl_matrix_new (nrow, 3);
1691 cpl_matrix * signal_matrix = cpl_matrix_new (nwave, nwave_ref);
1692 cpl_matrix * residual_matrix = cpl_matrix_new (nwave, nwave_ref);
1693 cpl_array * wave_reference_array = cpl_array_new (nwave_ref, CPL_TYPE_DOUBLE);
1694
1695 /*initialise arrays */
1696
1697 cpl_matrix_fill_column(vis_to_flux_matrix,1,2);
1698 for (cpl_size wave_ref = 0; wave_ref < nwave_ref; wave_ref+=1)
1699 {
1700 // make matrix todo
1701 double wave_value=1.95e-6+wave_ref*0.6e-6/((double) nwave_ref);
1702 cpl_array_set_double(wave_reference_array,wave_ref,wave_value);
1703 }
1704
1705 CPLCHECK_NUL ("Cannot initialize arrays for wavelength fit");
1706
1707
1708 for (cpl_size region = 0 ; region < n_region; region ++) {
1709
1710 int base=gravi_region_get_base (detector_table, region);
1711 char * data_x = GRAVI_DATA[region];
1712
1713 cpl_msg_info_overwritable (cpl_func, "Least square fitting of wavelength for region %s", data_x);
1714 // get data_flux_matrix
1715
1716 for (cpl_size row = 0; row < nrow; row ++ ) {
1717 cpl_array * flux_array= cpl_table_get_data_array(spectrum_table,data_x)[row];
1718 for (cpl_size wave = 0; wave < nwave; wave ++) {
1719 cpl_matrix_set (data_flux_matrix, row, wave, cpl_array_get(flux_array,wave, NULL));
1720 }
1721 }
1722
1723
1724 for (cpl_size wave_ref = 0; wave_ref < nwave_ref; wave_ref+=1)
1725 {
1726 // make matrix
1727 double wave_value=cpl_array_get(wave_reference_array,wave_ref,NULL);
1728
1729 for (cpl_size row = 0; row < nrow; row ++ ) {
1730 double opd = cpl_table_get (opd_table, "OPD", row*GRAVI_NBASE+base, NULL);
1731 double coherence_loss=1;
1732 if (fabs(opd) > 1e-9)
1733 {
1734 if (nwave <30)
1735 {
1736 coherence_loss=sin(opd*19500)/(opd*19500);
1737 }
1738 }
1739 cpl_matrix_set(vis_to_flux_matrix,row,0,cos(opd*6.28318530718/wave_value)*coherence_loss);
1740 cpl_matrix_set(vis_to_flux_matrix,row,1,sin(opd*6.28318530718/wave_value)*coherence_loss);
1741 }
1742
1743 cpl_matrix * coef_vis = cpl_matrix_solve_normal(vis_to_flux_matrix,data_flux_matrix); // coef_vis is 3xnwave
1744 cpl_matrix * data_flux_fit = cpl_matrix_product_create(vis_to_flux_matrix,coef_vis);
1745 cpl_matrix * residuals_fit = cpl_matrix_duplicate(data_flux_fit);
1746 cpl_matrix_subtract (residuals_fit,data_flux_matrix); //residuals_fit is nrowxnwave
1747
1748 for (cpl_size wave = 0; wave < nwave; wave ++ ) {
1749 cpl_matrix * temp_matrix = cpl_matrix_extract_column (data_flux_fit, wave);
1750 cpl_matrix_set(signal_matrix,wave,wave_ref,cpl_matrix_get_stdev(temp_matrix));
1751 FREE (cpl_matrix_delete, temp_matrix);
1752 cpl_matrix * temp_matrix2 = cpl_matrix_extract_column (residuals_fit, wave);
1753 cpl_matrix_set(residual_matrix,wave,wave_ref,cpl_matrix_get_stdev(temp_matrix2));
1754 FREE (cpl_matrix_delete, temp_matrix2);
1755 }
1756
1757 FREE (cpl_matrix_delete, coef_vis);
1758 FREE (cpl_matrix_delete, data_flux_fit);
1759 FREE (cpl_matrix_delete, residuals_fit);
1760
1761 CPLCHECK_NUL ("Cannot do Matrix inversion to calculate optimum wavelength");
1762
1763 }
1764
1765 // get minimum chi2 and amplitude signal for each wave
1766 cpl_size wave_ref=1;
1767 cpl_size discarded=1;
1768 for (cpl_size wave = 0; wave < nwave; wave ++ ) {
1769
1770 cpl_matrix * chi2_extract=cpl_matrix_extract_row(residual_matrix,wave);
1771
1772 cpl_matrix_get_minpos(chi2_extract,&discarded, &wave_ref );
1773
1774 double wave_value = cpl_array_get(wave_reference_array, wave_ref, NULL );
1775 double weight_value = cpl_matrix_get(signal_matrix, wave , wave_ref )/(0.1+cpl_matrix_get(residual_matrix, wave, wave_ref ));
1776
1777 cpl_array_set_double (wave_individual_array, wave, wave_value);
1778 cpl_array_set_double (weight_individual_array, wave, weight_value);
1779
1780 FREE (cpl_matrix_delete, chi2_extract);
1781
1782 }
1783
1784 /* Add column */
1785 cpl_table_new_column_array (wave_individual_table, data_x, CPL_TYPE_DOUBLE, nwave);
1786 cpl_table_set_array (wave_individual_table, data_x, 0, wave_individual_array);
1787 cpl_table_new_column_array (weight_individual_table, data_x, CPL_TYPE_DOUBLE, nwave);
1788 cpl_table_set_array (weight_individual_table, data_x, 0, weight_individual_array);
1789 cpl_table_new_column_array (wave_fitted_table, data_x, CPL_TYPE_DOUBLE, nwave);
1790 cpl_table_set_array (wave_fitted_table, data_x, 0, wave_individual_array);
1791 }
1792
1793 CPLCHECK_NUL ("Cannot get individual wavelength for each pixel");
1794
1795 FREE (cpl_array_delete ,wave_individual_array);
1796 FREE (cpl_array_delete ,weight_individual_array);
1797 FREE (cpl_array_delete ,wave_reference_array);
1798 FREE (cpl_matrix_delete ,data_flux_matrix);
1799 FREE (cpl_matrix_delete ,vis_to_flux_matrix);
1800 FREE (cpl_matrix_delete ,signal_matrix);
1801 FREE (cpl_matrix_delete ,residual_matrix);
1802
1803 cpl_msg_info (cpl_func, "Now fitting polynomials on wavelength channels");
1804
1805 cpl_matrix * coef_to_wave = cpl_matrix_new (n_region / npol,5);
1806 cpl_matrix * coef_to_wave_weight = cpl_matrix_new (n_region / npol,n_region / npol);
1807 cpl_matrix * wavelength = cpl_matrix_new(n_region / npol,1);
1808
1809 // set coordinates
1810 for (cpl_size region = 0 ; region < n_region/ npol; region ++)
1811 {
1812 double mean_region = region - (n_region/npol-1)*0.5;
1813 cpl_matrix_set (coef_to_wave, region, 0, 1);
1814 cpl_matrix_set (coef_to_wave, region, 1, mean_region);
1815 cpl_matrix_set (coef_to_wave, region, 2, mean_region*mean_region);
1816 cpl_matrix_set (coef_to_wave, region, 3, mean_region*mean_region*mean_region);
1817 cpl_matrix_set (coef_to_wave, region, 4, mean_region*mean_region*mean_region*mean_region);
1818 }
1819
1820 for (int pol = 0; pol < npol; pol++) {
1821
1822 cpl_msg_info (cpl_func, "Looping for polyfit now, with pol: %i",(int) pol);
1823
1824 for (cpl_size wave = 0; wave < nwave; wave ++) {
1825
1826 // get the data for a common wave row
1827 for (cpl_size region = 0 ; region < n_region/ npol; region ++) {
1828
1829 // Get the pointers to the table arrays
1830 char * data_x = GRAVI_DATA[region*npol+pol];
1831
1832
1833 const cpl_array * wave_array = cpl_table_get_array (wave_individual_table, data_x, 0);
1834 const cpl_array * weight_array = cpl_table_get_array (weight_individual_table, data_x, 0);
1835
1836
1837 cpl_matrix_set (wavelength, region, 0, cpl_array_get(wave_array,wave,NULL));
1838 double weight_value=cpl_array_get(weight_array,wave,NULL);
1839 cpl_matrix_set (coef_to_wave_weight, region, region, weight_value*weight_value);
1840
1841 }
1842
1843
1844 cpl_matrix * coef_to_wave2 = cpl_matrix_product_create(coef_to_wave_weight,coef_to_wave);
1845 cpl_matrix * wavelength2 = cpl_matrix_product_create(coef_to_wave_weight,wavelength);
1846
1847 // Fit a second order polynomial
1848 cpl_matrix * coeff = cpl_matrix_solve_normal(coef_to_wave2, wavelength2); // 5 x 1
1849 cpl_matrix * wavelength_fitted = cpl_matrix_product_create(coef_to_wave, coeff); // n_region / npol x 1
1850
1851 CPLCHECK_NUL ("Cannot do Matrix inversion to calculate optimum polynomial for wavelength");
1852
1853 // Write result to new wave table
1854 for (cpl_size region = 0 ; region < n_region/ npol; region ++) {
1855
1856 char * data_x = GRAVI_DATA[region*npol+pol];
1857 //cpl_msg_info (cpl_func, "Writing: region %i",region);
1858 cpl_array * wave_array = cpl_table_get_data_array (wave_fitted_table, data_x)[0];
1859
1860 cpl_array_set_double(wave_array,wave,cpl_matrix_get(wavelength_fitted,region,0));
1861
1862 }
1863
1864 FREE (cpl_matrix_delete ,coef_to_wave2);
1865 FREE (cpl_matrix_delete ,wavelength2);
1866 FREE (cpl_matrix_delete ,coeff);
1867 FREE (cpl_matrix_delete ,wavelength_fitted);
1868 }
1869
1870 }
1871 FREE (cpl_matrix_delete ,coef_to_wave);
1872 FREE (cpl_matrix_delete ,coef_to_wave_weight);
1873 FREE (cpl_matrix_delete ,wavelength);
1874
1875 CPLCHECK_NUL ("Cannot fit individual wavelength with 3rd order polynomial");
1876
1877 cpl_msg_info (cpl_func, "Correcting for wavelength error");
1878
1879
1880 /* Correct the computed wavelength from the dispersion */
1881 for (cpl_size region = 0 ; region < n_region; region ++)
1882 {
1883 char * data_x = GRAVI_DATA[region];
1884 cpl_array * wavelength_fitted = cpl_table_get_data_array (wave_fitted_table, data_x)[0];
1885 cpl_size nwave_fitted = cpl_array_get_size (wavelength_fitted);
1886 for (cpl_size wave = 0 ; wave < nwave_fitted ; wave ++ ) {
1887
1888 double result = cpl_array_get (wavelength_fitted, wave, NULL);
1889 double d_met = (result - LAMBDA_MET) / LAMBDA_MET;
1890 cpl_array_set (wavelength_fitted, wave, result * (n0 + n1*d_met + n2*d_met*d_met));
1891
1892 }
1893 for (cpl_size wave = nwave_fitted/2 ; wave < nwave_fitted-1 ; wave ++ ) {
1894 double result = cpl_array_get (wavelength_fitted, wave, NULL);
1895 double result2 = cpl_array_get (wavelength_fitted, wave+1, NULL);
1896 if (result2<result+2e-10) {
1897 result2=result+2e-10;
1898 cpl_array_set (wavelength_fitted, wave+1, result2);
1899 }
1900 }
1901 for (cpl_size wave = nwave_fitted/2 ; wave > 0 ; wave -- ) {
1902 double result = cpl_array_get (wavelength_fitted, wave, NULL);
1903 double result2 = cpl_array_get (wavelength_fitted, wave-1, NULL);
1904 if (result2>result-2e-10) {
1905 result2=result-2e-10;
1906 cpl_array_set (wavelength_fitted, wave-1, result2);
1907 }
1908 }
1909 }
1910
1911 CPLCHECK_NUL ("Error in correcting for dispersion");
1912
1913 gravi_msg_function_exit(1);
1914 return wave_fitted_table;
1915}
1916
1917
1918
1919
1920/*----------------------------------------------------------------------------*/
1927/*----------------------------------------------------------------------------*/
1928
1929cpl_error_code gravi_wave_correct_dispersion (cpl_table * wave_fibre,
1930 double n0, double n1, double n2)
1931{
1932 gravi_msg_function_start(1);
1933 cpl_ensure_code (wave_fibre, CPL_ERROR_NULL_INPUT);
1934
1935 char name[100];
1936
1937 /* Get the number of polarisation */
1938 int npol = cpl_table_has_column (wave_fibre, "BASE_12_P") ? 2 : 1;
1939
1940 /* Loop on columns in the table */
1941 for (int pol = 0; pol < npol; pol ++) {
1942 for (int base = 0; base < GRAVI_NBASE; base ++) {
1943
1944 /* Get data of this region */
1945 sprintf (name, "BASE_%s_%s", GRAVI_BASE_NAME[base], GRAVI_POLAR(pol, npol));
1946 cpl_array * wavelength = cpl_table_get_data_array (wave_fibre, name)[0];
1947 CPLCHECK_MSG ("Cannot get data");
1948
1949 /* Loop on wave */
1950 cpl_size nwave = cpl_array_get_size (wavelength);
1951 for (cpl_size wave = 0 ; wave < nwave ; wave ++ ) {
1952
1953 /* Correct the computed wavelength from the dispersion */
1954 double result = cpl_array_get (wavelength, wave, NULL);
1955 double d_met = (result - LAMBDA_MET) / LAMBDA_MET;
1956 cpl_array_set (wavelength, wave, result * (n0 + n1*d_met + n2*d_met*d_met));
1957
1958 }
1959 } /* End loop on base */
1960 } /* End loop on pol */
1961
1962 gravi_msg_function_exit(1);
1963 return CPL_ERROR_NONE;
1964}
1965
1966/*----------------------------------------------------------------------------*/
1972/*----------------------------------------------------------------------------*/
1973
1974cpl_error_code gravi_wave_correct_color (gravi_data * vis_data)
1975{
1976 gravi_msg_function_start(1);
1977 cpl_ensure_code (vis_data, CPL_ERROR_NULL_INPUT);
1978
1979 cpl_propertylist * primary_header = gravi_data_get_header (vis_data);
1980 cpl_propertylist * oiwave_header = NULL;
1981 cpl_table * oiwave_table = NULL;
1982
1983 /* For each type of data SC / FT */
1984 int ntype_data = 2;
1985 for (int type_data = 0; type_data < ntype_data ; type_data ++) {
1986
1987 /* Loop on polarisation */
1988 int npol = gravi_pfits_get_pola_num (primary_header, type_data);
1989 for (int pol = 0 ; pol<npol ; pol++) {
1990 oiwave_table = cpl_table_duplicate ( gravi_data_get_oi_wave ( vis_data, type_data, pol, npol ) );
1991 oiwave_header = cpl_propertylist_duplicate ( gravi_data_get_oi_wave_plist ( vis_data, type_data, pol, npol ) );
1992
1993 /* here you can do what you want on this duplicated oi_wave_table
1994 * to get OI_FLUX table :
1995 * cpl_table * oiflux_table = gravi_data_get_oi_flux(vis_data, type_data, pol, npol)
1996 * */
1997
1998
1999 /* save the new extension with name OI_WAVELENGTH_CORR */
2000 gravi_data_add_table(vis_data, oiwave_header, "OI_WAVELENGTH_CORR", oiwave_table);
2001
2002 CPLCHECK_MSG("Cannot apply color wave correction");
2003 }
2004 /* End loop on polarisation */
2005 }
2006 /* End loop on data_type */
2007
2008 gravi_msg_function_exit(1);
2009 return CPL_ERROR_NONE;
2010}
2011
2012
2013/*----------------------------------------------------------------------------*/
2025/*----------------------------------------------------------------------------*/
2026
2027cpl_imagelist * gravi_wave_test_image (cpl_table * wavedata_table,
2028 cpl_table * wavefibre_table,
2029 cpl_table * profile_table,
2030 cpl_table * detector_table)
2031{
2032 gravi_msg_function_start(1);
2033 cpl_ensure (wavedata_table, CPL_ERROR_NULL_INPUT, NULL);
2034 cpl_ensure (wavefibre_table, CPL_ERROR_NULL_INPUT, NULL);
2035 cpl_ensure (detector_table, CPL_ERROR_NULL_INPUT, NULL);
2036 cpl_ensure (profile_table, CPL_ERROR_NULL_INPUT, NULL);
2037
2038 int nv = 0;
2039 char name[100];
2040
2041 cpl_size n_region = cpl_table_get_nrow (detector_table);
2042 int npol = (n_region > 24 ? 2 : 1);
2043
2044 CPLCHECK_NUL ("Cannot get data");
2045
2046 /* Init the output images */
2047 cpl_size sizex = cpl_table_get_column_dimension (profile_table, "DATA1", 0);
2048 cpl_size sizey = cpl_table_get_column_dimension (profile_table, "DATA1", 1);
2049
2050 cpl_image * profilesum_image = cpl_image_new (sizex, sizey, CPL_TYPE_DOUBLE);
2051 cpl_image_fill_window (profilesum_image, 1, 1, sizex, sizey, 0.0);
2052
2053 cpl_image * wave_image = cpl_image_new (sizex, sizey, CPL_TYPE_DOUBLE);
2054 cpl_image_fill_window (wave_image, 1, 1, sizex, sizey, 0.0);
2055
2056 cpl_image * realwave_image = cpl_image_new (sizex, sizey, CPL_TYPE_DOUBLE);
2057 cpl_image_fill_window (realwave_image, 1, 1, sizex, sizey, 0.0);
2058
2059 CPLCHECK_NUL ("Cannot prepare output image");
2060
2061 /* Loop on regions */
2062 for (cpl_size region = 0 ; region < n_region; region ++) {
2063
2064 /* Load the profile of this region as an image */
2065 cpl_imagelist * profile_imglist = gravi_imagelist_wrap_column (profile_table, GRAVI_DATA[region]);
2066 cpl_image * profile_image = cpl_imagelist_get (profile_imglist, 0);
2067
2068 CPLCHECK_NUL ("Cannot get data");
2069
2070 /* Sum all profils of all regions */
2071 cpl_image_add (profilesum_image, profile_image);
2072
2073 /*
2074 * Define an image of each region with its WAVE_DATA
2075 */
2076 const cpl_array * wavelength;
2077 wavelength = cpl_table_get_array (wavedata_table, GRAVI_DATA[region], 0);
2078 CPLCHECK_NUL ("Cannot get data");
2079
2080 for (cpl_size x = 0; x < sizex; x ++){
2081 for (cpl_size y = 0; y < sizey; y ++){
2082 if (cpl_image_get (profile_image, x+1, y+1, &nv) > 0.01)
2083 cpl_image_set (wave_image, x+1, y+1,
2084 cpl_array_get (wavelength, x, NULL));
2085 }
2086 }
2087 CPLCHECK_NUL ("Cannot make image of wave_map");
2088
2089 /*
2090 * Define an image of each region with its WAVE_FIBER
2091 */
2092 int base = gravi_region_get_base (detector_table, region);
2093 int pol = gravi_region_get_pol (detector_table, region);
2094 sprintf (name, "BASE_%s_%s", GRAVI_BASE_NAME[base], GRAVI_POLAR(pol, npol));
2095 wavelength = cpl_table_get_array (wavefibre_table, name, 0);
2096 CPLCHECK_NUL ("Cannot get data");
2097
2098 for (cpl_size x = 0; x < sizex; x ++){
2099 for (cpl_size y = 0; y < sizey; y ++){
2100 if (cpl_image_get (profile_image, x+1, y+1, &nv) > 0.01)
2101 cpl_image_set (realwave_image, x+1, y+1,
2102 cpl_array_get (wavelength, x, NULL));
2103 }
2104 }
2105
2106 gravi_imagelist_unwrap_images (profile_imglist);
2107 CPLCHECK_NUL ("Cannot make image of wave_fibre");
2108 } /* End loop on regions */
2109
2110 /* Create an imagelist */
2111 cpl_imagelist * testwave_imglist = cpl_imagelist_new ();
2112 cpl_imagelist_set (testwave_imglist, wave_image, 0);
2113 cpl_imagelist_set (testwave_imglist, profilesum_image, 1);
2114 cpl_imagelist_set (testwave_imglist, realwave_image, 2);
2115
2116 gravi_msg_function_exit(1);
2117 return testwave_imglist;
2118}
2119
2120/*----------------------------------------------------------------------------*/
2130/*----------------------------------------------------------------------------*/
2131
2132cpl_error_code gravi_wave_qc (gravi_data * wave_map, gravi_data * profile_map)
2133{
2134 gravi_msg_function_start(1);
2135 cpl_ensure_code (wave_map, CPL_ERROR_NULL_INPUT);
2136
2137 char name[100];
2138
2139 cpl_propertylist * wave_header = gravi_data_get_header (wave_map);
2140
2141 /* Loop on extensions (thus SC/FT) */
2142 for (int type_data = 0; type_data < 2; type_data ++ ) {
2143
2144 /* Get WAVE_FIBRE and IMAGING_DETECTOR tables */
2145 cpl_table * detector_table = gravi_data_get_imaging_detector (wave_map, type_data);
2146 cpl_table * wave_data = gravi_data_get_wave_data (wave_map, type_data);
2147 cpl_size n_region = cpl_table_get_nrow (detector_table);
2148
2149 /* Get the full STARTX */
2150 cpl_propertylist * plist = gravi_data_get_wave_data_plist (wave_map, type_data);
2151 int fullstartx = gravi_pfits_get_fullstartx (plist);
2152
2153 CPLCHECK_MSG ("Cannot get data");
2154
2155 /*
2156 * Compute the min and max wave over all regions
2157 */
2158 double minwave = -1e10;
2159 double maxwave = 1e10;
2160
2161 for (int region = 0; region < n_region; region++) {
2162 const cpl_array * wavelength = cpl_table_get_array (wave_data, GRAVI_DATA[region], 0);
2163
2164 minwave = CPL_MAX (minwave, cpl_array_get_min (wavelength));
2165 maxwave = CPL_MIN (maxwave, cpl_array_get_max (wavelength));
2166 } /* End loop on regions */
2167
2168 cpl_msg_info (cpl_func,"%s = %g [m]", QC_MINWAVE(type_data), minwave);
2169 cpl_msg_info (cpl_func,"%s = %g [m]", QC_MAXWAVE(type_data), maxwave);
2170 cpl_propertylist_update_double (wave_header, QC_MINWAVE(type_data), minwave);
2171 cpl_propertylist_update_double (wave_header, QC_MAXWAVE(type_data), maxwave);
2172 cpl_propertylist_update_double (wave_header, QC_MINWAVE_UM(type_data), minwave * 1e6);
2173 cpl_propertylist_update_double (wave_header, QC_MAXWAVE_UM(type_data), maxwave * 1e6);
2174
2175 CPLCHECK_MSG ("Cannot compute minwave or maxwave");
2176
2177
2178 /*
2179 * Compute the pixel position of QC specified argon wavelength.
2180 * Having two lines allow to check the position and dispersion
2181 */
2182 int qc_reg[3] = {0,23,47};
2183 double qc_wave[2] = {2.099184e-6, 2.313952e-6};
2184
2185 /* Loop on regions */
2186 for (int reg = 0; reg < (n_region > 24 ? 3 : 2); reg++) {
2187 cpl_size region = qc_reg[reg];
2188
2189 const cpl_array * wavelength = cpl_table_get_array (wave_data, GRAVI_DATA[region], 0);
2190 cpl_size nwave = cpl_array_get_size (wavelength);
2191 const double * wave_tab = cpl_array_get_data_double_const (wavelength);
2192 CPLCHECK_MSG ("Cannot get wavelength data");
2193
2194 /* Loop on the two argon lines */
2195 for (int iqc = 0 ; iqc < 2 ; iqc++) {
2196 cpl_size l2 = 0;
2197 if ( wave_tab[0] < wave_tab[nwave-1]) {while (wave_tab[l2] < qc_wave[iqc]) l2 ++;}
2198 else {while (wave_tab[l2] > qc_wave[iqc]) l2 ++;}
2199
2200 if (l2-1 < 0 || l2 > nwave-1) {
2201 cpl_msg_error (cpl_func, "Cannot find the QC position for lbd=%g", qc_wave[iqc]);
2202 continue;
2203 }
2204
2205 /* Position on full detector */
2206 double qc_pos = 0.0;
2207 qc_pos = fullstartx + (l2-1) + (qc_wave[iqc] - wave_tab[l2-1]) / (wave_tab[l2] - wave_tab[l2-1]);
2208
2209 sprintf (name, "ESO QC REFWAVE%i", iqc+1);
2210 cpl_propertylist_update_double (wave_header, name, qc_wave[iqc]);
2211 cpl_propertylist_set_comment (wave_header, name, "[m] value of ref wave");
2212
2213 sprintf (name, "ESO QC REFPOS%i %s%lli", iqc+1, GRAVI_TYPE(type_data),region+1);
2214 cpl_propertylist_update_double (wave_header, name, qc_pos);
2215 cpl_propertylist_set_comment (wave_header, name, "[pix] position of ref wave");
2216
2217 cpl_msg_info (cpl_func, "%s = %f [pix] for %e [m]", name, qc_pos, qc_wave[iqc]);
2218
2219 CPLCHECK_MSG ("Cannot set QC");
2220 } /* End loop on the 2 argon lines */
2221 } /* End loop on 2 or 3 regions */
2222 } /* End loop on SC / FT */
2223
2224
2225 /*
2226 * Create the test image for SC (only used for debug)
2227 */
2228 cpl_table * wavedata_table = gravi_data_get_wave_data (wave_map, GRAVI_SC);
2229 cpl_table * wavefibre_table = gravi_data_get_wave_fibre (wave_map, GRAVI_SC);
2230 cpl_table * profile_table = gravi_data_get_table (profile_map, GRAVI_PROFILE_DATA_EXT);
2231 cpl_table * detector_table = gravi_data_get_imaging_detector (wave_map, GRAVI_SC);
2232
2233 cpl_imagelist * testwave_imglist;
2234 testwave_imglist = gravi_wave_test_image (wavedata_table, wavefibre_table,
2235 profile_table, detector_table);
2236 CPLCHECK_MSG ("Cannot compute TEST_WAVE");
2237
2238 /* Set TEST_WAVE in output */
2239 cpl_propertylist * plist = gravi_data_get_plist (profile_map, GRAVI_PROFILE_DATA_EXT);
2240 gravi_data_add_cube (wave_map, cpl_propertylist_duplicate (plist),
2241 "TEST_WAVE", testwave_imglist);
2242
2243 CPLCHECK_MSG ("Cannot set data");
2244
2245 gravi_msg_function_exit(1);
2246 return CPL_ERROR_NONE;
2247}
2248
2249/*----------------------------------------------------------------------------*/
2277/*----------------------------------------------------------------------------*/
2278
2279cpl_error_code gravi_compute_wave (gravi_data * wave_map,
2280 gravi_data * spectrum_data,
2281 int type_data, const cpl_parameterlist * parlist,
2282 gravi_data * wave_param)
2283{
2284 gravi_msg_function_start(1);
2285 cpl_ensure_code (wave_map, CPL_ERROR_NULL_INPUT);
2286 cpl_ensure_code (spectrum_data, CPL_ERROR_NULL_INPUT);
2287 cpl_ensure_code (type_data == GRAVI_SC || type_data == GRAVI_FT,
2288 CPL_ERROR_ILLEGAL_INPUT);
2289
2290 /* Get headers */
2291 cpl_propertylist * raw_header = gravi_data_get_header (spectrum_data);
2292 cpl_propertylist * wave_header = gravi_data_get_header (wave_map);
2293
2294 /* Dump full header of wave data */
2295 if (type_data == GRAVI_FT) {
2296 cpl_propertylist_append (wave_header, raw_header);
2297 }
2298
2299
2300 /* Copy IMAGING_DETECTOR in output WAVE */
2301 const char * extname = (type_data == GRAVI_SC) ? GRAVI_IMAGING_DETECTOR_SC_EXT : GRAVI_IMAGING_DETECTOR_FT_EXT;
2302 gravi_data_copy_ext (wave_map, spectrum_data, extname);
2303
2304 /*
2305 * Compute WAVE_FIBRE map
2306 */
2307 cpl_msg_info (cpl_func, "Compute WAVE_FIBRE for %s", GRAVI_TYPE(type_data));
2308
2309 cpl_table * spectrum_table = gravi_data_get_spectrum_data (spectrum_data, type_data);
2310 cpl_table * detector_table = gravi_data_get_imaging_detector (spectrum_data, type_data);
2311 cpl_table * opd_table = gravi_data_get_table (spectrum_data, type_data==GRAVI_SC ? "OPD_SC" : "OPD_FT");
2312
2313 cpl_table * wavefibre_table;
2314 wavefibre_table = gravi_wave_fibre (spectrum_table, detector_table, opd_table);
2315 CPLCHECK_MSG ("Cannot compute the WAVE_FIBRE");
2316
2317 /*
2318 * Correct WAVE_FIBRE from the dispersion model
2319 */
2320 cpl_msg_info (cpl_func, "Correct dispersion in WAVE_FIBRE of %s", GRAVI_TYPE(type_data));
2321
2322 /* Hardcoded values to correct for the fiber dispersion */
2323 double n0 = 1.0, n1 = -0.0165448, n2 = 0.00256002;
2324 cpl_msg_info (cpl_func,"Rescale wavelengths with dispersion (%g,%g,%g)",n0,n1,n2);
2325
2326 cpl_propertylist_update_string (wave_header, "ESO QC WAVE_CORR", "lbd*(N0+N1*(lbd-lbd0)/lbd0+N2*(lbd-lbd0)^2/lbd0^2)");
2327 cpl_propertylist_update_double (wave_header, "ESO QC WAVE_CORR N0", n0);
2328 cpl_propertylist_update_double (wave_header, "ESO QC WAVE_CORR N1", n1);
2329 cpl_propertylist_update_double (wave_header, "ESO QC WAVE_CORR N2", n2);
2330 CPLCHECK_MSG ("Cannot set Keywords");
2331
2332 gravi_wave_correct_dispersion (wavefibre_table, n0, n1, n2);
2333 CPLCHECK_MSG ("Cannot correct dispersion in wave");
2334
2335 /*
2336 * Fit the 2d dispersion WAVE_FIBRE into WAVE_DATA
2337 */
2338 cpl_msg_info (cpl_func,"Fit WAVE_DATA map for %s", GRAVI_TYPE(type_data));
2339
2340 /* Get the fullstartx */
2341 cpl_propertylist * spectrum_plist = gravi_data_get_spectrum_data_plist (spectrum_data, type_data);
2342 int fullstartx = gravi_pfits_get_fullstartx (spectrum_plist);
2343
2344 /* Interpolate table 2D */
2345 cpl_table * wavedata_table;
2346 int spatial_order=2; // default spatial order
2347 int spectral_order=3; // default spectral order
2348 double rms_residuals;
2349 if (type_data == GRAVI_FT && gravi_param_get_bool(parlist, "gravity.calib.force-wave-ft-equal")) {
2350 spatial_order = 0;
2351 cpl_msg_info (cpl_func, "Option force-waveFT-equal applied");
2352 }
2353 // Keep default value
2354 // if (type_data == GRAVI_SC) {
2355 // spectral_order = gravi_param_get_int(parlist, "gravity.calib.wave-spectral-order");
2356 // cpl_msg_info (cpl_func, "Option set_spectral order to %d", spectral_order);
2357 // }
2358
2359 wavedata_table = gravi_wave_fit_2d (wavefibre_table,
2360 detector_table,
2361 wave_param,
2362 fullstartx, spatial_order, spectral_order, &rms_residuals);
2363 cpl_propertylist_update_double (wave_header, QC_RMS_RESIDUALS(type_data), rms_residuals);
2364 cpl_propertylist_update_double (wave_header, QC_RMS_RESIDUALS_UM(type_data), rms_residuals * 1e6);
2365
2366
2367 double rms_fit=cpl_propertylist_get_double (raw_header, QC_PHASECHI2);
2368 cpl_propertylist_update_double (wave_header, QC_CHI2WAVE(type_data),rms_fit*1e9);
2369 cpl_propertylist_set_comment (wave_header, QC_CHI2WAVE(type_data),"[nm]rms a.SC-b.FT+c=MET");
2370
2371 CPLCHECK_MSG ("Cannot fit 2d data");
2372
2373 /*
2374 * Addtional wavelength calibration for LOW mode.
2375 * Default method: "New Wavelength interpolation made by sylvestre on January 30 2018"
2376 * New method: fitting using empirically measured bandpass.
2377 */
2378 if (type_data == GRAVI_SC && !strcmp (gravi_pfits_get_spec_res(raw_header), "LOW"))
2379 {
2380 char *wave_mode = cpl_strdup(cpl_parameter_get_string(
2381 cpl_parameterlist_find_const(parlist, "gravity.calib.wave-mode")));
2382 /* just in case: tolerate mixed/lower case arguments */
2383 for(int i = 0; wave_mode[i]; i++)
2384 wave_mode[i] = toupper(wave_mode[i]);
2385
2386 if (!strcmp(wave_mode, "INDIVIDUAL")) {
2387 cpl_msg_info (cpl_func, "Wave calibration using independent polynomial fit");
2388 cpl_table * wave_individual_table = cpl_table_new (1);
2389 cpl_table * weight_individual_table = cpl_table_new (1);
2390 cpl_table * wave_fitted_table = cpl_table_new (1);
2391
2392 gravi_wave_fit_individual (wave_individual_table,
2393 weight_individual_table,
2394 wave_fitted_table,
2395 opd_table,
2396 spectrum_table,
2397 detector_table,
2398 n0,n1,n2);
2399
2400 cpl_msg_info (cpl_func,"Add tables in wave_map");
2401
2402 gravi_data_add_table (wave_map, cpl_propertylist_duplicate (spectrum_plist),
2403 "WAVE_INDIV_SC", wave_individual_table);
2404 gravi_data_add_table (wave_map, cpl_propertylist_duplicate (spectrum_plist),
2405 "WAVE_WEIGHT_SC", weight_individual_table);
2406 gravi_data_add_table (wave_map, cpl_propertylist_duplicate (spectrum_plist),
2407 "WAVE_FITTED_SC", wavedata_table);
2408
2409 gravi_data_add_table (wave_map, cpl_propertylist_duplicate (spectrum_plist),
2410 GRAVI_WAVE_FIBRE_EXT(type_data), wavefibre_table);
2411 gravi_data_add_table (wave_map, cpl_propertylist_duplicate (spectrum_plist),
2412 GRAVI_WAVE_DATA_EXT(type_data), wave_fitted_table);
2413 } else if (!strcmp(wave_mode, "BANDPASS")) {
2414 cpl_msg_info (cpl_func, "Wave calibration using bandpass");
2415
2416 cpl_table *wave_bandpass_table = gravi_wave_fit_bandpass(
2417 spectrum_table, detector_table, opd_table);
2418
2419 cpl_msg_info (cpl_func,"Add tables in wave_map");
2420
2421 gravi_data_add_table (wave_map, cpl_propertylist_duplicate (spectrum_plist),
2422 "WAVE_FITTED_SC", wavedata_table);
2423
2424 gravi_data_add_table (wave_map, cpl_propertylist_duplicate (spectrum_plist),
2425 GRAVI_WAVE_FIBRE_EXT(type_data), wavefibre_table);
2426 gravi_data_add_table (wave_map, cpl_propertylist_duplicate (spectrum_plist),
2427 GRAVI_WAVE_DATA_EXT(type_data), wave_bandpass_table);
2428 }
2429
2430 cpl_free(wave_mode);
2431 } else {
2432 /*
2433 * Add the WAVE_FIBRE and WAVE_DATA table in the wave_map
2434 */
2435 cpl_msg_info (cpl_func,"Add WAVE_FIBRE and WAVE_DATA in wave_map");
2436
2437 gravi_data_add_table (wave_map, cpl_propertylist_duplicate (spectrum_plist),
2438 GRAVI_WAVE_FIBRE_EXT(type_data), wavefibre_table);
2439
2440 gravi_data_add_table (wave_map, cpl_propertylist_duplicate (spectrum_plist),
2441 GRAVI_WAVE_DATA_EXT(type_data), wavedata_table);
2442
2443 }
2444
2445 CPLCHECK_MSG ("Cannot set data");
2446
2447 gravi_msg_function_exit(1);
2448 return CPL_ERROR_NONE;
2449}
2450
2451
GRAVI_PRIMARY_HDR_EXT
#define GRAVI_PRIMARY_HDR_EXT
Definition gravi-test.h:212
gravi_calib.h
gravi_cpl.h
gravi_data.h
gravi_data
typedefCPL_BEGIN_DECLS struct _gravi_data_ gravi_data
Definition gravi_data.h:39
gravi_data_get_spectrum_data
#define gravi_data_get_spectrum_data(data, type)
Definition gravi_data.h:63
gravi_data_get_wave_data_plist
#define gravi_data_get_wave_data_plist(data, type)
Definition gravi_data.h:54
gravi_data_get_wave_data
#define gravi_data_get_wave_data(data, type)
Definition gravi_data.h:53
gravi_data_get_spectrum_data_plist
#define gravi_data_get_spectrum_data_plist(data, type)
Definition gravi_data.h:64
gravi_data_get_header
#define gravi_data_get_header(data)
Definition gravi_data.h:75
gravi_data_get_oi_wave
#define gravi_data_get_oi_wave(data, type, pol, npol)
Definition gravi_data.h:45
gravi_data_get_imaging_detector
#define gravi_data_get_imaging_detector(data, type)
Definition gravi_data.h:60
gravi_data_get_wave_fibre
#define gravi_data_get_wave_fibre(data, type)
Definition gravi_data.h:51
gravi_data_get_oi_wave_plist
#define gravi_data_get_oi_wave_plist(data, type, pol, npol)
Definition gravi_data.h:69
ntel
const cpl_size ntel
Definition gravi_demodulate.c:84
gravi_dfs.h
gravi_ellipse.h
USE_LINEAR_ELLIPSE
#define USE_LINEAR_ELLIPSE
Definition gravi_ellipse.h:35
gravi_metrology.h
wave_data
gravi_data * wave_data
Definition gravi_old.c:1994
cpl_msg_debug
cpl_msg_debug(cpl_func, "Spectra has <50 pixels -> don't flat")
plist
cpl_propertylist * plist
Definition gravi_old.c:2000
cpl_msg_info
cpl_msg_info(cpl_func, "Compute WAVE_SCAN for %s", GRAVI_TYPE(type_data))
cpl_propertylist_update_double
cpl_propertylist_update_double(header, "ESO QC MINWAVE SC", cpl_propertylist_get_double(plist, "ESO QC MINWAVE SC"))
wave_data_sc
cpl_table * wave_data_sc
Definition gravi_old.c:1997
gravi_vis_create_met_ft
gravi_vis_create_met_ft(phase_ft, vis_met)
gravi_pfits.h
GRAVI_OI_VIS_MET_EXT
#define GRAVI_OI_VIS_MET_EXT
Definition gravi_pfits.h:89
GRAVI_IMAGING_DETECTOR_FT_EXT
#define GRAVI_IMAGING_DETECTOR_FT_EXT
Definition gravi_pfits.h:50
GRAVI_SC
#define GRAVI_SC
Definition gravi_pfits.h:166
QC_MAXWAVE
#define QC_MAXWAVE(type)
Definition gravi_pfits.h:111
QC_MINWAVE_UM
#define QC_MINWAVE_UM(type)
Definition gravi_pfits.h:112
GRAVI_PROFILE_DATA_EXT
#define GRAVI_PROFILE_DATA_EXT
Definition gravi_pfits.h:80
GRAVI_WAVE_FIBRE_EXT
#define GRAVI_WAVE_FIBRE_EXT(type)
Definition gravi_pfits.h:70
QC_MAXWAVE_UM
#define QC_MAXWAVE_UM(type)
Definition gravi_pfits.h:113
OPD_COEFF_SIGN
#define OPD_COEFF_SIGN(type)
Definition gravi_pfits.h:116
QC_CHI2WAVE
#define QC_CHI2WAVE(type)
Definition gravi_pfits.h:109
GRAVI_WAVE_DATA_EXT
#define GRAVI_WAVE_DATA_EXT(type)
Definition gravi_pfits.h:67
QC_RMS_RESIDUALS_UM
#define QC_RMS_RESIDUALS_UM(type)
Definition gravi_pfits.h:115
GRAVI_TYPE
#define GRAVI_TYPE(type)
Definition gravi_pfits.h:168
QC_MINWAVE
#define QC_MINWAVE(type)
Definition gravi_pfits.h:110
LAMBDA_MET
#define LAMBDA_MET
Definition gravi_pfits.h:104
GRAVI_IMAGING_DETECTOR_SC_EXT
#define GRAVI_IMAGING_DETECTOR_SC_EXT
Definition gravi_pfits.h:49
GRAVI_FT
#define GRAVI_FT
Definition gravi_pfits.h:167
QC_PHASECHI2
#define QC_PHASECHI2
Definition gravi_pfits.h:108
QC_RMS_RESIDUALS
#define QC_RMS_RESIDUALS(type)
Definition gravi_pfits.h:114
gravi_preproc.h
gravi_signal.h
gravi_utils.h
GRAVI_POLAR
#define GRAVI_POLAR(pol, npol)
Definition gravi_utils.h:144
gravi_spectrum_get_npol
#define gravi_spectrum_get_npol(table)
Definition gravi_utils.h:92
gravi_msg_function_exit
#define gravi_msg_function_exit(flag)
Definition gravi_utils.h:85
FREE
#define FREE(function, variable)
Definition gravi_utils.h:69
CPLCHECK_NUL
#define CPLCHECK_NUL(msg)
Definition gravi_utils.h:48
gravi_msg_function_start
#define gravi_msg_function_start(flag)
Definition gravi_utils.h:84
CPLCHECK_MSG
#define CPLCHECK_MSG(msg)
Definition gravi_utils.h:45
GRAVI_LBD_FTSC
#define GRAVI_LBD_FTSC
Definition gravi_utils.h:116
FREELOOP
#define FREELOOP(function, variable, n)
Definition gravi_utils.h:72
GRAVI_NBASE
#define GRAVI_NBASE
Definition gravi_utils.h:105
gravi_wave.h
gravi_table_get_column_mean
double gravi_table_get_column_mean(cpl_table *table, const char *name, int base, int nbase)
Definition gravi_cpl.c:342
gravi_table_new_column
cpl_error_code gravi_table_new_column(cpl_table *table, const char *name, const char *unit, cpl_type type)
Definition gravi_cpl.c:1655
gravi_vector_unwrap_with_guess
cpl_error_code gravi_vector_unwrap_with_guess(cpl_vector *vector, cpl_vector *ref, double ref_to_phase)
Unwrap a phase vector following a guess vector. The difference is actually unwrap and shall thus be s...
Definition gravi_cpl.c:2718
gravi_imagelist_wrap_column
cpl_imagelist * gravi_imagelist_wrap_column(cpl_table *table_data, const char *data_x)
Wrap a column array of a table into an imagelist.
Definition gravi_cpl.c:1756
gravi_imagelist_unwrap_images
cpl_error_code gravi_imagelist_unwrap_images(cpl_imagelist *imglist)
Unwrap an imagelist an all its images.
Definition gravi_cpl.c:1726
gravi_table_get_vector
cpl_vector * gravi_table_get_vector(cpl_table *spectrum_data, cpl_size index, const char *regname)
Create a vector from the row index of the column regname.
Definition gravi_cpl.c:2234
gravi_table_add_scalar
cpl_error_code gravi_table_add_scalar(cpl_table *table, const char *name, int base, int nbase, double value)
Multply scalar or array column by scalar.
Definition gravi_cpl.c:2815
gravi_table_new_column_array
cpl_error_code gravi_table_new_column_array(cpl_table *table, const char *name, const char *unit, cpl_type type, cpl_size size)
Definition gravi_cpl.c:1677
gravi_data_add_cube
cpl_error_code gravi_data_add_cube(gravi_data *self, cpl_propertylist *plist, const char *extname, cpl_imagelist *imglist)
Add an IMAGE (imagelist) extension in gravi_data.
Definition gravi_data.c:2353
gravi_data_get_plist
cpl_propertylist * gravi_data_get_plist(gravi_data *self, const char *extname)
Get the propertylist from EXTNAME.
Definition gravi_data.c:2049
gravi_data_add_table
cpl_error_code gravi_data_add_table(gravi_data *self, cpl_propertylist *plist, const char *extname, cpl_table *table)
Add a BINTABLE extension in gravi_data.
Definition gravi_data.c:2289
gravi_data_get_table
cpl_table * gravi_data_get_table(gravi_data *self, const char *extname)
Return a pointer on a table extension by its EXTNAME.
Definition gravi_data.c:2096
gravi_data_copy_ext
cpl_error_code gravi_data_copy_ext(gravi_data *output, gravi_data *input, const char *name)
Copy extensions from one data to another.
Definition gravi_data.c:1690
gravi_param_get_bool
int gravi_param_get_bool(const cpl_parameterlist *parlist, const char *name)
Definition gravi_dfs.c:1532
gravi_ellipse_meanopd_create
cpl_vector * gravi_ellipse_meanopd_create(cpl_table *spectrum_table, cpl_table *detector_table, cpl_table **oiwave_tables, cpl_vector *guess_vector, int base)
Compute the OPD modulation of a baseline from spectrum.
Definition gravi_ellipse.c:351
gravi_ellipse_phase_create
cpl_vector * gravi_ellipse_phase_create(cpl_vector *vectCA, cpl_vector *vectDB, cpl_vector *envelope)
Compute the phase atan{X',Y'}, unwraped from first sample.
Definition gravi_ellipse.c:107
gravi_metrology_create
cpl_table * gravi_metrology_create(cpl_table *metrology_table, cpl_propertylist *header)
Create the VIS_MET table.
Definition gravi_metrology.c:1734
gravi_metrology_drs
cpl_error_code gravi_metrology_drs(cpl_table *metrology_table, cpl_table *vismet_table, cpl_propertylist *header, const cpl_parameterlist *parlist)
Fill the VIS_MET table with the DRS algorithm.
Definition gravi_metrology.c:1957
gravi_metrology_tac
cpl_error_code gravi_metrology_tac(cpl_table *metrology_table, cpl_table *vismet_table, cpl_propertylist *header)
Compute the metrology signal from TAC algorithm.
Definition gravi_metrology.c:2464
gravi_pfits_get_pola_num
int gravi_pfits_get_pola_num(const cpl_propertylist *plist, int type_data)
Definition gravi_pfits.c:263
gravi_pfits_get_fullstartx
int gravi_pfits_get_fullstartx(const cpl_propertylist *plist)
Definition gravi_pfits.c:76
gravi_pfits_add_check
cpl_error_code gravi_pfits_add_check(cpl_propertylist *header, const char *msg)
Add a QC.CHECK keyword to the header.
Definition gravi_pfits.c:1682
gravi_pfits_get_spec_res
const char * gravi_pfits_get_spec_res(const cpl_propertylist *plist)
Definition gravi_pfits.c:162
gravi_pfits_get_dit_sc
double gravi_pfits_get_dit_sc(const cpl_propertylist *plist)
Definition gravi_pfits.c:664
gravi_pfits_get_met_wavelength_mean
double gravi_pfits_get_met_wavelength_mean(const cpl_propertylist *plist, cpl_table *met_table)
Definition gravi_pfits.c:312
gravi_pfits_get_double_default
double gravi_pfits_get_double_default(const cpl_propertylist *plist, const char *name, double def)
Definition gravi_pfits.c:1626
gravi_vis_create_opdsc_ft
cpl_error_code gravi_vis_create_opdsc_ft(cpl_table *vis_FT, cpl_table *vis_SC, double dit_sc)
Compute the resampled SC signal for each FT DIT per base.
Definition gravi_signal.c:1843
gravi_region_get_pol
int gravi_region_get_pol(cpl_table *imaging_detector, int region)
Return the polarisation id of a region.
Definition gravi_utils.c:470
gravi_region_get_base
int gravi_region_get_base(cpl_table *imaging_detector, int region)
Return the base of a region.
Definition gravi_utils.c:386
GRAVI_BASE_NAME
char GRAVI_BASE_NAME[6][3]
Definition gravi_utils.c:57
gravi_spectrum_get_nwave
cpl_size gravi_spectrum_get_nwave(const cpl_table *table)
Definition gravi_utils.c:1038
GRAVI_DATA
char GRAVI_DATA[50][7]
Definition gravi_utils.c:71
gravi_spectrum_get_nregion
cpl_size gravi_spectrum_get_nregion(const cpl_table *table)
Definition gravi_utils.c:1043
gravi_compute_envelope
cpl_vector * gravi_compute_envelope(const cpl_vector *opd, int wave, int nwave)
Compute the envelope value.
Definition gravi_utils.c:1106
GRAVI_BASE_TEL
int GRAVI_BASE_TEL[GRAVI_NBASE][2]
Definition gravi_utils.c:56
gravi_get_region
int gravi_get_region(cpl_table *img_det, int base, char phase, int pol)
Find the region matching base, phase and pol.
Definition gravi_utils.c:609
gravi_region_get_base_sign
int gravi_region_get_base_sign(cpl_table *imaging_detector, int base)
Return the sign of a base by looking at the PORT order.
Definition gravi_utils.c:424
GRAVI_CLO_BASE
int GRAVI_CLO_BASE[4][3]
Definition gravi_utils.c:68
WAVE_TO_PLOT
#define WAVE_TO_PLOT
Definition gravi_wave.c:80
gravi_opds_calibration
cpl_table * gravi_opds_calibration(cpl_table *spectrum_table, cpl_table *detector_table, cpl_table *opdguess_table)
Compute the mean opd of each baseline from spectrum.
Definition gravi_wave.c:704
gravi_opds_fit_opdmet
cpl_vector * gravi_opds_fit_opdmet(cpl_table *opd_ft, double lbd_met)
Compute the absolute scaling coefficients of SC and FT OPDs.
Definition gravi_wave.c:465
gravi_wave_fibre
cpl_table * gravi_wave_fibre(cpl_table *spectrum_table, cpl_table *detector_table, cpl_table *opd_table)
Compute the wavelength of each channel for each 6 baselines.
Definition gravi_wave.c:995
gravi_wave_fit_2d
cpl_table * gravi_wave_fit_2d(cpl_table *wavefibre_table, cpl_table *detector_table, gravi_data *wave_param, cpl_size fullstartx, int spatial_order, int spectral_order, double *rms_residuals)
Compute the WAVE_DATA table (1 wavelength per region) from the WAVE_FIBRE table (1 wavelength per bas...
Definition gravi_wave.c:1185
PLOT_WAVE_PHASE_VS_OPD
#define PLOT_WAVE_PHASE_VS_OPD
Definition gravi_wave.c:79
CPL_ATTR_ALLOC
CPL_ATTR_ALLOC
Definition gravi_wave.c:1480
my_ft_discrete
cpl_array * my_ft_discrete(const cpl_vector *x, const cpl_vector *y, const cpl_vector *freqs)
Compute discrete fourier transform of input data.
Definition gravi_wave.c:1436
gravi_wave_compute_opds
cpl_error_code gravi_wave_compute_opds(gravi_data *spectrum_data, cpl_table *met_table, const cpl_parameterlist *parlist)
Recover the OPD modulation from a spectrum_data and vismet_table.
Definition gravi_wave.c:826
gravi_opds_correct_closures
cpl_error_code gravi_opds_correct_closures(cpl_table *opd_sc, const char *name)
Correct the input table to ensure constant closure phase.
Definition gravi_wave.c:360
gravi_compute_wave
cpl_error_code gravi_compute_wave(gravi_data *wave_map, gravi_data *spectrum_data, int type_data, const cpl_parameterlist *parlist, gravi_data *wave_param)
Create the WAVE calibration map.
Definition gravi_wave.c:2279
gravi_compute_argon_wave
cpl_table * gravi_compute_argon_wave(cpl_table *spectrum_table)
Compute a WAVE calibration from the ARGON data (SC only)
Definition gravi_wave.c:155
gravi_wave_fit_individual
cpl_table * gravi_wave_fit_individual(cpl_table *wave_individual_table, cpl_table *weight_individual_table, cpl_table *wave_fitted_table, cpl_table *opd_table, cpl_table *spectrum_table, cpl_table *detector_table, double n0, double n1, double n2)
Definition gravi_wave.c:1653
gravi_wave_qc
cpl_error_code gravi_wave_qc(gravi_data *wave_map, gravi_data *profile_map)
Compute the QC parameters of the WAVE product.
Definition gravi_wave.c:2132
gravi_wave_correct_dispersion
cpl_error_code gravi_wave_correct_dispersion(cpl_table *wave_fibre, double n0, double n1, double n2)
Correct the WAVE_FIBRE table from a harcoded dispersion model.
Definition gravi_wave.c:1929
gravi_wave_correct_color
cpl_error_code gravi_wave_correct_color(gravi_data *vis_data)
Create a OI_WAVELENGTH_CORR table with color corrected wavelength.
Definition gravi_wave.c:1974
gravi_wave_test_image
cpl_imagelist * gravi_wave_test_image(cpl_table *wavedata_table, cpl_table *wavefibre_table, cpl_table *profile_table, cpl_table *detector_table)
Compute the (useless) TEST_WAVE table from the WAVE_FIBRE and the PROFILE maps.
Definition gravi_wave.c:2027
gravi_wave_fit_bandpass
cpl_table * gravi_wave_fit_bandpass(cpl_table *spectrum_table, cpl_table *detector_table, cpl_table *opd_table)
Compute wavelength calibration using empirical bandpass.
Definition gravi_wave.c:1494
gravi_opds_compute_guess
cpl_table * gravi_opds_compute_guess(cpl_table *spectrumsc_table, cpl_table *opdft_table, cpl_table *met_table, double dit_sc, double lbd_met)
Compute a guess of the OPD modulation of SC from FT and MET.
Definition gravi_wave.c:572
Generated on Wed Jul 29 2026 08:27:42 for GRAVI Pipeline Reference Manual by doxygen 1.13.2