GRAVI Pipeline Reference Manual 1.10.1
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 <complex.h>
59
60#include "gravi_data.h"
61#include "gravi_dfs.h"
62#include "gravi_pfits.h"
63#include "gravi_cpl.h"
64
65#include "gravi_utils.h"
66#include "gravi_ellipse.h"
67#include "gravi_signal.h"
68
69#include "gravi_preproc.h"
70#include "gravi_calib.h"
71#include "gravi_wave.h"
72#include "gravi_metrology.h"
73
74/*----------------------------------------------------------------------------
75 DEBUG
76 -----------------------------------------------------------------------------*/
77
78#define PLOT_WAVE_PHASE_VS_OPD 0
79#define WAVE_TO_PLOT 5
80
81/*-----------------------------------------------------------------------------
82 Private prototypes
83 -----------------------------------------------------------------------------*/
84
85cpl_table * gravi_opds_compute_guess (cpl_table * spectrumsc_table,
86 cpl_table * opdft_table,
87 cpl_table * met_table,
88 double dit_sc,
89 double lbd_met);
90
91cpl_table * gravi_opds_calibration (cpl_table * spectrum_table,
92 cpl_table * detector_table,
93 cpl_table * opdguess_table);
94
95cpl_error_code gravi_opds_correct_closures (cpl_table * opd_sc,
96 const char *name);
97
98cpl_vector * gravi_opds_fit_opdmet (cpl_table *opd_ft, double lbd_met);
99
100cpl_table * gravi_wave_fibre (cpl_table * spectrum_table,
101 cpl_table * detector_table,
102 cpl_table * opd_table);
103
104cpl_table * gravi_wave_fit_2d (cpl_table * wavefibre_table,
105 cpl_table * detector_table,
106 gravi_data * wave_param,
107 cpl_size fullstartx,
108 int spatial_order,
109 int spectral_order,
110 double * rms_residuals);
111
112cpl_table * gravi_wave_fit_individual (cpl_table * wave_individual_table,
113 cpl_table * weight_individual_table,
114 cpl_table * wave_fitted_table,
115 cpl_table * opd_table,
116 cpl_table * spectrum_table,
117 cpl_table * detector_table,
118 double n0, double n1, double n2);
119
120cpl_error_code gravi_wave_correct_dispersion (cpl_table * wave_fibre,
121 double n0, double n1, double n2);
122
123cpl_imagelist * gravi_wave_test_image (cpl_table * wavedata_table,
124 cpl_table * wavefibre_table,
125 cpl_table * profile_table,
126 cpl_table * detector_table);
127
128
129
130
131/*-----------------------------------------------------------------------------
132 Functions code
133 -----------------------------------------------------------------------------*/
134
135/*----------------------------------------------------------------------------*/
146/*----------------------------------------------------------------------------*/
147
148cpl_table * gravi_compute_argon_wave (cpl_table * spectrum_table)
149{
150 gravi_msg_function_start(1);
151 cpl_ensure (spectrum_table, CPL_ERROR_NULL_INPUT, NULL);
152
153 /* Get the SC spectrums */
154 cpl_size nwave = cpl_table_get_column_depth (spectrum_table,"DATA1");
155 cpl_size nregion = gravi_spectrum_get_nregion (spectrum_table);
156 CPLCHECK_NUL ("Cannot get data");
157
158 cpl_msg_info (cpl_func, "spectrum_table = %lld",cpl_table_get_column_depth (spectrum_table, "DATA1"));
159
160 /* Prepare the outputs */
161 cpl_table * wave_data_sc = cpl_table_new (1);
162 CPLCHECK_NUL ("Cannot create output");
163
164 /* Init the list of lines*/
165 cpl_size nlines = 14;
166 double plines_str[] = {60,
167 160,
168 90,
169 1800,
170 50,
171 30,
172 30,
173 1500,
174 130,
175 510,
176 240,
177 169,
178 64,
179 73
180 };
181
182 double plines_lbd[] = {1.982291,
183 1.997118,
184 2.032256,
185 2.062186,
186 2.065277,
187 2.073922,
188 2.081672,
189 2.099184,
190 2.133871,
191 2.154009,
192 2.208321,
193 2.313952,
194 2.385154,
195 2.397306
196 };
197
198 cpl_vector * lines_lbd = cpl_vector_wrap (nlines,plines_lbd);
199 cpl_vector * lines_str = cpl_vector_wrap (nlines,plines_str);
200 cpl_bivector * lines = cpl_bivector_wrap_vectors (lines_lbd, lines_str);
201
202 /* Init the slit and lamp model */
203 cpl_wlcalib_slitmodel * model = cpl_wlcalib_slitmodel_new ();
204 cpl_wlcalib_slitmodel_set_wslit (model, 0.1);
205 cpl_wlcalib_slitmodel_set_wfwhm (model, 1.25);
206 cpl_wlcalib_slitmodel_set_threshold (model, 5.0);
207 cpl_wlcalib_slitmodel_set_catalog (model, lines);
208
209 /* Starting point of dispersion */
210 double values_hr[] = {1.97048,0.2228e-3,2.7728e-08,-4.72e-12};
211 double values_mr[] = {1.9697,0.2179e-2,2.5e-07,-5e-10};
212 double * values = nwave > 1000 ? values_hr : values_mr;
213
214 /* Fill the dispersion polynomial */
215 cpl_polynomial * dispersion0 = cpl_polynomial_new (1);
216 for (cpl_size order = 0 ; order < 4 ; order ++) {
217 cpl_polynomial_set_coeff (dispersion0, &order, values[order]);
218 }
219
220 /* QC parameters */
221 double minwave = -1e10;
222 double maxwave = +1e10;
223
224 cpl_vector * spectrum = cpl_vector_new (nwave);
225 cpl_vector * input_spectrum = cpl_vector_new (nwave);
226 cpl_vector * model_spectrum = cpl_vector_new (nwave);
227
228 /* Prepare output table */
229 cpl_table * fit_table = cpl_table_new (nregion);
230 gravi_table_new_column_array (fit_table, "DATA", "adu", CPL_TYPE_DOUBLE, nwave);
231 gravi_table_new_column_array (fit_table, "DATA_MODEL", "adu", CPL_TYPE_DOUBLE, nwave);
232 gravi_table_new_column_array (fit_table, "DATA_WAVE", "adu", CPL_TYPE_DOUBLE, nwave);
233
234 /* The starting point will be the fit of previous region */
235 cpl_polynomial * dispersion = cpl_polynomial_duplicate (dispersion0);
236
237 /* Loop on regions */
238 for (cpl_size reg = 0; reg < nregion ; reg ++ ) {
239 cpl_msg_info (cpl_func,"Fit region %lld over %lld", reg+1, nregion);
240
241 /* Copy the spectrum of this region into a vector */
242 const cpl_array * data = cpl_table_get_array (spectrum_table, GRAVI_DATA[reg], 0);
243 for (cpl_size wave=0; wave < nwave; wave ++) {
244 cpl_vector_set (spectrum, wave, log (CPL_MAX (cpl_array_get (data, wave, NULL), 1.0)));
245 cpl_vector_set (input_spectrum, wave, cpl_array_get (data, wave, NULL));
246 }
247
248 CPLCHECK_NUL ("Cannot prepare data");
249
250 /* Parameters: Explore a large hsize [pixel]
251 * only for the first region, since we keep
252 * the starting point */
253 int maxdeg = 3, nmaxima = 0, linelim = 99;
254 int maxite = 100 * maxdeg, maxfail = 10, maxcont = 10;
255 int hsize = (reg == 0 ? 40 : 5);
256 double pixtol = 1e-6, pixstep = 0.5, pxc = 0.0;
257
258 /* Use the log of spectrum to put
259 * less weight on the intensity */
260 irplib_polynomial_find_1d_from_correlation_all (dispersion, maxdeg, spectrum,
261 nmaxima, linelim,
262 (irplib_base_spectrum_model*)model,
263 &irplib_vector_fill_logline_spectrum,
264 pixtol, pixstep,
265 hsize, maxite, maxfail, maxcont,
266 CPL_FALSE,&pxc);
267 CPLCHECK_NUL ("Cannot fit polynomial");
268
269 /* Fill in the final fit */
270 irplib_vector_fill_line_spectrum (model_spectrum, dispersion, (void *)model);
271 cpl_polynomial_dump (dispersion, stdout);
272
273 /* Evaluate the polynomial for the output */
274 cpl_array * wave_data = cpl_array_new (nwave, CPL_TYPE_DOUBLE);
275 for (cpl_size pix = 0; pix < nwave ; pix++)
276 cpl_array_set (wave_data, pix, cpl_polynomial_eval_1d (dispersion, (double)pix, NULL) * 1e-6);
277
278 /* Check the limits for the QC */
279 minwave = CPL_MAX (cpl_array_get_min (wave_data), minwave);
280 maxwave = CPL_MIN (cpl_array_get_max (wave_data), maxwave);
281
282 /* Fill the WAVE_DATA_SC */
283 cpl_table_new_column_array (wave_data_sc, GRAVI_DATA[reg], CPL_TYPE_DOUBLE, nwave);
284 cpl_table_set_array (wave_data_sc, GRAVI_DATA[reg], 0, wave_data);
285 FREE (cpl_array_delete, wave_data);
286
287 /* Save the test tables */
288 cpl_array * tmp_array;
289 tmp_array = cpl_array_wrap_double (cpl_vector_get_data (input_spectrum), nwave);
290 cpl_table_set_array (fit_table, "DATA", reg, tmp_array);
291 cpl_array_unwrap (tmp_array);
292 tmp_array = cpl_array_wrap_double (cpl_vector_get_data (model_spectrum), nwave);
293 cpl_table_set_array (fit_table, "DATA_MODEL", reg, tmp_array);
294 cpl_array_unwrap (tmp_array);
295 tmp_array = cpl_array_new (nwave, CPL_TYPE_DOUBLE);
296 for (cpl_size wave = 0; wave < nwave; wave ++)
297 cpl_array_set (tmp_array, wave, cpl_polynomial_eval_1d (dispersion, (double)wave, NULL) * 1e-6);
298 cpl_table_set_array (fit_table, "DATA_WAVE", reg, tmp_array);
299 cpl_array_delete (tmp_array);
300
301 CPLCHECK_NUL ("Cannot set data");
302 } /* End loop on regions */
303
304 // cpl_table_save (fit_table, NULL, NULL, "fit_wave.fits", CPL_IO_CREATE);
305 FREE (cpl_table_delete, fit_table);
306
307 /* Delete */
308 FREE (cpl_vector_delete, spectrum);
309 FREE (cpl_vector_delete, input_spectrum);
310 FREE (cpl_vector_delete, model_spectrum);
311 FREE (cpl_polynomial_delete, dispersion0);
312 FREE (cpl_polynomial_delete, dispersion);
313
314 /* This memory desalocation does not work: */
315 // FREE (cpl_wlcalib_slitmodel_delete, model);
316
317 gravi_msg_function_exit(1);
318 return wave_data_sc;
319}
320
321/*----------------------------------------------------------------------------*/
351/*----------------------------------------------------------------------------*/
352
353cpl_error_code gravi_opds_correct_closures (cpl_table * phase_table,
354 const char *name)
355{
356 gravi_msg_function_start(1);
357 cpl_ensure_code (phase_table, CPL_ERROR_NULL_INPUT);
358 cpl_ensure_code (name, CPL_ERROR_NULL_INPUT);
359
360 int nclo = 4;
361 cpl_size nrow = cpl_table_get_nrow (phase_table) / GRAVI_NBASE;
362
363 /* Get the data */
364 double * phase = cpl_table_get_data_double (phase_table, name);
365 CPLCHECK_MSG ("Cannot get data");
366
367 /* Model and right-hand-side for the lineary system
368 * (unfilled matrix are 0.0) */
369 cpl_matrix * rhs_matrix = cpl_matrix_new (nrow * nclo, 1);
370 cpl_matrix * model_matrix = cpl_matrix_new (nrow * nclo, GRAVI_NBASE-1 + nclo);
371
372 for (cpl_size row = 0; row < nrow; row++) {
373 for (int clo = 0; clo < nclo; clo++) {
374 int b0 = GRAVI_CLO_BASE[clo][0];
375 int b1 = GRAVI_CLO_BASE[clo][1];
376 int b2 = GRAVI_CLO_BASE[clo][2];
377
378 /* Fill the rhs with measured closure */
379 double closure = phase[row*GRAVI_NBASE+b0] +
380 phase[row*GRAVI_NBASE+b1] -
381 phase[row*GRAVI_NBASE+b2];
382 cpl_matrix_set (rhs_matrix, row*nclo+clo, 0, closure);
383 CPLCHECK_MSG ("Cannot fill rhs_matrix");
384
385 /* Fill the k (unfilled are 0.0) */
386 cpl_matrix_set (model_matrix, row*nclo+clo, clo, 1.0);
387 CPLCHECK_MSG ("Cannot fill k in model_matrix");
388
389 /* Fill the f in of the model_matrix */
390 if (b0 < GRAVI_NBASE-1) cpl_matrix_set (model_matrix, row*nclo+clo, nclo+b0,
391 +phase[row*GRAVI_NBASE+b0]);
392 if (b1 < GRAVI_NBASE-1) cpl_matrix_set (model_matrix, row*nclo+clo, nclo+b1,
393 +phase[row*GRAVI_NBASE+b1]);
394 if (b2 < GRAVI_NBASE-1) cpl_matrix_set (model_matrix, row*nclo+clo, nclo+b2,
395 -phase[row*GRAVI_NBASE+b2]);
396 CPLCHECK_MSG ("Cannot fill phase in model_matrix");
397 }
398 } /* End loop on clo and rows */
399
400 /* Solve the system */
401 cpl_matrix * res_matrix = cpl_matrix_solve_normal (model_matrix, rhs_matrix);
402 CPLCHECK_MSG ("Cannot solve system");
403
404 /* Compute residuals */
405 cpl_matrix * residual_matrix = cpl_matrix_product_create (model_matrix, res_matrix);
406 cpl_matrix_subtract (residual_matrix, rhs_matrix);
407 double rms_fit = cpl_matrix_get_stdev (residual_matrix);
408
409 /* Correct baseline with (1-f). Last baseline is kept unchanged */
410 for (int base = 0; base < GRAVI_NBASE - 1; base ++) {
411 double f = cpl_matrix_get (res_matrix, nclo + base, 0);
412 cpl_msg_info (cpl_func,"correction factor f = 1 %+.20f", -1*f);
413 for (cpl_size row = 0; row < nrow; row++) phase[row*GRAVI_NBASE+base] *= 1 - f;
414 }
415
416 /* Dump the residual with their units */
417 const char * unit = cpl_table_get_column_unit (phase_table, name);
418 cpl_msg_info (cpl_func, "residual stdev = %.20g [%s]", rms_fit, unit);
419
420 /* Delete matrix */
421 FREE (cpl_matrix_delete, residual_matrix);
422 FREE (cpl_matrix_delete, res_matrix);
423 FREE (cpl_matrix_delete, model_matrix);
424 FREE (cpl_matrix_delete, rhs_matrix);
425
426 gravi_msg_function_exit(1);
427 return CPL_ERROR_NONE;
428}
429
430/*----------------------------------------------------------------------------*/
456/*----------------------------------------------------------------------------*/
457
458cpl_vector * gravi_opds_fit_opdmet (cpl_table * ft_table, double lbd_met)
459{
460 gravi_msg_function_start(1);
461 cpl_ensure (ft_table, CPL_ERROR_NULL_INPUT, NULL);
462
463
464 /* Get the number of acquisitions */
465 cpl_size nrow = cpl_table_get_nrow (ft_table) / GRAVI_NBASE;
466
467 /* Get the pointer to data */
468 double * opd_sc = cpl_table_get_data_double (ft_table, "OPD_SC");
469 double * opd_ft = cpl_table_get_data_double (ft_table, "OPD");
470 double * phase_met = cpl_table_get_data_double (ft_table, "PHASE_MET_FC");
471 CPLCHECK_NUL ("Cannot get data");
472
473 /* Number of valid rows */
474 cpl_size nrow_valid = 0;
475 for (cpl_size row = 0; row < nrow; row++) if (opd_sc[row*GRAVI_NBASE] != 0) nrow_valid++;
476 cpl_msg_info (cpl_func,"nrow_valid = %lld", nrow_valid);
477
478 cpl_ensure (nrow_valid > 100, CPL_ERROR_ILLEGAL_INPUT, NULL);
479
480 /* Model and right-hand-side for the lineary system
481 * (unfilled matrix are 0.0) */
482 cpl_matrix * rhs_matrix = cpl_matrix_new (nrow_valid*GRAVI_NBASE, 1);
483 cpl_matrix * model_matrix = cpl_matrix_new (nrow_valid*GRAVI_NBASE, GRAVI_NBASE + 2);
484
485 for (int base = 0; base < GRAVI_NBASE; base++) {
486 cpl_size row_valid = 0;
487
488 for (cpl_size row=0; row<nrow; row++) {
489 if (opd_sc[row*GRAVI_NBASE] == 0) continue;
490
491 int idv = row_valid * GRAVI_NBASE + base;
492 int id = row * GRAVI_NBASE + base;
493
494 /* Fill the OPD metrology */
495 cpl_matrix_set (rhs_matrix, idv, 0, phase_met[id] * lbd_met / CPL_MATH_2PI);
496 CPLCHECK_NUL ("Cannot set OPD");
497
498 /* Fill the model b and c */
499 cpl_matrix_set (model_matrix, idv, 0, 1*opd_sc[id]);
500 cpl_matrix_set (model_matrix, idv, 1, -1*opd_ft[id]);
501 CPLCHECK_NUL ("Cannot set SC or FT");
502
503 /* Fill the model Aij (unfilled matrix are 0.0) */
504 cpl_matrix_set (model_matrix, idv, 2 + base, 1.0);
505 CPLCHECK_NUL ("Cannot set the zero points");
506
507 row_valid++;
508 } /* End loop rows */
509 } /* End loop on base*/
510
511 /* Solve the system */
512 cpl_matrix * res_matrix = cpl_matrix_solve_normal (model_matrix, rhs_matrix);
513
514 /* Compute residuals */
515 cpl_matrix * residual_matrix = cpl_matrix_product_create (model_matrix, res_matrix);
516 cpl_matrix_subtract (residual_matrix, rhs_matrix);
517 double rms_fit = cpl_matrix_get_stdev (residual_matrix);
518 //cpl_plot_vector(NULL, NULL, NULL, cpl_vector_wrap(nrow_valid*GRAVI_NBASE,cpl_matrix_get_data(residual_matrix)));
519
520 /* Verbose */
521 cpl_msg_info (cpl_func, "coeff SC = %.20g ", cpl_matrix_get (res_matrix, 0, 0));
522 cpl_msg_info (cpl_func, "coeff FT = %.20g ", cpl_matrix_get (res_matrix, 1, 0));
523 cpl_msg_info (cpl_func, "residual stdev = %.20g [m]", rms_fit);
524
525 /* Set the Results */
526 cpl_vector * opd_coeff = cpl_vector_new(3);
527 cpl_vector_set (opd_coeff, GRAVI_SC, cpl_matrix_get (res_matrix, 0, 0));
528 cpl_vector_set (opd_coeff, GRAVI_FT, cpl_matrix_get (res_matrix, 1, 0));
529 cpl_vector_set (opd_coeff, 2, rms_fit);
530
531 /* Delete matrix */
532 FREE (cpl_matrix_delete, residual_matrix);
533 FREE (cpl_matrix_delete, res_matrix);
534 FREE (cpl_matrix_delete, model_matrix);
535 FREE (cpl_matrix_delete, rhs_matrix);
536
537 gravi_msg_function_exit(1);
538 return opd_coeff;
539}
540
541/*----------------------------------------------------------------------------*/
563/*----------------------------------------------------------------------------*/
564
565cpl_table * gravi_opds_compute_guess (cpl_table * spectrumsc_table,
566 cpl_table * ft_table,
567 cpl_table * vismet_table,
568 double dit_sc,
569 double lbd_met)
570{
571 gravi_msg_function_start(1);
572 cpl_ensure (spectrumsc_table, CPL_ERROR_NULL_INPUT, NULL);
573 cpl_ensure (ft_table, CPL_ERROR_NULL_INPUT, NULL);
574 cpl_ensure (vismet_table, CPL_ERROR_NULL_INPUT, NULL);
575 cpl_ensure (dit_sc>0, CPL_ERROR_ILLEGAL_INPUT, NULL);
576
577 int ntel = 4;
578
579 cpl_size nrow = cpl_table_get_nrow (spectrumsc_table);
580 cpl_size nrow_met = cpl_table_get_nrow (vismet_table) / ntel;
581 cpl_size nrow_ft = cpl_table_get_nrow (ft_table) / GRAVI_NBASE;
582 int * time_SC = cpl_table_get_data_int (spectrumsc_table, "TIME");
583 CPLCHECK_NUL ("Cannot get data");
584
585 /* Pointer on MET data */
586 int * time_MET = cpl_table_get_data_int (vismet_table, "TIME");
587 double * phase_MET = cpl_table_get_data_double (vismet_table, "PHASE_FC_DRS");
588 CPLCHECK_NUL ("Cannot get MET data");
589
590 /* Pointer on OPD data */
591 int * time_FT = cpl_table_get_data_int (ft_table, "TIME");
592 double * opd_FT = cpl_table_get_data_double (ft_table, "OPD");
593 CPLCHECK_NUL ("Cannot get FT data");
594
595 /* Create table */
596 cpl_table * guess_table = cpl_table_new (nrow * GRAVI_NBASE);
597 gravi_table_new_column (guess_table, "OPD", "m", CPL_TYPE_DOUBLE);
598
599 /* Loop on base */
600 for (int base = 0; base < GRAVI_NBASE; base++) {
601 int tel0 = GRAVI_BASE_TEL[base][0];
602 int tel1 = GRAVI_BASE_TEL[base][1];
603
604 /* Loop on the SC frames */
605 cpl_size row_met = 0, row_ft = 0;
606 for (cpl_size row = 0 ; row < nrow ; row++) {
607
608 /*
609 * Average the MET
610 */
611 int counter_met = 0;
612 double opd_met = 0.0;
613 while (time_MET[row_met*ntel] < (time_SC[row] + dit_sc/2.)) {
614 if ((time_MET[row_met*ntel] > (time_SC[row] - dit_sc/2.)) && (row_met < nrow_met)) {
615 opd_met += phase_MET[row_met*ntel+tel1] - phase_MET[row_met*ntel+tel0];
616 counter_met ++;
617 }
618
619 /* If not enough data to synchronize */
620 if (row_met > nrow_met - 2) {
621 cpl_msg_warning (cpl_func,"Not enough data to synchronise the MET with SC");
622 break;
623 }
624 row_met ++;
625 }
626 /* End sum the MET over the current SC DIT */
627
628 opd_met = opd_met * lbd_met / CPL_MATH_2PI / counter_met; // [m]
629
630 /*
631 * Average the FT
632 */
633 int counter_ft = 0;
634 double opd_ft = 0.0;
635 while (time_FT[row_ft*GRAVI_NBASE+base] < (time_SC[row] + dit_sc/2.)) {
636 if ((time_FT[row_ft*GRAVI_NBASE+base] > (time_SC[row] - dit_sc/2.)) && (row_ft < nrow_ft)) {
637 opd_ft += opd_FT[row_ft*GRAVI_NBASE+base];
638 counter_ft ++;
639 }
640
641 /* If not enough data to synchronize */
642 if (row_ft > nrow_ft - 2) {
643 cpl_msg_warning (cpl_func,"Not enough data to synchronise the FT with SC");
644 break;
645 }
646 row_ft ++;
647 }
648 /* End sum the FT over the current SC DIT */
649
650 opd_ft = opd_ft / counter_ft;
651
652 /* Set the total guess OPD as OPD_FT - OPD_MET */
653 cpl_table_set (guess_table, "OPD", row*GRAVI_NBASE+base, opd_ft - opd_met);
654 CPLCHECK_NUL ("Cannot compute opd guess");
655 } /* End loop on row */
656 } /* End loop on base */
657
658 /*
659 * Remove the mean of each base
660 */
661 cpl_msg_info (cpl_func, "Remove the mean from opdguess table");
662
663 for (int base = 0; base < GRAVI_NBASE; base++) {
664 double mean = gravi_table_get_column_mean (guess_table, "OPD", base, GRAVI_NBASE);
665 gravi_table_add_scalar (guess_table, "OPD", base, GRAVI_NBASE, -1*mean);
666 }
667
668 gravi_msg_function_exit(1);
669 return guess_table;
670}
671
672/*----------------------------------------------------------------------------*/
695/*----------------------------------------------------------------------------*/
696
697cpl_table * gravi_opds_calibration (cpl_table * spectrum_table,
698 cpl_table * detector_table,
699 cpl_table * guess_table)
700{
701 /* Verbose */
702 gravi_msg_function_start(1);
703 cpl_ensure (spectrum_table, CPL_ERROR_NULL_INPUT, NULL);
704 cpl_ensure (detector_table, CPL_ERROR_NULL_INPUT, NULL);
705
706 /* Verbose the phase algorithm selected */
707 if (USE_LINEAR_ELLIPSE) {
708 cpl_msg_info (cpl_func, "Solve linearly the ellipse X,Y");
709 }
710 else {
711 cpl_msg_info (cpl_func, "Solve non-linearly the ellipse X,Y");
712 }
713
714 /* Get the size of the vectors */
715 cpl_size nrow = cpl_table_get_nrow (spectrum_table);
716 cpl_size nwave = gravi_spectrum_get_nwave (spectrum_table);
717 cpl_size npol = gravi_spectrum_get_npol (spectrum_table);
718 CPLCHECK_NUL ("Cannot get size");
719
720 /* Create a fake OI_WAVE table, to average the different
721 * spectral channels [m] */
722
723 cpl_table ** oiwave_tables = cpl_calloc (npol, sizeof (cpl_table *));
724 for (int pol = 0; pol < npol; pol++) {
725 oiwave_tables[pol] = cpl_table_new (nwave);
726 cpl_table_new_column (oiwave_tables[pol], "EFF_WAVE", CPL_TYPE_DOUBLE);
727 for (cpl_size wave = 0; wave < nwave; wave++) {
728 double value = 1.97e-6 + (2.48e-6 - 1.97e-6) / nwave * wave;
729 cpl_table_set (oiwave_tables[pol], "EFF_WAVE", wave, value);
730 }
731 }
732
733 /* Create column in phase table */
734 cpl_table * output_table = cpl_table_new (GRAVI_NBASE * nrow);
735 gravi_table_new_column (output_table,"OPD", "m", CPL_TYPE_DOUBLE);
736 gravi_table_new_column (output_table,"TIME", "us", CPL_TYPE_INT);
737
738 /* Set the time */
739 for (cpl_size row = 0; row < nrow; row++) {
740 double value = cpl_table_get (spectrum_table, "TIME", row, NULL);
741 for (int base = 0; base < GRAVI_NBASE; base++)
742 cpl_table_set (output_table, "TIME", row*GRAVI_NBASE+base, value);
743 }
744
745 /* Loop on base */
746 for (int base = 0; base < GRAVI_NBASE; base ++) {
747
748 /* Check if a guess exists */
749 cpl_vector * opd_guess = NULL;
750 if (guess_table) {
751 opd_guess = cpl_vector_new (nrow);
752 for (cpl_size row = 0; row < nrow; row++) {
753 double value = cpl_table_get (guess_table, "OPD", row*GRAVI_NBASE+base, NULL);
754 cpl_vector_set (opd_guess, row, value);
755 }
756 }
757
758 /* Recover the mean OPD modulation (averaved over pol
759 * and channels) using the ellipses. In [m]. However since oiwave
760 * is not known... this is to a scaling factor correction. */
761 cpl_vector * mean_opd;
762 mean_opd = gravi_ellipse_meanopd_create (spectrum_table, detector_table,
763 oiwave_tables, opd_guess, base);
764 CPLCHECK_NUL ("Cannot compute opd");
765
766 /* Save the mean opd [m] */
767 for (cpl_size row = 0; row < nrow; row++) {
768 double value = cpl_vector_get (mean_opd, row);
769 cpl_table_set (output_table, "OPD", row*GRAVI_NBASE+base, value);
770 }
771 CPLCHECK_NUL ("Cannot set phase");
772
773 FREE (cpl_vector_delete, mean_opd);
774 FREE (cpl_vector_delete, opd_guess);
775
776 } /* End loop on base */
777
778 FREELOOP (cpl_table_delete, oiwave_tables, npol);
779
780 /* Verbose */
781 gravi_msg_function_exit(1);
782 return output_table;
783}
784
785/*----------------------------------------------------------------------------*/
817/*----------------------------------------------------------------------------*/
818
819cpl_error_code gravi_wave_compute_opds (gravi_data * spectrum_data,
820 cpl_table * met_table,
821 const cpl_parameterlist * parlist)
822{
823 gravi_msg_function_start(1);
824 cpl_ensure_code (spectrum_data, CPL_ERROR_NULL_INPUT);
825 cpl_ensure_code (met_table, CPL_ERROR_NULL_INPUT);
826
827 /* Get DIT */
828 cpl_propertylist * spectrum_header = gravi_data_get_header (spectrum_data);
829 double dit_sc = gravi_pfits_get_dit_sc (spectrum_header)*1e6; // [us]
830 CPLCHECK_MSG ("Cannot get DIT");
831
832 /* Get the input */
833 cpl_table * spectrumft_table = gravi_data_get_spectrum_data (spectrum_data, GRAVI_FT);
834 cpl_table * detectorft_table = gravi_data_get_imaging_detector (spectrum_data, GRAVI_FT);
835 cpl_table * spectrumsc_table = gravi_data_get_spectrum_data (spectrum_data, GRAVI_SC);
836 cpl_table * detectorsc_table = gravi_data_get_imaging_detector (spectrum_data, GRAVI_SC);
837 CPLCHECK_MSG ("Cannot get data");
838
839
840 /*
841 * Reduce the raw METROLOGY into VIS_MET
842 */
843 cpl_msg_info (cpl_func, "Compute the phase of MET_FC");
844
845 cpl_table * vismet_table = gravi_metrology_create (met_table, spectrum_header);
846 gravi_metrology_drs (met_table, vismet_table, spectrum_header, parlist);
847 /* also apply the TAC reduction */
848 gravi_metrology_tac(met_table, vismet_table, spectrum_header);
849
850 /* Compute the mean LBD_MET for this file */
851 double lbd_met = gravi_pfits_get_met_wavelength_mean (spectrum_header, met_table);
852
853 /*
854 * Compute the phase of SC and FT. For the SC, we use a guess
855 * of the opd modulation, computed from FT and MET, to unwrap.
856 */
857 cpl_msg_info (cpl_func, "Compute OPD of FT from ellipse");
858
859 cpl_table * ft_table;
860 ft_table = gravi_opds_calibration (spectrumft_table, detectorft_table, NULL);
861 gravi_opds_correct_closures (ft_table, "OPD");
862
863 cpl_msg_info (cpl_func, "Compute OPD of SC from ellipse");
864
865 cpl_table * guess_table;
866 guess_table = gravi_opds_compute_guess (spectrumsc_table, ft_table, vismet_table, dit_sc, lbd_met);
867
868 cpl_table * sc_table;
869 sc_table = gravi_opds_calibration (spectrumsc_table, detectorsc_table, guess_table);
870 gravi_opds_correct_closures (sc_table, "OPD");
871 FREE (cpl_table_delete, guess_table);
872
873 CPLCHECK_MSG ("Cannot calibrate phase");
874
875 /*
876 * Interpolate MET and SC at the sampling frequency of FT
877 * Results are saved in the ft_table table
878 */
879 cpl_msg_info (cpl_func, "Fit MET = a.SC - b.FT + c to get absolute modulation");
880
881 gravi_vis_create_met_ft (ft_table, vismet_table);
882 gravi_vis_create_opdsc_ft (ft_table, sc_table, dit_sc);
883
884 CPLCHECK_MSG ("Cannot resample SC or MET at the FT frequency");
885
886 /*
887 * Compute the scaling coefficients of OPDs by fitting:
888 * OPD_MET_ijt = a.OPD_SC_ijt - b.OPD_FT_ijt + c_ij
889 */
890 cpl_vector * coeff_vector = gravi_opds_fit_opdmet (ft_table, lbd_met);
891
892 CPLCHECK_MSG ("Cannot fit opdmet");
893
894 /* Set the CHI2 of the fit in the QC parameters */
895 cpl_propertylist_update_double (spectrum_header, QC_PHASECHI2,
896 cpl_vector_get (coeff_vector, 2));
897 cpl_propertylist_set_comment (spectrum_header, QC_PHASECHI2,
898 "chi2 of a.SC-b.FT+c=MET");
899
900 /* Add the OPD COEFF in QC parameters */
901 for (int type_data = 0; type_data < 2; type_data++) {
902 double tmp = cpl_vector_get (coeff_vector, type_data);
903 cpl_propertylist_update_float (spectrum_header, OPD_COEFF_SIGN(type_data), tmp);
904 cpl_propertylist_set_comment (spectrum_header, OPD_COEFF_SIGN(type_data), "wavelength correction");
905 }
906
907 /* Set a warning */
908 if (cpl_vector_get (coeff_vector, 2) > 1e-7) {
909 gravi_pfits_add_check (spectrum_header,"Residuals of fit MET=a.SC-b.FT+c are high");
910 }
911
912 if (cpl_vector_get (coeff_vector, 2) > 1.2e-7) {
913 cpl_msg_warning (cpl_func, "*************************************************");
914 cpl_msg_warning (cpl_func, "**** !!! residuals of the fit too high !!! ****");
915 cpl_msg_warning (cpl_func, "**** Residuals are:%7.0f nm ****",cpl_vector_get (coeff_vector, 2)*1e9);
916 cpl_msg_warning (cpl_func, "**** SC and RMN may be desynchronized ****");
917 cpl_msg_warning (cpl_func, "**** (or out of the envelope in LOW) ****");
918 cpl_msg_warning (cpl_func, "*************************************************");
919 }
920
921 CPLCHECK_MSG ("Cannot set QC parameter");
922
923 /*
924 * Correct opd from overall scaling by fitting MET
925 */
926 double coeff_sc = cpl_vector_get (coeff_vector, GRAVI_SC);
927 cpl_table_multiply_scalar (sc_table, "OPD", coeff_sc);
928
929 double coeff_ft = cpl_vector_get (coeff_vector, GRAVI_FT);
930 cpl_table_multiply_scalar (ft_table, "OPD", coeff_ft);
931
932 FREE (cpl_vector_delete, coeff_vector);
933 CPLCHECK_MSG ("Cannot correct OPDs from scaling coefficients");
934
935 /*
936 * Fill the output
937 */
938// if ((det_type == GRAVI_DET_SC || det_type == GRAVI_DET_ALL))
939// {
940 gravi_data_add_table (spectrum_data, NULL, "OPD_SC", sc_table);
941 CPLCHECK_MSG ("Cannot set OPD_SC table");
942// }
943// else
944// FREE (cpl_table_delete, sc_table);
945//
946//
947// if ((det_type == GRAVI_DET_FT || det_type == GRAVI_DET_ALL))
948// {
949 gravi_data_add_table (spectrum_data, NULL, "OPD_FT", ft_table);
950 CPLCHECK_MSG ("Cannot set OPD_FT table");
951// }
952// else
953// FREE (cpl_table_delete, ft_table);
954
955 gravi_data_add_table (spectrum_data, NULL, GRAVI_OI_VIS_MET_EXT, vismet_table);
956 CPLCHECK_MSG ("Cannot set OI_VIS_MET table");
957
958 gravi_msg_function_exit(1);
959 return CPL_ERROR_NONE;
960}
961
962/*----------------------------------------------------------------------------*/
986/*----------------------------------------------------------------------------*/
987
988cpl_table * gravi_wave_fibre (cpl_table * spectrum_table,
989 cpl_table * detector_table,
990 cpl_table * opd_table)
991{
992 gravi_msg_function_start(1);
993 cpl_ensure (spectrum_table, CPL_ERROR_NULL_INPUT, NULL);
994 cpl_ensure (detector_table, CPL_ERROR_NULL_INPUT, NULL);
995 cpl_ensure (opd_table, CPL_ERROR_NULL_INPUT, NULL);
996
997 char name[100];
998
999 /* Create the output table */
1000 cpl_table * wave_fibre = cpl_table_new (1);
1001
1002 /* Get the number of wavelength, region, polarisation... */
1003 cpl_size nwave = cpl_table_get_column_depth (spectrum_table, "DATA1");
1004 cpl_size n_region = cpl_table_get_nrow (detector_table);
1005 cpl_size nrow = cpl_table_get_nrow (spectrum_table);
1006
1007 int npol = (n_region > 24 ? 2 : 1);
1008
1009 /*
1010 * Calibration of each polarization and base
1011 */
1012
1013 for (int pol = 0; pol < npol; pol++) {
1014 for (int base = 0; base < GRAVI_NBASE; base ++) {
1015 cpl_msg_info (cpl_func, "Compute wave fibre for pol %i over %i, base %i over %i",
1016 pol+1, npol, base+1, GRAVI_NBASE);
1017
1018 /* Get the index of the ABCD. */
1019 int iA = gravi_get_region (detector_table, base, 'A', pol);
1020 int iB = gravi_get_region (detector_table, base, 'B', pol);
1021 int iC = gravi_get_region (detector_table, base, 'C', pol);
1022 int iD = gravi_get_region (detector_table, base, 'D', pol);
1023 if (iA<0 || iB<0 || iC<0 || iD<0){
1024 cpl_msg_warning (cpl_func, "Don't found the A, B, C or D !!!");
1025 }
1026
1027 /* Sign of this baseline */
1028 double phi_sign = gravi_region_get_base_sign (detector_table, base);
1029
1030 /* Get the OPD into various flarous */
1031 cpl_matrix * opd_matrix = cpl_matrix_new (1, nrow);
1032 cpl_vector * opd_vector = cpl_vector_new (nrow);
1033 for (cpl_size row = 0; row < nrow; row ++ ) {
1034 double value = cpl_table_get (opd_table, "OPD", row*GRAVI_NBASE+base, NULL);
1035 cpl_matrix_set (opd_matrix, 0, row, value);
1036 cpl_vector_set (opd_vector, row, value);
1037 }
1038 CPLCHECK_NUL ("Cannot extract the OPD");
1039
1040 /* Create array to fill the wavelenghts */
1041 cpl_array * wavelength = cpl_array_new (nwave, CPL_TYPE_DOUBLE);
1042 cpl_array * wavechi2 = cpl_array_new (nwave, CPL_TYPE_DOUBLE);
1043 cpl_array_fill_window (wavelength, 0, nwave, 0.0);
1044 cpl_array_fill_window (wavechi2, 0, nwave, 1e10);
1045
1046 /*
1047 * Loop on spectral channel
1048 */
1049 for (cpl_size wave = 0; wave < nwave; wave++) {
1050
1051 cpl_vector * vector_T = NULL, * vector_X, * vector_Y;
1052
1053 /* Define the vector_X = C - A */
1054 vector_X = gravi_table_get_vector (spectrum_table, wave, GRAVI_DATA[iC]);
1055 vector_T = gravi_table_get_vector (spectrum_table, wave, GRAVI_DATA[iA]);
1056 cpl_vector_subtract (vector_X, vector_T);
1057 FREE (cpl_vector_delete, vector_T);
1058
1059 /* Define the vector_Y = D - B */
1060 vector_Y = gravi_table_get_vector (spectrum_table, wave, GRAVI_DATA[iD]);
1061 vector_T = gravi_table_get_vector (spectrum_table, wave, GRAVI_DATA[iB]);
1062 cpl_vector_subtract (vector_Y, vector_T);
1063 FREE (cpl_vector_delete, vector_T);
1064
1065 CPLCHECK_NUL ("Cannot compute vector_X or Y");
1066
1067 /* Compute envelope from OPD for this channel */
1068 cpl_vector * envelope_vector = gravi_compute_envelope (opd_vector, wave, nwave);
1069
1070 /* Compute the phase from the ellipse */
1071 cpl_vector * phase;
1072 phase = gravi_ellipse_phase_create (vector_X, vector_Y,
1073 envelope_vector);
1074 FREE (cpl_vector_delete, vector_X);
1075 FREE (cpl_vector_delete, vector_Y);
1076 FREE (cpl_vector_delete, envelope_vector);
1077
1078 /* If the computation of the ellipse fails, we continue with next wave */
1079 if (phase == NULL) {
1080 cpl_msg_warning (cpl_func, "Cannot compute wave for %lld and base %d", wave, base);
1081 continue;
1082 }
1083
1084 /* Multiply by the sign */
1085 cpl_vector_multiply_scalar (phase, phi_sign);
1086
1087 /* Unwrap phase from the OPD */
1088 double lbd_channel = 1.95e-6 + (2.46e-6 - 1.95e-6) / nwave * wave;
1089 gravi_vector_unwrap_with_guess (phase, opd_vector, CPL_MATH_2PI / lbd_channel);
1090
1091 /* Fit the slope of the phase versus OPD gives the
1092 * wavelength of the spectral element */
1093 const cpl_size mindeg = 0, maxdeg = 1;
1094 cpl_polynomial * fit_slope = cpl_polynomial_new (1);
1095 cpl_polynomial_fit (fit_slope, opd_matrix, NULL, phase, NULL, CPL_FALSE, &mindeg, &maxdeg);
1096
1097 /* Compute residuals */
1098 cpl_vector * residuals = cpl_vector_new (nwave);
1099 double rechisq;
1100 cpl_vector_fill_polynomial_fit_residual (residuals, phase, NULL, fit_slope, opd_matrix, &rechisq);
1101
1102 if(PLOT_WAVE_PHASE_VS_OPD && (wave == WAVE_TO_PLOT))
1103 {
1104 char gnuplot_str[200];
1105 sprintf (gnuplot_str, "set title 'Wavelength (base %d)'; set xlabel 'Phase [rad]'; set ylabel 'OPD (m)';", base);
1106 cpl_plot_vector (gnuplot_str, NULL, NULL, opd_vector);
1107 sprintf (gnuplot_str, "set title 'Wavelength residuals (base %d)'; set xlabel 'Phase [rad]'; set ylabel 'OPD (m)';", base);
1108 cpl_plot_vector (gnuplot_str, NULL, NULL, residuals);
1109 CPLCHECK_NUL ("Cannot plot OPD versus phase");
1110 }
1111
1112 CPLCHECK_NUL ("Cannot fit the OPD and phase");
1113
1114 /* Get the slope */
1115 const cpl_size pow_slope = 1;
1116 double slope = cpl_polynomial_get_coeff (fit_slope, &pow_slope);
1117
1118 /* Check slope sign. Should be positive */
1119 if (slope < 0.0 && wave == 0) {
1120 cpl_msg_warning (cpl_func, "Negative wavelength!! "
1121 "Report to DRS team");
1122 }
1123
1124 /* Set the wavelength */
1125 cpl_array_set (wavechi2, wave, sqrt(rechisq));
1126 cpl_array_set (wavelength, wave, CPL_MATH_2PI / fabs(slope));
1127
1128 /* Free memory */
1129 cpl_vector_delete (phase);
1130 cpl_vector_delete (residuals);
1131 cpl_polynomial_delete (fit_slope);
1132
1133 } /* End loop on wave */
1134
1135 cpl_matrix_delete (opd_matrix);
1136 cpl_vector_delete (opd_vector);
1137
1138 /* Set the wavelength array in the output table */
1139 sprintf (name, "BASE_%s_%s", GRAVI_BASE_NAME[base], GRAVI_POLAR(pol,npol));
1140 cpl_table_new_column_array (wave_fibre, name, CPL_TYPE_DOUBLE, nwave);
1141 cpl_table_set_column_unit (wave_fibre, name, "m");
1142 cpl_table_set_array (wave_fibre, name, 0, wavelength);
1143 cpl_array_delete (wavelength);
1144
1145 /* Set the chi2 array in the output table */
1146 sprintf (name, "BASE_%s_%s_CHI2", 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_array (wave_fibre, name, 0, wavechi2);
1149 cpl_array_delete (wavechi2);
1150
1151 } /* End loop on base*/
1152 } /* End loop on polar */
1153
1154 gravi_msg_function_exit(1);
1155 return wave_fibre;
1156}
1157
1158/*----------------------------------------------------------------------------*/
1176/*----------------------------------------------------------------------------*/
1177
1178cpl_table * gravi_wave_fit_2d (cpl_table * wavefibre_table,
1179 cpl_table * detector_table,
1180 gravi_data * wave_param,
1181 cpl_size fullstartx,
1182 int spatial_order,
1183 int spectral_order,
1184 double * rms_residuals)
1185{
1186 gravi_msg_function_start(1);
1187 cpl_ensure (wavefibre_table, CPL_ERROR_NULL_INPUT, NULL);
1188 cpl_ensure (detector_table, CPL_ERROR_NULL_INPUT, NULL);
1189
1190 char name[100];
1191 *rms_residuals = 0;
1192
1193 /* Get numbers */
1194 cpl_size n_region = cpl_table_get_nrow (detector_table);
1195 int npol = n_region > 24 ? 2 : 1;
1196 cpl_size nwave = cpl_table_get_column_depth (wavefibre_table, npol > 1 ? "BASE_12_S" : "BASE_12_C");
1197 CPLCHECK_NUL ("Cannot get data");
1198
1199 /* get the wave param */
1200 cpl_propertylist * wave_param_plist = gravi_data_get_plist(wave_param,GRAVI_PRIMARY_HDR_EXT);
1201
1202 /* Odd index, for SC only */
1203 cpl_vector * odd_index = cpl_vector_new (nwave);
1204 for (int i = fullstartx; i < fullstartx + nwave; i++) {
1205 if (nwave > GRAVI_LBD_FTSC) cpl_vector_set (odd_index, i - fullstartx, ((i/64)%2 == 0) ? 0 : 1);
1206 else cpl_vector_set (odd_index, i - fullstartx, 0);
1207 }
1208
1209 CPLCHECK_NUL ("Cannot buid odd_index");
1210
1211 /* Save the 2D coeficient of each polar */
1212 cpl_polynomial ** coef_poly = cpl_calloc (npol, sizeof (cpl_polynomial*));
1213
1214 /* Loop on polarisation */
1215 for (int pol = 0; pol < npol; pol++) {
1216
1217 /* Prepare the 2D coordinates
1218 * and values to fit */
1219 cpl_vector * coord_X = cpl_vector_new (GRAVI_NBASE * nwave);
1220 cpl_vector * coord_Y = cpl_vector_new (GRAVI_NBASE * nwave);
1221
1222 cpl_vector * all_wavelength = cpl_vector_new (GRAVI_NBASE * nwave);
1223 cpl_vector * all_wavechi2 = cpl_vector_new (GRAVI_NBASE * nwave);
1224 cpl_vector * all_valid = cpl_vector_new (GRAVI_NBASE * nwave);
1225
1226 /*
1227 * Loop on base and wave to get all
1228 * wavelenght and coordinates
1229 */
1230 for (int base = 0; base < GRAVI_NBASE; base ++) {
1231
1232 /* Mean position of this baseline */
1233 int iA = gravi_get_region (detector_table, base, 'A', pol);
1234 int iB = gravi_get_region (detector_table, base, 'B', pol);
1235 int iC = gravi_get_region (detector_table, base, 'C', pol);
1236 int iD = gravi_get_region (detector_table, base, 'D', pol);
1237
1238 /* WAVE_FIBRE data */
1239 sprintf (name, "BASE_%s_%s", GRAVI_BASE_NAME[base], GRAVI_POLAR(pol,npol));
1240 cpl_array * wavelength = cpl_table_get_data_array (wavefibre_table, name)[0];
1241
1242 sprintf (name, "BASE_%s_%s_CHI2", GRAVI_BASE_NAME[base], GRAVI_POLAR(pol,npol));
1243 cpl_array * wavechi2 = cpl_table_get_data_array (wavefibre_table, name)[0];
1244
1245 /* Loop on wave */
1246 for (cpl_size wave = 0; wave < nwave; wave++) {
1247 cpl_size pos = base * nwave + wave;
1248
1249 /* Get the values */
1250 int nv = 0;
1251 double wave_value = cpl_array_get (wavelength, wave, &nv);
1252 double chi2_value = cpl_array_get (wavechi2, wave, &nv);
1253
1254 /* The FT accept all channel for the 2D fit */
1255 cpl_vector_set (all_valid, pos, 1);
1256
1257 if (nwave > 1000) {
1258 /* SC HIGH */
1259 if ((chi2_value > M_PI_4 ||
1260 wave_value < gravi_pfits_get_double_default(wave_param_plist, "ESO OIWAVE HIGH LBD MIN", 2.01e-6) ||
1261 wave_value > gravi_pfits_get_double_default(wave_param_plist, "ESO OIWAVE HIGH LBD MAX", 2.43e-6)))
1262 cpl_vector_set (all_valid, pos, 0);
1263 }
1264 else if (nwave > 100) {
1265 /* SC MEDIUM */
1266 if ((chi2_value > M_PI_4 ||
1267 wave_value < gravi_pfits_get_double_default(wave_param_plist, "ESO OIWAVE MED LBD MIN", 2.01e-6) ||
1268 wave_value > gravi_pfits_get_double_default(wave_param_plist, "ESO OIWAVE MED LBD MAX", 2.43e-6)))
1269 cpl_vector_set (all_valid, pos, 0);
1270 }
1271 else if (nwave > GRAVI_LBD_FTSC) {
1272 /* SC LOW */
1273 if ((chi2_value > M_PI_4 ||
1274 //wave_value < 2.01e-6 ||
1275 wave_value < gravi_pfits_get_double_default(wave_param_plist, "ESO OIWAVE LOW LBD MIN", 1.99e-6) ||
1276 wave_value > gravi_pfits_get_double_default(wave_param_plist, "ESO OIWAVE LOW LBD MAX", 2.5e-6) ||
1277 wave == 0 || wave == nwave-1))
1278 cpl_vector_set (all_valid, pos, 0);
1279 }
1280 else if (nwave == GRAVI_LBD_FTSC) {
1281 /* FT */
1282 if ((chi2_value > M_PI_4 ||
1283 //wave_value < 2.01e-6 ||
1284 wave_value < 1.99e-6 ||
1285 wave_value > 2.5e-6))
1286 cpl_vector_set (all_valid, pos, 0);
1287 }
1288
1289 /* Set the wavelength */
1290 cpl_vector_set (all_wavelength, pos, wave_value);
1291 cpl_vector_set (all_wavechi2, pos, chi2_value);
1292
1293 /* Set the X position as the mean of the 4 regions */
1294 cpl_vector_set (coord_X, pos, (double)(iA + iB + iC + iD) / 4.);
1295
1296 /* Set the Y position. Add a shift of 0.15 pixels
1297 * on odd output for SC detector */
1298 cpl_vector_set (coord_Y, pos, wave + cpl_vector_get (odd_index, wave)*0.15);
1299
1300 } /* End loop on wave */
1301 } /* End loop on base */
1302
1303 CPLCHECK_NUL ("Error get all wavelength");
1304
1305 /*
1306 * Reformat the valid point in vector and matrix
1307 */
1308 cpl_size nvalid = cpl_vector_get_sum (all_valid);
1309
1310 cpl_msg_info (cpl_func, "Remove %lld/%lld bad wave",
1311 GRAVI_NBASE*nwave - nvalid, GRAVI_NBASE * nwave);
1312
1313 cpl_vector * vector = cpl_vector_new (nvalid);
1314 cpl_matrix * matrix = cpl_matrix_new (2, nvalid);
1315
1316 for (cpl_size c = 0, i = 0 ; i < nwave * GRAVI_NBASE; i ++) {
1317 if (!cpl_vector_get (all_valid, i)) continue;
1318 cpl_vector_set (vector, c, cpl_vector_get (all_wavelength, i));
1319 cpl_matrix_set (matrix, 0, c, cpl_vector_get (coord_X, i));
1320 cpl_matrix_set (matrix, 1, c, cpl_vector_get (coord_Y, i));
1321 c++;
1322 }
1323
1324 CPLCHECK_NUL ("Error cropping all wavelength");
1325
1326 /*
1327 * Perform a 2D fit with a polynomial model
1328 * between the position and the wavelength = F(y, x)
1329 */
1330 cpl_size deg2d[2] = {2, 3};
1331 if ( (nwave < 20) && (nwave > 8) ) {deg2d[0] = 2; deg2d[1] = 7;} // FIXME Useless line ?
1332 deg2d[0] = spatial_order;
1333 deg2d[1] = spectral_order;
1334
1335 cpl_msg_info (cpl_func, "Fit a 2d polynomial {%lli..%lli} to the wavelengths map", deg2d[0], deg2d[1]);
1336
1337 cpl_polynomial * fit2d = cpl_polynomial_new (2);
1338 cpl_polynomial_fit (fit2d, matrix, NULL, vector, NULL, CPL_TRUE, NULL, deg2d);
1339 coef_poly[pol] = fit2d;
1340
1341 CPLCHECK_NUL ("Cannot fit 2D");
1342
1343 /*
1344 * Compute residuals
1345 */
1346 double rechisq = 0.0;
1347 cpl_vector * residuals = cpl_vector_new (nvalid);
1348 cpl_vector_fill_polynomial_fit_residual (residuals, vector, NULL, fit2d, matrix, &rechisq);
1349 *rms_residuals += cpl_vector_get_stdev(residuals)/npol;
1350 FREE (cpl_vector_delete, residuals);
1351 CPLCHECK_NUL ("Cannot compute residuals");
1352
1353
1354 FREE (cpl_matrix_delete, matrix);
1355 FREE (cpl_vector_delete, vector);
1356 FREE (cpl_vector_delete, all_wavelength);
1357 FREE (cpl_vector_delete, all_wavechi2);
1358 FREE (cpl_vector_delete, all_valid);
1359 FREE (cpl_vector_delete, coord_X);
1360 FREE (cpl_vector_delete, coord_Y);
1361
1362 } /* End loop on polarisation */
1363
1364
1365 /*
1366 * Create and fill the interpolated WAVE_DATA table
1367 */
1368
1369 cpl_table * wavedata_table = cpl_table_new (1);
1370 cpl_vector * pos = cpl_vector_new (2);
1371 cpl_array * value = cpl_array_new (nwave, CPL_TYPE_DOUBLE);
1372
1373 for (cpl_size region = 0 ; region < n_region; region ++) {
1374
1375 int pol = gravi_region_get_pol (detector_table, region);
1376 /* Loop on wave to evaluate the 2D polynome */
1377 for (cpl_size wave = 0; wave < nwave; wave ++) {
1378
1379 /* Evaluate */
1380 cpl_vector_set (pos, 0, region);
1381 cpl_vector_set (pos, 1, wave + cpl_vector_get (odd_index, wave)*0.15);
1382
1383 double result = cpl_polynomial_eval (coef_poly[pol], pos);
1384 cpl_array_set (value, wave, result);
1385
1386 }
1387 /* ensure cresent wavelength */
1388 double previous_wave = cpl_array_get(value, nwave/2, NULL);
1389 for (cpl_size wave_loop = nwave/2 ; wave_loop >= 0 ; wave_loop --){
1390 if (previous_wave < cpl_array_get(value, wave_loop, NULL))
1391 cpl_array_set(value, wave_loop, previous_wave);
1392 else previous_wave = cpl_array_get(value, wave_loop, NULL);
1393 }
1394
1395 previous_wave = cpl_array_get(value, nwave/2, NULL);
1396 for (cpl_size wave_loop = nwave/2 ; wave_loop < nwave ; 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 /* Add column */
1403 char * data_x = GRAVI_DATA[region];
1404 cpl_table_new_column_array (wavedata_table, data_x, CPL_TYPE_DOUBLE, nwave);
1405 cpl_table_set_array (wavedata_table, data_x, 0, value);
1406
1407 } /* End loop on regions */
1408
1409 /* Delete allocations */
1410 FREE (cpl_vector_delete, pos);
1411 FREE (cpl_array_delete, value);
1412 FREELOOP (cpl_polynomial_delete, coef_poly, npol);
1413 FREE (cpl_vector_delete, odd_index);
1414
1415 gravi_msg_function_exit(1);
1416 return wavedata_table;
1417}
1418
1419
1420
1421/*----------------------------------------------------------------------------*/
1422/*
1423 * @brief TBD
1424 */
1425/*----------------------------------------------------------------------------*/
1426
1427cpl_table * gravi_wave_fit_individual (cpl_table * wave_individual_table,
1428 cpl_table * weight_individual_table,
1429 cpl_table * wave_fitted_table,
1430 cpl_table * opd_table,
1431 cpl_table * spectrum_table,
1432 cpl_table * detector_table,
1433 double n0, double n1, double n2)
1434{
1435
1436 gravi_msg_function_start(1);
1437
1438 cpl_ensure (wave_individual_table, CPL_ERROR_NULL_INPUT, NULL);
1439 cpl_ensure (weight_individual_table, CPL_ERROR_NULL_INPUT, NULL);
1440 cpl_ensure (wave_fitted_table, CPL_ERROR_NULL_INPUT, NULL);
1441 cpl_ensure (opd_table, CPL_ERROR_NULL_INPUT, NULL);
1442 cpl_ensure (spectrum_table, CPL_ERROR_NULL_INPUT, NULL);
1443 cpl_ensure (detector_table, CPL_ERROR_NULL_INPUT, NULL);
1444
1445 /* Get the number of wavelength, region, polarisation... */
1446 cpl_size nwave = cpl_table_get_column_depth (spectrum_table, "DATA1");
1447 cpl_size n_region = cpl_table_get_nrow (detector_table);
1448 cpl_size nrow = cpl_table_get_nrow (spectrum_table);
1449 int npol = (n_region > 24 ? 2 : 1);
1450 cpl_size nwave_ref=3000;
1451 if (nwave<10) nwave_ref=600;
1452
1453 CPLCHECK_NUL ("Cannot buid odd_index");
1454
1455 /* Get OPD Table */
1456
1457
1458
1459 /* create temporary variables */
1460
1461 cpl_array * wave_individual_array = cpl_array_new (nwave, CPL_TYPE_DOUBLE);
1462 cpl_array * weight_individual_array = cpl_array_new (nwave, CPL_TYPE_DOUBLE);
1463 cpl_matrix * data_flux_matrix = cpl_matrix_new (nrow, nwave);
1464 cpl_matrix * vis_to_flux_matrix = cpl_matrix_new (nrow, 3);
1465 cpl_matrix * signal_matrix = cpl_matrix_new (nwave, nwave_ref);
1466 cpl_matrix * residual_matrix = cpl_matrix_new (nwave, nwave_ref);
1467 cpl_array * wave_reference_array = cpl_array_new (nwave_ref, CPL_TYPE_DOUBLE);
1468
1469 /*initialise arrays */
1470
1471 cpl_matrix_fill_column(vis_to_flux_matrix,1,2);
1472 for (cpl_size wave_ref = 0; wave_ref < nwave_ref; wave_ref+=1)
1473 {
1474 // make matrix todo
1475 double wave_value=1.95e-6+wave_ref*0.6e-6/((double) nwave_ref);
1476 cpl_array_set_double(wave_reference_array,wave_ref,wave_value);
1477 }
1478
1479 CPLCHECK_NUL ("Cannot initialize arrays for wavelength fit");
1480
1481
1482 for (cpl_size region = 0 ; region < n_region; region ++) {
1483
1484 int base=gravi_region_get_base (detector_table, region);
1485 char * data_x = GRAVI_DATA[region];
1486
1487 cpl_msg_info_overwritable (cpl_func, "Least square fitting of wavelength for region %s", data_x);
1488 // get data_flux_matrix
1489
1490 for (cpl_size row = 0; row < nrow; row ++ ) {
1491 cpl_array * flux_array= cpl_table_get_data_array(spectrum_table,data_x)[row];
1492 for (cpl_size wave = 0; wave < nwave; wave ++) {
1493 cpl_matrix_set (data_flux_matrix, row, wave, cpl_array_get(flux_array,wave, NULL));
1494 }
1495 }
1496
1497
1498 for (cpl_size wave_ref = 0; wave_ref < nwave_ref; wave_ref+=1)
1499 {
1500 // make matrix
1501 double wave_value=cpl_array_get(wave_reference_array,wave_ref,NULL);
1502
1503 for (cpl_size row = 0; row < nrow; row ++ ) {
1504 double opd = cpl_table_get (opd_table, "OPD", row*GRAVI_NBASE+base, NULL);
1505 double coherence_loss=1;
1506 if (fabs(opd) > 1e-9)
1507 {
1508 if (nwave <30)
1509 {
1510 coherence_loss=sin(opd*19500)/(opd*19500);
1511 }
1512 }
1513 cpl_matrix_set(vis_to_flux_matrix,row,0,cos(opd*6.28318530718/wave_value)*coherence_loss);
1514 cpl_matrix_set(vis_to_flux_matrix,row,1,sin(opd*6.28318530718/wave_value)*coherence_loss);
1515 }
1516
1517 cpl_matrix * coef_vis = cpl_matrix_solve_normal(vis_to_flux_matrix,data_flux_matrix); // coef_vis is 3xnwave
1518 cpl_matrix * data_flux_fit = cpl_matrix_product_create(vis_to_flux_matrix,coef_vis);
1519 cpl_matrix * residuals_fit = cpl_matrix_duplicate(data_flux_fit);
1520 cpl_matrix_subtract (residuals_fit,data_flux_matrix); //residuals_fit is nrowxnwave
1521
1522 for (cpl_size wave = 0; wave < nwave; wave ++ ) {
1523 cpl_matrix * temp_matrix = cpl_matrix_extract_column (data_flux_fit, wave);
1524 cpl_matrix_set(signal_matrix,wave,wave_ref,cpl_matrix_get_stdev(temp_matrix));
1525 FREE (cpl_matrix_delete, temp_matrix);
1526 cpl_matrix * temp_matrix2 = cpl_matrix_extract_column (residuals_fit, wave);
1527 cpl_matrix_set(residual_matrix,wave,wave_ref,cpl_matrix_get_stdev(temp_matrix2));
1528 FREE (cpl_matrix_delete, temp_matrix2);
1529 }
1530
1531 FREE (cpl_matrix_delete, coef_vis);
1532 FREE (cpl_matrix_delete, data_flux_fit);
1533 FREE (cpl_matrix_delete, residuals_fit);
1534
1535 CPLCHECK_NUL ("Cannot do Matrix inversion to calculate optimum wavelength");
1536
1537 }
1538
1539 // get minimum chi2 and amplitude signal for each wave
1540 cpl_size wave_ref=1;
1541 cpl_size discarded=1;
1542 for (cpl_size wave = 0; wave < nwave; wave ++ ) {
1543
1544 cpl_matrix * chi2_extract=cpl_matrix_extract_row(residual_matrix,wave);
1545
1546 cpl_matrix_get_minpos(chi2_extract,&discarded, &wave_ref );
1547
1548 double wave_value = cpl_array_get(wave_reference_array, wave_ref, NULL );
1549 double weight_value = cpl_matrix_get(signal_matrix, wave , wave_ref )/(0.1+cpl_matrix_get(residual_matrix, wave, wave_ref ));
1550
1551 cpl_array_set_double (wave_individual_array, wave, wave_value);
1552 cpl_array_set_double (weight_individual_array, wave, weight_value);
1553
1554 FREE (cpl_matrix_delete, chi2_extract);
1555
1556 }
1557
1558 /* Add column */
1559 cpl_table_new_column_array (wave_individual_table, data_x, CPL_TYPE_DOUBLE, nwave);
1560 cpl_table_set_array (wave_individual_table, data_x, 0, wave_individual_array);
1561 cpl_table_new_column_array (weight_individual_table, data_x, CPL_TYPE_DOUBLE, nwave);
1562 cpl_table_set_array (weight_individual_table, data_x, 0, weight_individual_array);
1563 cpl_table_new_column_array (wave_fitted_table, data_x, CPL_TYPE_DOUBLE, nwave);
1564 cpl_table_set_array (wave_fitted_table, data_x, 0, wave_individual_array);
1565 }
1566
1567 CPLCHECK_NUL ("Cannot get individual wavelength for each pixel");
1568
1569 FREE (cpl_array_delete ,wave_individual_array);
1570 FREE (cpl_array_delete ,weight_individual_array);
1571 FREE (cpl_array_delete ,wave_reference_array);
1572 FREE (cpl_matrix_delete ,data_flux_matrix);
1573 FREE (cpl_matrix_delete ,vis_to_flux_matrix);
1574 FREE (cpl_matrix_delete ,signal_matrix);
1575 FREE (cpl_matrix_delete ,residual_matrix);
1576
1577 cpl_msg_info (cpl_func, "Now fitting polynomials on wavelength channels");
1578
1579 cpl_matrix * coef_to_wave = cpl_matrix_new (n_region / npol,5);
1580 cpl_matrix * coef_to_wave_weight = cpl_matrix_new (n_region / npol,n_region / npol);
1581 cpl_matrix * wavelength = cpl_matrix_new(n_region / npol,1);
1582
1583 // set coordinates
1584 for (cpl_size region = 0 ; region < n_region/ npol; region ++)
1585 {
1586 double mean_region = region - (n_region/npol-1)*0.5;
1587 cpl_matrix_set (coef_to_wave, region, 0, 1);
1588 cpl_matrix_set (coef_to_wave, region, 1, mean_region);
1589 cpl_matrix_set (coef_to_wave, region, 2, mean_region*mean_region);
1590 cpl_matrix_set (coef_to_wave, region, 3, mean_region*mean_region*mean_region);
1591 cpl_matrix_set (coef_to_wave, region, 4, mean_region*mean_region*mean_region*mean_region);
1592 }
1593
1594 for (int pol = 0; pol < npol; pol++) {
1595
1596 cpl_msg_info (cpl_func, "Looping for polyfit now, with pol: %i",(int) pol);
1597
1598 for (cpl_size wave = 0; wave < nwave; wave ++) {
1599
1600 // get the data for a common wave row
1601 for (cpl_size region = 0 ; region < n_region/ npol; region ++) {
1602
1603 // Get the pointers to the table arrays
1604 char * data_x = GRAVI_DATA[region*npol+pol];
1605
1606
1607 const cpl_array * wave_array = cpl_table_get_array (wave_individual_table, data_x, 0);
1608 const cpl_array * weight_array = cpl_table_get_array (weight_individual_table, data_x, 0);
1609
1610
1611 cpl_matrix_set (wavelength, region, 0, cpl_array_get(wave_array,wave,NULL));
1612 double weight_value=cpl_array_get(weight_array,wave,NULL);
1613 cpl_matrix_set (coef_to_wave_weight, region, region, weight_value*weight_value);
1614
1615 }
1616
1617
1618 cpl_matrix * coef_to_wave2 = cpl_matrix_product_create(coef_to_wave_weight,coef_to_wave);
1619 cpl_matrix * wavelength2 = cpl_matrix_product_create(coef_to_wave_weight,wavelength);
1620
1621 // Fit a second order polynomial
1622 cpl_matrix * coeff = cpl_matrix_solve_normal(coef_to_wave2, wavelength2); // 5 x 1
1623 cpl_matrix * wavelength_fitted = cpl_matrix_product_create(coef_to_wave, coeff); // n_region / npol x 1
1624
1625 CPLCHECK_NUL ("Cannot do Matrix inversion to calculate optimum polynomial for wavelength");
1626
1627 // Write result to new wave table
1628 for (cpl_size region = 0 ; region < n_region/ npol; region ++) {
1629
1630 char * data_x = GRAVI_DATA[region*npol+pol];
1631 //cpl_msg_info (cpl_func, "Writing: region %i",region);
1632 cpl_array * wave_array = cpl_table_get_data_array (wave_fitted_table, data_x)[0];
1633
1634 cpl_array_set_double(wave_array,wave,cpl_matrix_get(wavelength_fitted,region,0));
1635
1636 }
1637
1638 FREE (cpl_matrix_delete ,coef_to_wave2);
1639 FREE (cpl_matrix_delete ,wavelength2);
1640 FREE (cpl_matrix_delete ,coeff);
1641 FREE (cpl_matrix_delete ,wavelength_fitted);
1642 }
1643
1644 }
1645 FREE (cpl_matrix_delete ,coef_to_wave);
1646 FREE (cpl_matrix_delete ,coef_to_wave_weight);
1647 FREE (cpl_matrix_delete ,wavelength);
1648
1649 CPLCHECK_NUL ("Cannot fit individual wavelength with 3rd order polynomial");
1650
1651 cpl_msg_info (cpl_func, "Correcting for wavelength error");
1652
1653
1654 /* Correct the computed wavelength from the dispersion */
1655 for (cpl_size region = 0 ; region < n_region; region ++)
1656 {
1657 char * data_x = GRAVI_DATA[region];
1658 cpl_array * wavelength_fitted = cpl_table_get_data_array (wave_fitted_table, data_x)[0];
1659 cpl_size nwave_fitted = cpl_array_get_size (wavelength_fitted);
1660 for (cpl_size wave = 0 ; wave < nwave_fitted ; wave ++ ) {
1661
1662 double result = cpl_array_get (wavelength_fitted, wave, NULL);
1663 double d_met = (result - LAMBDA_MET) / LAMBDA_MET;
1664 cpl_array_set (wavelength_fitted, wave, result * (n0 + n1*d_met + n2*d_met*d_met));
1665
1666 }
1667 for (cpl_size wave = nwave_fitted/2 ; wave < nwave_fitted-1 ; wave ++ ) {
1668 double result = cpl_array_get (wavelength_fitted, wave, NULL);
1669 double result2 = cpl_array_get (wavelength_fitted, wave+1, NULL);
1670 if (result2<result+2e-10) {
1671 result2=result+2e-10;
1672 cpl_array_set (wavelength_fitted, wave+1, result2);
1673 }
1674 }
1675 for (cpl_size wave = nwave_fitted/2 ; wave > 0 ; wave -- ) {
1676 double result = cpl_array_get (wavelength_fitted, wave, NULL);
1677 double result2 = cpl_array_get (wavelength_fitted, wave-1, NULL);
1678 if (result2>result-2e-10) {
1679 result2=result-2e-10;
1680 cpl_array_set (wavelength_fitted, wave-1, result2);
1681 }
1682 }
1683 }
1684
1685 CPLCHECK_NUL ("Error in correcting for dispersion");
1686
1687 gravi_msg_function_exit(1);
1688 return wave_fitted_table;
1689}
1690
1691
1692
1693
1694/*----------------------------------------------------------------------------*/
1701/*----------------------------------------------------------------------------*/
1702
1703cpl_error_code gravi_wave_correct_dispersion (cpl_table * wave_fibre,
1704 double n0, double n1, double n2)
1705{
1706 gravi_msg_function_start(1);
1707 cpl_ensure_code (wave_fibre, CPL_ERROR_NULL_INPUT);
1708
1709 char name[100];
1710
1711 /* Get the number of polarisation */
1712 int npol = cpl_table_has_column (wave_fibre, "BASE_12_P") ? 2 : 1;
1713
1714 /* Loop on columns in the table */
1715 for (int pol = 0; pol < npol; pol ++) {
1716 for (int base = 0; base < GRAVI_NBASE; base ++) {
1717
1718 /* Get data of this region */
1719 sprintf (name, "BASE_%s_%s", GRAVI_BASE_NAME[base], GRAVI_POLAR(pol, npol));
1720 cpl_array * wavelength = cpl_table_get_data_array (wave_fibre, name)[0];
1721 CPLCHECK_MSG ("Cannot get data");
1722
1723 /* Loop on wave */
1724 cpl_size nwave = cpl_array_get_size (wavelength);
1725 for (cpl_size wave = 0 ; wave < nwave ; wave ++ ) {
1726
1727 /* Correct the computed wavelength from the dispersion */
1728 double result = cpl_array_get (wavelength, wave, NULL);
1729 double d_met = (result - LAMBDA_MET) / LAMBDA_MET;
1730 cpl_array_set (wavelength, wave, result * (n0 + n1*d_met + n2*d_met*d_met));
1731
1732 }
1733 } /* End loop on base */
1734 } /* End loop on pol */
1735
1736 gravi_msg_function_exit(1);
1737 return CPL_ERROR_NONE;
1738}
1739
1740/*----------------------------------------------------------------------------*/
1746/*----------------------------------------------------------------------------*/
1747
1748cpl_error_code gravi_wave_correct_color (gravi_data * vis_data)
1749{
1750 gravi_msg_function_start(1);
1751 cpl_ensure_code (vis_data, CPL_ERROR_NULL_INPUT);
1752
1753 cpl_propertylist * primary_header = gravi_data_get_header (vis_data);
1754 cpl_propertylist * oiwave_header = NULL;
1755 cpl_table * oiwave_table = NULL;
1756
1757 /* For each type of data SC / FT */
1758 int ntype_data = 2;
1759 for (int type_data = 0; type_data < ntype_data ; type_data ++) {
1760
1761 /* Loop on polarisation */
1762 int npol = gravi_pfits_get_pola_num (primary_header, type_data);
1763 for (int pol = 0 ; pol<npol ; pol++) {
1764 oiwave_table = cpl_table_duplicate ( gravi_data_get_oi_wave ( vis_data, type_data, pol, npol ) );
1765 oiwave_header = cpl_propertylist_duplicate ( gravi_data_get_oi_wave_plist ( vis_data, type_data, pol, npol ) );
1766
1767 /* here you can do what you want on this duplicated oi_wave_table
1768 * to get OI_FLUX table :
1769 * cpl_table * oiflux_table = gravi_data_get_oi_flux(vis_data, type_data, pol, npol)
1770 * */
1771
1772
1773 /* save the new extension with name OI_WAVELENGTH_CORR */
1774 gravi_data_add_table(vis_data, oiwave_header, "OI_WAVELENGTH_CORR", oiwave_table);
1775
1776 CPLCHECK_MSG("Cannot apply color wave correction");
1777 }
1778 /* End loop on polarisation */
1779 }
1780 /* End loop on data_type */
1781
1782 gravi_msg_function_exit(1);
1783 return CPL_ERROR_NONE;
1784}
1785
1786
1787/*----------------------------------------------------------------------------*/
1799/*----------------------------------------------------------------------------*/
1800
1801cpl_imagelist * gravi_wave_test_image (cpl_table * wavedata_table,
1802 cpl_table * wavefibre_table,
1803 cpl_table * profile_table,
1804 cpl_table * detector_table)
1805{
1806 gravi_msg_function_start(1);
1807 cpl_ensure (wavedata_table, CPL_ERROR_NULL_INPUT, NULL);
1808 cpl_ensure (wavefibre_table, CPL_ERROR_NULL_INPUT, NULL);
1809 cpl_ensure (detector_table, CPL_ERROR_NULL_INPUT, NULL);
1810 cpl_ensure (profile_table, CPL_ERROR_NULL_INPUT, NULL);
1811
1812 int nv = 0;
1813 char name[100];
1814
1815 cpl_size n_region = cpl_table_get_nrow (detector_table);
1816 int npol = (n_region > 24 ? 2 : 1);
1817
1818 CPLCHECK_NUL ("Cannot get data");
1819
1820 /* Init the output images */
1821 cpl_size sizex = cpl_table_get_column_dimension (profile_table, "DATA1", 0);
1822 cpl_size sizey = cpl_table_get_column_dimension (profile_table, "DATA1", 1);
1823
1824 cpl_image * profilesum_image = cpl_image_new (sizex, sizey, CPL_TYPE_DOUBLE);
1825 cpl_image_fill_window (profilesum_image, 1, 1, sizex, sizey, 0.0);
1826
1827 cpl_image * wave_image = cpl_image_new (sizex, sizey, CPL_TYPE_DOUBLE);
1828 cpl_image_fill_window (wave_image, 1, 1, sizex, sizey, 0.0);
1829
1830 cpl_image * realwave_image = cpl_image_new (sizex, sizey, CPL_TYPE_DOUBLE);
1831 cpl_image_fill_window (realwave_image, 1, 1, sizex, sizey, 0.0);
1832
1833 CPLCHECK_NUL ("Cannot prepare output image");
1834
1835 /* Loop on regions */
1836 for (cpl_size region = 0 ; region < n_region; region ++) {
1837
1838 /* Load the profile of this region as an image */
1839 cpl_imagelist * profile_imglist = gravi_imagelist_wrap_column (profile_table, GRAVI_DATA[region]);
1840 cpl_image * profile_image = cpl_imagelist_get (profile_imglist, 0);
1841
1842 CPLCHECK_NUL ("Cannot get data");
1843
1844 /* Sum all profils of all regions */
1845 cpl_image_add (profilesum_image, profile_image);
1846
1847 /*
1848 * Define an image of each region with its WAVE_DATA
1849 */
1850 const cpl_array * wavelength;
1851 wavelength = cpl_table_get_array (wavedata_table, GRAVI_DATA[region], 0);
1852 CPLCHECK_NUL ("Cannot get data");
1853
1854 for (cpl_size x = 0; x < sizex; x ++){
1855 for (cpl_size y = 0; y < sizey; y ++){
1856 if (cpl_image_get (profile_image, x+1, y+1, &nv) > 0.01)
1857 cpl_image_set (wave_image, x+1, y+1,
1858 cpl_array_get (wavelength, x, NULL));
1859 }
1860 }
1861 CPLCHECK_NUL ("Cannot make image of wave_map");
1862
1863 /*
1864 * Define an image of each region with its WAVE_FIBER
1865 */
1866 int base = gravi_region_get_base (detector_table, region);
1867 int pol = gravi_region_get_pol (detector_table, region);
1868 sprintf (name, "BASE_%s_%s", GRAVI_BASE_NAME[base], GRAVI_POLAR(pol, npol));
1869 wavelength = cpl_table_get_array (wavefibre_table, name, 0);
1870 CPLCHECK_NUL ("Cannot get data");
1871
1872 for (cpl_size x = 0; x < sizex; x ++){
1873 for (cpl_size y = 0; y < sizey; y ++){
1874 if (cpl_image_get (profile_image, x+1, y+1, &nv) > 0.01)
1875 cpl_image_set (realwave_image, x+1, y+1,
1876 cpl_array_get (wavelength, x, NULL));
1877 }
1878 }
1879
1880 gravi_imagelist_unwrap_images (profile_imglist);
1881 CPLCHECK_NUL ("Cannot make image of wave_fibre");
1882 } /* End loop on regions */
1883
1884 /* Create an imagelist */
1885 cpl_imagelist * testwave_imglist = cpl_imagelist_new ();
1886 cpl_imagelist_set (testwave_imglist, wave_image, 0);
1887 cpl_imagelist_set (testwave_imglist, profilesum_image, 1);
1888 cpl_imagelist_set (testwave_imglist, realwave_image, 2);
1889
1890 gravi_msg_function_exit(1);
1891 return testwave_imglist;
1892}
1893
1894/*----------------------------------------------------------------------------*/
1904/*----------------------------------------------------------------------------*/
1905
1906cpl_error_code gravi_wave_qc (gravi_data * wave_map, gravi_data * profile_map)
1907{
1908 gravi_msg_function_start(1);
1909 cpl_ensure_code (wave_map, CPL_ERROR_NULL_INPUT);
1910
1911 char name[100];
1912
1913 cpl_propertylist * wave_header = gravi_data_get_header (wave_map);
1914
1915 /* Loop on extensions (thus SC/FT) */
1916 for (int type_data = 0; type_data < 2; type_data ++ ) {
1917
1918 /* Get WAVE_FIBRE and IMAGING_DETECTOR tables */
1919 cpl_table * detector_table = gravi_data_get_imaging_detector (wave_map, type_data);
1920 cpl_table * wave_data = gravi_data_get_wave_data (wave_map, type_data);
1921 cpl_size n_region = cpl_table_get_nrow (detector_table);
1922
1923 /* Get the full STARTX */
1924 cpl_propertylist * plist = gravi_data_get_wave_data_plist (wave_map, type_data);
1925 int fullstartx = gravi_pfits_get_fullstartx (plist);
1926
1927 CPLCHECK_MSG ("Cannot get data");
1928
1929 /*
1930 * Compute the min and max wave over all regions
1931 */
1932 double minwave = -1e10;
1933 double maxwave = 1e10;
1934
1935 for (int region = 0; region < n_region; region++) {
1936 const cpl_array * wavelength = cpl_table_get_array (wave_data, GRAVI_DATA[region], 0);
1937
1938 minwave = CPL_MAX (minwave, cpl_array_get_min (wavelength));
1939 maxwave = CPL_MIN (maxwave, cpl_array_get_max (wavelength));
1940 } /* End loop on regions */
1941
1942 cpl_msg_info (cpl_func,"%s = %g [m]", QC_MINWAVE(type_data), minwave);
1943 cpl_msg_info (cpl_func,"%s = %g [m]", QC_MAXWAVE(type_data), maxwave);
1944 cpl_propertylist_update_double (wave_header, QC_MINWAVE(type_data), minwave);
1945 cpl_propertylist_update_double (wave_header, QC_MAXWAVE(type_data), maxwave);
1946 cpl_propertylist_update_double (wave_header, QC_MINWAVE_UM(type_data), minwave * 1e6);
1947 cpl_propertylist_update_double (wave_header, QC_MAXWAVE_UM(type_data), maxwave * 1e6);
1948
1949 CPLCHECK_MSG ("Cannot compute minwave or maxwave");
1950
1951
1952 /*
1953 * Compute the pixel position of QC specified argon wavelength.
1954 * Having two lines allow to check the position and dispersion
1955 */
1956 int qc_reg[3] = {0,23,47};
1957 double qc_wave[2] = {2.099184e-6, 2.313952e-6};
1958
1959 /* Loop on regions */
1960 for (int reg = 0; reg < (n_region > 24 ? 3 : 2); reg++) {
1961 cpl_size region = qc_reg[reg];
1962
1963 const cpl_array * wavelength = cpl_table_get_array (wave_data, GRAVI_DATA[region], 0);
1964 cpl_size nwave = cpl_array_get_size (wavelength);
1965 const double * wave_tab = cpl_array_get_data_double_const (wavelength);
1966 CPLCHECK_MSG ("Cannot get wavelength data");
1967
1968 /* Loop on the two argon lines */
1969 for (int iqc = 0 ; iqc < 2 ; iqc++) {
1970 cpl_size l2 = 0;
1971 if ( wave_tab[0] < wave_tab[nwave-1]) {while (wave_tab[l2] < qc_wave[iqc]) l2 ++;}
1972 else {while (wave_tab[l2] > qc_wave[iqc]) l2 ++;}
1973
1974 if (l2-1 < 0 || l2 > nwave-1) {
1975 cpl_msg_error (cpl_func, "Cannot find the QC position for lbd=%g", qc_wave[iqc]);
1976 continue;
1977 }
1978
1979 /* Position on full detector */
1980 double qc_pos = 0.0;
1981 qc_pos = fullstartx + (l2-1) + (qc_wave[iqc] - wave_tab[l2-1]) / (wave_tab[l2] - wave_tab[l2-1]);
1982
1983 sprintf (name, "ESO QC REFWAVE%i", iqc+1);
1984 cpl_propertylist_update_double (wave_header, name, qc_wave[iqc]);
1985 cpl_propertylist_set_comment (wave_header, name, "[m] value of ref wave");
1986
1987 sprintf (name, "ESO QC REFPOS%i %s%lli", iqc+1, GRAVI_TYPE(type_data),region+1);
1988 cpl_propertylist_update_double (wave_header, name, qc_pos);
1989 cpl_propertylist_set_comment (wave_header, name, "[pix] position of ref wave");
1990
1991 cpl_msg_info (cpl_func, "%s = %f [pix] for %e [m]", name, qc_pos, qc_wave[iqc]);
1992
1993 CPLCHECK_MSG ("Cannot set QC");
1994 } /* End loop on the 2 argon lines */
1995 } /* End loop on 2 or 3 regions */
1996 } /* End loop on SC / FT */
1997
1998
1999 /*
2000 * Create the test image for SC (only used for debug)
2001 */
2002 cpl_table * wavedata_table = gravi_data_get_wave_data (wave_map, GRAVI_SC);
2003 cpl_table * wavefibre_table = gravi_data_get_wave_fibre (wave_map, GRAVI_SC);
2004 cpl_table * profile_table = gravi_data_get_table (profile_map, GRAVI_PROFILE_DATA_EXT);
2005 cpl_table * detector_table = gravi_data_get_imaging_detector (wave_map, GRAVI_SC);
2006
2007 cpl_imagelist * testwave_imglist;
2008 testwave_imglist = gravi_wave_test_image (wavedata_table, wavefibre_table,
2009 profile_table, detector_table);
2010 CPLCHECK_MSG ("Cannot compute TEST_WAVE");
2011
2012 /* Set TEST_WAVE in output */
2013 cpl_propertylist * plist = gravi_data_get_plist (profile_map, GRAVI_PROFILE_DATA_EXT);
2014 gravi_data_add_cube (wave_map, cpl_propertylist_duplicate (plist),
2015 "TEST_WAVE", testwave_imglist);
2016
2017 CPLCHECK_MSG ("Cannot set data");
2018
2019 gravi_msg_function_exit(1);
2020 return CPL_ERROR_NONE;
2021}
2022
2023/*----------------------------------------------------------------------------*/
2051/*----------------------------------------------------------------------------*/
2052
2053cpl_error_code gravi_compute_wave (gravi_data * wave_map,
2054 gravi_data * spectrum_data,
2055 int type_data, const cpl_parameterlist * parlist,
2056 gravi_data * wave_param)
2057{
2058 gravi_msg_function_start(1);
2059 cpl_ensure_code (wave_map, CPL_ERROR_NULL_INPUT);
2060 cpl_ensure_code (spectrum_data, CPL_ERROR_NULL_INPUT);
2061 cpl_ensure_code (type_data == GRAVI_SC || type_data == GRAVI_FT,
2062 CPL_ERROR_ILLEGAL_INPUT);
2063
2064 /* Get headers */
2065 cpl_propertylist * raw_header = gravi_data_get_header (spectrum_data);
2066 cpl_propertylist * wave_header = gravi_data_get_header (wave_map);
2067
2068 /* Dump full header of wave data */
2069 if (type_data == GRAVI_FT) {
2070 cpl_propertylist_append (wave_header, raw_header);
2071 }
2072
2073
2074 /* Copy IMAGING_DETECTOR in output WAVE */
2075 const char * extname = (type_data == GRAVI_SC) ? GRAVI_IMAGING_DETECTOR_SC_EXT : GRAVI_IMAGING_DETECTOR_FT_EXT;
2076 gravi_data_copy_ext (wave_map, spectrum_data, extname);
2077
2078 /*
2079 * Compute WAVE_FIBRE map
2080 */
2081 cpl_msg_info (cpl_func, "Compute WAVE_FIBRE for %s", GRAVI_TYPE(type_data));
2082
2083 cpl_table * spectrum_table = gravi_data_get_spectrum_data (spectrum_data, type_data);
2084 cpl_table * detector_table = gravi_data_get_imaging_detector (spectrum_data, type_data);
2085 cpl_table * opd_table = gravi_data_get_table (spectrum_data, type_data==GRAVI_SC ? "OPD_SC" : "OPD_FT");
2086
2087 cpl_table * wavefibre_table;
2088 wavefibre_table = gravi_wave_fibre (spectrum_table, detector_table, opd_table);
2089 CPLCHECK_MSG ("Cannot compute the WAVE_FIBRE");
2090
2091 /*
2092 * Correct WAVE_FIBRE from the dispersion model
2093 */
2094 cpl_msg_info (cpl_func, "Correct dispersion in WAVE_FIBRE of %s", GRAVI_TYPE(type_data));
2095
2096 /* Hardcoded values to correct for the fiber dispersion */
2097 double n0 = 1.0, n1 = -0.0165448, n2 = 0.00256002;
2098 cpl_msg_info (cpl_func,"Rescale wavelengths with dispersion (%g,%g,%g)",n0,n1,n2);
2099
2100 cpl_propertylist_update_string (wave_header, "ESO QC WAVE_CORR", "lbd*(N0+N1*(lbd-lbd0)/lbd0+N2*(lbd-lbd0)^2/lbd0^2)");
2101 cpl_propertylist_update_double (wave_header, "ESO QC WAVE_CORR N0", n0);
2102 cpl_propertylist_update_double (wave_header, "ESO QC WAVE_CORR N1", n1);
2103 cpl_propertylist_update_double (wave_header, "ESO QC WAVE_CORR N2", n2);
2104 CPLCHECK_MSG ("Cannot set Keywords");
2105
2106 gravi_wave_correct_dispersion (wavefibre_table, n0, n1, n2);
2107 CPLCHECK_MSG ("Cannot correct dispersion in wave");
2108
2109 /*
2110 * Fit the 2d dispersion WAVE_FIBRE into WAVE_DATA
2111 */
2112 cpl_msg_info (cpl_func,"Fit WAVE_DATA map for %s", GRAVI_TYPE(type_data));
2113
2114 /* Get the fullstartx */
2115 cpl_propertylist * spectrum_plist = gravi_data_get_spectrum_data_plist (spectrum_data, type_data);
2116 int fullstartx = gravi_pfits_get_fullstartx (spectrum_plist);
2117
2118 /* Interpolate table 2D */
2119 cpl_table * wavedata_table;
2120 int spatial_order=2; // default spatial order
2121 int spectral_order=3; // default spectral order
2122 double rms_residuals;
2123 if (type_data == GRAVI_FT && gravi_param_get_bool(parlist, "gravity.calib.force-wave-ft-equal")) {
2124 spatial_order = 0;
2125 cpl_msg_info (cpl_func, "Option force-waveFT-equal applied");
2126 }
2127 // Keep default value
2128 // if (type_data == GRAVI_SC) {
2129 // spectral_order = gravi_param_get_int(parlist, "gravity.calib.wave-spectral-order");
2130 // cpl_msg_info (cpl_func, "Option set_spectral order to %d", spectral_order);
2131 // }
2132
2133 wavedata_table = gravi_wave_fit_2d (wavefibre_table,
2134 detector_table,
2135 wave_param,
2136 fullstartx, spatial_order, spectral_order, &rms_residuals);
2137 cpl_propertylist_update_double (wave_header, QC_RMS_RESIDUALS(type_data), rms_residuals);
2138 cpl_propertylist_update_double (wave_header, QC_RMS_RESIDUALS_UM(type_data), rms_residuals * 1e6);
2139
2140
2141 double rms_fit=cpl_propertylist_get_double (raw_header, QC_PHASECHI2);
2142 cpl_propertylist_update_double (wave_header, QC_CHI2WAVE(type_data),rms_fit*1e9);
2143 cpl_propertylist_set_comment (wave_header, QC_CHI2WAVE(type_data),"[nm]rms a.SC-b.FT+c=MET");
2144
2145 CPLCHECK_MSG ("Cannot fit 2d data");
2146
2147 /*
2148 * New Wavelength interpolation made by sylvestre on January 30 2018
2149 */
2150
2151 if (type_data == GRAVI_SC && !strcmp (gravi_pfits_get_spec_res(raw_header), "LOW"))
2152 {
2153 cpl_msg_info (cpl_func, "Additional Wavelength Fit");
2154 cpl_table * wave_individual_table = cpl_table_new (1);
2155 cpl_table * weight_individual_table = cpl_table_new (1);
2156 cpl_table * wave_fitted_table = cpl_table_new (1);
2157
2158 gravi_wave_fit_individual (wave_individual_table,
2159 weight_individual_table,
2160 wave_fitted_table,
2161 opd_table,
2162 spectrum_table,
2163 detector_table,
2164 n0,n1,n2);
2165
2166 cpl_msg_info (cpl_func,"Add tables in wave_map");
2167
2168 gravi_data_add_table (wave_map, cpl_propertylist_duplicate (spectrum_plist),
2169 "WAVE_INDIV_SC", wave_individual_table);
2170 gravi_data_add_table (wave_map, cpl_propertylist_duplicate (spectrum_plist),
2171 "WAVE_WEIGHT_SC", weight_individual_table);
2172 gravi_data_add_table (wave_map, cpl_propertylist_duplicate (spectrum_plist),
2173 "WAVE_FITTED_SC", wavedata_table);
2174
2175 gravi_data_add_table (wave_map, cpl_propertylist_duplicate (spectrum_plist),
2176 GRAVI_WAVE_FIBRE_EXT(type_data), wavefibre_table);
2177 gravi_data_add_table (wave_map, cpl_propertylist_duplicate (spectrum_plist),
2178 GRAVI_WAVE_DATA_EXT(type_data), wave_fitted_table);
2179 } else {
2180 /*
2181 * Add the WAVE_FIBRE and WAVE_DATA table in the wave_map
2182 */
2183 cpl_msg_info (cpl_func,"Add WAVE_FIBRE and WAVE_DATA in wave_map");
2184
2185 gravi_data_add_table (wave_map, cpl_propertylist_duplicate (spectrum_plist),
2186 GRAVI_WAVE_FIBRE_EXT(type_data), wavefibre_table);
2187
2188 gravi_data_add_table (wave_map, cpl_propertylist_duplicate (spectrum_plist),
2189 GRAVI_WAVE_DATA_EXT(type_data), wavedata_table);
2190
2191 }
2192
2193
2194
2195
2196 CPLCHECK_MSG ("Cannot set data");
2197
2198 gravi_msg_function_exit(1);
2199 return CPL_ERROR_NONE;
2200}
2201
2202
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
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:88
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:165
QC_MAXWAVE
#define QC_MAXWAVE(type)
Definition gravi_pfits.h:110
QC_MINWAVE_UM
#define QC_MINWAVE_UM(type)
Definition gravi_pfits.h:111
GRAVI_PROFILE_DATA_EXT
#define GRAVI_PROFILE_DATA_EXT
Definition gravi_pfits.h:79
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:112
OPD_COEFF_SIGN
#define OPD_COEFF_SIGN(type)
Definition gravi_pfits.h:115
QC_CHI2WAVE
#define QC_CHI2WAVE(type)
Definition gravi_pfits.h:108
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:114
GRAVI_TYPE
#define GRAVI_TYPE(type)
Definition gravi_pfits.h:167
QC_MINWAVE
#define QC_MINWAVE(type)
Definition gravi_pfits.h:109
LAMBDA_MET
#define LAMBDA_MET
Definition gravi_pfits.h:103
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:166
QC_PHASECHI2
#define QC_PHASECHI2
Definition gravi_pfits.h:107
QC_RMS_RESIDUALS
#define QC_RMS_RESIDUALS(type)
Definition gravi_pfits.h:113
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:1534
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:1665
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:1609
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:445
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:361
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:1013
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:1018
gravi_compute_envelope
cpl_vector * gravi_compute_envelope(const cpl_vector *opd, int wave, int nwave)
Compute the envelope value.
Definition gravi_utils.c:1081
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:584
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:399
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:79
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:697
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:458
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:988
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:1178
PLOT_WAVE_PHASE_VS_OPD
#define PLOT_WAVE_PHASE_VS_OPD
Definition gravi_wave.c:78
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:819
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:353
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:2053
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:148
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:1427
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:1906
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:1703
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:1748
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:1801
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:565
Generated on Tue Apr 28 2026 07:47:55 for GRAVI Pipeline Reference Manual by doxygen 1.13.2