GRAVI Pipeline Reference Manual 1.10.1
Loading...
Searching...
No Matches
gravi_p2vmred.c
Go to the documentation of this file.
1/* $Id: gravi_tf.c,v 1.10 2014/11/12 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
32
33/*
34 * History :
35 * 11/01/2019 Fix Warnings , unused parameter : npol_ft, npol_sc, p2vmreducedFlag ,init_type_data, end_type_data
36 */
37
38/*-----------------------------------------------------------------------------
39 Includes
40 -----------------------------------------------------------------------------*/
41
42#ifdef HAVE_CONFIG_H
43#include <config.h>
44#endif
45
46#include <cpl.h>
47#include <string.h>
48#include <stdio.h>
49#include <math.h>
50#include <time.h>
51#include <complex.h>
52
53#include "gravi_data.h"
54#include "gravi_dfs.h"
55#include "gravi_pfits.h"
56#include "gravi_cpl.h"
57
58#include "gravi_utils.h"
59
60#include "gravi_p2vmred.h"
61#include "gravi_vis.h"
62
63/*-----------------------------------------------------------------------------
64 Private prototypes
65 -----------------------------------------------------------------------------*/
66
67cpl_table * gravi_create_oitarget_table (const cpl_propertylist *header,
68 const char * mode);
69
70cpl_table * gravi_create_oiarray_table (const cpl_table * array_geometry,
71 int is_cal);
72
73/*-----------------------------------------------------------------------------
74 Function code
75 -----------------------------------------------------------------------------*/
76
77/*----------------------------------------------------------------------------*/
88/*----------------------------------------------------------------------------*/
89
90cpl_table * gravi_create_oitarget_table (const cpl_propertylist *header,
91 const char * mode)
92{
93 /* Message and timer */
94 gravi_msg_function_start(1);
95 cpl_ensure (header, CPL_ERROR_NULL_INPUT, NULL);
96 cpl_ensure (mode, CPL_ERROR_NULL_INPUT, NULL);
97
98 /* Check if this observation has a reference object in FT and a science object in SC */
99 int n_target = 1;
100 if (!(strcmp (mode, "gravi_dual"))) n_target = 2;
101 else if (!(strcmp (mode, "gravi_single"))) n_target = 1;
102 else {cpl_error_set_message(cpl_func, CPL_ERROR_NULL_INPUT, "Invalid mode (bug)");return NULL;}
103
104 /* Create the table */
105 cpl_table * oi_target = cpl_table_new (n_target);
106
107 /* TARGET_ID */
108 cpl_table_new_column (oi_target, "TARGET_ID", CPL_TYPE_INT);
109 cpl_table_set_column_savetype (oi_target, "TARGET_ID", CPL_TYPE_SHORT);
110 /* TARGET */
111 cpl_table_new_column (oi_target, "TARGET", CPL_TYPE_STRING);
112 /* EQUINOX */
113 cpl_table_new_column (oi_target, "EQUINOX", CPL_TYPE_FLOAT);
114 cpl_table_set_column_unit (oi_target, "EQUINOX", "yr");
115 /* RAEP0 and DECEP0 */
116 cpl_table_new_column (oi_target, "RAEP0", CPL_TYPE_DOUBLE);
117 cpl_table_set_column_unit (oi_target, "RAEP0", "deg");
118 cpl_table_new_column (oi_target, "DECEP0", CPL_TYPE_DOUBLE);
119 cpl_table_set_column_unit (oi_target, "DECEP0", "deg");
120 /* RA_ERR and DEC_ERR */
121 cpl_table_new_column (oi_target, "RA_ERR", CPL_TYPE_DOUBLE);
122 cpl_table_set_column_unit (oi_target, "RA_ERR", "deg");
123 cpl_table_new_column (oi_target, "DEC_ERR", CPL_TYPE_DOUBLE);
124 cpl_table_set_column_unit (oi_target, "DEC_ERR", "deg");
125 /* SYSVEL, VELTYP, VELDEF */
126 cpl_table_new_column (oi_target, "SYSVEL", CPL_TYPE_DOUBLE);
127 cpl_table_set_column_unit (oi_target, "SYSVEL", "m/s");
128 cpl_table_new_column (oi_target, "VELTYP", CPL_TYPE_STRING);
129 cpl_table_new_column (oi_target, "VELDEF", CPL_TYPE_STRING);
130 /* PMRA, PMDEC and associated errors */
131 cpl_table_new_column (oi_target, "PMRA", CPL_TYPE_DOUBLE);
132 cpl_table_set_column_unit (oi_target, "PMRA", "deg/yr");
133 cpl_table_new_column (oi_target, "PMDEC", CPL_TYPE_DOUBLE);
134 cpl_table_set_column_unit (oi_target, "PMDEC", "deg/yr");
135 cpl_table_new_column (oi_target, "PMRA_ERR", CPL_TYPE_DOUBLE);
136 cpl_table_set_column_unit (oi_target, "PMRA_ERR", "deg/yr");
137 cpl_table_new_column (oi_target, "PMDEC_ERR", CPL_TYPE_DOUBLE);
138 cpl_table_set_column_unit (oi_target, "PMDEC_ERR", "deg/yr");
139 /* PARALLAX and PARA_ERR */
140 cpl_table_new_column (oi_target, "PARALLAX", CPL_TYPE_FLOAT);
141 cpl_table_set_column_unit (oi_target, "PARALLAX", "deg");
142 cpl_table_new_column (oi_target, "PARA_ERR", CPL_TYPE_FLOAT);
143 cpl_table_set_column_unit (oi_target, "PARA_ERR", "deg");
144 /* SPECTYP */
145 cpl_table_new_column (oi_target, "SPECTYP", CPL_TYPE_STRING);
146 CPLCHECK_NUL("Cannot create columns in OI_TARGET");
147
148 /* loop on target to fill the columns*/
149 for (int ti=0; ti<n_target; ti++) {
150
151 /* TARGET_ID=1 is for FT and TARGET_ID=2 is for SC (if any) */
152 cpl_table_set_int (oi_target, "TARGET_ID", ti, ti+1);
153
154 /* TARGET name, forced to match 16 chars */
155 const char *name = (ti?gravi_pfits_get_sobj(header):gravi_pfits_get_robj(header));
156 gravi_table_set_string_fixlen (oi_target, "TARGET", ti, name, 16);
157
158 CPLCHECK_NUL("Cannot set Target");
159
160 /* At the moment, we can't fill the following columns. We need to put something for
161 * standard compliance. Some third party software fail if those columns are missing.
162 * Should be checked whether better values are available. */
163
164 /* EQUINOX, RAEP0 and DECEP0 at mean equinox, for FT or SC */
165 float equinox = (float)gravi_pfits_get_double_silentdefault (header, "EQUINOX", 0.);
166 double raep = (ti?gravi_pfits_get_sobj_raep (header):gravi_pfits_get_robj_raep (header));
167 double decp = (ti?gravi_pfits_get_sobj_decep (header):gravi_pfits_get_robj_decep (header));
168 if (ti==0) gravi_dump_the_boss (raep, decp);
169
170 double ra_err = 0.0;
171 double dec_err = 0.0;
172 raep *= CPL_MATH_DEG_RAD; // [deg]
173 decp *= CPL_MATH_DEG_RAD; // [deg]
174 CPLCHECK_NUL("Cannot get EQUINOX, RA and DEC");
175
176 /* Some verbose */
177 cpl_msg_info (cpl_func,"Found target '%s' with ra=%fd and dec=%fd at equinox %.2f", name, raep, decp, equinox);
178
179 /* SYSVEL, VELTYP, VELDEF */
180 double sysvel = gravi_pfits_get_double_silentdefault (header, "ESO FT ROBJ RADVEL", 0.);
181 const char * veltyp = "UNKNOWN ";
182 const char * veldef = "OPTICAL ";
183 CPLCHECK_NUL("Cannot get SYSVEL, VELTYPE...");
184
185 /* PMRA, PMDEC and associated errors in deg/years
186 * These quantities are in as/years in the HEADER */
187 double pmra = gravi_pfits_get_pmra (header) / 3600.0; // [deg/year]
188 double pmdec = gravi_pfits_get_pmdec (header) / 3600.0; // [deg/year]
189 double pmra_err = 0.0;
190 double pmdec_err = 0.0;
191 CPLCHECK_NUL("Cannot get PMRA...");
192
193 /* PARALLAX and PARA_ERR */
194 float parallax = gravi_pfits_get_plx (header); // [as]
195 float para_err = 0.0;
196 CPLCHECK_NUL("Cannot get PARRALAX...");
197
198 /* SPECTYP */
199 const char * spectyp = "UNKNOWN ";
200
201 /* Fill common columns */
202 cpl_table_set_float (oi_target, "EQUINOX", ti, equinox);
203 cpl_table_set_double (oi_target, "RAEP0", ti, raep);
204 cpl_table_set_double (oi_target, "DECEP0", ti, decp);
205 cpl_table_set_double (oi_target, "RA_ERR", ti, ra_err);
206 cpl_table_set_double (oi_target, "DEC_ERR", ti, dec_err);
207 cpl_table_set_double (oi_target, "SYSVEL", ti, sysvel);
208 cpl_table_set_string (oi_target, "VELTYP", ti, veltyp);
209 cpl_table_set_string (oi_target, "VELDEF", ti, veldef);
210 cpl_table_set_double (oi_target, "PMRA", ti, pmra);
211 cpl_table_set_double (oi_target, "PMDEC", ti, pmdec);
212 cpl_table_set_double (oi_target, "PMRA_ERR", ti, pmra_err);
213 cpl_table_set_double (oi_target, "PMDEC_ERR",ti, pmdec_err);
214 cpl_table_set_float (oi_target, "PARALLAX", ti, parallax);
215 cpl_table_set_float (oi_target, "PARA_ERR", ti, para_err);
216 cpl_table_set_string (oi_target, "SPECTYP", ti, spectyp);
217 }
218 /* End loop on target_id */
219
220 gravi_msg_function_exit(1);
221 return oi_target;
222}
223
224/*----------------------------------------------------------------------------*/
236/*----------------------------------------------------------------------------*/
237
238cpl_table * gravi_create_oiarray_table (const cpl_table * array_geometry,
239 int is_cal)
240{
241 gravi_msg_function_start(1);
242 cpl_ensure (array_geometry, CPL_ERROR_NULL_INPUT, NULL);
243
244 cpl_table * oi_array = cpl_table_duplicate (array_geometry);
245
246 /* Rename column STA to STAXYZ */
247 if (cpl_table_has_column (oi_array, "STA") ) {
248 cpl_table_name_column (oi_array, "STA", "STAXYZ");
249 }
250 CPLCHECK_NUL ("Cannot rename STA into STAXYZ");
251
252 /* Make sure there are 4 rows (one per beam) */
253 cpl_size n_row = cpl_table_get_nrow (oi_array);
254 if (n_row < 4) {
255
256 if (is_cal) cpl_msg_info (cpl_func, "ARRAY_GEOMETRY has %lld rows instead of 4", n_row);
257 else cpl_msg_warning (cpl_func, "ARRAY_GEOMETRY has %lld rows instead of 4", n_row);
258 cpl_table_set_size (oi_array, 4);
259
260 cpl_array * zero_array = gravi_array_init_double (3,0.0);
261
262 /* Fill these new rows with default values. */
263 for ( int row=n_row; row<4; row++ ) {
264 char name16[17];
265 cpl_table_set_int (oi_array, "STA_INDEX", row, 100+row);
266 sprintf (name16, "S%i", row);
267 cpl_table_set_string (oi_array, "STA_NAME", row, name16);
268 sprintf (name16, "T%i", row);
269 cpl_table_set_string (oi_array, "TEL_NAME", row, name16);
270 cpl_table_set_array (oi_array, "STAXYZ", row, zero_array);
271 cpl_table_set_float (oi_array, "DIAMETER", row, 0.0);
272 cpl_table_set_int (oi_array, "MNTSTA", row, 0);
273 if (is_cal) cpl_msg_info (cpl_func,"Add STA_INDEX %i with TEL_NAME=%s in OI_ARRAY",100+row,name16);
274 else cpl_msg_warning (cpl_func,"Add STA_INDEX %i with TEL_NAME=%s in OI_ARRAY",100+row,name16);
275 CPLCHECK_NUL ("Cannot insert rows");
276 }
277
278 FREE (cpl_array_delete, zero_array);
279
280 }/* End case there are missing rows */
281
282 /* Update TEL_NAME and STA_NAME to use 16 char for OIFITS compliance */
283 for ( int row=0; row<cpl_table_get_nrow (oi_array); row++ ) {
284 const char * tel_name = cpl_table_get_string (oi_array, "TEL_NAME", row);
285 gravi_table_set_string_fixlen (oi_array, "TEL_NAME", row, tel_name, 16);
286 const char * sta_name = cpl_table_get_string (oi_array, "STA_NAME", row);
287 gravi_table_set_string_fixlen (oi_array, "STA_NAME", row, sta_name, 16);
288 }
289 CPLCHECK_NUL ("Cannot force 16 chars in OI_ARRAY");
290
291 gravi_msg_function_exit(1);
292 return oi_array;
293}
294
295/*----------------------------------------------------------------------------*/
328/*----------------------------------------------------------------------------*/
329
330gravi_data * gravi_compute_p2vmred (gravi_data * preproc_data, gravi_data * p2vm_map,
331 const char * mode, const cpl_parameterlist * parlist)
332{
333 /* Message and timer */
334 gravi_msg_function_start(1);
335 cpl_ensure (preproc_data, CPL_ERROR_NULL_INPUT, NULL);
336 cpl_ensure (p2vm_map, CPL_ERROR_NULL_INPUT, NULL);
337 cpl_ensure (mode, CPL_ERROR_NULL_INPUT, NULL);
338 cpl_ensure (parlist, CPL_ERROR_NULL_INPUT, NULL);
339
340 char qc_name[90];
341 int ntel = 4, nclo = 4 ;
342 int nv; /* npol_ft, npol_sc; */
343
344 /* Timers */
345 clock_t timer_ft = 0, timer_sc = 0, timer_other = 0, timer_start = 0;
346 timer_start = clock();
347
348 /* Global properties */
349 cpl_propertylist * preproc_header = gravi_data_get_header (preproc_data);
350
351 /* Flag if the p2vmreduced-file is requested */
352 /* int p2vmreducedFlag = 0; */
353 if ( (cpl_parameterlist_find_const (parlist, "gravity.dfs.p2vmred-file") != NULL &&
354 gravi_param_get_bool (parlist,"gravity.dfs.p2vmred-file")) ||
355 (cpl_parameterlist_find_const (parlist, "gravity.dfs.astro-file") != NULL &&
356 gravi_param_get_bool (parlist,"gravity.dfs.astro-file")) ) {
357 cpl_msg_info (cpl_func, "p2vmreduced-file is requested, will add computation time.");
358 /* p2vmreducedFlag = 1; */
359 } else {
360 cpl_msg_info (cpl_func, "p2vmreduced-file is not requested.");
361 }
362
363 CPLCHECK_NUL("Cannot check the parameter");
364
365 /* Get this flag if internal calib or on-sky */
366 int is_cal = gravi_pfits_is_calib (preproc_header);
367
368 /* Build the output gravi_data. */
369 gravi_data * p2vmred_data = gravi_data_new (0);
370 cpl_propertylist * p2vmred_header = gravi_data_get_header (p2vmred_data);
371
372 /* Dump all header from RAW product */
373 cpl_propertylist_append (p2vmred_header, preproc_header);
374
375 /* Duplicate necessary tables in product */
376 gravi_data_copy_ext (p2vmred_data, preproc_data, GRAVI_OI_WAVELENGTH_EXT);
377
378 /* Create OI_TARGET from header and add it to product */
379 cpl_table * oi_target = gravi_create_oitarget_table (preproc_header, mode);
380 gravi_data_add_table (p2vmred_data, NULL, GRAVI_OI_TARGET_EXT, oi_target);
381 CPLCHECK_NUL("Cannot create OI_TARGET");
382
383 /* Duplicate ARRAY_GEOMETRY as OI_ARRAY */
384 cpl_table * array_geo = gravi_data_get_table (preproc_data, GRAVI_ARRAY_GEOMETRY_EXT);
385 CPLCHECK_NUL("Cannot get ARRAY_GEOMETRY");
386 cpl_table * oi_array = gravi_create_oiarray_table (array_geo, is_cal);
387 CPLCHECK_NUL("Cannot convert ARRAY_GEOMETRY into OI_ARRAY");
388
389 /* Build the header of OI_ARRAY from the header of OI_ARRAY */
390 cpl_propertylist * oiarray_plist = gravi_data_get_plist (preproc_data, GRAVI_ARRAY_GEOMETRY_EXT);
391 oiarray_plist = cpl_propertylist_duplicate (oiarray_plist);
392
393 /* Set OI_ARRAY in output data */
394 gravi_data_add_table (p2vmred_data, oiarray_plist,
395 GRAVI_OI_ARRAY_EXT, oi_array);
396 CPLCHECK_NUL ("Cannot add OI_ARRAY");
397
398 /* Update the keyword ARRAYX/Y/Z to ensure double */
399 cpl_msg_info (cpl_func, "Update the ARRAYX/Y/Z keywords as double precision");
400 gravi_pfits_ensure_double (oiarray_plist, "ARRAYX");
401 gravi_pfits_ensure_double (oiarray_plist, "ARRAYY");
402 gravi_pfits_ensure_double (oiarray_plist, "ARRAYZ");
403
404 /* Read OPTICAL_TRAIN */
405 cpl_table * optical_train_table = gravi_data_get_table (preproc_data, GRAVI_OPTICAL_TRAIN_EXT);
406
407 /* Loop on lab input corresponding to GRAVITY to make sure all beams exists
408 * If they don't, they are filled with fake value to macth the OI_ARRAY */
409 int labInput[4] = {GRAVI_LABINPUT_1,GRAVI_LABINPUT_2,GRAVI_LABINPUT_3,GRAVI_LABINPUT_4};
410 for (int lab = 0 ; lab < 4 ; lab++) {
411
412 /* Search if this LABINIT already exist */
413 cpl_table_unselect_all (optical_train_table);
414 if ( cpl_table_or_selected_int (optical_train_table, "VALUE2", CPL_EQUAL_TO, labInput[lab]) == 1 ) {
415 continue;
416 }
417
418 /* If not, add a rows and fill it with default values */
419 char name16[17];
420 int nrow = cpl_table_get_nrow (optical_train_table);
421 cpl_table_set_size (optical_train_table, nrow + 1);
422 cpl_table_set_int (optical_train_table, "VALUE2", nrow, labInput[lab]);
423 sprintf (name16, "T%i", nrow);
424 cpl_table_set_string (optical_train_table, "TEL_NAME", nrow, name16);
425 if (is_cal) cpl_msg_info (cpl_func,"Add LABINPUT %i with TEL_NAME=%s in OPTICAL_TRAIN",labInput[lab],name16);
426 else cpl_msg_warning (cpl_func,"Add LABINPUT %i with TEL_NAME=%s in OPTICAL_TRAIN",labInput[lab],name16);
427 CPLCHECK_NUL("Cannot insert rows");
428 } /* End loop in labInput */
429
430 /*
431 * For each type of data SC / FT
432 */
433 /* int init_type_data = 1;
434 int end_type_data = 1;
435
436 if(det_type == GRAVI_DET_SC || det_type == GRAVI_DET_ALL)
437 init_type_data = 0;
438 if(det_type == GRAVI_DET_FT || det_type == GRAVI_DET_ALL)
439 end_type_data = 1;
440 */
441 for (int type_data = 0; type_data < 2; type_data ++){
442
443 /* Check if this data is present */
444 if (!gravi_data_has_extension (preproc_data, GRAVI_SPECTRUM_DATA_EXT(type_data))) {
445 cpl_msg_info (cpl_func,"No data for %s, skip it", GRAVI_TYPE(type_data));
446 continue;
447 }
448
449 /* Get the table and property list utilities */
450 cpl_table * p2vm_table = gravi_data_get_p2vm_data (p2vm_map, type_data);
451 cpl_table * spectrum_table = gravi_data_get_spectrum_data (preproc_data, type_data);
452 cpl_table * detector_table = gravi_data_get_imaging_detector (preproc_data, type_data);
453 cpl_table * detector_p2vm = gravi_data_get_imaging_detector (p2vm_map, type_data);
454 CPLCHECK_NUL ("Cannot get data");
455
456 /* Verify the data and P2VM are conformable */
457 cpl_ensure (cpl_table_get_nrow (detector_table) ==
458 cpl_table_get_nrow (detector_p2vm),
459 CPL_ERROR_ILLEGAL_INPUT, NULL);
460
461 /* Read nrow and nregion */
462 cpl_size nrow = cpl_table_get_nrow (spectrum_table);
463 int n_detregion = cpl_table_get_nrow (detector_table);
464
465 /* Guess the number of polarisation by looking at the number of region
466 * Save the number of polarisation for further use */
467 int npol = (n_detregion>24 ? 2 : 1);
468 /*if (type_data == GRAVI_FT) npol_ft = npol;
469 if (type_data == GRAVI_SC) npol_sc = npol; */
470
471 /* FIXME: Here we assume the two polarisation (if any)
472 * have the same number of channels ? */
473 cpl_size nwave = gravi_spectrum_get_nwave (spectrum_table);
474
475 /*
476 * Loop on polarisations
477 */
478 for (int p = 0; p < npol; p ++) {
479
480 cpl_msg_info(cpl_func, "Visibility extraction: polarisation %d over %d for %s",
481 p+1, npol, type_data==GRAVI_FT?"FT":"SC");
482
483 timer_other += (clock() - timer_start);
484 timer_start = clock();
485
486 /*
487 *
488 * Create output tables and set them in output data
489 *
490 */
491 cpl_msg_info(cpl_func, "Create and set output tables...");
492
493 /* Create table */
494 cpl_table * oi_vis = gravi_table_oi_create (nwave, nrow, GRAVI_OI_VIS_EXT);
495 cpl_table * oi_t3 = gravi_table_oi_create (nwave, nrow, GRAVI_OI_T3_EXT);
496 cpl_table * oi_flux = gravi_table_oi_create (nwave, nrow, GRAVI_OI_FLUX_EXT);
497
498 /* Delete useless columns in this product */
499 cpl_table_erase_column (oi_vis, "VISAMP");
500 cpl_table_erase_column (oi_vis, "VISPHI");
501 cpl_table_erase_column (oi_vis, "VISAMPERR");
502 cpl_table_erase_column (oi_vis, "VISPHIERR");
503 cpl_table_erase_column (oi_vis, "RVIS");
504 cpl_table_erase_column (oi_vis, "RVISERR");
505 cpl_table_erase_column (oi_vis, "IVIS");
506 cpl_table_erase_column (oi_vis, "IVISERR");
507
508 cpl_table_erase_column (oi_t3, "T3AMP");
509 cpl_table_erase_column (oi_t3, "T3PHI");
510 cpl_table_erase_column (oi_t3, "T3AMPERR");
511 cpl_table_erase_column (oi_t3, "T3PHIERR");
512 cpl_table_erase_column (oi_t3, "U1COORD");
513 cpl_table_erase_column (oi_t3, "V1COORD");
514 cpl_table_erase_column (oi_t3, "U2COORD");
515 cpl_table_erase_column (oi_t3, "V2COORD");
516
517 /* Set table in p2vmred */
518 cpl_propertylist * oi_plist = cpl_propertylist_new ();
519 cpl_propertylist_copy_property (oi_plist, oiarray_plist, "ARRNAME");
520 cpl_propertylist_copy_property (oi_plist, preproc_header, "DATE-OBS");
521 cpl_propertylist_update_int (oi_plist, "OI_REVN", 1);
522 cpl_propertylist_update_int (oi_plist, "NWAVE", nwave);
523 cpl_propertylist_update_string (oi_plist, "INSNAME", GRAVI_INSNAME(type_data,p,npol));
524 cpl_propertylist_update_int (oi_plist, "EXTVER", GRAVI_EXTVER(type_data,p,npol));
525
526 gravi_data_add_table (p2vmred_data, oi_plist,
527 GRAVI_OI_VIS_EXT, oi_vis);
528
529 oi_plist = cpl_propertylist_duplicate (oi_plist);
530
531 gravi_data_add_table (p2vmred_data, oi_plist,
532 GRAVI_OI_T3_EXT, oi_t3);
533
534 oi_plist = cpl_propertylist_duplicate (oi_plist);
535
536 gravi_data_add_table (p2vmred_data, oi_plist,
537 GRAVI_OI_FLUX_EXT, oi_flux);
538
539 cpl_msg_info(cpl_func, "Total time to create tables: %.4f s", (double)(clock()-timer_start)/(double)CLOCKS_PER_SEC );
540
541 /*
542 *
543 * Compute the P2VM and the V2PM for all wavelengths
544 *
545 */
546 timer_other += (clock() - timer_start);
547 timer_start = clock();
548
549 cpl_msg_info (cpl_func, "Compute the invers V2P2M -> P2VM matrix...");
550
551 /* Construction of output tables */
552 cpl_matrix ** p2vm = cpl_calloc (nwave,sizeof (cpl_matrix*));
553 cpl_matrix ** v2pm = cpl_calloc (nwave,sizeof (cpl_matrix*));
554
555 /* Get the list of regions having this pol */
556 int all_region[48];
557 int n_region = 0;
558 for (int detreg = 0; detreg < n_detregion; detreg++) {
559 if (gravi_region_get_pol (detector_p2vm, detreg) != p) continue;
560 all_region[n_region] = detreg;
561 n_region++;
562 }
563
564 /* Get direct pointer to arrays */
565 const cpl_array ** trans, ** coh, ** phase;
566 trans = cpl_table_get_data_array_const (p2vm_table, TRANSMISSION);
567 coh = cpl_table_get_data_array_const (p2vm_table, COHERENCE);
568 phase = cpl_table_get_data_array_const (p2vm_table, PHASE);
569 CPLCHECK_NUL ("Cannot load the p2vm data");
570
571 /* Allocate memory for direct access to P2VM */
572 double ** pP2VM = cpl_malloc (sizeof(double*) * nwave);
573 double ** pV2PM = cpl_malloc (sizeof(double*) * nwave);
574
575 /* Loop on wave */
576 for (cpl_size wave = 0; wave < nwave; wave ++) {
577
578 /* Construction of the v2pm matrix */
579 v2pm[wave] = cpl_matrix_new (n_region, 16);
580
581 /* Loop on region */
582 for (cpl_size region = 0; region < n_region; region++) {
583 int detregion = all_region[region];
584
585 /* Replace the four first columns of the v2pm by
586 * the transmission */
587 for (int i = 0; i < 4; i++){
588 cpl_matrix_set (v2pm[wave], region, i,
589 cpl_array_get (trans[detregion], wave + i * nwave, &nv));
590 }
591
592 /* Replace the twelve other columns of the v2pm by
593 * the real part and the imaginary part of the coherence
594 * and the phase */
595 for (int i = 0; i < 6; i++) {
596 cpl_matrix_set (v2pm[wave], region, i + 4,
597 cpl_array_get (coh[detregion], wave + i * nwave, &nv) *
598 cos(cpl_array_get (phase[detregion], wave + i * nwave, &nv)));
599 cpl_matrix_set (v2pm[wave], region, i + 10,
600 cpl_array_get (coh[detregion], wave + i * nwave, &nv) *
601 sin(cpl_array_get (phase[detregion], wave + i * nwave, &nv)));
602 }
603 }
604 CPLCHECK_NUL("Cannot fill the V2PM");
605
606 /* Ensure the V2PM is flux conservative */
607 cpl_matrix_multiply_scalar (v2pm[wave], 1./n_region);
608
609 /* Compute the matrix inversion of the v2pm using the
610 * singular value decomposition method */
611 p2vm[wave] = gravi_matrix_invertSV_create (v2pm[wave]);
612
613 /* Keep a pointer to the data */
614 pP2VM[wave] = cpl_matrix_get_data (p2vm[wave]);
615 pV2PM[wave] = cpl_matrix_get_data (v2pm[wave]);
616
617 CPLCHECK_NUL ("Cannot invers V2PM");
618 }
619 /* End loop on wavelength */
620
621 cpl_msg_info (cpl_func, "Total time to invers matrix: %.4f s", (double)(clock()-timer_start)/(double)CLOCKS_PER_SEC );
622
623 /*
624 * Fill the STA_INDEX and TIME.
625 */
626 cpl_msg_info (cpl_func, "Fill the STA_INDEX and TIME.");
627 timer_other += (clock() - timer_start);
628 timer_start = clock();
629
630 /* Wrap the TIME column into an array. TIME is in [us] from
631 * the PCR.ACQ.START (start of RMN recording) */
632 cpl_array * times = cpl_array_wrap_int (cpl_table_get_data_int (spectrum_table, "TIME"),
633 cpl_table_get_nrow (spectrum_table));
634 cpl_ensure (times, CPL_ERROR_ILLEGAL_INPUT, NULL);
635
636 /* Compute the MJDs for each row, by adding PCR.ACQ.START
637 * (in MJD) to the TIME column */
638 cpl_array * mjds = cpl_array_cast (times, CPL_TYPE_DOUBLE);
639 cpl_array_divide_scalar (mjds, 86400.E6);
640 double mjd0 = gravi_convert_to_mjd (gravi_pfits_get_start_prcacq (preproc_header));
641 cpl_array_add_scalar (mjds, mjd0);
642
643 /* We keep the TIME column of the P2VMRED in [us] from
644 * the PCR.ACQ.START, in order to correlate with RMN tables.
645 * Thus this TIME column is *not* at the OIFITS standart */
646 cpl_table_set_column_unit (oi_vis, "TIME", "us");
647 cpl_table_set_column_unit (oi_t3, "TIME", "us");
648 cpl_table_set_column_unit (oi_flux, "TIME", "us");
649
650 /* Fill STA_INDEX and times for OI_VIS */
651 cpl_array * sta_index = cpl_array_new (2, CPL_TYPE_INT);
652
653 for (int base=0; base < GRAVI_NBASE; ++base) {
654 /* Build sta_index */
655 int sta0 = gravi_sta_index(GRAVI_BASE_TEL[base][0]+1, optical_train_table, oi_array);
656 int sta1 = gravi_sta_index(GRAVI_BASE_TEL[base][1]+1, optical_train_table, oi_array);
657 cpl_array_set_int (sta_index, 0, sta0);
658 cpl_array_set_int (sta_index, 1, sta1);
659 CPLCHECK_NUL ("Cannot find the sta_index");
660
661 /* loop on rows */
662 for (int row=0; row < nrow; ++row) {
663 int idx = row * GRAVI_NBASE + base;
664 cpl_table_set_array (oi_vis, "STA_INDEX", idx, sta_index);
665 cpl_table_set (oi_vis, "TIME", idx, cpl_array_get (times, row, &nv));
666 cpl_table_set (oi_vis, "MJD", idx, cpl_array_get (mjds, row, &nv));
667 }
668 } /* End loop on base */
669
670 cpl_array_delete (sta_index);
671 CPLCHECK_NUL ("Cannot fill sta_index or time in OI_VIS");
672
673 /* Fill STA_INDEX and times for OI_T3 */
674 sta_index = cpl_array_new (3, CPL_TYPE_INT);
675
676 for (int closure=0; closure < nclo; ++closure) {
677 /* Build sta_index */
678 int sta0 = gravi_sta_index(GRAVI_CLO_TEL[closure][0]+1, optical_train_table, oi_array);
679 int sta1 = gravi_sta_index(GRAVI_CLO_TEL[closure][1]+1, optical_train_table, oi_array);
680 int sta2 = gravi_sta_index(GRAVI_CLO_TEL[closure][2]+1, optical_train_table, oi_array);
681 cpl_array_set_int (sta_index, 0, sta0);
682 cpl_array_set_int (sta_index, 1, sta1);
683 cpl_array_set_int (sta_index, 2, sta2);
684 CPLCHECK_NUL ("Cannot find the sta_index");
685
686 /* loop on rows */
687 for (int row=0; row < nrow; ++row) {
688 int idx = row * nclo + closure;
689 cpl_table_set_array (oi_t3, "STA_INDEX", idx, sta_index);
690 cpl_table_set (oi_t3, "TIME", idx, cpl_array_get (times, row, &nv));
691 cpl_table_set (oi_t3, "MJD", idx, cpl_array_get (mjds, row, &nv));
692 }
693 } /* End loop on closures */
694
695 cpl_array_delete (sta_index);
696 CPLCHECK_NUL ("Cannot fill sta_index or time in OI_T3");
697
698 /* Fill STA_INDEX and times for OI_FLUX */
699 for (int tel = 0; tel < ntel; tel++){
700 int sta0 = gravi_sta_index(tel+1, optical_train_table, oi_array);
701 for (int row=0; row < nrow; ++row) {
702 int idx = row * ntel + tel;
703 cpl_table_set_int (oi_flux, "STA_INDEX", idx, sta0);
704 cpl_table_set (oi_flux, "TIME", idx, cpl_array_get (times, row, &nv));
705 cpl_table_set (oi_flux, "MJD", idx, cpl_array_get (mjds, row, &nv));
706 }
707 } /* End loop on tel */
708
709 cpl_array_delete (mjds);
710 cpl_array_unwrap (times);
711 CPLCHECK_NUL ("Cannot fill sta_index or time in OI_FLUX");
712
713 /* Fill the exposure time */
714 double exptime = gravi_pfits_get_dit (gravi_data_get_header (preproc_data), type_data);
715 cpl_table_fill_column_window (oi_vis, "INT_TIME", 0, nrow * GRAVI_NBASE, exptime);
716 cpl_table_fill_column_window (oi_t3, "INT_TIME", 0, nrow * nclo, exptime);
717 cpl_table_fill_column_window (oi_flux, "INT_TIME", 0, nrow * ntel, exptime);
718 CPLCHECK_NUL ("Cannot fill exptime");
719
720 /* target_id is 1 unless type_data==SC and mode==dual_field
721 * Same definition applies when building the OI_TARGET */
722 int target_id = (!strcmp(mode, "gravi_dual") && (type_data==GRAVI_SC) )?2:1;
723 cpl_table_fill_column_window_int (oi_vis, "TARGET_ID", 0, nrow * GRAVI_NBASE, target_id);
724 cpl_table_fill_column_window_int (oi_t3, "TARGET_ID", 0, nrow * nclo, target_id);
725 cpl_table_fill_column_window_int (oi_flux, "TARGET_ID", 0, nrow * ntel, target_id);
726 CPLCHECK_NUL ("Cannot fill target_id");
727
728 cpl_msg_info (cpl_func, "Total time to fill STA_INDEX and TIME: %.4f s",
729 (double)(clock()-timer_start)/(double)CLOCKS_PER_SEC );
730
731 /*
732 *
733 * Extract visibility with p2vm
734 *
735 */
736
737 cpl_msg_info (cpl_func,"Apply p2vm to the data...");
738 timer_other += (clock() - timer_start);
739 timer_start = clock();
740
741 /* Get the pointers to the input data */
742 double ** pReg = cpl_malloc (n_region * sizeof(double*));
743 double ** pErr = cpl_malloc (n_region * sizeof(double*));
744 cpl_array *** pRegArr = cpl_malloc (n_region * sizeof(cpl_array**));
745 cpl_array *** pErrArr = cpl_malloc (n_region * sizeof(cpl_array**));
746 for (int reg = 0; reg < n_region; reg++) {
747 pRegArr[reg] = cpl_table_get_data_array (spectrum_table, GRAVI_DATA[all_region[reg]]);
748 pErrArr[reg] = cpl_table_get_data_array (spectrum_table, GRAVI_DATAERR[all_region[reg]]);
749 }
750 CPLCHECK_NUL ("Cannot get data");
751
752 /* Get the pointers to the output data */
753 cpl_array** tFlux = cpl_table_get_data_array (oi_flux, "FLUX");
754 cpl_array** tFluxErr = cpl_table_get_data_array (oi_flux, "FLUXERR");
755 cpl_array** tVis = cpl_table_get_data_array (oi_vis, "VISDATA");
756 cpl_array** tVisErr = cpl_table_get_data_array (oi_vis, "VISERR");
757 CPLCHECK_NUL ("Cannot get data");
758
759 /* Create column for chi2 */
760 gravi_table_new_column_array (oi_flux, "CHI2", NULL, CPL_TYPE_DOUBLE, nwave);
761 cpl_array** tChi2 = cpl_table_get_data_array (oi_flux, "CHI2");
762 int ndof = n_region - 16;
763
764 /* Temporary matrix output memory */
765 double* pOut = cpl_malloc (16 * nwave * sizeof(double));
766 double* pOutVar = cpl_malloc (16 * nwave * sizeof(double));
767 double* pChi2 = cpl_malloc (nwave * sizeof(double));
768
769 /* Quantities to test flux conservation */
770 double full_flux_reg = 0.0, full_flux_tel = 0.0;
771
772 cpl_msg_debug (cpl_func, "Matrix multiplication");
773
774 /* Loop on the frames */
775 for (cpl_size row = 0; row < nrow; row ++) {
776
777 /* Get pointers to the input data of this row */
778 for (int reg = 0; reg < n_region; reg++) {
779 pReg[reg] = cpl_array_get_data_double (pRegArr[reg][row]);
780 pErr[reg] = cpl_array_get_data_double (pErrArr[reg][row]);
781 }
782
783 /* Loop on wavelength */
784 for (cpl_size wave = 0 ; wave < nwave ; wave++ ) {
785
786 /* Integration the input flux of the row */
787 for (int reg = 0; reg < n_region; reg++) full_flux_reg += pReg[reg][wave];
788
789 /* Matrix multiplication, loop on outputs and regions
790 * We neglect the input correlation and don't compute
791 * the output correlations. FIXME: we need to compute
792 * the output covariance of {R,I} */
793 for (int out=0; out<16; out++) {
794 pOut[out*nwave+wave] = 0.0;
795 pOutVar[out*nwave+wave] = 0.0;
796 for (cpl_size reg = 0; reg < n_region; reg++) {
797 pOut[out*nwave+wave] += pReg[reg][wave] * pP2VM[wave][out*n_region+reg];
798 pOutVar[out*nwave+wave] += gravi_pow2 (pErr[reg][wave] * pP2VM[wave][out*n_region+reg]);
799 }
800 } /* End outputs and regions */
801
802 /* We compute the reduced chi2. We loop on regions
803 * to recompute the expected value, and accumulate the chi2 */
804 pChi2[wave] = 0.0;
805 for (cpl_size reg = 0; reg < n_region; reg++) {
806 double value = 0.0;
807 for (int out=0; out<16; out++)
808 value += pV2PM[wave][reg*16+out] * pOut[out*nwave+wave];
809 pChi2[wave] += gravi_pow2 ((value-pReg[reg][wave]) / pErr[reg][wave]) / ndof;
810 }
811
812 /* Integration the output flux of the row */
813 for (int tel = 0; tel < 4; tel++) full_flux_tel += pOut[tel*nwave+wave];
814
815 } /* End loop on wavelengths */
816
817 /* Set CHI2 (wrap is the fastest to create an array).
818 Same for all telescope since all-together */
819 for (int tel = 0; tel < ntel; tel++){
820 double * data = cpl_malloc (nwave * sizeof(double));
821 for (cpl_size wave = 0 ; wave < nwave ; wave++ )
822 data[wave] = pChi2[wave];
823 tChi2[row*ntel+tel] = cpl_array_wrap_double (data, nwave);
824 }
825
826 /* Set FLUX (wrap is the fastest to create an array) */
827 for (int tel = 0; tel < ntel; tel++){
828 double * data = cpl_malloc (nwave * sizeof(double));
829 for (cpl_size wave = 0 ; wave < nwave ; wave++ )
830 data[wave] = pOut[tel*nwave+wave];
831 tFlux[row*ntel+tel] = cpl_array_wrap_double (data, nwave);
832 }
833
834 /* Set FLUXERR */
835 for (int tel = 0; tel < ntel; tel++){
836 double * data = cpl_malloc (nwave * sizeof(double));
837 for (cpl_size wave = 0 ; wave < nwave ; wave++ )
838 data[wave] = sqrt (pOutVar[tel*nwave+wave]);
839 tFluxErr[row*ntel+tel] = cpl_array_wrap_double (data, nwave);
840 }
841
842 /* Set VISDATA */
843 for (int base = 0; base < GRAVI_NBASE; base++){
844 double complex * data = cpl_malloc (nwave * sizeof(double complex));
845 for (cpl_size wave = 0 ; wave < nwave ; wave++ )
846 data[wave] = (double complex)( pOut[(base+4)*nwave+wave] + 1.*I * pOut[(base+10)*nwave+wave] );
847 tVis[row*GRAVI_NBASE+base] = cpl_array_wrap_double_complex (data, nwave);
848 }
849
850 /* Set VISDATAERR */
851 for (int base = 0; base < GRAVI_NBASE; base++){
852 double complex * data = cpl_malloc (nwave * sizeof(double complex));
853 for (cpl_size wave = 0 ; wave < nwave ; wave++ )
854 data[wave] = (double complex)( sqrt(pOutVar[(base+4)*nwave+wave]) + 1.*I * sqrt(pOutVar[(base+10)*nwave+wave]) );
855 tVisErr[row*GRAVI_NBASE+base] = cpl_array_wrap_double_complex (data, nwave);
856 }
857
858 /* Free the PREPROC data to save memory
859 * (ex: max pointer 34065038 versus 48947400) */
860 for (cpl_size reg = 0; reg < n_region; reg++) {
861 FREE (cpl_array_delete, pRegArr[reg][row]);
862 FREE (cpl_array_delete, pErrArr[reg][row]);
863 }
864
865
866 }/* End loop on frames */
867
868 /* Deallocation of variables */
869 cpl_msg_debug (cpl_func, "Free pointers to data");
870 FREE (cpl_free, pReg);
871 FREE (cpl_free, pErr);
872 FREE (cpl_free, pRegArr);
873 FREE (cpl_free, pErrArr);
874 FREE (cpl_free, pOut);
875 FREE (cpl_free, pOutVar);
876 FREE (cpl_free, pChi2);
877 FREE (cpl_free, pP2VM);
878 FREE (cpl_free, pV2PM);
879 FREELOOP (cpl_matrix_delete, p2vm, nwave);
880 FREELOOP (cpl_matrix_delete, v2pm, nwave);
881
882 /* Check how "flux conservative" is the P2VM */
883 cpl_msg_info (cpl_func, "Total flux in TELs: %.2f [e], in REGIONs:%.2f [e] (ratio=%.5f)",
884 full_flux_tel,full_flux_reg,full_flux_tel/full_flux_reg);
885
886 sprintf (qc_name, "ESO QC TRANS P2VM %s",GRAVI_TYPE(type_data));
887 cpl_propertylist_update_double (p2vmred_header, qc_name, full_flux_tel/full_flux_reg);
888 cpl_propertylist_set_comment (p2vmred_header, qc_name, "[e/e] at P2VM extraction");
889
890 /* Count time */
891 cpl_msg_info (cpl_func, "Total time to apply matrix: %.4f s", (double)(clock()-timer_start)/(double)CLOCKS_PER_SEC );
892
893 if( type_data==GRAVI_FT )
894 timer_ft += (clock() - timer_start);
895 else
896 timer_sc += (clock() - timer_start);
897 timer_start = clock();
898 }
899 /* End loop on polarisation */
900 }
901 /* Loop on data type SC / FT */
902
903 /*
904 * Add QC for lambda met wavelength
905 */
906 if (gravi_data_has_extension(preproc_data, GRAVI_METROLOGY_EXT)){
907 double lambda_met = gravi_pfits_get_met_wavelength_mean(preproc_header, gravi_data_get_table(preproc_data, GRAVI_METROLOGY_EXT));
908 sprintf (qc_name, "ESO QC MET LAMBDA MEAN");
909 cpl_propertylist_update_double (p2vmred_header, qc_name, lambda_met);
910 cpl_propertylist_set_comment (p2vmred_header, qc_name, "[m] Effective mean metrology wavelength");
911 }
912
913
914 /*
915 * Print timing
916 */
917 timer_other += (clock() - timer_start);
918 cpl_msg_info (cpl_func, "Total time for FT: %10.4f s", (double)(timer_ft)/(double)CLOCKS_PER_SEC);
919 cpl_msg_info (cpl_func, "Total time for SC: %10.4f s", (double)(timer_sc)/(double)CLOCKS_PER_SEC);
920 cpl_msg_info (cpl_func, "Total time for OTHER: %7.4f s", (double)(timer_other)/(double)CLOCKS_PER_SEC);
921
922 /* Message and timer */
923 gravi_msg_function_exit(1);
924 return p2vmred_data;
925}
926
927/*----------------------------------------------------------------------------*/
939/*----------------------------------------------------------------------------*/
940
941cpl_error_code gravi_compute_opdc_state (gravi_data * p2vmred_data)
942{
943 /* Message and timer */
944 gravi_msg_function_start(1);
945 cpl_ensure_code (p2vmred_data, CPL_ERROR_NULL_INPUT);
946
947 int ntel = 4;
948 char qc_name[90];
949
950 /* Get necessary data */
951 cpl_propertylist * header = gravi_data_get_header (p2vmred_data);
952 int npol_ft = gravi_pfits_get_pola_num (header, GRAVI_FT);
953 cpl_table * oi_vis = gravi_data_get_oi_vis (p2vmred_data, GRAVI_FT, 0, npol_ft);
954 cpl_table * oi_flux = gravi_data_get_oi_flux (p2vmred_data, GRAVI_FT, 0, npol_ft);
955 cpl_table * opdc = gravi_data_get_table (p2vmred_data, GRAVI_OPDC_EXT);
956 cpl_size nrow_ft = cpl_table_get_nrow (oi_vis) / GRAVI_NBASE;
957 cpl_size nrow_opdc = cpl_table_get_nrow (opdc);
958 CPLCHECK_MSG ("Cannot get data");
959
960 /* Create the target phase for each baseline in the OI_VIS table */
961 gravi_table_new_column (oi_vis, "TARGET_PHASE", "rad", CPL_TYPE_DOUBLE);
962
963 /* Create the OPDC state for each telescope in the OI_FLUX table */
964 cpl_table_new_column (oi_flux, "STATE", CPL_TYPE_INT);
965 cpl_table_fill_column_window (oi_flux, "STATE", 0, nrow_ft * ntel, -1);
966
967 /* Create the OPDC state for each baseline in the OI_VIS table */
968 cpl_table_new_column (oi_vis, "STATE", CPL_TYPE_INT);
969 cpl_table_fill_column_window (oi_vis, "STATE", 0, nrow_ft * GRAVI_NBASE, -1);
970
971 /* Create the Global OPDC state in the OI_VIS table */
972 cpl_table_new_column (oi_vis, "OPDC_STATE", CPL_TYPE_INT);
973 cpl_table_fill_column_window (oi_vis, "OPDC_STATE", 0, nrow_ft * GRAVI_NBASE, -1);
974
975 if (nrow_opdc < nrow_ft)
976 cpl_msg_warning (cpl_func,"Missing FT or OPDC data: nrow_ft - nrow_opdc = %lli", nrow_ft-nrow_opdc);
977
978
979 /* Set the global OPDC state */
980 int * time_opdc = cpl_table_get_data_int (opdc, "TIME");
981 double * time_ft = cpl_table_get_data_double (oi_vis, "TIME");
982 int * global_state_opdc = cpl_table_get_data_int (opdc, "STATE");
983 int * global_state = cpl_table_get_data_int (oi_vis, "OPDC_STATE");
984 CPLCHECK_MSG ("Cannot get data");
985
986 /* Loop on FT rows */
987 for (cpl_size row_opdc=0, row_ft=0 ; row_ft<nrow_ft ; row_ft++) {
988
989 /* Check bounds or find the OPDC sample just following the current FT
990 * FIXME: We should use the closesd OPDC sample in the past, not future */
991 if ( (time_ft[row_ft*GRAVI_NBASE] < time_opdc[0]) || (time_ft[row_ft*GRAVI_NBASE] > time_opdc[nrow_opdc-1]) ) continue;
992 while ( time_ft[row_ft*GRAVI_NBASE] > time_opdc[row_opdc] ) row_opdc ++;
993
994 /* Set the global OPDC state */
995 for (int base = 0; base < GRAVI_NBASE; base++)
996 global_state[row_ft*GRAVI_NBASE+base] = global_state_opdc[row_opdc];
997 }
998
999 /* BASELINE_STATE was not in the original data of the instrument */
1000 if ( cpl_table_has_column (opdc,"BASELINE_STATE") ) {
1001
1002 int * steps_opdc = cpl_table_get_data_int (opdc, "STEPS");
1003 int * state_opdc = cpl_table_get_data_int (opdc, "BASELINE_STATE");
1004 int * base_state = cpl_table_get_data_int (oi_vis, "STATE");
1005 int * tel_state = cpl_table_get_data_int (oi_flux, "STATE");
1006 double * base_steps = cpl_table_get_data_double (oi_vis, "TARGET_PHASE");
1007 CPLCHECK_MSG ("Cannot get data");
1008
1009 /* Loop on FT rows */
1010 for (cpl_size row_opdc=0, row_ft=0 ; row_ft<nrow_ft ; row_ft++) {
1011
1012 /* Check bounds or find the OPDC sample just following the current FT
1013 * FIXME: We should use the closesd OPDC sample in the past, not future */
1014 if ( (time_ft[row_ft*GRAVI_NBASE] < time_opdc[0]) || (time_ft[row_ft*GRAVI_NBASE] > time_opdc[nrow_opdc-1]) ) continue;
1015 while ( time_ft[row_ft*GRAVI_NBASE] > time_opdc[row_opdc] ) row_opdc ++;
1016
1017 /* Get the flag */
1018 int state_flag = state_opdc[row_opdc];
1019
1020 /* Disentangle the state of each telescope:
1021 * The 4 beams are the first 4 bits */
1022 for (int tel = 0; tel < ntel; tel++)
1023 tel_state[row_ft*ntel+tel] = gravi_bit_get (state_flag, tel);
1024
1025 /* Disentangle the state of each baseline:
1026 * The 6 baselines are the bits from 5 to 10 */
1027 for (int base = 0; base < GRAVI_NBASE; base++)
1028 base_state[row_ft*GRAVI_NBASE+base] = gravi_bit_get (state_flag, ntel+base);
1029
1030 /* Bit 11 (>>10) is the Kalman state */
1031 // kalman = gravi_bit_get (state_opdc[row_opdc], ntel+GRAVI_NBASE);
1032
1033 /* Remaing bits are the off-load system */
1034 // offload = state_opdc[row_opdc] >> 11;
1035
1036 /* Use the closing triangles to recover the tracking on some baselines */
1037 for (cpl_size base = 0; base < GRAVI_NBASE; base++) {
1038 base_state[row_ft*GRAVI_NBASE+base] = CPL_MAX( base_state[row_ft*GRAVI_NBASE+base],
1039 base_state[row_ft*GRAVI_NBASE+GRAVI_TRI_BASE[base][0][0]] *
1040 base_state[row_ft*GRAVI_NBASE+GRAVI_TRI_BASE[base][0][1]] );
1041 base_state[row_ft*GRAVI_NBASE+base] = CPL_MAX( base_state[row_ft*GRAVI_NBASE+base],
1042 base_state[row_ft*GRAVI_NBASE+GRAVI_TRI_BASE[base][1][0]] *
1043 base_state[row_ft*GRAVI_NBASE+GRAVI_TRI_BASE[base][1][1]] );
1044 }
1045
1046 /* Disentangle the piezo steps, in units of [pi/8 rad]
1047 * Each beam is coded over 4 bits (16 values over the circle)
1048 * *BUT* the last telescope is coded in wrong memory place !
1049 * beam0 = 0-3, beam1 = 4-7, beam2 = 8-11, beam3 = 16-19 */
1050 int tmp_arr[4];
1051 for (int tel = 0; tel < ntel; tel++) {
1052 int pos = (tel<3) ? 4*tel : 4*(tel+1);
1053 tmp_arr[tel] = 15 & ((steps_opdc[row_opdc]) >> pos);
1054 }
1055
1056 /* Compute the FT target phase of each baseline, in [rad] */
1057 for (int base=0; base<GRAVI_NBASE; base++) {
1058 base_steps[row_ft*GRAVI_NBASE+base] = (tmp_arr[GRAVI_BASE_TEL[base][1]] - tmp_arr[GRAVI_BASE_TEL[base][0]]) * CPL_MATH_PI / 8.0;
1059 }
1060
1061 } /* End loop on FT rows */
1062
1063 } else {
1064 cpl_msg_warning (cpl_func,"No column BASELINE_STATE in OPDC... old data ?");
1065 cpl_msg_warning (cpl_func,"The STATE flags are set to 1 (valid) although the information is unavailable");
1066 cpl_table_fill_column_window (oi_flux, "STATE", 0, nrow_ft * ntel, 1);
1067 cpl_table_fill_column_window (oi_vis, "STATE", 0, nrow_ft * GRAVI_NBASE, 1);
1068 }
1069
1070 /* Duplicate in the second polarisation if any */
1071 if ( npol_ft>1 ) {
1072 cpl_table * tmp;
1073 cpl_msg_debug (cpl_func, "Duplicate the FT tracking state for 2sd polarisation");
1074
1075 tmp = gravi_data_get_oi_vis (p2vmred_data, GRAVI_FT, 1, npol_ft);
1076 cpl_table_duplicate_column (tmp, "TARGET_PHASE", oi_vis, "TARGET_PHASE");
1077 cpl_table_duplicate_column (tmp, "STATE", oi_vis, "STATE");
1078 cpl_table_duplicate_column (tmp, "OPDC_STATE", oi_vis, "OPDC_STATE");
1079
1080 tmp = gravi_data_get_oi_flux (p2vmred_data, GRAVI_FT, 1, npol_ft);
1081 cpl_table_duplicate_column (tmp, "STATE", oi_flux, "STATE");
1082 }
1083
1084
1085 /*
1086 * QC: Compute the phase RMS per base when in tracking state
1087 */
1088
1089 cpl_array **visdata = cpl_table_get_data_array (oi_vis,"VISDATA");
1090 int * state = cpl_table_get_data_int (oi_vis,"STATE");
1091 double * target_phase = cpl_table_get_data_double (oi_vis,"TARGET_PHASE");
1092 CPLCHECK_MSG ("Cannot get data");
1093
1094 double complex * tmp_cpx = cpl_malloc (sizeof(double complex)*nrow_ft);
1095
1096 /* Loop on base */
1097 for (int base = 0; base < GRAVI_NBASE; base ++) {
1098
1099 /* Compute real-time phase and its mean (as phasors) */
1100 double complex mean_cpx = 0.0 + I * 0.0;
1101 for (cpl_size row_ft=0 ; row_ft<nrow_ft ; row_ft++) {
1102 if (state[row_ft*GRAVI_NBASE+base] < 1) continue;
1103 tmp_cpx[row_ft] = cpl_array_get_mean_complex (visdata[row_ft*GRAVI_NBASE+base]) * cexp (-I*target_phase[row_ft*GRAVI_NBASE+base]);
1104 mean_cpx += tmp_cpx[row_ft];
1105 }
1106
1107 /* Compute <phi**2> and the number of valid point */
1108 double sum = 0.0000001, sum2 = 0.0;
1109 for (cpl_size row_ft=0 ; row_ft<nrow_ft ; row_ft++) {
1110 if (state[row_ft*GRAVI_NBASE+base] < 1) continue;
1111 sum2 += pow (carg (tmp_cpx[row_ft] * conj(mean_cpx)), 2);
1112 sum += 1.0;
1113 }
1114
1115 /* Write the QC parameter for this base */
1116 sprintf (qc_name, "ESO QC PHASE_FT%d%d RMS",GRAVI_BASE_TEL[base][0]+1, GRAVI_BASE_TEL[base][1]+1);
1117 cpl_propertylist_update_double (header, qc_name, sqrt (sum2 / sum));
1118 cpl_propertylist_set_comment (header, qc_name, "[rad] residuals when tracking");
1119
1120 /* Write the QC parameter for this base */
1121 sprintf (qc_name, "ESO QC TRACKING_RATIO_FT%d%d",GRAVI_BASE_TEL[base][0]+1, GRAVI_BASE_TEL[base][1]+1);
1122 cpl_propertylist_update_int (header, qc_name, sum*100.0/nrow_ft);
1123 cpl_propertylist_set_comment (header, qc_name, "[%] ratio of time with tracking");
1124 }
1125 /* End loop on base */
1126
1127 FREE (cpl_free,tmp_cpx);
1128
1129 /*
1130 * QC: Compute fraction of time traking
1131 */
1132
1133 int* tab_state = cpl_table_get_data_int (gravi_data_get_table (p2vmred_data, GRAVI_OPDC_EXT), "STATE");
1134 CPLCHECK_MSG ("Cannot get data");
1135
1136 long tracking = 0;
1137 for (cpl_size row_opdc=0; row_opdc < nrow_opdc; row_opdc++) {
1138 if (tab_state[row_opdc] == 3 || tab_state[row_opdc] == 2) tracking ++;
1139 }
1140
1141 sprintf (qc_name, "ESO QC TRACKING_RATIO");
1142 cpl_propertylist_update_int (header, qc_name, tracking*100/nrow_ft);
1143 cpl_propertylist_set_comment (header, qc_name, "[%] ratio of time with full FT tracking");
1144
1145 CPLCHECK_MSG ("Cannot put QC parameters");
1146
1147 /* Message and timer */
1148 gravi_msg_function_exit(1);
1149 return CPL_ERROR_NONE;
1150}
1151
1152
1153/*----------------------------------------------------------------------------*/
1163/*----------------------------------------------------------------------------*/
1164
1165cpl_error_code gravi_compute_tau0 (gravi_data * data)
1166{
1167 /* verbose */
1168 gravi_msg_function_start(1);
1169 cpl_ensure_code (data, CPL_ERROR_NULL_INPUT);
1170
1171 double gain = 16.8; // [rad/V]
1172 char qc_name[90];
1173
1174 /* Get the OPDC table */
1175 cpl_propertylist * header = gravi_data_get_header (data);
1176 cpl_table * opdc = gravi_data_get_table (data, GRAVI_OPDC_EXT);
1177 cpl_size nrow = cpl_table_get_nrow (opdc);
1178
1179 CPLCHECK_MSG("Cannot load the OPDC table");
1180
1181 /* Compute the delays to explore */
1182 cpl_size max_delay = 200;
1183
1184 /* Time in microseconds [us] and offset in [V] */
1185 int *time = cpl_table_get_data_int (opdc, "TIME");
1186 CPLCHECK_MSG("Cannot load the TIME column");
1187
1188 float ** piezo = gravi_table_get_data_array_float (opdc, "PIEZO_DL_OFFSET");
1189 CPLCHECK_MSG("Cannot load the PIEZO_OFFSET columns");
1190
1191 /* Loop on base */
1192 for (int base=0; base<GRAVI_NBASE; base++) {
1193 int t1 = GRAVI_BASE_TEL[base][0];
1194 int t2 = GRAVI_BASE_TEL[base][1];
1195
1196 /* Init the QC parameters */
1197 sprintf (qc_name, "ESO QC TAU0 OPDC%d%d", t1+1,t2+1);
1198 cpl_propertylist_update_double (header, qc_name, GRAVI_NAN_DOUBLE);
1199 cpl_propertylist_set_comment (header, qc_name, "[s] tau0 for variance of 1 rad2");
1200
1201
1202 /* Loop on delays */
1203 for (cpl_size delay=0; delay < max_delay; delay++) {
1204
1205 /* srtf(d) = < (opd(t) - opd(t+d))^2> in [rad^2] */
1206 double strf = 0.0;
1207 for (cpl_size row=0; row < nrow - max_delay - 1 ; row++) {
1208 float diff = (piezo[row][t1]-piezo[row][t2]) - (piezo[row+delay][t1]-piezo[row+delay][t2]);
1209 strf += (double)(diff * diff);
1210 }
1211 strf = strf * gain * gain / (nrow - max_delay - 1);
1212
1213 /* If the variance is larger than 1rad2,
1214 * we found the tau0 */
1215 if ( strf > 1.0 ) {
1216 double tau0 = (time[delay+1] - time[1]) * 1e-6;
1217 cpl_msg_info (cpl_func,"Compute %s = %.2f [ms]", qc_name, 1e3 * tau0);
1218 cpl_propertylist_update_double (header, qc_name, tau0);
1219 break;
1220 }
1221
1222 } /* End loop on delay */
1223
1224 } /* End loop on base */
1225
1226 /* Free the list of arrays */
1227 FREE (cpl_free, piezo);
1228
1229 gravi_msg_function_exit(1);
1230 return CPL_ERROR_NONE;
1231}
1232
1233/*----------------------------------------------------------------------------*/
1243/*----------------------------------------------------------------------------*/
1244
1245cpl_error_code gravi_compute_qc_injection (gravi_data * data)
1246{
1247 gravi_msg_function_start(1);
1248 cpl_ensure_code (data, CPL_ERROR_NULL_INPUT);
1249
1250 char qc_name[100];
1251 double p05 = 0.0, p95 = 0.0;
1252
1253 cpl_propertylist * header = gravi_data_get_header (data);
1254 int npol = gravi_pfits_get_pola_num (header, GRAVI_FT);
1255 int ntel = 4, hw = 2;
1256
1257 /* Create the kernel to smooth the flux series */
1258 cpl_msg_info (cpl_func, "Smooth flux sery over %i samples", hw);
1259
1260 /* How many rows do we have per telescope */
1261 cpl_table * table = gravi_data_get_oi_flux (data, GRAVI_FT, 0, npol);
1262 cpl_size nrow = cpl_table_get_nrow (table) / ntel;
1263
1264
1265 /* For each telescope */
1266 for (int tel = 0; tel<ntel; tel++) {
1267
1268 /* Build a vector of total flux */
1269 cpl_vector * raw_flux = cpl_vector_new (nrow);
1270 cpl_vector_fill (raw_flux, 0.0);
1271
1272 for (int pol = 0; pol<npol; pol++) {
1273
1274 /* Get table */
1275 cpl_table * table = gravi_data_get_oi_flux (data, GRAVI_FT, pol, npol);
1276 const cpl_array ** flux_array = cpl_table_get_data_array_const (table, "FLUX");
1277
1278 for (int row = 0; row<nrow; row++)
1279 cpl_vector_set (raw_flux, row, cpl_vector_get (raw_flux, row) +
1280 cpl_array_get_mean (flux_array[tel+row*ntel]));
1281 }
1282
1283 /* Convolve to filter noise, delete raw flux */
1284 cpl_vector * flux;
1285 flux = cpl_vector_filter_lowpass_create (raw_flux, CPL_LOWPASS_GAUSSIAN, hw);
1286 cpl_vector_delete (raw_flux);
1287
1288 /* Define positive */
1289 for (int row = 0; row<nrow; row++)
1290 cpl_vector_set (flux, row, CPL_MAX (0.0, cpl_vector_get (flux, row)));
1291
1292 /* Sort the flux array and normalise */
1293 cpl_vector_sort (flux, CPL_SORT_ASCENDING);
1294 /* Do not divide if flux equal to 0*/
1295 if (cpl_vector_get_max (flux) > 0)
1296 cpl_vector_divide_scalar (flux, cpl_vector_get_max (flux) / 100);
1297
1298 /* Histogram as a string, values are 0-100 */
1299 char qc_value[100];
1300 sprintf (qc_value,"%.0f,%.0f,%.0f,%.0f,%.0f,%.0f,%.0f,%.0f,%.0f",
1301 cpl_vector_get (flux, (cpl_size)(0.1*(nrow-1))),
1302 cpl_vector_get (flux, (cpl_size)(0.2*(nrow-1))),
1303 cpl_vector_get (flux, (cpl_size)(0.3*(nrow-1))),
1304 cpl_vector_get (flux, (cpl_size)(0.4*(nrow-1))),
1305 cpl_vector_get (flux, (cpl_size)(0.5*(nrow-1))),
1306 cpl_vector_get (flux, (cpl_size)(0.6*(nrow-1))),
1307 cpl_vector_get (flux, (cpl_size)(0.7*(nrow-1))),
1308 cpl_vector_get (flux, (cpl_size)(0.8*(nrow-1))),
1309 cpl_vector_get (flux, (cpl_size)(0.9*(nrow-1))));
1310
1311 /* Create the QC entry in the FITS header */
1312 cpl_msg_info (cpl_func, "Flux Histo: %s", qc_value);
1313 sprintf (qc_name, "ESO QC FLUX_FT%d HISTO", tel+1);
1314 cpl_propertylist_update_string (header, qc_name, qc_value);
1315 cpl_propertylist_set_comment (header, qc_name, "decile in percent");
1316
1317 /* Compute the 5 percentile and 95 percentile */
1318 p05 = cpl_vector_get (flux, (cpl_size)(0.05*(nrow-1)));
1319 p95 = cpl_vector_get (flux, (cpl_size)(0.95*(nrow-1)));
1320 /* remove error if p95 equal to 0*/
1321 if (p95<1) p95=1;
1322
1323 /* Create the QC entry in the FITS header */
1324 sprintf (qc_name, "ESO QC FLUX_FT%d P05P95", tel+1);
1325 cpl_propertylist_update_double (header, qc_name, p05/p95);
1326 cpl_propertylist_set_comment (header, qc_name, "injected flux 5 percentile over 95 percentile");
1327
1328 /* Deleve vector */
1329 cpl_vector_delete (flux);
1330 }
1331
1332 gravi_msg_function_exit(1);
1333 return CPL_ERROR_NONE;
1334}
1335
1336
1337/*----------------------------------------------------------------------------*/
1347/*----------------------------------------------------------------------------*/
1348
1349cpl_error_code gravi_compute_qc_ft_opd_estimator (gravi_data * p2vmred_data)
1350{
1351 gravi_msg_function_start(1);
1352 cpl_ensure_code (p2vmred_data, CPL_ERROR_NULL_INPUT);
1353
1354 char qc_name[100];
1355 int nv;
1356
1357 /* Get necessary data */
1358 cpl_propertylist * header = gravi_data_get_header (p2vmred_data);
1359 int npol = gravi_pfits_get_pola_num (header, GRAVI_FT);
1360 cpl_table * opdc = gravi_data_get_table (p2vmred_data, GRAVI_OPDC_EXT);
1361 int * time_opdc = cpl_table_get_data_int (opdc, "TIME");
1362 if (cpl_table_has_column(opdc, "OPD") == 1){
1363 cpl_array **opd_data_opdc = cpl_table_get_data_array (opdc,"OPD");
1364 cpl_size nrow = cpl_table_get_nrow (opdc);
1365 CPLCHECK_MSG ("Cannot get data");
1366
1367 /* initialize arrays */
1368 cpl_array * opd_res_array = cpl_array_new(nrow, CPL_TYPE_DOUBLE_COMPLEX);
1369 cpl_array * opd_array = cpl_array_new(nrow, CPL_TYPE_DOUBLE);
1370 cpl_array * snr_array = cpl_array_new(nrow, CPL_TYPE_DOUBLE);
1371
1372 /* do it for both polarisations (if needed) */
1373 for (int pol = 0; pol<npol; pol++) {
1374 cpl_table * oi_vis = gravi_data_get_oi_vis (p2vmred_data, GRAVI_FT, pol, npol);
1375 double * time_oivis = cpl_table_get_data_double (oi_vis, "TIME");
1376 cpl_array **visdata = cpl_table_get_data_array (oi_vis,"VISDATA");
1377 cpl_array **visdata_err = cpl_table_get_data_array (oi_vis,"VISERR");
1378
1379 cpl_size nrow_oivis = cpl_table_get_nrow (oi_vis);
1380 cpl_size nwave= cpl_array_get_size (visdata[0]);
1381
1382
1383 /* For each baseline */
1384 for (int base = 0; base<GRAVI_NBASE; base++) {
1385 cpl_size row_oivis = base;
1386 /* For each row */
1387 for (cpl_size row = 0; row<nrow; row++)
1388 {
1389
1390 /* FIXME: get rid of the 0.0011/2 (half step) */
1391 while (( fabs (time_oivis[row_oivis+GRAVI_NBASE]-time_opdc[row]+0.0011/2*1e6 )< fabs (time_oivis[row_oivis]-time_opdc[row]+0.0011/2*1e6 ) ) && (row_oivis<nrow_oivis-2*GRAVI_NBASE))
1392 {
1393 row_oivis+=GRAVI_NBASE;
1394 }
1395
1396 double complex visdata_mean=cpl_array_get_mean_complex (visdata[row_oivis]);
1397 cpl_array_set_double_complex(opd_res_array,row, visdata_mean*cexp ( I * cpl_array_get_float (opd_data_opdc[row], base, &nv) ));
1398 cpl_array_set_double(opd_array, row, carg(visdata_mean));
1399
1400
1401 /* get SNR */
1402 double signal =0.0;
1403 double noise =0.0;
1404 for (cpl_size wave = 0 ; wave < nwave ; wave ++)
1405 {
1406 signal+=cabs(cpl_array_get_double_complex (visdata[row_oivis],wave, &nv));
1407 double err=cabs(cpl_array_get_double_complex (visdata_err[row_oivis],wave, &nv));
1408 noise+=err*err;
1409 }
1410 noise=sqrt(noise+1);
1411 if (signal > noise)
1412 cpl_array_set_double(snr_array,row, signal /noise -1 );
1413 else
1414 cpl_array_set_double(snr_array,row, 0 );
1415
1416 }
1417 CPLCHECK_MSG ("Cannot compute new array for cross-referencing OPD and FT OPD");
1418
1419 double complex opd_mean = cpl_array_get_mean_complex (opd_res_array);
1420 cpl_array_multiply_scalar_complex(opd_res_array,conj(opd_mean));
1421
1422 double c_o=0,s_o=0;
1423 double a_num=0,b_num=0,c_num=0,d_num=0,a_den=1,b_den=1,c_den=1,d_den=1;
1424 double opd_std_num=0,opd_std_den=1;
1425 double opd, opd_ref, snr_ref;
1426
1427 for (cpl_size row = 0; row<nrow; row++)
1428 {
1429 opd=cpl_array_get_double(opd_array,row,&nv);
1430 opd_ref=carg(cpl_array_get_double_complex(opd_res_array,row,&nv));
1431 snr_ref=cpl_array_get_double(snr_array,row,&nv);
1432 c_o=cos(opd);
1433 a_num+=snr_ref*opd_ref*c_o;
1434 a_den+=snr_ref*c_o*c_o;
1435 s_o=sin(opd);
1436 b_num+=snr_ref*opd_ref*s_o;
1437 b_den+=snr_ref*s_o*s_o;
1438 c_o=c_o*c_o;
1439 c_num+=snr_ref*opd_ref*c_o;
1440 c_den+=snr_ref*c_o*c_o;
1441 s_o=s_o*s_o;
1442 d_num+=snr_ref*opd_ref*s_o;
1443 d_den+=snr_ref*s_o*s_o;
1444
1445 /* calculate the standard deviation: numerator and denominator*/
1446 opd_std_num+=snr_ref*opd_ref*snr_ref*opd_ref;
1447 opd_std_den+=snr_ref;
1448
1449 }
1450
1451 CPLCHECK_MSG ("Cannot calculate variance of linearity");
1452
1453 double a=a_num/a_den;
1454 double b=b_num/b_den;
1455 double c=c_num/c_den;
1456 double d=d_num/d_den;
1457 /* calculate standard deviation */
1458 double opd_std=sqrt(0.1+opd_std_num)/opd_std_den;
1459
1460 /* calculate the average standard error of cos, sin, sin2 and cos2 */
1461 double res_car=a*a+b*b+c*c+d*d-16*opd_std*opd_std;
1462 if (res_car > 0) res_car=sqrt(res_car);
1463 else res_car=0;
1464
1465 /* Create the QC entry in the FITS header*/
1466 if (npol == 2) sprintf (qc_name, "ESO QC LIN_FT P%d_B%d", pol,base+1);
1467 if (npol == 1) sprintf (qc_name, "ESO QC LIN_FT P%d_B%d", 3,base+1);
1468 cpl_propertylist_update_double (header, qc_name, res_car);
1469 cpl_propertylist_set_comment (header, qc_name, "FT nonlinearity biases [rad]");
1470
1471 CPLCHECK_MSG ("Cannot store QCs");
1472
1473 }
1474 }
1475
1476 cpl_array_delete (opd_array);
1477 cpl_array_delete (opd_res_array);
1478 cpl_array_delete (snr_array);
1479 CPLCHECK_MSG ("Cannot close arrays");
1480 }
1481 else{
1482 cpl_msg_info(cpl_func, "Cannot found the OPD column (probably old data). Skip computation of QC LIN_FT");
1483 }
1484
1485 gravi_msg_function_exit(1);
1486 return CPL_ERROR_NONE;
1487}
1488
1489
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_oi_flux
#define gravi_data_get_oi_flux(data, type, pol, npol)
Definition gravi_data.h:49
gravi_data_get_header
#define gravi_data_get_header(data)
Definition gravi_data.h:75
gravi_data_get_p2vm_data
#define gravi_data_get_p2vm_data(data, type)
Definition gravi_data.h:57
gravi_data_get_imaging_detector
#define gravi_data_get_imaging_detector(data, type)
Definition gravi_data.h:60
gravi_data_get_oi_vis
#define gravi_data_get_oi_vis(data, type, pol, npol)
Definition gravi_data.h:46
ntel
const cpl_size ntel
Definition gravi_demodulate.c:84
gravi_dfs.h
cpl_msg_debug
cpl_msg_debug(cpl_func, "Spectra has <50 pixels -> don't flat")
header
cpl_propertylist * header
Definition gravi_old.c:2004
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"))
gravi_p2vmred.h
gravi_pfits.h
GRAVI_INSNAME
#define GRAVI_INSNAME(type, pol, npol)
Definition gravi_pfits.h:198
GRAVI_OI_TARGET_EXT
#define GRAVI_OI_TARGET_EXT
Definition gravi_pfits.h:86
GRAVI_OI_ARRAY_EXT
#define GRAVI_OI_ARRAY_EXT
Definition gravi_pfits.h:83
GRAVI_OPDC_EXT
#define GRAVI_OPDC_EXT
Definition gravi_pfits.h:62
GRAVI_NAN_DOUBLE
#define GRAVI_NAN_DOUBLE
Definition gravi_pfits.h:97
GRAVI_SC
#define GRAVI_SC
Definition gravi_pfits.h:165
GRAVI_OI_VIS_EXT
#define GRAVI_OI_VIS_EXT
Definition gravi_pfits.h:92
GRAVI_OI_FLUX_EXT
#define GRAVI_OI_FLUX_EXT
Definition gravi_pfits.h:93
gravi_pfits_get_dit
#define gravi_pfits_get_dit(plist, type)
Definition gravi_pfits.h:213
GRAVI_EXTVER
#define GRAVI_EXTVER(type, pol, npol)
Definition gravi_pfits.h:204
gravi_pfits_get_sobj_decep
#define gravi_pfits_get_sobj_decep(plist)
Definition gravi_pfits.h:206
GRAVI_OI_T3_EXT
#define GRAVI_OI_T3_EXT
Definition gravi_pfits.h:94
GRAVI_TYPE
#define GRAVI_TYPE(type)
Definition gravi_pfits.h:167
gravi_pfits_get_sobj_raep
#define gravi_pfits_get_sobj_raep(plist)
Definition gravi_pfits.h:207
GRAVI_SPECTRUM_DATA_EXT
#define GRAVI_SPECTRUM_DATA_EXT(type)
Definition gravi_pfits.h:55
GRAVI_ARRAY_GEOMETRY_EXT
#define GRAVI_ARRAY_GEOMETRY_EXT
Definition gravi_pfits.h:84
GRAVI_METROLOGY_EXT
#define GRAVI_METROLOGY_EXT
Definition gravi_pfits.h:60
GRAVI_FT
#define GRAVI_FT
Definition gravi_pfits.h:166
GRAVI_OI_WAVELENGTH_EXT
#define GRAVI_OI_WAVELENGTH_EXT
Definition gravi_pfits.h:91
GRAVI_OPTICAL_TRAIN_EXT
#define GRAVI_OPTICAL_TRAIN_EXT
Definition gravi_pfits.h:85
gravi_utils.h
TRANSMISSION
#define TRANSMISSION
Definition gravi_utils.h:150
gravi_bit_get
#define gravi_bit_get(number, pos)
Definition gravi_utils.h:38
GRAVI_LABINPUT_4
#define GRAVI_LABINPUT_4
Definition gravi_utils.h:127
COHERENCE
#define COHERENCE
Definition gravi_utils.h:149
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
GRAVI_LABINPUT_3
#define GRAVI_LABINPUT_3
Definition gravi_utils.h:126
CPLCHECK_NUL
#define CPLCHECK_NUL(msg)
Definition gravi_utils.h:48
GRAVI_LABINPUT_1
#define GRAVI_LABINPUT_1
Definition gravi_utils.h:124
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
PHASE
#define PHASE
Definition gravi_utils.h:151
GRAVI_LABINPUT_2
#define GRAVI_LABINPUT_2
Definition gravi_utils.h:125
gravi_pow2
#define gravi_pow2(data)
Definition gravi_utils.h:81
FREELOOP
#define FREELOOP(function, variable, n)
Definition gravi_utils.h:72
GRAVI_NBASE
#define GRAVI_NBASE
Definition gravi_utils.h:105
gravi_vis.h
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_matrix_invertSV_create
cpl_matrix * gravi_matrix_invertSV_create(cpl_matrix *a_in)
Invers a matrix with singular value decomposition.
Definition gravi_cpl.c:2920
gravi_array_init_double
cpl_array * gravi_array_init_double(long n, double value)
Definition gravi_cpl.c:596
gravi_table_set_string_fixlen
cpl_error_code gravi_table_set_string_fixlen(cpl_table *table, const char *name, int row, const char *value, int len)
Set string in table, ensuring fixed length (right space padding). see cpl_table_set_string.
Definition gravi_cpl.c:768
gravi_table_get_data_array_float
float ** gravi_table_get_data_array_float(cpl_table *table, const char *name)
Definition gravi_cpl.c:496
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_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_new
gravi_data * gravi_data_new(int nb_ext)
Create an empty gravi_data.
Definition gravi_data.c:110
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_has_extension
int gravi_data_has_extension(gravi_data *raw_calib, const char *ext_name)
Check if data has extension with given EXTNAME.
Definition gravi_data.c:1808
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_compute_qc_injection
cpl_error_code gravi_compute_qc_injection(gravi_data *data)
Compute the QC for the injection stability.
Definition gravi_p2vmred.c:1245
gravi_compute_tau0
cpl_error_code gravi_compute_tau0(gravi_data *data)
Compute the QC TAU0 parameter.
Definition gravi_p2vmred.c:1165
gravi_compute_opdc_state
cpl_error_code gravi_compute_opdc_state(gravi_data *p2vmred_data)
Compute the real-time tracking state from OPDC.
Definition gravi_p2vmred.c:941
gravi_compute_p2vmred
gravi_data * gravi_compute_p2vmred(gravi_data *preproc_data, gravi_data *p2vm_map, const char *mode, const cpl_parameterlist *parlist)
Converts preprocessed data into coherent fluxes using the P2VM.
Definition gravi_p2vmred.c:330
gravi_create_oitarget_table
cpl_table * gravi_create_oitarget_table(const cpl_propertylist *header, const char *mode)
Create the OI_TARGET table from the main header.
Definition gravi_p2vmred.c:90
gravi_create_oiarray_table
cpl_table * gravi_create_oiarray_table(const cpl_table *array_geometry, int is_cal)
Create the OI_ARRAY table from the ARRAY_GEOMETRY.
Definition gravi_p2vmred.c:238
gravi_compute_qc_ft_opd_estimator
cpl_error_code gravi_compute_qc_ft_opd_estimator(gravi_data *p2vmred_data)
Compute the QC for the FT linearity.
Definition gravi_p2vmred.c:1349
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_sobj
const char * gravi_pfits_get_sobj(const cpl_propertylist *plist)
Definition gravi_pfits.c:441
gravi_convert_to_mjd
double gravi_convert_to_mjd(const char *start)
Definition gravi_pfits.c:1188
gravi_pfits_get_pmra
double gravi_pfits_get_pmra(const cpl_propertylist *plist)
Definition gravi_pfits.c:119
gravi_pfits_ensure_double
cpl_error_code gravi_pfits_ensure_double(cpl_propertylist *self, const char *name)
Definition gravi_pfits.c:1518
gravi_pfits_get_plx
double gravi_pfits_get_plx(const cpl_propertylist *plist)
Definition gravi_pfits.c:133
gravi_pfits_is_calib
int gravi_pfits_is_calib(const cpl_propertylist *plist)
Definition gravi_pfits.c:1452
gravi_pfits_get_robj_raep
double gravi_pfits_get_robj_raep(const cpl_propertylist *plist)
Definition gravi_pfits.c:1331
gravi_pfits_get_robj
const char * gravi_pfits_get_robj(const cpl_propertylist *plist)
Definition gravi_pfits.c:499
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_start_prcacq
const char * gravi_pfits_get_start_prcacq(const cpl_propertylist *plist)
Definition gravi_pfits.c:647
gravi_pfits_get_pmdec
double gravi_pfits_get_pmdec(const cpl_propertylist *plist)
Definition gravi_pfits.c:126
gravi_pfits_get_double_silentdefault
double gravi_pfits_get_double_silentdefault(const cpl_propertylist *plist, const char *name, double def)
Definition gravi_pfits.c:1634
gravi_pfits_get_robj_decep
double gravi_pfits_get_robj_decep(const cpl_propertylist *plist)
Definition gravi_pfits.c:1383
gravi_dump_the_boss
cpl_error_code gravi_dump_the_boss(double ra, double dec)
Definition gravi_utils.c:1217
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_CLO_TEL
int GRAVI_CLO_TEL[4][3]
Definition gravi_utils.c:64
gravi_spectrum_get_nwave
cpl_size gravi_spectrum_get_nwave(const cpl_table *table)
Definition gravi_utils.c:1013
GRAVI_DATAERR
char GRAVI_DATAERR[50][10]
Definition gravi_utils.c:79
GRAVI_DATA
char GRAVI_DATA[50][7]
Definition gravi_utils.c:71
gravi_sta_index
short gravi_sta_index(int gravi_input, cpl_table *optical_train_table, cpl_table *array_geometry_table)
Retrieve STA_INDEX corresponding to a given input.
Definition gravi_utils.c:828
GRAVI_BASE_TEL
int GRAVI_BASE_TEL[GRAVI_NBASE][2]
Definition gravi_utils.c:56
gravi_table_oi_create
cpl_table * gravi_table_oi_create(int nwave, int nrow, const char *oi_name)
Create the oi table (oi_vis, oi_vis2, oi_t3)
Definition gravi_utils.c:153
GRAVI_TRI_BASE
int GRAVI_TRI_BASE[GRAVI_NBASE][2][2]
Definition gravi_utils.c:60
Generated on Tue Apr 28 2026 07:47:55 for GRAVI Pipeline Reference Manual by doxygen 1.13.2