GRAVI Pipeline Reference Manual 1.10.1
Loading...
Searching...
No Matches
gravi_metrology.c
Go to the documentation of this file.
1/* $Id: gravi_vis.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
33
34/*
35 *
36 * History :
37 * 08/10/2019 two places : correct bug for AT values : fc_focus_at
38 * 12/11/2018 add gravi_data *static_param_data to gravi_metrology_telfc
39 * add gravi_data *static_param_data to gravi_metrology_get_fc_focus
40 * 10/01/2019 remove int old_flag_telescope[4][4][2]; int old_flag_fiber_coupler[4][2];
41 */
42
43/*-----------------------------------------------------------------------------
44 Includes
45 -----------------------------------------------------------------------------*/
46
47#ifdef HAVE_CONFIG_H
48#include <config.h>
49#endif
50
51#include <cpl.h>
52#include <string.h>
53#include <stdio.h>
54#include <math.h>
55#include <time.h>
56
57#include "gravi_data.h"
58#include "gravi_pfits.h"
59#include "gravi_dfs.h"
60
61#include "gravi_cpl.h"
62#include "gravi_utils.h"
63#include "gravi_ellipse.h"
64#include "gravi_signal.h"
65#include "gravi_eop.h"
66
67#include "gravi_metrology.h"
68
69/*-----------------------------------------------------------------------------
70 Private prototypes
71 -----------------------------------------------------------------------------*/
72
73
74cpl_error_code gravi_metrology_tac (cpl_table * metrology_table,
75 cpl_table * vismet_table,
76 cpl_propertylist * header);
77
78cpl_error_code gravi_metrology_telfc (cpl_table * metrology_table,
79 cpl_table * vismet_table,
80 gravi_data *static_param_data,
81 cpl_propertylist * header,
82 const cpl_parameterlist * parlist);
83
84cpl_error_code gravi_metrology_acq (cpl_table * visacq_table,
85 cpl_table * vismet_table,
86 double delay,
87 cpl_propertylist * header);
88
89cpl_error_code gravi_metrology_update_receiverpos (cpl_propertylist * header,
90 cpl_table *receiver_table);
91
92double gravi_metrology_get_posx (cpl_propertylist * header,
93 int tel, int diode);
94
95double gravi_metrology_get_posy (cpl_propertylist * header,
96 int tel, int diode);
97
98double gravi_metrology_get_fc_focus (cpl_propertylist * header, int gv, gravi_data *static_param_data);
99
100double gravi_metrology_get_fc_shift (cpl_propertylist * header, int gv, gravi_data *static_param_data);
101
102cpl_error_code gravi_metrology_get_astig (cpl_propertylist * header, int gv,
103 double * ampltiude, double * angle,
104 double * radius, gravi_data *static_param_data);
105
106long gravi_round (double number);
107long gravi_round (double number)
108{
109 return (number >= 0) ? (long)(number + 0.5) : (long)(number - 0.5);
110}
111
112/*-----------------------------------------------------------------------------
113 TAC definitions
114 -----------------------------------------------------------------------------*/
115
116#define DEBUG_LEVEL 2
117
118#define PI 3.14159265359
119#define TWOPI 6.28318530718
120
121#define MAX_DATA 1000000
122#define MET_MAX_HISTORY 100
123
124#define FLAG_FLUX 1 /* error L on gui */
125#define FLAG_SNR 2 /* error K on gui */
126#define FLAG_UNWRAP 4 /* error J on gui */
127#define FLAG_SPEED 8 /* error I on gui */
128#define FLAG_JUMP 16 /* error H on gui */
129#define FLAG_OVERLOAD 32 /* error G on gui */
130#define FLAG_UNLOCKED 64 /* error F on gui */
131#define FLAG_LASER_MV 128 /* error E on gui */
132#define FLAG_LASER_MW 256 /* error D on gui */
133#define FLAG_LASER_WAVE 512 /* error C on gui */
134#define FLAG_VOLT 1024 /* error B on gui */
135#define FLAG_PHASOR0 2048 /* error A on gui */
136#define FLAG_LOST_PHASE 4096 /* error T on gui */
137
138#define FT 0
139#define SC 1
140#define BOTH 2
141
142#define SIN 0
143#define COS 1
144
145#define NCOL 80
146
147/* Structure definitions */
148
149typedef struct {
150
151 double min_allowed_flux_telescope[4][4][2]; /* minimum allowed flux per telescope and diode and side */
152 double min_allowed_flux_fiber_coupler[4][2]; /* minimum allowed flux per fiber coupler and side */
153 double min_allowed_snr_telescope[4][4][2]; /* minimum allowed SNR per telescope and diode and side */
154 double min_allowed_snr_fiber_coupler[4][2]; /* minimum allowed SNR per fiber coupler and side */
155 double min_allowed_snr_diff[4][4][2]; /* minimum allowed SNR per telescope and diode and side */
156 long number_to_average; /* how many cycles used for smoothing signal */
157 long number_for_rms; /* how many cycles used for calculation of rms */
158 long number_for_speed; /* how many cycles used for calculation of speed */
159 double nominal_wavelength; /* nominal wavelength of laser in nm */
160 double max_allowed_phase_difference; /* in rad */
161 double max_allowed_phase_speed; /* in rad per phase step */
162 double opl_base; /* proportionality constant between radian and meters */
163 double offset_volt_telescope[4][4][2][2]; /* offset voltage per telescope and diode and side and sin/cos */
164 double offset_volt_fiber_coupler[4][2][2]; /* offset voltage per fiber coupler and side and sin/cos */
165 double rms_flux_telescope[4][4][2]; /* rms of voltage per telescope and diode and side */
166 double rms_flux_fiber_coupler[4][2]; /* rms of voltage per fiber coupler and side */
167 double lockin_constant_telescope[4][4][2][2]; /* multiplicative constant for lockin voltage per telescope and diode and side and sin/cos */
168 double lockin_constant_fiber_coupler[4][2][2]; /* multiplicative constant for lockin voltage per fiber coupler and side and sin/cos */
169 double max_allowed_voltage;
170 double min_allowed_voltage;
171 long number_to_smooth_for_telescope; /* how many cycles are smoothed for difference Tel - FC phasor */
172 long sign_of_phase;
173 long check_for_jumps;
174 long calc_phase_speed;
175 double decrement_factor_speed; /* how much the weight of each speed value in buffer is decreased after a cycle */
176 double norm_speed; /* normalization factor for weighted speed adding */
177 double sigma_clip_speed; /* at what sigma level the speed predictor is clipped */
178
179} structTacConfiguration;
180
181
182typedef struct {
183 long sample_number; /* Just a TAC loop counter : 0, 1, 2, ... */
184 long buffer_idx_avg;
185
186 double volt_lockin_telescope[4][4][2][2]; /* per telescope and diode and side and sin/cos */
187 double sum_volt_lockin_telescope[4][4][2][2]; /* per telescope and diode and side and sin/cos */
188 double buffer_volt_lockin_telescope[4][4][2][2][MET_MAX_HISTORY]; /* per telescope and diode and side and sin/cos */
189 double flux_telescope[4][4][2];
190 double sum_flux_telescope[4][4][2]; /* 4 x 4 x FT/SC */
191 double sum_sq_flux_telescope[4][4][2]; /* 4 x 4 x FT/SC */
192 double buffer_flux_telescope[4][4][2][MET_MAX_HISTORY]; /* 4 x 4 x FT/SC */
193 double buffer_sq_flux_telescope[4][4][2][MET_MAX_HISTORY]; /* 4 x 4 x FT/SC */
194 double rms_flux_telescope[4][4][2]; /* 4 x 4 x FT/SC */
195 double snr_flux_telescope[4][4][2]; /* 4 x 4 x FT/SC */
196 double buffer_delta_phasor_telescope[4][4][2][2][MET_MAX_HISTORY];
197 double delta_phasor_telescope[4][4][2][2];
198 double best_estimate_phasor_telescope[4][4][2][2];
199 double best_estimate_flux_telescope[4][4][2];
200
201 double buffer_delta_phase_telescope[4][4][2][MET_MAX_HISTORY];
202 double sum_delta_phase_telescope[4][4][2];
203 double buffer_sq_delta_phase_telescope[4][4][2][MET_MAX_HISTORY];
204 double sum_sq_delta_phase_telescope[4][4][2];
205
206 double sum_speed_telescope[4][4][2];
207 double buffer_speed_telescope[4][4][2][MET_MAX_HISTORY];
208
209 int flag_volt_telescope[4][4][2];
210 int flag_flux_telescope[4][4][2]; /* 4 x 4 x FT/SC */
211 int flag_snr_telescope[4][4][2]; /* 4 x 4 x FT/SC */
212
213 double volt_lockin_fiber_coupler[4][2][2]; /* per fiber coupler and side and sin/cos */
214 double sum_volt_lockin_fiber_coupler[4][2][2]; /* per fiber coupler and side and sin/cos */
215 double buffer_volt_lockin_fiber_coupler[4][2][2][MET_MAX_HISTORY]; /* per fiber coupler and side and sin/cos */
216 double flux_fiber_coupler[4][2];
217 double sum_flux_fiber_coupler[4][2];
218 double sum_sq_flux_fiber_coupler[4][2];
219 double buffer_flux_fiber_coupler[4][2][MET_MAX_HISTORY];
220 double buffer_sq_flux_fiber_coupler[4][2][MET_MAX_HISTORY];
221 double rms_flux_fiber_coupler[4][2];
222 double snr_flux_fiber_coupler[4][2];
223 double buffer_delta_phasor_fiber_coupler[4][2][2][MET_MAX_HISTORY];
224 double delta_phasor_fiber_coupler[4][2][2];
225 double best_estimate_phasor_fiber_coupler[4][2][2];
226 double best_estimate_flux_fiber_coupler[4][2];
227 double buffer_best_estimate_phasor_fiber_coupler[4][2][2][MET_MAX_HISTORY];
228
229 double buffer_delta_phase_fiber_coupler[4][2][MET_MAX_HISTORY];
230 double sum_delta_phase_fiber_coupler[4][2];
231 double buffer_sq_delta_phase_fiber_coupler[4][2][MET_MAX_HISTORY];
232 double sum_sq_delta_phase_fiber_coupler[4][2];
233
234 double sum_speed_fiber_coupler[4][2];
235 double buffer_speed_fiber_coupler[4][2][MET_MAX_HISTORY];
236 double buffer_used_speed_fiber_coupler[4][2][MET_MAX_HISTORY];
237
238 int flag_volt_fiber_coupler[4][2];
239 int flag_flux_fiber_coupler[4][2]; /* 4 fiber couplers, sin/cos */
240 int flag_snr_fiber_coupler[4][2]; /* 4 fiber couplers, sin/cos */
241
242 double raw_phase_telescope[4][4][2]; /* 4 x 4 (tel, diode), FT/SC */
243 double prev_phase_telescope[4][4][2];
244 double unwrapped_phase_telescope[4][4][2];
245 int flag_unwrap_telescope[4][4][2];
246 int flag_jump_telescope[4][4][2];
247 int flag_speed_telescope[4][4][2];
248 int flag_phasor0_telescope[4][4][2];
249 double delta_phase_telescope[4][4]; /* FT - SC */
250 double mean_phase_telescope[4]; /* mean over tge 4 diodes */
251
252 double raw_phase_fiber_coupler[4][2]; /* 4 fiber couplers , FT/SC */
253 double prev_phase_fiber_coupler[4][2];
254 double unwrapped_phase_fiber_coupler[4][2];
255 double buffer_unwrapped_phase_fiber_coupler[4][2][MET_MAX_HISTORY];
256 int flag_unwrap_fiber_coupler[4][2];
257 int flag_jump_fiber_coupler[4][2];
258 int flag_speed_fiber_coupler[4][2];
259 int flag_phasor0_fiber_coupler[4][2];
260 double delta_phase_fiber_coupler[4]; /* FT - SC */
261
262 double raw_phase_diff[4][4][2]; /* difference: tel 1/2/3/4 minus FC */
263 double prev_phase_diff[4][4][2];
264 double unwrapped_phase_diff[4][4][2];
265 int flag_unwrap_diff[4][4][2];
266 int flag_phasor0_diff[4][4][2];
267 int flag_snr_diff[4][4][2];
268 int flag_flux_diff[4][4][2];
269
270 double buffer_phasor_diff[4][4][2][2][MET_MAX_HISTORY]; /* in volt, tel diode, side, sin/cos, index */
271 double sum_phasor_diff[4][4][2][2]; /* in volt, tel diode, side, sin/cos */
272 double flux_diff[4][4][2]; /* phasor length in volt, tel, diode, side */
273 double sum_flux_diff[4][4][2]; /* 4 x 4 x FT/SC */
274 double sum_sq_flux_diff[4][4][2]; /* 4 x 4 x FT/SC */
275 double buffer_flux_diff[4][4][2][MET_MAX_HISTORY]; /* 4 x 4 x FT/SC */
276 double buffer_sq_flux_diff[4][4][2][MET_MAX_HISTORY]; /* 4 x 4 x FT/SC */
277 double rms_flux_diff[4][4][2]; /* 4 x 4 x FT/SC */
278 double snr_flux_diff[4][4][2]; /* 4 x 4 x FT/SC */
279
280 int total_flag_telescope[4][4][2]; /* 4 x 4 (tel, diode), FT/SC */
281 int total_flag_fiber_coupler[4][2]; /* 4 fiber couplers , FT/SC */
282
283 double opl_telescope[4]; /* diff FT-SC in meter */
284 double opl_telescope_diode[4][4]; /* diff FT-SC in meter */
285 double opl_fiber_coupler[4]; /* diff FT-SC in meter */
286 double opl_zero_telescope[4]; /* zero point in meter */
287 double opl_zero_telescope_diode[4][4]; /* zero point in meter */
288 double opl_zero_fiber_coupler[4]; /* zero point meter */
289
290 double start_phase_fiber_coupler[4][2]; /* in rad */
291 double start_phase_telescope[4][4][2]; /* in rad */
292
293 long freeze_count_fiber_coupler[4][2];
294 double freeze_phasor_fiber_coupler[4][2][2];
295 double freeze_speed_fiber_coupler[4][2];
296 double freeze_unwrapped_phase_fiber_coupler[4][2];
297 double used_unwrapped_phase_fiber_coupler[4][2];
298
299 structTacConfiguration * tacConfiguration;
300
301} structTacData;
302
303/*-----------------------------------------------------------------------------
304 TAC prototypes
305 -----------------------------------------------------------------------------*/
306
307double metrology_sq(double x);
308double myAtan(double x, double y, int* flag);
309
310structTacConfiguration * metrology_makeDefaultTacConfiguration (double lambda_met);
311
312structTacData * metrology_makeDefaultTacData (double lambda_met);
313
314int metrology_algorithm(structTacData * tacData);
315
316int metrology_unwrap(double raw_phase, double previous_phase, double max_allowed_phase_diff,
317 double *unwrapped_phase);
318
319int metrology_read_voltages(structTacData * tacData, double * volt);
320
321
322/*-----------------------------------------------------------------------------
323 TAC function code
324 -----------------------------------------------------------------------------*/
325
326structTacConfiguration * metrology_makeDefaultTacConfiguration(double lambda_met)
327{
328 structTacConfiguration *defaultTacConfiguration;
329
330 defaultTacConfiguration = cpl_malloc((size_t) sizeof(structTacConfiguration));
331
332 long tel, diode, side, comp, i;
333
334 for (tel = 0; tel < 4; tel++) {
335 for(side = FT; side <= SC; side++) {
336 for (diode = 0; diode < 4; diode++) {
337 for (comp = SIN; comp <= COS; comp++) {
338 defaultTacConfiguration->offset_volt_telescope[tel][diode][side][comp] = 0.0;
339 defaultTacConfiguration->lockin_constant_telescope[tel][diode][side][comp] = 1.0;
340 }
341 defaultTacConfiguration->min_allowed_flux_telescope[tel][diode][side] = 0.01;
342 defaultTacConfiguration->min_allowed_snr_telescope[tel][diode][side] = 2.0;
343 defaultTacConfiguration->min_allowed_snr_diff[tel][diode][side] = 2.0;
344 defaultTacConfiguration->rms_flux_telescope[tel][diode][side] = 0.005;
345 }
346 for (comp = SIN; comp <= COS; comp++) {
347 defaultTacConfiguration->offset_volt_fiber_coupler[tel][side][comp] = 0.00;
348 defaultTacConfiguration->lockin_constant_fiber_coupler[tel][side][comp] = 1.0;
349 }
350 defaultTacConfiguration->min_allowed_flux_fiber_coupler[tel][side] = 0.01;
351 defaultTacConfiguration->min_allowed_snr_fiber_coupler[tel][side] = 2.0;
352 defaultTacConfiguration->rms_flux_fiber_coupler[tel][side] = 0.004;
353 }
354 }
355
356 defaultTacConfiguration->min_allowed_voltage = -11.0;
357 defaultTacConfiguration->max_allowed_voltage = 11.0;
358
359 defaultTacConfiguration->number_to_average = 20;
360 defaultTacConfiguration->number_to_average = 20;
361 defaultTacConfiguration->number_for_rms = 50;
362
363 defaultTacConfiguration->max_allowed_phase_difference = TWOPI/3;
364 defaultTacConfiguration->max_allowed_phase_speed = TWOPI/4;
365
366 defaultTacConfiguration->nominal_wavelength = lambda_met * 1e9;
367 defaultTacConfiguration->opl_base = defaultTacConfiguration->nominal_wavelength / 1e9 / TWOPI;
368
369 defaultTacConfiguration->number_to_smooth_for_telescope = 50;
370 defaultTacConfiguration->sign_of_phase = 1;
371
372 defaultTacConfiguration->check_for_jumps = 1;
373 defaultTacConfiguration->calc_phase_speed = 1;
374
375 defaultTacConfiguration->number_for_speed = 3;
376 defaultTacConfiguration->decrement_factor_speed = 0.8;
377 defaultTacConfiguration->sigma_clip_speed = 0.0;
378
379 defaultTacConfiguration->norm_speed = 0;
380 for(i=0; i<defaultTacConfiguration->number_for_speed; i++) {
381 defaultTacConfiguration->norm_speed += pow(defaultTacConfiguration->decrement_factor_speed,i);
382 }
383
384 return defaultTacConfiguration;
385}
386
387structTacData * metrology_makeDefaultTacData(double lambda_met)
388{
389 long tel, diode, idx, side, comp;
390 structTacData *defaultTacData;
391
392 defaultTacData = cpl_malloc((size_t) sizeof(structTacData));
393
394 defaultTacData->tacConfiguration = metrology_makeDefaultTacConfiguration(lambda_met);
395
396 defaultTacData->sample_number = 0;
397
398 for (tel = 0; tel < 4; tel++) {
399 for(side = FT; side <= SC; side++) {
400 for (diode = 0; diode < 4; diode++) {
401 for (comp = SIN; comp <= COS; comp++) {
402 defaultTacData->volt_lockin_telescope[tel][diode][side][comp] = 0.0;
403 defaultTacData->sum_volt_lockin_telescope[tel][diode][side][comp] = 0.0;
404 defaultTacData->delta_phasor_telescope[tel][diode][side][comp] = 0.0;
405 defaultTacData->best_estimate_phasor_telescope[tel][diode][side][comp] = 0.0;
406 defaultTacData->sum_phasor_diff[tel][diode][side][comp] = 0.0;
407 }
408 defaultTacData->flux_telescope[tel][diode][side] = 0.0;
409 defaultTacData->rms_flux_telescope[tel][diode][side] = 0.0;
410 defaultTacData->snr_flux_telescope[tel][diode][side] = 0.0;
411 defaultTacData->sum_flux_telescope[tel][diode][side] = 0.0;
412 defaultTacData->sum_sq_flux_telescope[tel][diode][side] = 0.0;
413 defaultTacData->flag_flux_telescope[tel][diode][side] = 0;
414 defaultTacData->flag_snr_telescope[tel][diode][side] = 0;
415 defaultTacData->raw_phase_telescope[tel][diode][side] = 0.0;
416 defaultTacData->prev_phase_telescope[tel][diode][side] = 0.0;
417 defaultTacData->unwrapped_phase_telescope[tel][diode][side] = 0.0;
418 defaultTacData->flag_unwrap_telescope[tel][diode][side] = 0;
419 defaultTacData->flag_speed_telescope[tel][diode][side] = 0;
420 defaultTacData->flag_jump_telescope[tel][diode][side] = 0;
421 defaultTacData->flag_phasor0_telescope[tel][diode][side] = 0;
422 defaultTacData->total_flag_telescope[tel][diode][side] = 0;
423 defaultTacData->flag_volt_telescope[tel][diode][side] = 0;
424 defaultTacData->start_phase_telescope[tel][diode][side] = 0.0;
425 defaultTacData->best_estimate_flux_telescope[tel][diode][side] = 0.0;
426 defaultTacData->sum_speed_telescope[tel][diode][side] = 0.0;
427
428 for (idx = 0; idx < MET_MAX_HISTORY; idx++) {
429 for (comp = SIN; comp <= COS; comp++) {
430 defaultTacData->buffer_volt_lockin_telescope[tel][diode][side][comp][idx] = 0.0;
431 defaultTacData->buffer_delta_phasor_telescope[tel][diode][side][comp][idx] = 0.0;
432 defaultTacData->buffer_phasor_diff[tel][diode][side][comp][idx] = 0.0;
433 }
434 defaultTacData->buffer_delta_phase_telescope[tel][diode][side][idx] = 0.0;
435 defaultTacData->buffer_sq_delta_phase_telescope[tel][diode][side][idx] = 0.0;
436 defaultTacData->buffer_flux_telescope[tel][diode][side][idx] = 0.0;
437 defaultTacData->buffer_sq_flux_telescope[tel][diode][side][idx] = 0.0;
438 defaultTacData->buffer_speed_telescope[tel][diode][side][idx] = 0.0;
439 defaultTacData->buffer_flux_diff[tel][diode][side][idx] = 0.0;
440 defaultTacData->buffer_sq_flux_diff[tel][diode][side][idx] = 0.0;
441 }
442 defaultTacData->flux_diff[tel][diode][side] = 0.0;
443 defaultTacData->rms_flux_diff[tel][diode][side] = 0.0;
444 defaultTacData->snr_flux_diff[tel][diode][side] = 0.0;
445 defaultTacData->sum_flux_diff[tel][diode][side] = 0.0;
446 defaultTacData->sum_sq_flux_diff[tel][diode][side] = 0.0;
447
448 defaultTacData->sum_delta_phase_telescope[tel][diode][side] = 0.0;
449 defaultTacData->sum_sq_delta_phase_telescope[tel][diode][side] = 0.0;
450
451 defaultTacData->raw_phase_diff[tel][diode][side] = 0.0;
452 defaultTacData->prev_phase_diff[tel][diode][side] = 0.0;
453 defaultTacData->unwrapped_phase_diff[tel][diode][side] = 0.0;
454 defaultTacData->flag_unwrap_diff[tel][diode][side] = 0;
455 defaultTacData->flag_phasor0_diff[tel][diode][side] = 0;
456 defaultTacData->flag_snr_diff[tel][diode][side] = 0;
457 } /* end loop diode */
458 for (comp = SIN; comp <= COS; comp++) {
459 defaultTacData->volt_lockin_fiber_coupler[tel][side][comp] = 0.0;
460 defaultTacData->sum_volt_lockin_fiber_coupler[tel][side][comp] = 0.0;
461 defaultTacData->delta_phasor_fiber_coupler[tel][side][comp] = 0.0;
462 defaultTacData->best_estimate_phasor_fiber_coupler[tel][side][comp] = 0.0;
463 defaultTacData->freeze_phasor_fiber_coupler[tel][side][comp] = 0.0;
464 }
465 defaultTacData->flux_fiber_coupler[tel][side] = 0.0;
466 defaultTacData->rms_flux_fiber_coupler[tel][side] = 0.0;
467 defaultTacData->snr_flux_fiber_coupler[tel][side] = 0.0;
468 defaultTacData->sum_flux_fiber_coupler[tel][side] = 0.0;
469 defaultTacData->sum_sq_flux_fiber_coupler[tel][side] = 0.0;
470 defaultTacData->flag_flux_fiber_coupler[tel][side] = 0;
471 defaultTacData->flag_snr_fiber_coupler[tel][side] = 0;
472 defaultTacData->raw_phase_fiber_coupler[tel][side] = 0.0;
473 defaultTacData->prev_phase_fiber_coupler[tel][side] = 0.0;
474 defaultTacData->unwrapped_phase_fiber_coupler[tel][side] = 0.0;
475 defaultTacData->flag_unwrap_fiber_coupler[tel][side] = 0;
476 defaultTacData->flag_jump_fiber_coupler[tel][side] = 0;
477 defaultTacData->flag_speed_fiber_coupler[tel][side] = 0;
478 defaultTacData->flag_phasor0_fiber_coupler[tel][side] = 0;
479 defaultTacData->total_flag_fiber_coupler[tel][side] = 0;
480 defaultTacData->flag_volt_fiber_coupler[tel][side] = 0;
481 defaultTacData->start_phase_fiber_coupler[tel][side] = 0.0;
482 defaultTacData->best_estimate_flux_fiber_coupler[tel][side] = 0.0;
483 defaultTacData->sum_speed_fiber_coupler[tel][side] = 0.0;
484 defaultTacData->freeze_count_fiber_coupler[tel][side] = 0;
485 defaultTacData->freeze_speed_fiber_coupler[tel][side] = 0.0;
486 defaultTacData->freeze_unwrapped_phase_fiber_coupler[tel][side] = 0.0;
487 defaultTacData->used_unwrapped_phase_fiber_coupler[tel][side] = 0.0;
488
489 for (idx = 0; idx < MET_MAX_HISTORY; idx++) {
490 for (comp = SIN; comp <= COS; comp++) {
491 defaultTacData->buffer_volt_lockin_fiber_coupler[tel][side][comp][idx] = 0.0;
492 defaultTacData->buffer_delta_phasor_fiber_coupler[tel][side][comp][idx] = 0.0;
493 defaultTacData->buffer_best_estimate_phasor_fiber_coupler[tel][side][comp][idx] = 0.0;
494 }
495 defaultTacData->buffer_delta_phase_fiber_coupler[tel][side][idx] = 0.0;
496 defaultTacData->buffer_sq_delta_phase_fiber_coupler[tel][side][idx] = 0.0;
497 defaultTacData->buffer_flux_fiber_coupler[tel][side][idx] = 0.0;
498 defaultTacData->buffer_sq_flux_fiber_coupler[tel][side][idx] = 0.0;
499 defaultTacData->buffer_speed_fiber_coupler[tel][side][idx] = 0.0;
500 defaultTacData->buffer_used_speed_fiber_coupler[tel][side][idx] = 0.0;
501 defaultTacData->buffer_unwrapped_phase_fiber_coupler[tel][side][idx] = 0.0;
502 }
503 defaultTacData->sum_delta_phase_fiber_coupler[tel][side] = 0.0;
504 defaultTacData->sum_sq_delta_phase_fiber_coupler[tel][side] = 0.0;
505 } /* end loop over side */
506 defaultTacData->mean_phase_telescope[tel] = 0;
507
508 defaultTacData->delta_phase_fiber_coupler[tel] = 0.0;
509 for (diode = 0; diode < 4; diode++) {
510 defaultTacData->delta_phase_telescope[tel][diode] = 0.0;
511 defaultTacData->opl_telescope_diode[tel][diode] = 0.0;
512 defaultTacData->opl_zero_telescope_diode[tel][diode] = 0.0;
513 }
514
515 defaultTacData->opl_fiber_coupler[tel] = 0.0;
516 defaultTacData->opl_zero_fiber_coupler[tel] = 0.0;
517 defaultTacData->opl_telescope[tel] = 0.0;
518 defaultTacData->opl_zero_telescope[tel] = 0.0;
519
520 } /* end loop over tel */
521
522 return defaultTacData;
523
524}
525
526int metrology_read_voltages(structTacData * tacData, double * volt)
527{
528 int err = 0;
529 int tel, diode, side, comp;
530 int idx = 0;
531 long buffer_idx_avg = tacData->buffer_idx_avg;
532 double sqflux, flux;
533 structTacConfiguration * tacConfiguration = tacData->tacConfiguration;
534 long buffer_idx_rms = ((tacData->sample_number - 1) % tacConfiguration->number_for_rms);
535
536
537 for (side = FT; side <= SC; side++) {
538 for (tel = 0; tel < 4; tel++) {
539 for (diode = 0; diode < 4; diode++) {
540 for (comp = SIN; comp <= COS; comp++) {
541 tacData->volt_lockin_telescope[tel][diode][side][comp] = volt[idx];
542 idx++;
543 }
544 }
545 }
546 }
547 for (side = FT; side <= SC; side++) {
548 for (tel = 0; tel < 4; tel++) {
549 for (comp = SIN; comp <= COS; comp++) {
550 tacData->volt_lockin_fiber_coupler[tel][side][comp] = volt[idx];
551 idx++;
552 }
553 }
554 }
555
556 for (tel = 0; tel < 4; tel++) {
557 for (side = FT; side <= SC; side++) {
558 for (diode = 0; diode < 4; diode++) {
559 tacData->flag_volt_telescope[tel][diode][side] = 0;
560 if( tacData->volt_lockin_telescope[tel][diode][side][SIN] > tacConfiguration->max_allowed_voltage ||
561 tacData->volt_lockin_telescope[tel][diode][side][SIN] < tacConfiguration->min_allowed_voltage ||
562 tacData->volt_lockin_telescope[tel][diode][side][COS] > tacConfiguration->max_allowed_voltage ||
563 tacData->volt_lockin_telescope[tel][diode][side][COS] < tacConfiguration->min_allowed_voltage) {
564 tacData->flag_volt_telescope[tel][diode][side] = FLAG_VOLT;
565 }
566 tacData->volt_lockin_telescope[tel][diode][side][SIN] -=
567 tacConfiguration->offset_volt_telescope[tel][diode][side][SIN];
568 tacData->volt_lockin_telescope[tel][diode][side][SIN] *=
569 tacConfiguration->lockin_constant_telescope[tel][diode][side][SIN];
570
571 tacData->volt_lockin_telescope[tel][diode][side][COS] -=
572 tacConfiguration->offset_volt_telescope[tel][diode][side][COS];
573 tacData->volt_lockin_telescope[tel][diode][side][COS] *=
574 tacConfiguration->lockin_constant_telescope[tel][diode][side][COS];
575
576 }
577 tacData->flag_volt_fiber_coupler[tel][side] = 0;
578 if( tacData->volt_lockin_fiber_coupler[tel][side][SIN] > tacConfiguration->max_allowed_voltage ||
579 tacData->volt_lockin_fiber_coupler[tel][side][SIN] < tacConfiguration->min_allowed_voltage ||
580 tacData->volt_lockin_fiber_coupler[tel][side][COS] > tacConfiguration->max_allowed_voltage ||
581 tacData->volt_lockin_fiber_coupler[tel][side][COS] < tacConfiguration->min_allowed_voltage) {
582 tacData->flag_volt_fiber_coupler[tel][side] = FLAG_VOLT;
583 }
584 tacData->volt_lockin_fiber_coupler[tel][side][SIN] -=
585 tacConfiguration->offset_volt_fiber_coupler[tel][side][SIN];
586 tacData->volt_lockin_fiber_coupler[tel][side][SIN] *=
587 tacConfiguration->lockin_constant_fiber_coupler[tel][side][SIN];
588 tacData->volt_lockin_fiber_coupler[tel][side][COS] -=
589 tacConfiguration->offset_volt_fiber_coupler[tel][side][COS];
590 tacData->volt_lockin_fiber_coupler[tel][side][COS] *=
591 tacConfiguration->lockin_constant_fiber_coupler[tel][side][COS];
592 }
593 }
594 for (tel = 0; tel < 4; tel++) {
595 for (side = FT; side <= SC; side++) {
596 for (diode = 0; diode < 4; diode++) {
597
598 tacData->buffer_volt_lockin_telescope[tel][diode][side][SIN][buffer_idx_avg] = tacData->volt_lockin_telescope[tel][diode][side][SIN];
599 tacData->buffer_volt_lockin_telescope[tel][diode][side][COS][buffer_idx_avg] = tacData->volt_lockin_telescope[tel][diode][side][COS];
600 sqflux = metrology_sq(tacData->volt_lockin_telescope[tel][diode][side][SIN]) + metrology_sq(tacData->volt_lockin_telescope[tel][diode][side][COS]);
601 flux=sqrt(sqflux);
602 tacData->flux_telescope[tel][diode][side] = flux;
603 tacData->sum_flux_telescope[tel][diode][side] -= tacData->buffer_flux_telescope[tel][diode][side][buffer_idx_rms];
604 tacData->sum_sq_flux_telescope[tel][diode][side] -= tacData->buffer_sq_flux_telescope[tel][diode][side][buffer_idx_rms];
605 tacData->buffer_flux_telescope[tel][diode][side][buffer_idx_rms] = flux;
606 tacData->buffer_sq_flux_telescope[tel][diode][side][buffer_idx_rms] = sqflux;
607 tacData->sum_flux_telescope[tel][diode][side] += flux;
608 tacData->sum_sq_flux_telescope[tel][diode][side] += sqflux;
609
610 } /* end of loop over diodes */
611
612 tacData->buffer_volt_lockin_fiber_coupler[tel][side][SIN][buffer_idx_avg] = tacData->volt_lockin_fiber_coupler[tel][side][SIN];
613 tacData->buffer_volt_lockin_fiber_coupler[tel][side][COS][buffer_idx_avg] = tacData->volt_lockin_fiber_coupler[tel][side][COS];
614 sqflux = metrology_sq(tacData->volt_lockin_fiber_coupler[tel][side][SIN]) + metrology_sq(tacData->volt_lockin_fiber_coupler[tel][side][COS]);
615 flux=sqrt(sqflux);
616 tacData->flux_fiber_coupler[tel][side] = flux;
617 tacData->sum_flux_fiber_coupler[tel][side] -= tacData->buffer_flux_fiber_coupler[tel][side][buffer_idx_rms];
618 tacData->sum_sq_flux_fiber_coupler[tel][side] -= tacData->buffer_sq_flux_fiber_coupler[tel][side][buffer_idx_rms];
619 tacData->buffer_flux_fiber_coupler[tel][side][buffer_idx_rms] = flux;
620 tacData->buffer_sq_flux_fiber_coupler[tel][side][buffer_idx_rms] = sqflux;
621 tacData->sum_flux_fiber_coupler[tel][side] += flux;
622 tacData->sum_sq_flux_fiber_coupler[tel][side] += sqflux;
623
624 } /* end side loop */
625 } /* end telescope loop */
626
627 return err;
628}
629
630int metrology_unwrap(double raw_phase, double previous_phase, double max_allowed_phase_diff,
631 double *unwrapped_phase) {
632
633 double previous_phase_mod, twopicount, diff;
634 int flag = 0;
635
636 twopicount = floor(previous_phase / TWOPI); /* the number of full circles */
637 previous_phase_mod = previous_phase - TWOPI * twopicount; /* subtract the number of full circles to get the actual fraction of circle */
638
639 if (previous_phase_mod > PI) { /* if the fraction is bigger the PI subtract 2PI to get the rest */
640 previous_phase_mod -= TWOPI;
641 }
642
643 diff = raw_phase - previous_phase_mod;
644
645 if (diff > PI) {
646 diff -= TWOPI;
647 } else if (diff < (-1.0 * PI)) {
648 diff += TWOPI;
649 }
650
651 if (fabs(diff) > max_allowed_phase_diff) {
652 flag = FLAG_UNWRAP;
653 }
654
655 *unwrapped_phase = (previous_phase + diff);
656
657 return flag;
658}
659
660double metrology_sq(double x) {
661 return x * x;
662}
663
664double myAtan(double x, double y, int* flag)
665{
666 if(x==0 && y == 0) {
667 *flag = FLAG_PHASOR0;
668 return 0.;
669 }
670 *flag = 0;
671 return atan2(y,x);
672}
673
674
675int metrology_algorithm(structTacData * tacData)
676{
677 int err = 0;
678 long tel, diode, side;
679 long i, idx;
680 /* int old_flag_telescope[4][4][2];
681 int old_flag_fiber_coupler[4][2]; */
682 double x0, y0, x1, y1, xc, yc, delta_phase, sum_cos, sum_sin, new_phase, sq1, rms, tmp, speed;
683
684 int flag, comb_flag, sum_flag;
685 double sqflux, flux;
686 double cos_delta_phase, sin_delta_phase;
687 double *x0p, *y0p;
688
689 structTacConfiguration *tacConfiguration = tacData->tacConfiguration;
690 long sample_number = tacData->sample_number;
691
692 long buffer_idx_smooth_tel =((sample_number - 1) % tacConfiguration->number_to_smooth_for_telescope);
693 long buffer_idx_avg =((sample_number - 1) % tacConfiguration->number_to_average);
694 long prev_idx_avg = buffer_idx_avg - 1;
695 if(prev_idx_avg < 0) {
696 prev_idx_avg += tacConfiguration->number_to_average;
697 }
698
699 long buffer_idx_rms = ((sample_number - 1) % tacConfiguration->number_for_rms);
700 long buffer_idx_speed = ((sample_number - 1) % tacConfiguration->number_for_speed);
701
702 long prev_idx_speed = buffer_idx_speed - 1;
703 if(prev_idx_speed < 0) {
704 prev_idx_speed += tacConfiguration->number_for_speed;
705 }
706
707
708 /* check if phasor did not move too far */
709
710 if(sample_number > tacConfiguration->number_to_average) {
711 for (tel = 0; tel < 4; tel++) {
712 for (side = FT; side <= SC; side++) {
713 tacData->flag_jump_fiber_coupler[tel][side] = 0;
714 if((tacData->flux_fiber_coupler[tel][side] > tacConfiguration->min_allowed_flux_fiber_coupler[tel][side]) &&
715 (tacData->flux_fiber_coupler[tel][side] / tacConfiguration->rms_flux_fiber_coupler[tel][side] >
716 tacConfiguration->min_allowed_snr_fiber_coupler[tel][side]) &&
717 (tacData->snr_flux_fiber_coupler[tel][side] > tacConfiguration->min_allowed_snr_fiber_coupler[tel][side]) &&
718 (tacConfiguration->check_for_jumps == 1) ) {
719 x0 = tacData->volt_lockin_fiber_coupler[tel][side][COS];
720 y0 = tacData->volt_lockin_fiber_coupler[tel][side][SIN];
721 x1 = tacData->best_estimate_phasor_fiber_coupler[tel][side][COS];
722 y1 = tacData->best_estimate_phasor_fiber_coupler[tel][side][SIN];
723 xc = x1 * x0 + y1 * y0;
724 yc = x1 * y0 - y1 * x0;
725 if(fabs(myAtan(xc,yc,&flag)) > tacConfiguration->max_allowed_phase_difference) {
726 tacData->flag_jump_fiber_coupler[tel][side] = FLAG_JUMP;
727 x1 = tacData->buffer_volt_lockin_fiber_coupler[tel][side][COS][prev_idx_avg];
728 y1 = tacData->buffer_volt_lockin_fiber_coupler[tel][side][SIN][prev_idx_avg];
729 sq1 = sqrt(x1*x1+y1*y1);
730 if((sq1 > tacConfiguration->min_allowed_flux_fiber_coupler[tel][side]) &&
731 (sq1 / tacConfiguration->rms_flux_fiber_coupler[tel][side] >
732 tacConfiguration->min_allowed_snr_fiber_coupler[tel][side]) &&
733 (sq1 / tacData->rms_flux_fiber_coupler[tel][side] >
734 tacConfiguration->min_allowed_snr_fiber_coupler[tel][side]) ) {
735 xc = x1 * x0 + y1 * y0;
736 yc = x1 * y0 - y1 * x0;
737 if(fabs(myAtan(xc,yc,&flag)) < tacConfiguration->max_allowed_phase_difference) {
738 tacData->flag_jump_fiber_coupler[tel][side] = 0;
739 }
740 }
741 }
742 }
743 if(tacConfiguration->number_to_smooth_for_telescope < 0) {
744 for(diode = 0; diode < 4; diode++) {
745 if((tacData->flux_telescope[tel][diode][side] > tacConfiguration->min_allowed_flux_telescope[tel][diode][side]) &&
746 (tacData->flux_telescope[tel][diode][side] / tacConfiguration->rms_flux_telescope[tel][diode][side] >
747 tacConfiguration->min_allowed_snr_telescope[tel][diode][side]) &&
748 (tacData->snr_flux_telescope[tel][diode][side] > tacConfiguration->min_allowed_snr_telescope[tel][diode][side]) &&
749 (tacConfiguration->check_for_jumps == 1) ) {
750 tacData->flag_jump_telescope[tel][diode][side] = 0;
751 x0 = tacData->volt_lockin_telescope[tel][diode][side][COS];
752 y0 = tacData->volt_lockin_telescope[tel][diode][side][SIN];
753 x1 = tacData->best_estimate_phasor_telescope[tel][diode][side][COS];
754 y1 = tacData->best_estimate_phasor_telescope[tel][diode][side][SIN];
755 xc = x1 * x0 + y1 * y0;
756 yc = x1 * y0 - y1 * x0;
757
758 if(fabs(myAtan(xc,yc,&flag)) > tacConfiguration->max_allowed_phase_difference) {
759 tacData->flag_jump_telescope[tel][diode][side] = FLAG_JUMP;
760 x1 = tacData->buffer_volt_lockin_telescope[tel][diode][side][COS][prev_idx_avg];
761 y1 = tacData->buffer_volt_lockin_telescope[tel][diode][side][SIN][prev_idx_avg];
762 sq1 = sqrt(x1*x1+y1*y1);
763 if((sq1 > tacConfiguration->min_allowed_flux_telescope[tel][diode][side]) &&
764 (sq1 / tacConfiguration->rms_flux_telescope[tel][diode][side] >
765 tacConfiguration->min_allowed_snr_telescope[tel][diode][side]) &&
766 (sq1 / tacData->rms_flux_telescope[tel][diode][side] >
767 tacConfiguration->min_allowed_snr_telescope[tel][diode][side])) {
768 xc = x1 * x0 + y1 * y0;
769 yc = x1 * y0 - y1 * x0;
770 if(fabs(myAtan(xc,yc,&flag)) < tacConfiguration->max_allowed_phase_difference) {
771 tacData->flag_jump_telescope[tel][diode][side] = 0;
772 }
773 }
774 }
775 }
776 }
777 } /* end if num_smooth_tel < 0 */
778 }
779 } /* end loop tel */
780 }
781
782 /* raw phase calculation */
783
784 for (tel = 0; tel < 4; tel++) {
785 for (side = FT; side <= SC; side++) {
786 x0 = tacData->volt_lockin_fiber_coupler[tel][side][COS];
787 y0 = tacData->volt_lockin_fiber_coupler[tel][side][SIN];
788 x0 += tacData->best_estimate_phasor_fiber_coupler[tel][side][COS];
789 y0 += tacData->best_estimate_phasor_fiber_coupler[tel][side][SIN];
790 tacData->raw_phase_fiber_coupler[tel][side] = myAtan(x0,y0,&(tacData->flag_phasor0_fiber_coupler[tel][side]));
791
792 if(tacConfiguration->number_to_smooth_for_telescope < 0) {
793 for(diode = 0; diode < 4; diode++) {
794 x0 = tacData->volt_lockin_telescope[tel][diode][side][COS];
795 y0 = tacData->volt_lockin_telescope[tel][diode][side][SIN];
796 x0 += tacData->best_estimate_phasor_telescope[tel][diode][side][COS];
797 y0 += tacData->best_estimate_phasor_telescope[tel][diode][side][SIN];
798 tacData->raw_phase_telescope[tel][diode][side] = myAtan(x0,y0,&(tacData->flag_phasor0_telescope[tel][diode][side]));
799 }
800 }
801 }
802 }
803
804 /* unwrapping of phase */
805
806 for (tel = 0; tel < 4; tel++) {
807 for (side = FT; side <= SC; side++) {
808 if(sample_number == 1) {
809 tacData->prev_phase_fiber_coupler[tel][side] = tacData->raw_phase_fiber_coupler[tel][side];
810 } else {
811 tacData->prev_phase_fiber_coupler[tel][side] = tacData->unwrapped_phase_fiber_coupler[tel][side];
812 }
813
814 tacData->flag_unwrap_fiber_coupler[tel][side] = metrology_unwrap(tacData->raw_phase_fiber_coupler[tel][side],
815 tacData->prev_phase_fiber_coupler[tel][side], tacConfiguration->max_allowed_phase_difference,
816 &tacData->unwrapped_phase_fiber_coupler[tel][side]);
817
818 tacData->buffer_unwrapped_phase_fiber_coupler[tel][side][buffer_idx_avg] = tacData->unwrapped_phase_fiber_coupler[tel][side];
819
820 if(sample_number < tacConfiguration->number_to_average) {
821 tacData->flag_unwrap_fiber_coupler[tel][side] = 0;
822 }
823
824 if(tacConfiguration->number_to_smooth_for_telescope < 0) {
825 for(diode = 0; diode < 4; diode++) {
826 if(sample_number == 1) {
827 tacData->prev_phase_telescope[tel][diode][side] = tacData->raw_phase_telescope[tel][diode][side];
828 } else {
829 tacData->prev_phase_telescope[tel][diode][side] = tacData->unwrapped_phase_telescope[tel][diode][side];
830 }
831
832 tacData->flag_unwrap_telescope[tel][diode][side] = metrology_unwrap(tacData->raw_phase_telescope[tel][diode][side],
833 tacData->prev_phase_telescope[tel][diode][side], tacConfiguration->max_allowed_phase_difference,
834 &tacData->unwrapped_phase_telescope[tel][diode][side]);
835
836 if(sample_number < tacConfiguration->number_to_average) {
837 tacData->flag_unwrap_telescope[tel][diode][side] = 0;
838 }
839 }
840 }
841 }
842 }
843
844 /* calculate next best estimate for phasor */
845
846 for (tel = 0; tel < 4; tel++) {
847 for (side = FT; side <= SC; side++) {
848 if(tacConfiguration->number_to_average > 1) {
849 tacData->delta_phasor_fiber_coupler[tel][side][COS] -=
850 tacData->buffer_delta_phasor_fiber_coupler[tel][side][COS][buffer_idx_avg];
851 tacData->delta_phasor_fiber_coupler[tel][side][SIN] -=
852 tacData->buffer_delta_phasor_fiber_coupler[tel][side][SIN][buffer_idx_avg];
853
854 x0 = tacData->buffer_volt_lockin_fiber_coupler[tel][side][COS][prev_idx_avg];
855 y0 = tacData->buffer_volt_lockin_fiber_coupler[tel][side][SIN][prev_idx_avg];
856 x1 = tacData->volt_lockin_fiber_coupler[tel][side][COS];
857 y1 = tacData->volt_lockin_fiber_coupler[tel][side][SIN];
858 xc = x1 * x0 + y1 * y0;
859 yc = x1 * y0 - y1 * x0;
860 tacData->buffer_delta_phasor_fiber_coupler[tel][side][COS][buffer_idx_avg] = xc;
861 tacData->buffer_delta_phasor_fiber_coupler[tel][side][SIN][buffer_idx_avg] = yc;
862
863 tacData->delta_phasor_fiber_coupler[tel][side][COS] += xc;
864 tacData->delta_phasor_fiber_coupler[tel][side][SIN] += yc;
865
866 tacData->sum_delta_phase_fiber_coupler[tel][side] -= tacData->buffer_delta_phase_fiber_coupler[tel][side][buffer_idx_avg];
867 tacData->sum_sq_delta_phase_fiber_coupler[tel][side] -= tacData->buffer_sq_delta_phase_fiber_coupler[tel][side][buffer_idx_avg];
868
869 tacData->buffer_delta_phase_fiber_coupler[tel][side][buffer_idx_avg] = myAtan(yc, xc, &flag);
870 tacData->buffer_sq_delta_phase_fiber_coupler[tel][side][buffer_idx_avg] = metrology_sq(tacData->buffer_delta_phase_fiber_coupler[tel][side][buffer_idx_avg]);
871
872 tacData->sum_delta_phase_fiber_coupler[tel][side] += tacData->buffer_delta_phase_fiber_coupler[tel][side][buffer_idx_avg];
873 tacData->sum_sq_delta_phase_fiber_coupler[tel][side] += tacData->buffer_sq_delta_phase_fiber_coupler[tel][side][buffer_idx_avg];
874
875 rms = sqrt(tacData->sum_sq_delta_phase_fiber_coupler[tel][side]/tacConfiguration->number_to_average - metrology_sq(tacData->sum_delta_phase_fiber_coupler[tel][side]/tacConfiguration->number_to_average));
876
877 tacData->sum_speed_fiber_coupler[tel][side] -= tacData->buffer_speed_fiber_coupler[tel][side][buffer_idx_speed] * pow(tacConfiguration->decrement_factor_speed, tacConfiguration->number_for_speed - 1);
878
879 if(tacData->delta_phasor_fiber_coupler[tel][side][COS] == 0 && tacData->delta_phasor_fiber_coupler[tel][side][SIN] == 0) {
880 speed = tacData->buffer_speed_fiber_coupler[tel][side][prev_idx_speed];
881 } else {
882 speed = -myAtan(tacData->delta_phasor_fiber_coupler[tel][side][COS],tacData->delta_phasor_fiber_coupler[tel][side][SIN],&flag);
883 }
884
885 if(tacConfiguration->calc_phase_speed == 0) speed = 0.0;
886
887 if(tacConfiguration->sigma_clip_speed > 0) {
888 if(fabs(speed) * tacConfiguration->number_to_average < tacConfiguration->sigma_clip_speed * rms) speed = 0.0;
889 }
890
891 tacData->buffer_speed_fiber_coupler[tel][side][buffer_idx_speed] = speed;
892 tacData->sum_speed_fiber_coupler[tel][side] *= tacConfiguration->decrement_factor_speed;
893 tacData->sum_speed_fiber_coupler[tel][side] += speed;
894
895 delta_phase = tacData->sum_speed_fiber_coupler[tel][side] / tacConfiguration->norm_speed;
896
897 } else {
898 delta_phase = 0.;
899 }
900
901 tacData->buffer_used_speed_fiber_coupler[tel][side][buffer_idx_avg] = delta_phase;
902
903 tacData->flag_speed_fiber_coupler[tel][side] = 0;
904 if(sample_number > tacConfiguration->number_to_average) {
905 if(fabs(delta_phase) > tacConfiguration->max_allowed_phase_speed) {
906 tacData->flag_speed_fiber_coupler[tel][side] = FLAG_SPEED;
907 }
908 }
909
910 cos_delta_phase = cos(delta_phase);
911 sin_delta_phase = sin(delta_phase);
912
913 x1 = 1;
914 y1 = 0;
915
916 idx = buffer_idx_avg;
917 sum_cos = 0;
918 sum_sin = 0;
919 x0p = tacData->buffer_volt_lockin_fiber_coupler[tel][side][COS];
920 y0p = tacData->buffer_volt_lockin_fiber_coupler[tel][side][SIN];
921
922 for(i = 0; i < tacConfiguration->number_to_average; i++) {
923 x0 = x0p[idx];
924 y0 = y0p[idx];
925 sum_cos += x0 * x1 - y0 * y1;
926 sum_sin += x0 * y1 + x1 * y0;
927 idx -= 1;
928 if(idx < 0) {
929 idx += tacConfiguration->number_to_average;
930 }
931 tmp = x1 * cos_delta_phase - y1 * sin_delta_phase;
932 y1 = x1 * sin_delta_phase + y1 * cos_delta_phase;
933 x1 = tmp;
934 }
935
936 tacData->best_estimate_flux_fiber_coupler[tel][side] = sqrt(metrology_sq(sum_cos)+metrology_sq(sum_sin));
937 new_phase = myAtan(sum_cos, sum_sin, &flag) + delta_phase;
938
939 tacData->best_estimate_phasor_fiber_coupler[tel][side][COS] = tacData->best_estimate_flux_fiber_coupler[tel][side] * cos(new_phase);
940 tacData->best_estimate_phasor_fiber_coupler[tel][side][SIN] = tacData->best_estimate_flux_fiber_coupler[tel][side] * sin(new_phase);
941 tacData->best_estimate_flux_fiber_coupler[tel][side] /= tacConfiguration->number_to_average;
942 tacData->buffer_best_estimate_phasor_fiber_coupler[tel][side][COS][buffer_idx_avg] = tacData->best_estimate_phasor_fiber_coupler[tel][side][COS]/ tacConfiguration->number_to_average;
943 tacData->buffer_best_estimate_phasor_fiber_coupler[tel][side][SIN][buffer_idx_avg] = tacData->best_estimate_phasor_fiber_coupler[tel][side][SIN]/ tacConfiguration->number_to_average;
944
945 tacData->flag_flux_fiber_coupler[tel][side] = 0;
946 tacData->flag_snr_fiber_coupler[tel][side] = 0;
947 if (sample_number > tacConfiguration->number_to_average) {
948 if (tacData->best_estimate_flux_fiber_coupler[tel][side] < tacConfiguration->min_allowed_flux_fiber_coupler[tel][side]) {
949 tacData->flag_flux_fiber_coupler[tel][side] = FLAG_FLUX;
950 }
951 if (tacData->best_estimate_flux_fiber_coupler[tel][side] / tacConfiguration->rms_flux_fiber_coupler[tel][side] <
952 tacConfiguration->min_allowed_snr_fiber_coupler[tel][side]) {
953 tacData->flag_flux_fiber_coupler[tel][side] = FLAG_FLUX;
954 }
955 }
956
957 if(tacConfiguration->number_to_smooth_for_telescope < 0) {
958 for(diode = 0; diode < 4; diode++) {
959 if(tacConfiguration->number_to_average > 1) {
960 tacData->delta_phasor_telescope[tel][diode][side][COS] -=
961 tacData->buffer_delta_phasor_telescope[tel][diode][side][COS][buffer_idx_avg];
962 tacData->delta_phasor_telescope[tel][diode][side][SIN] -=
963 tacData->buffer_delta_phasor_telescope[tel][diode][side][SIN][buffer_idx_avg];
964
965 x0 = tacData->buffer_volt_lockin_telescope[tel][diode][side][COS][prev_idx_avg];
966 y0 = tacData->buffer_volt_lockin_telescope[tel][diode][side][SIN][prev_idx_avg];
967 x1 = tacData->volt_lockin_telescope[tel][diode][side][COS];
968 y1 = tacData->volt_lockin_telescope[tel][diode][side][SIN];
969 xc = x1 * x0 + y1 * y0;
970 yc = x1 * y0 - y1 * x0;
971 tacData->buffer_delta_phasor_telescope[tel][diode][side][COS][buffer_idx_avg] = xc;
972 tacData->buffer_delta_phasor_telescope[tel][diode][side][SIN][buffer_idx_avg] = yc;
973
974 tacData->delta_phasor_telescope[tel][diode][side][COS] += xc;
975 tacData->delta_phasor_telescope[tel][diode][side][SIN] += yc;
976
977 tacData->sum_delta_phase_telescope[tel][diode][side] -= tacData->buffer_delta_phase_telescope[tel][diode][side][buffer_idx_avg];
978 tacData->sum_sq_delta_phase_telescope[tel][diode][side] -= tacData->buffer_sq_delta_phase_telescope[tel][diode][side][buffer_idx_avg];
979
980 tacData->buffer_delta_phase_telescope[tel][diode][side][buffer_idx_avg] = myAtan(yc, xc, &flag);
981 tacData->buffer_sq_delta_phase_telescope[tel][diode][side][buffer_idx_avg] = metrology_sq(tacData->buffer_delta_phase_telescope[tel][diode][side][buffer_idx_avg]);
982
983 tacData->sum_delta_phase_telescope[tel][diode][side] += tacData->buffer_delta_phase_telescope[tel][diode][side][buffer_idx_avg];
984 tacData->sum_sq_delta_phase_telescope[tel][diode][side] += tacData->buffer_sq_delta_phase_telescope[tel][diode][side][buffer_idx_avg];
985
986 rms = sqrt(tacData->sum_sq_delta_phase_telescope[tel][diode][side]/tacConfiguration->number_to_average - metrology_sq(tacData->sum_delta_phase_telescope[tel][diode][side]/tacConfiguration->number_to_average));
987
988 tacData->sum_speed_telescope[tel][diode][side] -= tacData->buffer_speed_telescope[tel][diode][side][buffer_idx_speed] * pow(tacConfiguration->decrement_factor_speed, tacConfiguration->number_for_speed - 1);
989
990 if(tacData->delta_phasor_telescope[tel][diode][side][COS] == 0 && tacData->delta_phasor_telescope[tel][diode][side][SIN] == 0) {
991 speed = tacData->buffer_speed_telescope[tel][diode][side][prev_idx_speed];
992 } else {
993 speed = -myAtan(tacData->delta_phasor_telescope[tel][diode][side][COS],tacData->delta_phasor_telescope[tel][diode][side][SIN],&flag);
994 }
995
996 if(tacConfiguration->calc_phase_speed == 0) speed = 0.0;
997
998 if(tacConfiguration->sigma_clip_speed > 0) {
999 if(fabs(speed) * tacConfiguration->number_to_average < tacConfiguration->sigma_clip_speed * rms) speed = 0.0;
1000 }
1001
1002 tacData->buffer_speed_telescope[tel][diode][side][buffer_idx_speed] = speed;
1003 tacData->sum_speed_telescope[tel][diode][side] *= tacConfiguration->decrement_factor_speed;
1004 tacData->sum_speed_telescope[tel][diode][side] += speed;
1005
1006 delta_phase = tacData->sum_speed_telescope[tel][diode][side] / tacConfiguration->norm_speed;
1007
1008 } else {
1009 delta_phase = 0.;
1010 }
1011
1012 tacData->flag_speed_telescope[tel][diode][side] = 0;
1013 if(sample_number > tacConfiguration->number_to_average) {
1014 if(fabs(delta_phase) > tacConfiguration->max_allowed_phase_difference) {
1015 tacData->flag_speed_telescope[tel][diode][side] = FLAG_SPEED;
1016 }
1017 }
1018
1019 idx = buffer_idx_avg;
1020 sum_cos = 0;
1021 sum_sin = 0;
1022
1023 cos_delta_phase = cos(delta_phase);
1024 sin_delta_phase = sin(delta_phase);
1025
1026 x1 = 1;
1027 y1 = 0;
1028 x0p = tacData->buffer_volt_lockin_telescope[tel][diode][side][COS];
1029 y0p = tacData->buffer_volt_lockin_telescope[tel][diode][side][SIN];
1030
1031 for(i = 0; i < tacConfiguration->number_to_average; i++) {
1032 x0 = x0p[idx];
1033 y0 = y0p[idx];
1034 sum_cos += x0 * x1 - y0 * y1;
1035 sum_sin += x0 * y1 + x1 * y0;
1036 idx -= 1;
1037 if(idx < 0) {
1038 idx += tacConfiguration->number_to_average;
1039 }
1040 tmp = x1 * cos_delta_phase - y1 * sin_delta_phase;
1041 y1 = x1 * sin_delta_phase + y1 * cos_delta_phase;
1042 x1 = tmp;
1043 }
1044
1045 tacData->best_estimate_flux_telescope[tel][diode][side] = sqrt(metrology_sq(sum_cos)+metrology_sq(sum_sin));
1046 new_phase = myAtan(sum_cos, sum_sin, &flag) + delta_phase;
1047 tacData->best_estimate_phasor_telescope[tel][diode][side][COS] = tacData->best_estimate_flux_telescope[tel][diode][side] * cos(new_phase);
1048 tacData->best_estimate_phasor_telescope[tel][diode][side][SIN] = tacData->best_estimate_flux_telescope[tel][diode][side] * sin(new_phase);
1049 tacData->best_estimate_flux_telescope[tel][diode][side] /= tacConfiguration->number_to_average;
1050
1051 tacData->flag_flux_telescope[tel][diode][side] = 0;
1052 tacData->flag_snr_telescope[tel][diode][side] = 0;
1053 if (sample_number > tacConfiguration->number_to_average) {
1054 if (tacData->best_estimate_flux_telescope[tel][diode][side] < tacConfiguration->min_allowed_flux_telescope[tel][diode][side]) {
1055 tacData->flag_flux_telescope[tel][diode][side] = FLAG_FLUX;
1056 }
1057 if (tacData->best_estimate_flux_telescope[tel][diode][side] / tacConfiguration->rms_flux_telescope[tel][diode][side] <
1058 tacConfiguration->min_allowed_snr_telescope[tel][diode][side]) {
1059 tacData->flag_flux_telescope[tel][diode][side] = FLAG_FLUX;
1060 }
1061 }
1062 } /* end loop diode */
1063 } /* end if nsmooth_tel < 0*/
1064 } /* end loop side */
1065 } /* end loop telescope */
1066
1067 /* calclation of rms of flux and SNR check */
1068
1069 for (tel = 0; tel < 4; tel++) {
1070 for (side = FT; side <= SC; side++) {
1071 if(tacConfiguration->number_to_smooth_for_telescope < 0) {
1072 for(diode = 0; diode < 4; diode++) {
1073 tacData->rms_flux_telescope[tel][diode][side] =
1074 sqrt(tacData->sum_sq_flux_telescope[tel][diode][side]/tacConfiguration->number_for_rms
1075 - metrology_sq(tacData->sum_flux_telescope[tel][diode][side]/tacConfiguration->number_for_rms));
1076 tacData->snr_flux_telescope[tel][diode][side] = tacData->best_estimate_flux_telescope[tel][diode][side] /
1077 tacData->rms_flux_telescope[tel][diode][side];
1078 if (sample_number > tacConfiguration->number_for_rms) {
1079 if (tacData->snr_flux_telescope[tel][diode][side] <
1080 tacConfiguration->min_allowed_snr_telescope[tel][diode][side]) {
1081 tacData->flag_snr_telescope[tel][diode][side] = FLAG_SNR;
1082 }
1083 }
1084 }
1085 }
1086 tacData->rms_flux_fiber_coupler[tel][side] =
1087 sqrt(tacData->sum_sq_flux_fiber_coupler[tel][side]/tacConfiguration->number_for_rms
1088 - metrology_sq(tacData->sum_flux_fiber_coupler[tel][side]/tacConfiguration->number_for_rms));
1089 tacData->snr_flux_fiber_coupler[tel][side] = tacData->best_estimate_flux_fiber_coupler[tel][side] /
1090 tacData->rms_flux_fiber_coupler[tel][side];
1091
1092 if (sample_number > tacConfiguration->number_for_rms) {
1093 if (tacData->snr_flux_fiber_coupler[tel][side] <
1094 tacConfiguration->min_allowed_snr_fiber_coupler[tel][side]) {
1095 tacData->flag_snr_fiber_coupler[tel][side] = FLAG_SNR;
1096 }
1097 }
1098 }
1099 }
1100
1101 /* phase calculation for difference phasor Tel - FC */
1102
1103 if(tacConfiguration->number_to_smooth_for_telescope > 0) {
1104 for (tel = 0; tel < 4; tel++) {
1105 for (side = FT; side <= SC; side++) {
1106 x0 = tacData->volt_lockin_fiber_coupler[tel][side][COS];
1107 y0 = tacData->volt_lockin_fiber_coupler[tel][side][SIN];
1108 for (diode = 0; diode < 4; diode++) {
1109 x1 = tacData->volt_lockin_telescope[tel][diode][side][COS];
1110 y1 = tacData->volt_lockin_telescope[tel][diode][side][SIN];
1111 xc = x0 * x1 + y0 * y1;
1112 yc = x0 * y1 - y0 * x1;
1113
1114 /* maintain a buffer of x and y values to be able to average */
1115 tacData->sum_phasor_diff[tel][diode][side][COS] -= tacData->buffer_phasor_diff[tel][diode][side][COS][buffer_idx_smooth_tel];
1116 tacData->sum_phasor_diff[tel][diode][side][SIN] -= tacData->buffer_phasor_diff[tel][diode][side][SIN][buffer_idx_smooth_tel];
1117 tacData->buffer_phasor_diff[tel][diode][side][COS][buffer_idx_smooth_tel] = xc;
1118 tacData->buffer_phasor_diff[tel][diode][side][SIN][buffer_idx_smooth_tel] = yc;
1119 tacData->sum_phasor_diff[tel][diode][side][COS] += xc;
1120 tacData->sum_phasor_diff[tel][diode][side][SIN] += yc;
1121 /* the flux of the difference is calculated as the sum, not the average ! */
1122 sqflux = metrology_sq(tacData->sum_phasor_diff[tel][diode][side][COS])+
1123 metrology_sq(tacData->sum_phasor_diff[tel][diode][side][SIN]);
1124 flux = sqrt(sqflux);
1125 tacData->flux_diff[tel][diode][side] = flux;
1126
1127 tacData->raw_phase_diff[tel][diode][side] = myAtan(tacData->sum_phasor_diff[tel][diode][side][COS],
1128 tacData->sum_phasor_diff[tel][diode][side][SIN],
1129 &(tacData->flag_phasor0_diff[tel][diode][side]));
1130
1131 tacData->sum_flux_diff[tel][diode][side] -= tacData->buffer_flux_diff[tel][diode][side][buffer_idx_rms];
1132 tacData->sum_sq_flux_diff[tel][diode][side] -= tacData->buffer_sq_flux_diff[tel][diode][side][buffer_idx_rms];
1133 tacData->buffer_flux_diff[tel][diode][side][buffer_idx_rms] = flux;
1134 tacData->buffer_sq_flux_diff[tel][diode][side][buffer_idx_rms] = sqflux;
1135 tacData->sum_flux_diff[tel][diode][side] += flux;
1136 tacData->sum_sq_flux_diff[tel][diode][side] += sqflux;
1137
1138 tacData->rms_flux_diff[tel][diode][side] =
1139 sqrt(tacData->sum_sq_flux_diff[tel][diode][side]/tacConfiguration->number_for_rms
1140 - metrology_sq(tacData->sum_flux_diff[tel][diode][side]/tacConfiguration->number_for_rms));
1141 tacData->snr_flux_diff[tel][diode][side] = tacData->flux_diff[tel][diode][side] /
1142 tacData->rms_flux_diff[tel][diode][side];
1143
1144 tacData->flag_flux_diff[tel][diode][side] = 0;
1145 tacData->flag_snr_diff[tel][diode][side] = 0;
1146 if (sample_number > tacConfiguration->number_to_smooth_for_telescope) {
1147 if (tacData->flux_diff[tel][diode][side] < tacConfiguration->min_allowed_flux_telescope[tel][diode][side]) {
1148 tacData->flag_flux_diff[tel][diode][side] = FLAG_FLUX;
1149 }
1150 if (tacData->flux_diff[tel][diode][side] / tacConfiguration->rms_flux_telescope[tel][diode][side] <
1151 tacConfiguration->min_allowed_snr_diff[tel][diode][side]) {
1152 tacData->flag_flux_diff[tel][diode][side] = FLAG_FLUX;
1153 }
1154 if(sample_number > tacConfiguration->number_for_rms) {
1155 if (tacData->snr_flux_diff[tel][diode][side] <
1156 tacConfiguration->min_allowed_snr_diff[tel][diode][side]) {
1157 tacData->flag_snr_diff[tel][diode][side] = FLAG_SNR;
1158 }
1159 }
1160 }
1161 } /* end loop diode */
1162 } /* end loop side */
1163 } /* end loop tel */
1164
1165 /* unwrapping of difference phase */
1166 for (tel = 0; tel < 4; tel++) {
1167 for (side = FT; side <= SC; side++) {
1168 for (diode = 0; diode < 4; diode++) {
1169
1170 if(sample_number <= tacConfiguration->number_to_smooth_for_telescope) {
1171 tacData->prev_phase_diff[tel][diode][side] = tacData->raw_phase_diff[tel][diode][side];
1172 } else {
1173 tacData->prev_phase_diff[tel][diode][side] = tacData->unwrapped_phase_diff[tel][diode][side];
1174 }
1175
1176 tacData->flag_unwrap_diff[tel][diode][side] = metrology_unwrap(tacData->raw_phase_diff[tel][diode][side],
1177 tacData->prev_phase_diff[tel][diode][side], tacConfiguration->max_allowed_phase_difference,
1178 &tacData->unwrapped_phase_diff[tel][diode][side]);
1179
1180 } /* end loop diode */
1181 } /* end loop side */
1182 } /* end loop tel */
1183 } /* end if num_smooth_tel > 0 */
1184
1185
1186 /* filling the telescope flags and phase from difference phasor */
1187 if(tacConfiguration->number_to_smooth_for_telescope > 0) {
1188 for (tel = 0; tel < 4; tel++) {
1189 for (side = FT; side <= SC; side++) {
1190 for (diode = 0; diode < 4; diode++) {
1191 tacData->flag_unwrap_telescope[tel][diode][side] = tacData->flag_unwrap_diff[tel][diode][side];
1192 tacData->flag_flux_telescope[tel][diode][side] = tacData->flag_flux_diff[tel][diode][side];
1193 tacData->flag_snr_telescope[tel][diode][side] = tacData->flag_snr_diff[tel][diode][side];
1194 tacData->flag_phasor0_telescope[tel][diode][side] = tacData->flag_phasor0_diff[tel][diode][side];
1195 tacData->unwrapped_phase_telescope[tel][diode][side] = tacData->unwrapped_phase_fiber_coupler[tel][side] + tacData->unwrapped_phase_diff[tel][diode][side];
1196 }
1197 }
1198 }
1199 } /* end if num_smooth_tel > 0 */
1200
1201
1202
1203 /* Calculate the combined flags */
1204 for (tel = 0; tel < 4; tel++) {
1205 for (side = FT; side <= SC; side++) {
1206 comb_flag = 0;
1207 sum_flag = 0;
1208 for (diode = 0; diode < 4; diode++) {
1209 /* old_flag_telescope[tel][diode][side] = tacData->total_flag_telescope[tel][diode][side]; */
1210 tacData->total_flag_telescope[tel][diode][side] =
1211 tacData->flag_flux_telescope[tel][diode][side]
1212 + tacData->flag_unwrap_telescope[tel][diode][side]
1213 + tacData->flag_jump_telescope[tel][diode][side]
1214 + tacData->flag_speed_telescope[tel][diode][side]
1215 + tacData->flag_phasor0_telescope[tel][diode][side]
1216 + tacData->flag_snr_telescope[tel][diode][side]
1217 + tacData->flag_volt_telescope[tel][diode][side];
1218 comb_flag |= tacData->total_flag_telescope[tel][diode][side];
1219 if( tacData->total_flag_telescope[tel][diode][side] != 0) sum_flag++;
1220 } /* end loop diode */
1221
1222 /* old_flag_fiber_coupler[tel][side] = tacData->total_flag_fiber_coupler[tel][side]; */
1223 tacData->total_flag_fiber_coupler[tel][side] =
1224 tacData->flag_flux_fiber_coupler[tel][side]
1225 + tacData->flag_unwrap_fiber_coupler[tel][side]
1226 + tacData->flag_jump_fiber_coupler[tel][side]
1227 + tacData->flag_speed_fiber_coupler[tel][side]
1228 + tacData->flag_phasor0_fiber_coupler[tel][side]
1229 + tacData->flag_snr_fiber_coupler[tel][side]
1230 + tacData->flag_volt_fiber_coupler[tel][side];
1231 } /* end loop side */
1232 } /* end loop tel */
1233
1234 /* If there was an error, freeze + speed predict */
1235 for (tel = 0; tel < 4; tel++) {
1236 for (side = FT; side <= SC; side++) {
1237 if(tacData->total_flag_fiber_coupler[tel][side] != 0) {
1238 tacData->freeze_count_fiber_coupler[tel][side] += 1;
1239
1240 if(tacData->freeze_count_fiber_coupler[tel][side] == 1) {
1241
1242
1243 tacData->freeze_phasor_fiber_coupler[tel][side][COS] = tacData->buffer_best_estimate_phasor_fiber_coupler[tel][side][COS][prev_idx_avg];
1244 tacData->freeze_phasor_fiber_coupler[tel][side][SIN] = tacData->buffer_best_estimate_phasor_fiber_coupler[tel][side][SIN][prev_idx_avg];
1245 tacData->freeze_speed_fiber_coupler[tel][side] = tacData->buffer_used_speed_fiber_coupler[tel][side][prev_idx_avg];
1246
1247 if(tacConfiguration->calc_phase_speed == 0) tacData->freeze_speed_fiber_coupler[tel][side] = 0.0;
1248 if(tacData->flag_speed_fiber_coupler[tel][side] !=0) tacData->freeze_speed_fiber_coupler[tel][side] = 0.0;
1249 }
1250
1251 new_phase = myAtan(tacData->freeze_phasor_fiber_coupler[tel][side][COS], tacData->freeze_phasor_fiber_coupler[tel][side][SIN], &(tacData->flag_phasor0_fiber_coupler[tel][side])) + tacData->freeze_speed_fiber_coupler[tel][side];
1252 tmp = sqrt(metrology_sq(tacData->freeze_phasor_fiber_coupler[tel][side][COS])+metrology_sq(tacData->freeze_phasor_fiber_coupler[tel][side][SIN]));
1253 tacData->freeze_phasor_fiber_coupler[tel][side][COS] = tmp * cos(new_phase);
1254 tacData->freeze_phasor_fiber_coupler[tel][side][SIN] = tmp * sin(new_phase);
1255
1256 tmp = myAtan(tacData->freeze_phasor_fiber_coupler[tel][side][COS], tacData->freeze_phasor_fiber_coupler[tel][side][SIN], &(tacData->flag_phasor0_fiber_coupler[tel][side]));
1257
1258 } else {
1259 /* this should be the nominal case, no error occurred and the counter is 0 */
1260
1261 tacData->freeze_count_fiber_coupler[tel][side] = 0;
1262
1263 tmp = myAtan(cos(tacData->unwrapped_phase_fiber_coupler[tel][side]), sin(tacData->unwrapped_phase_fiber_coupler[tel][side]), &(tacData-> flag_phasor0_fiber_coupler[tel][side]));
1264
1265 if(tacData->flag_unwrap_fiber_coupler[tel][side] !=0 ) {
1266 printf("After freezing: speed predictor too far from where phase is, pred: %f, phase: %f\n",
1267 tacData->freeze_unwrapped_phase_fiber_coupler[tel][side],tacData->unwrapped_phase_fiber_coupler[tel][side]);
1268 }
1269 }
1270
1271 if( metrology_unwrap(tmp,tacData->used_unwrapped_phase_fiber_coupler[tel][side],
1272 tacConfiguration->max_allowed_phase_difference,
1273 &tacData->used_unwrapped_phase_fiber_coupler[tel][side]) != 0) {
1274
1275 if(sample_number >= tacConfiguration->number_to_average) {
1276
1277 tacData->total_flag_fiber_coupler[tel][side] += FLAG_LOST_PHASE ;
1278
1279 }
1280 }
1281
1282 } /* loop side */
1283 } /* loop tel */
1284
1285
1286 /* difference FT - SC; and calculate mean per telescope */
1287
1288 for (tel = 0; tel < 4; tel++) {
1289 tacData->mean_phase_telescope[tel] = 0.0;
1290 for (diode = 0; diode < 4; diode++) {
1291 tacData->delta_phase_telescope[tel][diode] = tacConfiguration->sign_of_phase * (tacData->unwrapped_phase_telescope[tel][diode][SC] - tacData->unwrapped_phase_telescope[tel][diode][FT]);
1292 tacData->mean_phase_telescope[tel] += tacData->delta_phase_telescope[tel][diode];
1293 }
1294 tacData->delta_phase_fiber_coupler[tel] = tacConfiguration->sign_of_phase * (tacData->used_unwrapped_phase_fiber_coupler[tel][SC] - tacData->used_unwrapped_phase_fiber_coupler[tel][FT]);
1295 tacData->mean_phase_telescope[tel] /= 4.0;
1296 }
1297
1298
1299 for (tel = 0; tel < 4; tel++) {
1300 for (diode = 0; diode < 4; diode++) {
1301 if ( tacData->total_flag_telescope[tel][diode][FT] == 0 && tacData->total_flag_telescope[tel][diode][SC] == 0) {
1302 tacData->opl_telescope_diode[tel][diode] = tacConfiguration->opl_base * tacData->delta_phase_telescope[tel][diode] - tacData->opl_zero_telescope_diode[tel][diode];
1303 }
1304 }
1305 tacData->opl_fiber_coupler[tel] = tacConfiguration->opl_base * tacData->delta_phase_fiber_coupler[tel]
1306 - tacData->opl_zero_fiber_coupler[tel];
1307
1308 }
1309
1310 return err;
1311}
1312
1313/*-----------------------------------------------------------------------------
1314 Functions code
1315 -----------------------------------------------------------------------------*/
1316
1317/*----------------------------------------------------------------------------*/
1331/*----------------------------------------------------------------------------*/
1332
1333cpl_error_code gravi_metrology_update_receiverpos (cpl_propertylist * header,
1334 cpl_table *receiver_table)
1335{
1336 gravi_msg_function_start(1);
1337 cpl_ensure_code (header, CPL_ERROR_NULL_INPUT);
1338 cpl_ensure_code (receiver_table, CPL_ERROR_NULL_INPUT);
1339
1340 /* Loop on telescope */
1341 for (int tel=0; tel<4; tel++) {
1342
1343 const char * telname = gravi_conf_get_telname (tel, header);
1344 if (telname == NULL) {
1345 cpl_msg_warning (cpl_func,"Cannot update receiver position for tel %i", tel);
1346 continue;
1347 }
1348
1349 /* Get row */
1350 cpl_size row;
1351 cpl_size nrow = cpl_table_get_nrow(receiver_table);
1352
1353 for (row = 0; row<nrow; row++) {
1354 if (!strcmp (telname, cpl_table_get_string (receiver_table, "TEL_NAME", row) )) break;
1355 }
1356 cpl_ensure_code (row<nrow, CPL_ERROR_ILLEGAL_INPUT);
1357
1358 /* Copy in header */
1359 for (int diode=0; diode<4; diode++) {
1360 char name[100];
1361
1362 /* Set in header */
1363 double posx = gravi_table_get_value (receiver_table,"RECX",row,diode);
1364 sprintf (name, "ESO MET %s REC%iX", telname, diode+1);
1365 cpl_propertylist_update_double (header, name, posx);
1366
1367 /* Set in header */
1368 double posy = gravi_table_get_value (receiver_table,"RECY",row,diode);
1369 sprintf (name, "ESO MET %s REC%iY", telname, diode+1);
1370 cpl_propertylist_update_double (header, name, posy);
1371
1372 cpl_msg_info (cpl_func, "Update diode %i of %s: x=%.3f, y=%.3f", diode+1, telname, posx, posy);
1373 }
1374 }
1375
1376 gravi_msg_function_exit(1);
1377 return CPL_ERROR_NONE;
1378}
1379
1380/*----------------------------------------------------------------------------*/
1390/*----------------------------------------------------------------------------*/
1391
1392double gravi_metrology_get_posx (cpl_propertylist * header,
1393 int tel, int diode)
1394{
1395 gravi_msg_function_start(0);
1396 cpl_ensure (header, CPL_ERROR_NULL_INPUT, 0);
1397
1398 /* Get telname */
1399 const char * telname = gravi_conf_get_telname (tel, header);
1400
1401 if (telname == NULL) {
1402 cpl_msg_warning (cpl_func,"Cannot read receiver x-position for tel %i (set 0.0)", tel);
1403 return 0.0;
1404 }
1405
1406 /* Read from header */
1407 char name[100];
1408 sprintf (name, "ESO MET %s REC%iX", telname, diode+1);
1409 double pos = cpl_propertylist_get_double (header, name)*1e-3;
1410
1411 gravi_msg_function_exit(0);
1412 return pos;
1413}
1414
1415double gravi_metrology_get_posy (cpl_propertylist * header,
1416 int tel, int diode)
1417{
1418 gravi_msg_function_start(0);
1419 cpl_ensure (header, CPL_ERROR_NULL_INPUT, 0);
1420
1421 /* Get telname */
1422 const char * telname = gravi_conf_get_telname (tel, header);
1423
1424 if (telname == NULL) {
1425 cpl_msg_warning (cpl_func,"Cannot read receiver y-position for tel %i", tel);
1426 return 0.0;
1427 }
1428
1429 /* Read from header */
1430 char name[100];
1431 sprintf (name, "ESO MET %s REC%iY", telname, diode+1);
1432 double pos = cpl_propertylist_get_double (header, name)*1e-3;
1433
1434 gravi_msg_function_exit(0);
1435 return pos;
1436}
1437
1438/*----------------------------------------------------------------------------*/
1447/*----------------------------------------------------------------------------*/
1448double gravi_metrology_get_fc_focus (cpl_propertylist * header, int gv, gravi_data *static_param_data)
1449{
1450 gravi_msg_function_start(0);
1451 cpl_ensure (header, CPL_ERROR_NULL_INPUT, 0);
1452// cpl_ensure (static_param_data, CPL_ERROR_NULL_INPUT, 0);
1453
1454 /* Identify telescope name of requested GV input */
1455 const char * telname = gravi_conf_get_telname (gv, header);
1456 if (telname == NULL) {
1457 cpl_msg_warning (cpl_func,"Cannot read fiber coupler focus offset for GV%i", gv+1);
1458 return 0.0;
1459 }
1460
1461 /* Loading defaults which will be used if no other choice with warning */
1462 double defocus=0.0;
1463 double defocus_at_default[4] = {-75.0, -100.0, 25.0, -75.0}; // Measured 2017-11-17
1464 double defocus_ut_default[4] = {-50.0, -125.0, -150.0, -175.0}; // Measured 2018-02-15
1465
1466 char column_name[80];
1467 int get_default=1;
1468
1469 /* Assemble column name
1470 * EKW 07/10/2019 : else shall be _at */
1471 if (telname[0] == 'U') strcpy(column_name,"fc_focus_ut");
1472 else strcpy(column_name,"fc_focus_at");
1473 if (gravi_pfits_get_axis (header) == MODE_ONAXIS) strcat(column_name,"_onaxis");
1474
1475 /* read value from calibration file */
1476 if (static_param_data)
1477 {
1478 cpl_table * static_param_table = gravi_data_get_table (static_param_data, "FOCUSPAR");
1479 if ( cpl_table_has_column(static_param_table , column_name) )
1480 {
1481 defocus = cpl_table_get_data_double (static_param_table, column_name)[gv];
1482 get_default=0;
1483 } else {
1484 cpl_msg_warning(cpl_func,"Cannot get the column %s from the calibration file",column_name);
1485 }
1486 } else {
1487 cpl_msg_error (cpl_func,"Cannot find static calibration file, using hard-coded values for fiber coupler focus shift for GV%i!", gv+1);
1488 }
1489
1490 /* if could not get the value for whatever reason, default to hardcoded values */
1491 if (get_default == 1)
1492 {
1493 cpl_msg_warning(cpl_func,"Using the default values for %s for GV%i",column_name,gv+1);
1494 if (telname[0] == 'U') {
1495 defocus=defocus_ut_default[gv];
1496 } else {
1497 defocus=defocus_at_default[gv];
1498 }
1499 }
1500
1501 cpl_msg_info(cpl_func,"value for %s on GV%i : %e",column_name,gv+1,defocus);
1502
1503 gravi_msg_function_exit(0);
1504 return defocus*1e-9; // Convert from [nm] in header to [m] in pipeline
1505}
1506
1507/*----------------------------------------------------------------------------*/
1518/*----------------------------------------------------------------------------*/
1519double gravi_metrology_get_fc_shift (cpl_propertylist * header, int gv, gravi_data *static_param_data)
1520{
1521 gravi_msg_function_start(0);
1522 cpl_ensure (header, CPL_ERROR_NULL_INPUT, 0);
1523 //cpl_ensure (static_param_data, CPL_ERROR_NULL_INPUT, 0);
1524
1525 /* Identify telescope name of requested GV input */
1526 const char * telname = gravi_conf_get_telname (gv, header);
1527 if (telname == NULL) {
1528 cpl_msg_warning (cpl_func,"Cannot read fiber coupler shift offset for GV%i", gv+1);
1529 return 0.0;
1530 }
1531
1532 /* Loading defaults which will be used if no other choice with warning */
1533 double shift=0.0;
1534 double shift_ut_default[4] = {-450.0, -350.0, -50.0, -525.0}; // Measured 2018-02-15
1535 double shift_at_default[4] = {-100.0*1.8/8.0, 50.0*1.8/8.0, 150.0*1.8/8.0, -100.0*1.8/8.0}; // Measured 2017-11-17
1536
1537 char column_name[80];
1538 int get_default=1;
1539
1540 /* Assemble column name
1541 * 08/10/2019 ug fix else fc_focus_shift_at */
1542 /* FE: "fc_focus_shift" is very misleading, because it is not related to "focus" at all, but only to
1543 the shift between fiber coupler pickup and pupil tracker reference spot. I don't know how the misleading
1544 name entered the fits header convention, and would propose to change at some point */
1545 if (telname[0] == 'U') strcpy(column_name,"fc_focus_shift_ut");
1546 else strcpy(column_name,"fc_focus_shift_at");
1547 if (gravi_pfits_get_axis (header) == MODE_ONAXIS) strcat(column_name,"_onaxis");
1548
1549 /* read value from calibration file */
1550 if (static_param_data)
1551 {
1552 cpl_table * static_param_table = gravi_data_get_table (static_param_data, "FOCUSPAR");
1553 if ( cpl_table_has_column(static_param_table , column_name) )
1554 {
1555 shift = cpl_table_get_data_double (static_param_table, column_name)[gv];
1556 get_default=0;
1557 } else {
1558 cpl_msg_warning(cpl_func,"Cannot get the column %s from the calibration file",column_name);
1559 }
1560 } else {
1561 cpl_msg_error (cpl_func,"Cannot find static calibration file, using hard-coded values for fiber coupler focus shift for GV%i!", gv+1);
1562 }
1563
1564 /* if could not get the value for whatever reason, default to hardcoded values */
1565 if (get_default == 1)
1566 {
1567 cpl_msg_warning(cpl_func,"Using the default values for %s for GV%i",column_name,gv+1);
1568 if (telname[0] == 'U') {
1569 shift=shift_ut_default[gv];
1570 } else {
1571 shift=shift_at_default[gv];
1572 }
1573 }
1574
1575 cpl_msg_info(cpl_func,"value for %s on GV%i : %e",column_name,gv+1,shift);
1576
1577 gravi_msg_function_exit(0);
1578 return shift*1e-9; // Convert from [nm/arcsec] in header to [m/arcsec] in pipeline
1579}
1580
1581/*----------------------------------------------------------------------------*/
1592/*----------------------------------------------------------------------------*/
1593cpl_error_code gravi_metrology_get_astig (cpl_propertylist * header, int gv,
1594 double * amplitude, double * angle,
1595 double * radius, gravi_data *static_param_data)
1596{
1597 gravi_msg_function_start(0);
1598 cpl_ensure (header, CPL_ERROR_NULL_INPUT, 0);
1599 cpl_ensure (static_param_data, CPL_ERROR_NULL_INPUT, 0);
1600
1601 /* load default "do nothing" values */
1602 *amplitude = 0.0;
1603 *angle = 0.0;
1604 *radius = 1.0;
1605
1606 /* variables to store the astigmatism values */
1607 double amplitude_header; // [nm]
1608 double angle_header; // [deg]]
1609 double amplitude_static; // [nm]
1610 double angle_static; // [deg]]
1611
1612 /* Identify telescope name of requested GV input */
1613 const char * telname = gravi_conf_get_telname (gv, header);
1614 if (telname == NULL) {
1615 cpl_msg_error (cpl_func,"Cannot read the astigmatism offset for GV%i", gv+1);
1616 return CPL_ERROR_NONE;
1617 }
1618
1619 /* Assemble header keywords */
1620 char name_amp[100], name_ang[100];
1621 if (telname[0] == 'U') {
1622 sprintf (name_amp, "ESO MET GV%i UT ASTIG AMP", gv+1);
1623 sprintf (name_ang, "ESO MET GV%i UT ASTIG ANG", gv+1);
1624 *radius = 4.0; // radius of an UT [m]
1625 } else {
1626 sprintf (name_amp, "ESO MET GV%i AT ASTIG AMP", gv+1);
1627 sprintf (name_ang, "ESO MET GV%i AT ASTIG ANG", gv+1);
1628 *radius = 0.9; // radius if an AT [m]
1629 }
1630 if (gravi_pfits_get_axis (header) == MODE_ONAXIS)
1631 {
1632 strcat(name_amp," ONA");
1633 strcat(name_ang," ONA");
1634 } else {
1635 strcat(name_amp," OFA");
1636 strcat(name_ang," OFA");
1637 }
1638
1639 /* Assemble columns names, to read the static file */
1640 char column_amp_name[80];
1641 if (telname[0] == 'U') strcpy(column_amp_name,"astig_amp_ut");
1642 else strcpy(column_amp_name,"astig_amp_at");
1643 if (gravi_pfits_get_axis (header) == MODE_ONAXIS) strcat(column_amp_name,"_onaxis");
1644 char column_ang_name[80];
1645 if (telname[0] == 'U') strcpy(column_ang_name,"astig_ang_ut");
1646 else strcpy(column_ang_name,"astig_ang_at");
1647 if (gravi_pfits_get_axis (header) == MODE_ONAXIS) strcat(column_ang_name,"_onaxis");
1648
1649 /* If static file exist, read them, otherwise read value from header */
1650 int got_static=0;
1651 if (static_param_data)
1652 {
1653 if (gravi_data_has_extension (static_param_data, "ASTIGPAR"))
1654 {
1655 cpl_table * static_param_table = gravi_data_get_table (static_param_data, "ASTIGPAR");
1656 if ( cpl_table_has_column(static_param_table , column_amp_name) & cpl_table_has_column(static_param_table , column_ang_name) )
1657 {
1658 amplitude_static = cpl_table_get_data_double (static_param_table, column_amp_name)[gv];
1659 angle_static = cpl_table_get_data_double (static_param_table, column_ang_name)[gv];
1660 got_static=1;
1661 cpl_msg_info(cpl_func,"Astigmatism (from static file): GV%i, Amplitude=%.2f nm, Angle=%.2f deg",
1662 gv+1,amplitude_static,angle_static);
1663 } else {
1664 cpl_msg_warning(cpl_func,"Cannot get the column %s or %s from the static calibration file",column_amp_name, column_ang_name);
1665 }
1666 } else {
1667 cpl_msg_warning(cpl_func,"Cannot get the sub-fit ASTIGPAR from the static calibration file");
1668 }
1669 } else {
1670 cpl_msg_warning (cpl_func,"Cannot find the static calibration file");
1671 }
1672
1673 /* If static file exist, read them, otherwise read value from header */
1674 int got_header=0;
1675 if (cpl_propertylist_has(header, name_amp)
1676 && cpl_propertylist_has(header, name_ang)) {
1677 amplitude_header = cpl_propertylist_get_double (header, name_amp);
1678 angle_header = cpl_propertylist_get_double (header, name_ang);
1679 got_header=1;
1680 cpl_msg_info(cpl_func,"Astigmatism (from header data): GV%i, Amplitude=%.2f nm, Angle=%.2f deg",
1681 gv+1,amplitude_header,angle_header);
1682 } else {
1683 cpl_msg_warning (cpl_func,"Cannot get astigmatism parameters from header");
1684 }
1685
1686 /* Here, the code decide which value to take */
1687 /* Also print the value in the log */
1688 if (got_static==1)
1689 {
1690 *amplitude = amplitude_static;
1691 *angle = angle_static;
1692 cpl_msg_info(cpl_func,"Using astigmatism parameters from static file");
1693 } else if (got_header==1)
1694 {
1695 *amplitude = amplitude_header;
1696 *angle = angle_header;
1697 cpl_msg_info(cpl_func,"Using astigmatism parameters from header");
1698 } else
1699 {
1700 cpl_msg_warning (cpl_func,"The static calibration file is missing astigmatism values");
1701 return CPL_ERROR_NONE;
1702 }
1703
1704 /* The astigmatism values are now converted with respect to the acquisition camera North Angle */
1705 if (gravi_pfits_get_axis (header) == MODE_ONAXIS)
1706 {
1707 *angle = *angle + 141;
1708 } else {
1709 *angle = *angle + 231;
1710 }
1711
1712
1713 gravi_msg_function_exit(0);
1714 return CPL_ERROR_NONE;
1715}
1716
1717/*----------------------------------------------------------------------------*/
1732/*----------------------------------------------------------------------------*/
1733
1734cpl_table * gravi_metrology_create (cpl_table * metrology_table,
1735 cpl_propertylist * header)
1736{
1737 gravi_msg_function_start(1);
1738 cpl_ensure (metrology_table, CPL_ERROR_NULL_INPUT, NULL);
1739 cpl_ensure (header, CPL_ERROR_NULL_INPUT, NULL);
1740
1741 /* Read MJD of PRC */
1742 double mjd0 = gravi_convert_to_mjd (gravi_pfits_get_start_prcacq (header));
1743
1744 /* Create the output table for VIS_MET */
1745 int ntel = 4;
1746 cpl_size nrow = cpl_table_get_nrow (metrology_table);
1747 cpl_table * vismet_table = cpl_table_new (nrow * ntel);
1748
1749 /* Create the TIME column */
1750 cpl_table_new_column (vismet_table, "TIME", CPL_TYPE_INT);
1751 cpl_table_set_column_unit (vismet_table, "TIME", "us");
1752
1753 /* Create the MJD column */
1754 cpl_table_new_column (vismet_table, "MJD", CPL_TYPE_DOUBLE);
1755 cpl_table_set_column_unit (vismet_table, "MJD", "d");
1756
1757 /* Fill the TIME and MJD column */
1758 for (cpl_size row = 0; row < nrow; row++) {
1759 int time_met = cpl_table_get_int (metrology_table, "TIME", row, NULL);
1760 double mjd_met = time_met / 86400.E6 + mjd0;
1761 for (cpl_size tel = 0; tel < ntel; tel++) {
1762 cpl_table_set_int (vismet_table, "TIME", row*ntel+tel, time_met);
1763 cpl_table_set_double (vismet_table, "MJD", row*ntel+tel, mjd_met);
1764 }
1765 }
1766
1767 /* Return */
1768 gravi_msg_function_exit(1);
1769 return vismet_table;
1770}
1771
1772/*----------------------------------------------------------------------------*/
1785/*----------------------------------------------------------------------------*/
1786
1787cpl_error_code gravi_metrology_acq (cpl_table * visacq_table,
1788 cpl_table * vismet_table,
1789 double delay,
1790 cpl_propertylist * header)
1791{
1792 gravi_msg_function_start(1);
1793 cpl_ensure_code (visacq_table, CPL_ERROR_NULL_INPUT);
1794 cpl_ensure_code (vismet_table, CPL_ERROR_NULL_INPUT);
1795 cpl_ensure_code (header, CPL_ERROR_NULL_INPUT);
1796
1797 cpl_msg_info (cpl_func,"Use acq-correction-delay = %.3f [s]",delay);
1798
1799 /* Get size */
1800 int ntel = 4;
1801 cpl_size nrow_met = cpl_table_get_nrow (vismet_table) / ntel;
1802 cpl_size nrow_acq = cpl_table_get_nrow (visacq_table) / ntel;
1803
1804 /* Create a temporary table with the time shifted
1805 * by the delta in [s] */
1806 cpl_table * visacq_tmp = cpl_table_new (nrow_acq * ntel);
1807
1808 cpl_table_duplicate_column (visacq_tmp, "TIME", visacq_table, "TIME");
1809 cpl_table_add_scalar (visacq_tmp, "TIME", 1e6 * delay);
1810
1811 /* Copy necessary data */
1812 cpl_table_duplicate_column (visacq_tmp, "OPD_PUPIL",
1813 visacq_table, "OPD_PUPIL");
1814 cpl_table_duplicate_column (visacq_tmp, "PUPIL_NSPOT",
1815 visacq_table, "PUPIL_NSPOT");
1816
1817 cpl_table_duplicate_column (visacq_tmp, "FIELD_FIBER_DX",
1818 visacq_table, "FIELD_FIBER_DX");
1819 cpl_table_duplicate_column (visacq_tmp, "FIELD_FIBER_DY",
1820 visacq_table, "FIELD_FIBER_DY");
1821
1822 cpl_table_duplicate_column (visacq_tmp, "PUPIL_U",
1823 visacq_table, "PUPIL_U");
1824 cpl_table_duplicate_column (visacq_tmp, "PUPIL_V",
1825 visacq_table, "PUPIL_V");
1826
1827 /* Get the ACQ DIT in [us] */
1828 double dit_acq = gravi_pfits_get_dit_acqcam (header) * 1e6;
1829 cpl_msg_info (cpl_func,"dit_acq = %g [us]", dit_acq);
1830
1831 /* Create SYNC information */
1832 gravi_signal_create_sync (visacq_tmp, ntel, dit_acq,
1833 vismet_table, ntel, "MET");
1834
1835 /* Create column in output table */
1836 gravi_table_new_column (vismet_table, "OPD_PUPIL", "m", CPL_TYPE_DOUBLE);
1837 double * opd_met = cpl_table_get_data_double (vismet_table, "OPD_PUPIL");
1838
1839 gravi_table_new_column (vismet_table, "FIELD_FIBER_DX", "pix", CPL_TYPE_DOUBLE);
1840 double * fdx_met = cpl_table_get_data_double (vismet_table, "FIELD_FIBER_DX");
1841
1842 gravi_table_new_column (vismet_table, "FIELD_FIBER_DY", "pix", CPL_TYPE_DOUBLE);
1843 double * fdy_met = cpl_table_get_data_double (vismet_table, "FIELD_FIBER_DY");
1844
1845 gravi_table_new_column (vismet_table, "PUPIL_U", "m", CPL_TYPE_DOUBLE);
1846 double * pupil_u = cpl_table_get_data_double (vismet_table, "PUPIL_U");
1847 gravi_table_new_column (vismet_table, "PUPIL_V", "m", CPL_TYPE_DOUBLE);
1848 double * pupil_v = cpl_table_get_data_double (vismet_table, "PUPIL_V");
1849
1850 /* Get data from input table */
1851 double * opd_acq = cpl_table_get_data_double (visacq_tmp, "OPD_PUPIL");
1852 double * fdx_acq = cpl_table_get_data_double (visacq_tmp, "FIELD_FIBER_DX");
1853 double * fdy_acq = cpl_table_get_data_double (visacq_tmp, "FIELD_FIBER_DY");
1854 int * nspot = cpl_table_get_data_int (visacq_tmp, "PUPIL_NSPOT");
1855 int * first = cpl_table_get_data_int (visacq_tmp, "FIRST_MET");
1856 int * last = cpl_table_get_data_int (visacq_tmp, "LAST_MET");
1857 double * pupil_u_acq = cpl_table_get_data_double (visacq_tmp, "PUPIL_U");
1858 double * pupil_v_acq = cpl_table_get_data_double (visacq_tmp, "PUPIL_V");
1859
1860 CPLCHECK_MSG ("Cannot load data");
1861
1862 /* Loop on beam */
1863 for (cpl_size tel = 0; tel < ntel; tel++) {
1864
1865 /* Loop on ACQ rows with undetected spot */
1866 if (gravi_pfits_get_mjd (header) < 57876.5) {
1867 cpl_msg_info (cpl_func, "Compute OPD_PUPIL for blink ACQ frames");
1868 for (cpl_size row = 1; row < nrow_acq-1; row++) {
1869 if (nspot[row*ntel+tel] == 0 &&
1870 nspot[(row-1)*ntel+tel] > 0 &&
1871 nspot[(row+1)*ntel+tel] > 0) {
1872 opd_acq[row*ntel+tel] = 0.5* (opd_acq[(row-1)*ntel+tel] + opd_acq[(row+1)*ntel+tel]);
1873 pupil_u_acq[row*ntel+tel] = 0.5* (pupil_u_acq[(row-1)*ntel+tel] + pupil_u_acq[(row+1)*ntel+tel]);
1874 pupil_v_acq[row*ntel+tel] = 0.5* (pupil_v_acq[(row-1)*ntel+tel] + pupil_v_acq[(row+1)*ntel+tel]);
1875 }
1876 }
1877 }
1878
1879 /* Loop on ACQ rows, fill the corresponding MET rows */
1880 for (cpl_size row = 0; row < nrow_acq; row++) {
1881 for (cpl_size row_met = first[row*ntel+tel]; row_met < last[row*ntel+tel]; row_met++) {
1882 opd_met[row_met*ntel+tel] = opd_acq[row*ntel+tel];
1883 fdx_met[row_met*ntel+tel] = fdx_acq[row*ntel+tel];
1884 fdy_met[row_met*ntel+tel] = fdy_acq[row*ntel+tel];
1885 pupil_u[row_met*ntel+tel] = pupil_u_acq[row*ntel+tel];
1886 pupil_v[row_met*ntel+tel] = pupil_v_acq[row*ntel+tel];
1887 }
1888 }
1889
1890 /* Loop on MET rows, to fill the empty by the closet futur value */
1891 double opd = opd_met[nrow_acq*ntel+tel];
1892 double fdx = fdx_met[nrow_acq*ntel+tel];
1893 double fdy = fdy_met[nrow_acq*ntel+tel];
1894 double pupu = pupil_u[nrow_acq*ntel+tel];
1895 double pupv = pupil_v[nrow_acq*ntel+tel];
1896
1897 for (cpl_size row = nrow_met-1; row >= 0; row--) {
1898 if (opd_met[row*ntel+tel] != 0 ) opd = opd_met[row*ntel+tel];
1899 else opd_met[row*ntel+tel] = opd;
1900 if (fdx_met[row*ntel+tel] != 0 ) fdx = fdx_met[row*ntel+tel];
1901 else fdx_met[row*ntel+tel] = fdx;
1902 if (fdy_met[row*ntel+tel] != 0 ) fdy = fdy_met[row*ntel+tel];
1903 else fdy_met[row*ntel+tel] = fdy;
1904 if (pupil_u[row*ntel+tel] != 0 ) pupu = pupil_u[row*ntel+tel];
1905 else pupil_u[row*ntel+tel] = pupu;
1906 if (pupil_v[row*ntel+tel] != 0 ) pupv = pupil_v[row*ntel+tel];
1907 else pupil_v[row*ntel+tel] = pupv;
1908 }
1909
1910 /* Loop on MET rows, to fill the empty by the closet past value */
1911 opd = opd_met[0*ntel+tel];
1912 fdx = fdx_met[0*ntel+tel];
1913 fdy = fdy_met[0*ntel+tel];
1914 pupu = pupil_u[0*ntel+tel];
1915 pupv = pupil_v[0*ntel+tel];
1916
1917 for (cpl_size row = 0; row < nrow_met; row++) {
1918 if (opd_met[row*ntel+tel] != 0) opd = opd_met[row*ntel+tel];
1919 else opd_met[row*ntel+tel] = opd;
1920 if (fdx_met[row*ntel+tel] != 0) fdx = fdx_met[row*ntel+tel];
1921 else fdx_met[row*ntel+tel] = fdx;
1922 if (fdy_met[row*ntel+tel] != 0) fdy = fdy_met[row*ntel+tel];
1923 else fdy_met[row*ntel+tel] = fdy;
1924 if (pupil_u[row*ntel+tel] != 0) pupu = pupil_u[row*ntel+tel];
1925 else pupil_u[row*ntel+tel] = pupu;
1926 if (pupil_v[row*ntel+tel] != 0) pupv = pupil_v[row*ntel+tel];
1927 else pupil_v[row*ntel+tel] = pupv;
1928 }
1929
1930 }/* End loop on beam */
1931
1932 /* Free the tmp table */
1933 FREE (cpl_table_delete, visacq_tmp);
1934
1935 /* Return */
1936 gravi_msg_function_exit(1);
1937 return CPL_ERROR_NONE;
1938}
1939
1940/*----------------------------------------------------------------------------*/
1955/*----------------------------------------------------------------------------*/
1956
1957cpl_error_code gravi_metrology_drs (cpl_table * metrology_table,
1958 cpl_table * vismet_table,
1959 cpl_propertylist * header,
1960 const cpl_parameterlist * parlist)
1961{
1962 gravi_msg_function_start(1);
1963 cpl_ensure_code (metrology_table, CPL_ERROR_NULL_INPUT);
1964 cpl_ensure_code (header, CPL_ERROR_NULL_INPUT);
1965
1966 int ntel = 4, ndiode = 4;
1967 char name[100];
1968
1969 double timer1_start=0;
1970 double timer2_start=0;
1971 double rate1=0;
1972 double rate2=0;
1973 double repeat1=0;
1974 double repeat2=0;
1975
1976 double preswitch_delay = gravi_param_get_int(parlist, "gravity.metrology.preswitch-delay");
1977 double postswitch_delay = gravi_param_get_int(parlist, "gravity.metrology.postswitch-delay");
1978
1979 if(gravi_pfits_get_met_mode(header)==MET_FAINT_HEADER)
1980 {
1981 cpl_vector * faint_params = gravi_pfits_get_met_faint_params(header);
1982
1983 rate1 = 1e6*cpl_vector_get(faint_params,0);
1984 repeat1 = cpl_vector_get(faint_params,1);
1985 rate2 = 1e6*cpl_vector_get(faint_params,3);
1986 repeat2 = cpl_vector_get(faint_params,4);
1987
1988 cpl_msg_debug(cpl_func,"ESO INS ANLO3 RATE1: %g",rate1);
1989 cpl_msg_debug(cpl_func,"ESO INS ANLO3 REPEAT1: %g",repeat1);
1990 cpl_msg_debug(cpl_func,"ESO INS ANLO3 TIMER1: %g",cpl_vector_get(faint_params,2));
1991 cpl_msg_debug(cpl_func,"ESO INS ANLO3 RATE2: %g",rate2);
1992 cpl_msg_debug(cpl_func,"ESO INS ANLO3 REPEAT2: %g",repeat2);
1993 cpl_msg_debug(cpl_func,"ESO INS ANLO3 TIMER2: %g",cpl_vector_get(faint_params,5));
1994
1995
1996
1997 /* Get the TIME of the header reference phase */
1998 /*const char * date = gravi_pfits_get_met_ph (header);*/
1999 /* Convert to UNIX time */
2000 double time_ref = 86400*(gravi_pfits_get_mjd(header) - gravi_convert_to_mjd ("1970-01-01T00:00:00"));
2001 /* Timer start times relative to metrology table time */
2002 timer1_start = 1e6*(cpl_vector_get(faint_params,2) - time_ref);
2003 timer2_start = 1e6*(cpl_vector_get(faint_params,5) - time_ref);
2004
2005
2006 cpl_msg_debug(cpl_func,"FAINT TIMER1 phase ref: %g",time_ref);
2007 cpl_msg_debug(cpl_func,"FAINT TIMER1 start: %g",timer1_start);
2008 cpl_msg_debug(cpl_func,"FAINT TIMER2 start: %g",timer2_start);
2009 /*int * time_met = cpl_table_get_data_int (metrology_table, "TIME");*/
2010
2011 cpl_vector_delete(faint_params);
2012
2013 }
2014
2015
2016 /* Parameters */
2017 int ind_sintel_FT[4][4]={{0,2,4,6},
2018 {8,10,12,14},
2019 {16,18,20,22},
2020 {24,26,28,30}};
2021 int ind_costel_FT[4][4]={{1,3,5,7},
2022 {9,11,13,15},
2023 {17,19,21,23},
2024 {25,27,29,31}};
2025 int ind_sintel_SC[4][4]={{32,34,36,38},
2026 {40,42,44,46},
2027 {48,50,52,54},
2028 {56,58,60,62}};
2029 int ind_costel_SC[4][4]={{33,35,37,39},
2030 {41,43,45,47},
2031 {49,51,53,55},
2032 {57,59,61,63}};
2033 int ind_sinfc_FT[4]={64,66,68,70};
2034 int ind_cosfc_FT[4]={65,67,69,71};
2035 int ind_sinfc_SC[4]={72,74,76,78};
2036 int ind_cosfc_SC[4]={73,75,77,79};
2037
2038 cpl_array ** raw_met=cpl_table_get_data_array(metrology_table,"VOLT");
2039 cpl_size nbrow_met = cpl_table_get_nrow (metrology_table);
2040 CPLCHECK_MSG ("Cannot load metrology data");
2041
2042 /* First: checking whether the metrology is in faint or not
2043 if metrology is in faint, read the FLAG column in the metrology table to known when it is faint or bright
2044 if metrology table does not have a FLAG column, use Felix recipy instead.
2045 Then store when the metrology is in faint into a vector
2046 Also flag as bad data the moment when the metrology changes in brightness
2047 */
2048
2049 cpl_array * bright_array = cpl_array_new(nbrow_met,CPL_TYPE_DOUBLE);
2050 cpl_array_fill_window (bright_array, 0, nbrow_met, 1);
2051
2052 cpl_array * good_met_array = cpl_array_new(nbrow_met,CPL_TYPE_DOUBLE);
2053 cpl_array_fill_window (good_met_array, 0, nbrow_met, 1);
2054
2055 int Nbefore=(int) (preswitch_delay*0.5);
2056 int Nafter=(int) (postswitch_delay*0.5);
2057
2058 if(gravi_pfits_get_met_mode(header)==MET_FAINT_HEADER)
2059 {
2060 if (cpl_table_has_column (metrology_table,"FLAG"))
2061 {
2062 cpl_msg_info(cpl_func,"Now reading the metrology flag to check the bright/faint status");
2063 for (cpl_size row = 0; row < nbrow_met; row ++)
2064 {
2065 int flag = cpl_table_get_int (metrology_table, "FLAG", row, NULL);
2066 if (flag == MET_FAINT_FLAG)
2067 cpl_array_set_double(bright_array,row,0.0);
2068 }
2069 } else {
2070 cpl_msg_warning(cpl_func,"Using theoretical information to check the bright/faint status");
2071 for (cpl_size row = 0; row < nbrow_met; row ++)
2072 {
2073 int time_met = cpl_table_get_int (metrology_table, "TIME", row, NULL);
2074 for (int repeat_counter = 0; repeat_counter < repeat1; repeat_counter ++)
2075 {
2076 if ( time_met > timer1_start+rate1*repeat_counter &&
2077 time_met < timer2_start+rate2*repeat_counter)
2078 cpl_array_set_double(bright_array,row,0.0);
2079 }
2080 }
2081 }
2082
2083 cpl_msg_info(cpl_func,"Removing part of the metrology that correspond to brigthnesse change ...");
2084 cpl_msg_info(cpl_func,"this means removing %.1fms before and %.1fms after (%i,%i DITs)",preswitch_delay,postswitch_delay,Nbefore,Nafter);
2085
2086 for (cpl_size row = 0; row < nbrow_met-1; row ++)
2087 {
2088 if (((cpl_array_get(bright_array,row, NULL) > 0.1)&&(cpl_array_get(bright_array,row+1, NULL) < 0.1))||
2089 ((cpl_array_get(bright_array,row, NULL) < 0.1)&&(cpl_array_get(bright_array,row+1, NULL) > 0.1)))
2090 {
2091 if (row-Nbefore <0)
2092 cpl_array_fill_window (good_met_array, 0, row+Nafter, 0);
2093 else if (row+Nafter >=nbrow_met)
2094 cpl_array_fill_window (good_met_array, row-Nbefore, nbrow_met-1-row+Nbefore, 0);
2095 else
2096 cpl_array_fill_window (good_met_array, row-Nbefore, Nbefore+Nafter, 0);
2097 }
2098 }
2099
2100 cpl_array_multiply (bright_array, good_met_array);
2101
2102 }
2103
2104 /*
2105 * Perform smoothing of the voltage values
2106 * First the Fiber coupler metrology
2107 * By 7 DITs (14ms)
2108 */
2109
2110 int DIT_smooth=3;
2111 int DIT_smooth_faint = DIT_smooth + Nbefore + Nafter;
2112
2113 cpl_msg_info (cpl_func,"Smoothing volts of the FC diodes by %d metrology DITS",DIT_smooth*2+1);
2114
2115 for (cpl_size diode = 64; diode < 80; diode ++) {
2116
2117 cpl_array * volt_array = cpl_array_new(nbrow_met,CPL_TYPE_DOUBLE);
2118
2119 for (cpl_size row = 0; row < nbrow_met; row ++)
2120 cpl_array_set_double(volt_array,row,cpl_array_get (raw_met[row], diode, NULL) );
2121
2122 cpl_array_multiply (volt_array, good_met_array);
2123
2124 cpl_array * volts_smooth_array=gravi_array_smooth (volt_array, DIT_smooth);
2125
2126 cpl_array * volts_smooth_array_faint=gravi_array_smooth (volt_array, DIT_smooth_faint);
2127
2128 CPLCHECK_MSG ("Cannot smooth the metrology data");
2129
2130
2131 for (cpl_size row = 0; row < nbrow_met; row ++)
2132 {
2133 if (cpl_array_get(bright_array,row, NULL) < 0.1)
2134 cpl_array_set_float(raw_met[row],diode,cpl_array_get_double (volts_smooth_array_faint, row, NULL) );
2135 else
2136 cpl_array_set_float(raw_met[row],diode,cpl_array_get_double (volts_smooth_array, row, NULL) );
2137 }
2138
2139
2140 cpl_array_delete(volt_array);
2141 cpl_array_delete(volts_smooth_array);
2142 cpl_array_delete(volts_smooth_array_faint);
2143
2144 }
2145
2146 cpl_msg_info (cpl_func,"Smoothing volts of the FC diodes done");
2147
2148 /*
2149 * Perform smoothing of the voltage values
2150 * By 41 DITS in off-axis (82ms), and by 81 DITS in on-axis (162ms)
2151 */
2152
2153 DIT_smooth=20;
2154 if ((gravi_pfits_get_axis (header) == MODE_ONAXIS)&&(gravi_pfits_get_met_mode(header)!=MET_FAINT_HEADER)) DIT_smooth=40;
2155
2156 cpl_msg_info (cpl_func,"Smoothing volts of the TEL diodes by %d metrology DITS",DIT_smooth*2+1);
2157
2158 for (cpl_size diode = 0; diode < 64; diode ++) {
2159
2160 cpl_array * volt_array = cpl_array_new(nbrow_met,CPL_TYPE_DOUBLE);
2161
2162 for (cpl_size row = 0; row < nbrow_met; row ++)
2163 cpl_array_set_double(volt_array,row,cpl_array_get (raw_met[row], diode, NULL) );
2164
2165 cpl_array * volts_smooth_array=gravi_array_smooth (volt_array, DIT_smooth);
2166
2167 CPLCHECK_MSG ("Cannot smooth the metrology data");
2168
2169 for (cpl_size row = 0; row < nbrow_met; row ++)
2170 cpl_array_set_float(raw_met[row],diode,cpl_array_get_double (volts_smooth_array, row, NULL) );
2171
2172
2173 cpl_array_delete(volt_array);
2174 cpl_array_delete(volts_smooth_array);
2175
2176 }
2177
2178 CPLCHECK_MSG ("Cannot smooth the metrology data");
2179 cpl_msg_info (cpl_func,"Smoothing volts of the TEL diodes done");
2180
2181
2182
2183 /*
2184 * Create the vismet_table
2185 */
2186
2187 cpl_msg_info (cpl_func,"Fill the OI_VIS_MET table with the DRS algorithm");
2188
2189
2190 /* Create the output table for VIS_MET
2191 * PHASE_TEL and PHASE_FC are defined as FT-SC */
2192 cpl_table_new_column (vismet_table, "PHASE_FC_DRS", CPL_TYPE_DOUBLE);
2193 cpl_table_set_column_unit (vismet_table, "PHASE_FC_DRS", "rad");
2194 cpl_table_new_column (vismet_table, "PHASE_FCFT_DRS", CPL_TYPE_DOUBLE);
2195 cpl_table_set_column_unit (vismet_table, "PHASE_FCFT_DRS", "rad");
2196 cpl_table_new_column (vismet_table, "PHASE_FCSC_DRS", CPL_TYPE_DOUBLE);
2197 cpl_table_set_column_unit (vismet_table, "PHASE_FCSC_DRS", "rad");
2198 cpl_table_new_column_array (vismet_table, "PHASE_TEL_DRS", CPL_TYPE_DOUBLE, ndiode);
2199 cpl_table_set_column_unit (vismet_table, "PHASE_TEL_DRS", "rad");
2200 cpl_table_new_column_array (vismet_table, "AMPLITUDE_TEL", CPL_TYPE_DOUBLE, ndiode);
2201 cpl_table_set_column_unit (vismet_table, "AMPLITUDE_TEL", "Volts");
2202 cpl_table_new_column (vismet_table, "FLAG_DRS", CPL_TYPE_INT);
2203 cpl_table_set_column_unit (vismet_table, "FLAG_DRS", "0/1");
2204
2205
2206
2207 /* get the reference DIT for phasing metrology */
2208
2209 const char * date = gravi_pfits_get_met_ph(header);
2210 const char * acq_date = gravi_pfits_get_start_prcacq(header);
2211 CPLCHECK_MSG ("Cannot get dates");
2212
2213 double met_mjd = 86400*1e6*(gravi_convert_to_mjd(date) -
2214 gravi_convert_to_mjd(acq_date));
2215 CPLCHECK_MSG ("Cannot convert dates");
2216
2217 /* Loop on row of the metrology */
2218 int met_date_row = -1;
2219 for (cpl_size row = 1; row < nbrow_met; row ++) {
2220 /* Look for the reference row */
2221 int time_met = cpl_table_get_int (metrology_table, "TIME", row, NULL);
2222 int time_met_minus = cpl_table_get_int (metrology_table, "TIME", row-1, NULL);
2223 CPLCHECK_MSG("Cannot get Time for Metrology phase date");
2224 if ((time_met > met_mjd) && (time_met_minus < met_mjd)) met_date_row = row;
2225 }
2226 cpl_msg_info(cpl_func,"Found DIT corresponding to metrology phase date at number %i",met_date_row);
2227 /* Store metrology DIT into header */
2228
2229 char qc_name[100];
2230 sprintf (qc_name, "ESO QC MET REF ROW");
2231 cpl_msg_info (cpl_func, "%s = %i", qc_name, met_date_row);
2232 cpl_propertylist_update_int (header, qc_name,met_date_row);
2233 cpl_propertylist_set_comment (header, qc_name, "met row at unwraped phase ref");
2234
2235 /* IF error of timing */
2236 if (met_date_row == -1){
2237 cpl_error_set_message (cpl_func, CPL_ERROR_ILLEGAL_INPUT,
2238 "The metrology phase date is not within"
2239 " the boundaries of RMN acquisition");
2240 return CPL_ERROR_ILLEGAL_INPUT;
2241 }
2242
2243 /*
2244 * Compute the metrology phase on the Fiber Coupler (FC), WITHOUT a temporal
2245 * smoothing over n_filter samples
2246 */
2247
2248 for (int tel = 0; tel < ntel; tel++) {
2249
2250 /* make an array of double complex */
2251
2252 cpl_array * phasor_array_ft = cpl_array_new(nbrow_met,CPL_TYPE_DOUBLE_COMPLEX);
2253 cpl_array * phasor_array_sc = cpl_array_new(nbrow_met,CPL_TYPE_DOUBLE_COMPLEX);
2254
2255 for (cpl_size row = 0; row < nbrow_met; row ++)
2256 {
2257 double complex phasor_ft =
2258 cpl_array_get (raw_met[row],ind_cosfc_FT[tel], NULL)
2259 + I*cpl_array_get (raw_met[row],ind_sinfc_FT[tel], NULL);
2260 double complex phasor_sc =
2261 cpl_array_get (raw_met[row],ind_cosfc_SC[tel], NULL)
2262 + I*cpl_array_get (raw_met[row],ind_sinfc_SC[tel], NULL);
2263
2264 cpl_array_set_double_complex(phasor_array_ft,row,phasor_ft);
2265 cpl_array_set_double_complex(phasor_array_sc,row,phasor_sc);
2266 }
2267 CPLCHECK_MSG("Error at computation of phasors");
2268
2269 cpl_array * phase_array_ft = cpl_array_duplicate (phasor_array_ft);
2270 cpl_array * phase_array_sc = cpl_array_duplicate (phasor_array_sc);
2271 cpl_array_arg(phase_array_ft);
2272 cpl_array_arg(phase_array_sc);
2273 gravi_array_phase_unwrap (phase_array_ft);
2274 gravi_array_phase_unwrap (phase_array_sc);
2275
2276 /* get the phase at correct time to compare with reference phase */
2277 double phase_reference_ft = cpl_array_get_double(phase_array_ft, met_date_row, NULL);
2278 double phase_reference_sc = cpl_array_get_double(phase_array_sc, met_date_row, NULL);
2279
2280 /* get the phase from rtc */
2281 sprintf (name, "ESO OCS MET PH_FC%d_FT", tel+1);
2282 double k_phase_ft = cpl_propertylist_get_double (header, name) - phase_reference_ft;
2283 sprintf (name, "ESO ADD MET PH_FC%d_FT", tel+1);
2284 if (cpl_propertylist_has(header, name))
2285 {
2286 /* if keyword in header, add 2 pi value to metrology FT -- dev use only */
2287 k_phase_ft += cpl_propertylist_get_int (header, name) * TWOPI;
2288 cpl_msg_warning( cpl_func, "Adding to FC FT metrology %i 2pi rad from header", cpl_propertylist_get_int (header, name));
2289 }
2290
2291 sprintf (name, "ESO OCS MET PH_FC%d_SC", tel+1);
2292 double k_phase_sc = cpl_propertylist_get_double (header, name) - phase_reference_sc;
2293 sprintf (name, "ESO ADD MET PH_FC%d_SC", tel+1);
2294 if (cpl_propertylist_has(header, name))
2295 {
2296 /* if keyword in header, add 2 pi value to metrology SC -- dev use only */
2297 k_phase_sc += cpl_propertylist_get_int (header, name) * TWOPI;
2298 cpl_msg_warning( cpl_func, "Adding to FC SC metrology %i 2pi rad from header", cpl_propertylist_get_int (header, name));
2299 }
2300
2301 /* get the 2 pi offset from rtc */
2302 k_phase_ft = floor(k_phase_ft/(2*M_PI)+0.5)*2*M_PI;
2303 k_phase_sc = floor(k_phase_sc/(2*M_PI)+0.5)*2*M_PI;
2304 cpl_array_add_scalar (phase_array_ft, k_phase_ft);
2305 cpl_array_add_scalar (phase_array_sc, k_phase_sc);
2306
2307 /* store data in gravi table */
2308 for (cpl_size row = 0; row < nbrow_met; row ++)
2309 {
2310 cpl_table_set (vismet_table, "PHASE_FC_DRS", row*ntel+tel, cpl_array_get_double (phase_array_ft, row, NULL)
2311 - cpl_array_get_double (phase_array_sc, row, NULL) );
2312 cpl_table_set (vismet_table, "PHASE_FCFT_DRS", row*ntel+tel, cpl_array_get_double (phase_array_ft, row, NULL));
2313 cpl_table_set (vismet_table, "PHASE_FCSC_DRS", row*ntel+tel, cpl_array_get_double (phase_array_sc, row, NULL));
2314 }
2315
2316 cpl_array_delete(phasor_array_ft);
2317 cpl_array_delete(phasor_array_sc);
2318 cpl_array_delete(phase_array_ft);
2319 cpl_array_delete(phase_array_sc);
2320 CPLCHECK_MSG("Error while computing the metrology phase at FC");
2321 }
2322 /* End loop on tel */
2323
2324 cpl_msg_info(cpl_func,"Finished computing phase (stored in PHASE_FC_DRS)");
2325
2326 /*
2327 * Compute the metrology phase at telescope, with a temporal
2328 * smoothing over twice the DIT_smoothed used on the raw signal
2329 */
2330
2331 for (int tel = 0; tel < ntel; tel++){
2332
2333 /* prepare the arrays to store the tel data */
2334 cpl_array ** tel_phase = cpl_malloc (nbrow_met * sizeof(cpl_array *)) ;
2335 for (cpl_size row = 0; row < nbrow_met; row ++)
2336 tel_phase[row]=cpl_array_new (ndiode, CPL_TYPE_DOUBLE);
2337 cpl_array ** tel_amplitude = cpl_malloc (nbrow_met * sizeof(cpl_array *)) ;
2338 for (cpl_size row = 0; row < nbrow_met; row ++)
2339 tel_amplitude[row]=cpl_array_new (ndiode, CPL_TYPE_DOUBLE);
2340
2341 for (int diode = 0; diode < ndiode; diode++){
2342
2343 /* make an array of double complex */
2344 cpl_array * phasor_array = cpl_array_new(nbrow_met,CPL_TYPE_DOUBLE_COMPLEX);
2345
2346 for (cpl_size row = 0; row < nbrow_met; row ++)
2347 {
2348 double complex phasor_ft =
2349 cpl_array_get(raw_met[row],ind_costel_FT[tel][diode], NULL)
2350 +I*cpl_array_get(raw_met[row],ind_sintel_FT[tel][diode], NULL);
2351 double complex conj_phasor_sc =
2352 cpl_array_get(raw_met[row],ind_costel_SC[tel][diode], NULL)
2353 -I*cpl_array_get(raw_met[row],ind_sintel_SC[tel][diode], NULL);
2354
2355 cpl_array_set_double_complex(phasor_array,row,phasor_ft*conj_phasor_sc);
2356 }
2357 CPLCHECK_MSG("Error at computation of phasors");
2358
2359 /* smoothing diode phase by a factor 2 with respect to before*/
2360 cpl_array * abs_array = gravi_array_smooth (phasor_array, DIT_smooth*2);
2361 cpl_array * phase_array = cpl_array_duplicate(abs_array);
2362
2363 cpl_array_arg(phase_array);
2364 gravi_array_phase_unwrap (phase_array);
2365
2366 cpl_array_abs(abs_array);
2367
2368 /* get the phase at correct time to compare with reference phase */
2369 double phase_reference = cpl_array_get_double(phase_array, met_date_row, NULL);
2370
2371 /* get the phase from rtc */
2372 sprintf (name, "ESO OCS MET PH_T%d_D%d_FT", tel+1,diode+1);
2373 double phase_rtc = cpl_propertylist_get_double (header, name);
2374 sprintf (name, "ESO OCS MET PH_T%d_D%d_SC", tel+1,diode+1);
2375 phase_rtc = phase_rtc - cpl_propertylist_get_double (header, name);
2376 double k_phase = phase_rtc - phase_reference;
2377 /* get the 2 pi offset from rtc */
2378 k_phase = floor(k_phase/(2*M_PI)+0.5)*2*M_PI;
2379 cpl_array_add_scalar (phase_array, k_phase);
2380 CPLCHECK_MSG("could not get met phase offset from headers");
2381
2382 for (cpl_size row = 0; row < nbrow_met; row ++)
2383 {
2384 cpl_array_set_double(tel_phase[row],diode,cpl_array_get_double (phase_array, row, NULL));
2385 cpl_array_set_double(tel_amplitude[row],diode,cpl_array_get_double (abs_array, row, NULL));
2386 }
2387
2388 cpl_array_delete(phasor_array);
2389 cpl_array_delete(phase_array);
2390 cpl_array_delete(abs_array);
2391 CPLCHECK_MSG("Computing the metrology phase for TEL diodes failed");
2392 } /* End loop on diodes */
2393 /* store into the table*/
2394 for (cpl_size row = 0; row < nbrow_met; row ++){
2395 cpl_table_set_array (vismet_table, "PHASE_TEL_DRS", row*ntel+tel, tel_phase[row]);
2396 cpl_table_set_array (vismet_table, "AMPLITUDE_TEL", row*ntel+tel, tel_amplitude[row]);
2397 cpl_table_set (vismet_table, "FLAG_DRS", row*ntel+tel, cpl_array_get (bright_array, row, NULL));
2398 cpl_array_delete(tel_phase[row]);
2399 cpl_array_delete(tel_amplitude[row]);
2400 }
2401 FREE (cpl_free, tel_phase);
2402 FREE (cpl_free, tel_amplitude);
2403 } /* End loop on tel */
2404
2405 cpl_msg_info(cpl_func,"Finished computing phase (stored in PHASE_TEL_DRS)");
2406
2407 /*
2408 * Compute the TEL-FC complex phasor = Vtel * conj (Vfc)
2409 * without any smoothing
2410 */
2411 cpl_msg_info (cpl_func,"Compute PHASOR_TELFC_DRS");
2412
2413 cpl_table_new_column_array (vismet_table, "PHASOR_TELFC_DRS", CPL_TYPE_DOUBLE_COMPLEX, ndiode);
2414 cpl_table_set_column_unit (vismet_table, "PHASOR_TELFC_DRS", "V^4");
2415 cpl_array ** phasor_telfc = cpl_table_get_data_array (vismet_table, "PHASOR_TELFC_DRS");
2416
2417 for (cpl_size row = 0; row < nbrow_met; row++) {
2418 for (int tel = 0; tel < ntel; tel++) {
2419 phasor_telfc[row*ntel+tel] = cpl_array_new (ndiode, CPL_TYPE_DOUBLE_COMPLEX);
2420
2421 double complex V_fc;
2422 V_fc = (cpl_array_get (raw_met[row],ind_cosfc_FT[tel], NULL) +
2423 cpl_array_get (raw_met[row],ind_sinfc_FT[tel], NULL) * I) *
2424 (cpl_array_get (raw_met[row],ind_cosfc_SC[tel], NULL) -
2425 cpl_array_get (raw_met[row],ind_sinfc_SC[tel], NULL) * I);
2426 CPLCHECK_MSG ("Cannot compute V_fc");
2427
2428 for (int diode = 0; diode < 4; diode++) {
2429 double complex V_tel;
2430 V_tel = (cpl_array_get (raw_met[row],ind_costel_FT[tel][diode], NULL) +
2431 cpl_array_get (raw_met[row],ind_sintel_FT[tel][diode], NULL) * I) *
2432 (cpl_array_get (raw_met[row],ind_costel_SC[tel][diode], NULL) -
2433 cpl_array_get (raw_met[row],ind_sintel_SC[tel][diode], NULL) * I);
2434 cpl_array_set_complex (phasor_telfc[row*ntel+tel], diode,
2435 V_tel * conj (V_fc) );
2436 CPLCHECK_MSG ("Cannot set V_tel * conj (V_fc)");
2437 } /* End loop on diode */
2438
2439 } /* End loop on tel */
2440 } /* End loop on rows */
2441
2442 cpl_array_delete(bright_array);
2443 cpl_array_delete(good_met_array);
2444
2445 gravi_msg_function_exit(1);
2446 return CPL_ERROR_NONE;
2447}
2448
2449
2450
2451/*----------------------------------------------------------------------------*/
2463/*----------------------------------------------------------------------------*/
2464cpl_error_code gravi_metrology_tac (cpl_table * metrology_table,
2465 cpl_table * vismet_table,
2466 cpl_propertylist * header)
2467{
2468 gravi_msg_function_start(1);
2469 cpl_ensure_code (metrology_table, CPL_ERROR_NULL_INPUT);
2470 cpl_ensure_code (header, CPL_ERROR_NULL_INPUT);
2471 cpl_ensure_code (vismet_table, CPL_ERROR_NULL_INPUT);
2472
2473 char card_ft[100];
2474 char card_sc[100];
2475 int ndiode = 4, ntel = 4;
2476 cpl_size nrow_met = cpl_table_get_nrow (metrology_table);
2477
2478 /*
2479 * Copy data to ensure they are of type DOUBLE
2480 * volts[nrow][ndiode]
2481 */
2482
2483 double** volts = cpl_malloc (nrow_met * sizeof(double*));
2484 cpl_array ** volts_array = cpl_table_get_data_array (metrology_table,"VOLT");
2485 for (cpl_size row = 0; row < nrow_met; row ++) {
2486 volts[row] = cpl_malloc (80 * sizeof(double));
2487 for (cpl_size diode = 0; diode < 80; diode ++) {
2488 volts[row][diode] = cpl_array_get (volts_array[row], diode, NULL);
2489 }
2490 }
2491
2492 CPLCHECK_MSG ("Cannot load metrology data");
2493
2494 /*
2495 * Allocate memory for FC the output
2496 * opd_fc[nrow*ntel]
2497 */
2498
2499 gravi_table_new_column (vismet_table, "FLAG_FC", NULL, CPL_TYPE_INT);
2500 int * flag_fc = cpl_table_get_data_int (vismet_table, "FLAG_FC");
2501
2502 /* allocate opd_fc space but not yet in vismet_table */
2503 double * opd_fc = cpl_malloc (sizeof(double) * nrow_met * ntel);
2504
2505 gravi_table_new_column (vismet_table, "PHASE_FC_TAC", "rad", CPL_TYPE_DOUBLE);
2506 double * phase_fc = cpl_table_get_data_double (vismet_table, "PHASE_FC_TAC");
2507
2508 gravi_table_new_column (vismet_table, "VAMP_FC_FT", "V", CPL_TYPE_DOUBLE);
2509 double * coher_fc_ft = cpl_table_get_data_double (vismet_table, "VAMP_FC_FT");
2510
2511 gravi_table_new_column (vismet_table, "VAMP_FC_SC", "V", CPL_TYPE_DOUBLE);
2512 double * coher_fc_sc = cpl_table_get_data_double (vismet_table, "VAMP_FC_SC");
2513
2514 /*
2515 * Allocate memory for TEL output
2516 * opd_tel[nrow*ntel][ndiode]
2517 */
2518
2519 gravi_table_new_column_array (vismet_table, "FLAG_TEL", NULL, CPL_TYPE_INT, ndiode);
2520 cpl_array ** flag_tel_array = cpl_table_get_data_array (vismet_table,"FLAG_TEL");
2521
2522 gravi_table_new_column_array (vismet_table, "PHASE_TEL_TAC", "rad", CPL_TYPE_DOUBLE, ndiode);
2523 cpl_array ** phase_tel_array = cpl_table_get_data_array (vismet_table, "PHASE_TEL_TAC");
2524
2525 gravi_table_new_column_array (vismet_table, "VAMP_TEL_FT", "V", CPL_TYPE_DOUBLE, ndiode);
2526 cpl_array ** coher_tel_ft_array = cpl_table_get_data_array (vismet_table,"VAMP_TEL_FT");
2527
2528 gravi_table_new_column_array (vismet_table, "VAMP_TEL_SC", "V", CPL_TYPE_DOUBLE, ndiode);
2529 cpl_array ** coher_tel_sc_array = cpl_table_get_data_array (vismet_table,"VAMP_TEL_SC");
2530
2531
2532 /* Wrap the newly computed phi_tel into the cpl_table. Wrapping make
2533 * the data 'valid' without having to copy them. Thus efficient */
2534 double ** opd_tel = cpl_malloc (sizeof(double*) * nrow_met * ntel);
2535 double ** phase_tel = cpl_malloc (sizeof(double*) * nrow_met * ntel);
2536 int ** flag_tel = cpl_malloc (sizeof(int*) * nrow_met * ntel);
2537 double ** coher_tel_ft = cpl_malloc (sizeof(double*) * nrow_met * ntel);
2538 double ** coher_tel_sc = cpl_malloc (sizeof(double*) * nrow_met * ntel);
2539
2540 for (cpl_size row = 0; row < nrow_met * ntel; row ++) {
2541 flag_tel[row] = cpl_malloc (sizeof(int) * ndiode);
2542 flag_tel_array[row] = cpl_array_wrap_int (flag_tel[row], ndiode);
2543
2544 opd_tel[row] = cpl_malloc (sizeof(double) * ndiode);
2545
2546 phase_tel[row] = cpl_malloc (sizeof(double) * ndiode);
2547 phase_tel_array[row] = cpl_array_wrap_double (phase_tel[row], ndiode);
2548
2549 coher_tel_ft[row] = cpl_malloc (sizeof(double) * ndiode);
2550 coher_tel_ft_array[row] = cpl_array_wrap_double (coher_tel_ft[row], ndiode);
2551
2552 coher_tel_sc[row] = cpl_malloc (sizeof(double) * ndiode);
2553 coher_tel_sc_array[row] = cpl_array_wrap_double (coher_tel_sc[row], ndiode);
2554 }
2555
2556 CPLCHECK_MSG ("Cannot allocate output memory");
2557
2558
2559
2560 /* Init the TAC algorithm -- this allocate memory */
2561 structTacData * tacData;
2562 structTacConfiguration * tacConfiguration;
2563
2564 /* get the laser wavelength data */
2565 double lambda_met_mean = gravi_pfits_get_met_wavelength_mean(header, metrology_table);
2566 cpl_msg_info (cpl_func,"Lambda met mean :%f nm", lambda_met_mean*1e9);
2567
2568 tacData = metrology_makeDefaultTacData(lambda_met_mean);
2569 tacConfiguration = tacData->tacConfiguration;
2570
2571 /* Loop on time sample to run the TAC algorithm */
2572 for(long sample_number = 1; sample_number <= nrow_met; sample_number++) {
2573
2574 /* Run the TAC algorithm */
2575 tacData->sample_number = sample_number;
2576 tacData->buffer_idx_avg = ((sample_number - 1) % tacConfiguration->number_to_average);
2577 metrology_read_voltages(tacData, volts[sample_number-1]);
2578 metrology_algorithm(tacData);
2579
2580 /* Fill the output arrays as D1T1, D2T1, D3T1, D4T1, D1T2, D2T2... */
2581 int idx_ft = 0;
2582 int idx_sc = 32;
2583 for (int idx = 0, tel = 0; tel < ntel; tel++) {
2584 cpl_size nmet = (sample_number-1)*ntel + tel;
2585
2586 for (int diode = 0; diode < ndiode; diode++) {
2587 /* TAC computation */
2588 opd_tel[nmet][diode] = tacData->opl_telescope_diode[tel][diode];
2589 phase_tel[nmet][diode] = - tacData->opl_telescope_diode[tel][diode] * CPL_MATH_2PI / lambda_met_mean;
2590 flag_tel[nmet][diode] = tacData->total_flag_telescope[tel][diode][FT] | tacData->total_flag_telescope[tel][diode][SC];
2591
2592 /* Volt amplitude */
2593 coher_tel_ft[nmet][diode] = sqrt (volts[sample_number-1][idx_ft+idx*2] * volts[sample_number-1][idx_ft+idx*2] +
2594 volts[sample_number-1][idx_ft+idx*2+1] * volts[sample_number-1][idx_ft+idx*2+1]);
2595 coher_tel_sc[nmet][diode] = sqrt (volts[sample_number-1][idx_sc+idx*2] * volts[sample_number-1][idx_sc+idx*2] +
2596 volts[sample_number-1][idx_sc+idx*2+1] * volts[sample_number-1][idx_sc+idx*2+1]);
2597 idx++;
2598 }
2599 }
2600
2601 /* Fill the output arrays as T1, T2, T3, T4... */
2602 idx_ft = 64;
2603 idx_sc = 72;
2604 for (int idx = 0, tel = 0; tel < ntel; tel++) {
2605 cpl_size nmet = (sample_number-1)*ntel + tel;
2606
2607 /* TAC computation */
2608 opd_fc[nmet] = tacData->opl_fiber_coupler[tel];
2609 phase_fc[nmet] = - tacData->opl_fiber_coupler[tel] * CPL_MATH_2PI / lambda_met_mean;
2610 flag_fc[nmet] = tacData->total_flag_fiber_coupler[tel][FT] | tacData->total_flag_fiber_coupler[tel][SC];
2611 /* Volt amplitude */
2612 coher_fc_ft[nmet] = sqrt (volts[sample_number-1][idx_ft+idx*2] * volts[sample_number-1][idx_ft+idx*2] +
2613 volts[sample_number-1][idx_ft+idx*2+1] * volts[sample_number-1][idx_ft+idx*2+1]);
2614 coher_fc_sc[nmet] = sqrt (volts[sample_number-1][idx_sc+idx*2] * volts[sample_number-1][idx_sc+idx*2] +
2615 volts[sample_number-1][idx_sc+idx*2+1] * volts[sample_number-1][idx_sc+idx*2+1]);
2616 idx++;
2617 }
2618
2619 }
2620 /* End loop on time sample */
2621
2622
2623 /* Get the TIME of the header reference phase */
2624 const char * date = gravi_pfits_get_met_ph (header);
2625 const char * acq_date = gravi_pfits_get_start_prcacq (header);
2626 double time_ref = 86400*1e6*(gravi_convert_to_mjd (date) - gravi_convert_to_mjd (acq_date));
2627 int * time_met = cpl_table_get_data_int (metrology_table, "TIME");
2628
2629 cpl_size row_ref = -1;
2630 for (cpl_size row = 1; row < nrow_met; row++)
2631 if (time_met[row]>time_ref && time_met[row-1]<time_ref) {
2632 row_ref = row; break;
2633 }
2634
2635 /* Check if we have found the metrology reference date inside the file */
2636 if (row_ref == -1) {
2637 cpl_error_set_message(cpl_func, CPL_ERROR_ILLEGAL_INPUT,
2638 "The metrology phase date is not within"
2639 " the boundaries of RMN acquisition");
2640 return CPL_ERROR_NULL_INPUT;
2641 } else {
2642 cpl_msg_info (cpl_func,"The metrology phase date is OK (row %lld over %lld).", row_ref, nrow_met);
2643 }
2644
2645 /* Re-install the absolute phase for the FC diodes */
2646 for (int tel = 0; tel < ntel; tel++) {
2647 sprintf (card_ft, "ESO OCS MET PH_FC%d_FT", tel+1);
2648 sprintf (card_sc, "ESO OCS MET PH_FC%d_SC", tel+1);
2649 double opd_ref = lambda_met_mean / CPL_MATH_2PI *
2650 (cpl_propertylist_get_double (header, card_sc) - cpl_propertylist_get_double (header, card_ft)) -
2651 opd_fc[row_ref*ntel+tel];
2652 double opd_ref_int = gravi_round (opd_ref / lambda_met_mean) * lambda_met_mean;
2653 double phase_ref_int = gravi_round (opd_ref / lambda_met_mean) * CPL_MATH_2PI;
2654 for (cpl_size row = 0; row < nrow_met; row++) {
2655 opd_fc[row*ntel+tel] += opd_ref_int;
2656 phase_fc[row*ntel+tel] -= phase_ref_int;
2657 }
2658 }
2659 /* End loop on tel and diodes */
2660
2661 /* Re-install the absolute phase for the TEL diodes */
2662 for (int tel = 0; tel < ntel; tel++) {
2663 for (int diode = 0; diode < ndiode; diode++) {
2664 sprintf (card_ft,"ESO OCS MET PH_T%d_D%d_FT", tel+1, diode+1);
2665 sprintf (card_sc,"ESO OCS MET PH_T%d_D%d_SC", tel+1, diode+1);
2666 double opd_ref = lambda_met_mean / CPL_MATH_2PI *
2667 (cpl_propertylist_get_double (header, card_sc) - cpl_propertylist_get_double (header, card_ft)) -
2668 opd_tel[row_ref*ntel+tel][diode];
2669 double opd_ref_int = gravi_round (opd_ref / lambda_met_mean) * lambda_met_mean;
2670 double phase_ref_int = gravi_round (opd_ref / lambda_met_mean) * CPL_MATH_2PI;
2671 for (cpl_size row = 0; row < nrow_met; row++) {
2672 opd_tel[row*ntel+tel][diode] += opd_ref_int;
2673 phase_tel[row*ntel+tel][diode] -= phase_ref_int;
2674 }
2675 }
2676 }
2677
2678 CPLCHECK_MSG ("Cannot fill result of metrology computation");
2679
2680 /* Free the pointer to pointer to data */
2681 FREELOOP (cpl_free, volts, nrow_met);
2682 FREELOOP (cpl_free, opd_tel, nrow_met*ntel);
2683 FREE (cpl_free, opd_fc);
2684 FREE (cpl_free, phase_tel);
2685 FREE (cpl_free, flag_tel);
2686 FREE (cpl_free, coher_tel_ft);
2687 FREE (cpl_free, coher_tel_sc);
2688
2689 /* Free the TAC data */
2690 FREE (cpl_free, tacConfiguration);
2691 FREE (cpl_free, tacData);
2692
2693 gravi_msg_function_exit(1);
2694 return CPL_ERROR_NONE;
2695}
2696
2697
2698
2699/*----------------------------------------------------------------------------*/
2713/*----------------------------------------------------------------------------*/
2714cpl_error_code gravi_metrology_telfc (cpl_table * metrology_table,
2715 cpl_table * vismet_table,
2716 gravi_data * static_param_data,
2717 cpl_propertylist * header,
2718 const cpl_parameterlist * parlist)
2719{
2720 gravi_msg_function_start(1);
2721
2722 /* get size of arrays */
2723 cpl_size ndiode = 4, ntel = 4;
2724 cpl_size nrow_met = cpl_table_get_nrow (metrology_table);
2725 char qc_name[100];
2726
2727 /* get the options */
2728 int use_faint_met = gravi_param_get_bool (parlist, "gravity.metrology.use-faint-met");
2729 int use_fiber_dxy = gravi_param_get_bool (parlist, "gravity.metrology.use-fiber-dxy");
2730 int use_met_rtc = gravi_param_get_bool (parlist, "gravity.metrology.use-met-rtc");
2731
2732 /* loading phase data
2733 * The choice is yours here. Either you use PHASE_FC_DRS and PHASE_TEL_DRS
2734 * Or you can use PHASE_FC_TAC and PHASE_TEL_TAC
2735 * DSR version is prefered by me
2736 */
2737 double * phase_fc;
2738 double * phase_fc_2;
2739 double ** phase_tel;
2740 double ** amplitude_tel;
2741 int * met_flag;
2742
2743 if (use_met_rtc == 0) {
2744 cpl_msg_info (cpl_func,"Using DRS metrology algorithm");
2745
2746 phase_fc = cpl_table_get_data_double (vismet_table, "PHASE_FC_DRS");
2747 CPLCHECK_MSG ("Cannot get PHASE_FC_DRS from vismet_table");
2748
2749 phase_fc_2 = cpl_table_get_data_double (vismet_table, "PHASE_FC_TAC");
2750 CPLCHECK_MSG ("Cannot get PHASE_FC_TAC from vismet_table");
2751
2752 phase_tel = gravi_table_get_data_array_double (vismet_table, "PHASE_TEL_DRS");
2753 CPLCHECK_MSG ("Cannot get PHASE_TEL_DRS from vismet_table");
2754 }
2755 else{
2756 cpl_msg_info (cpl_func,"Using RTC (TAC) metrology algorithm");
2757
2758 phase_fc = cpl_table_get_data_double (vismet_table, "PHASE_FC_TAC");
2759 CPLCHECK_MSG ("Cannot get PHASE_FC_DRS from vismet_table");
2760
2761 phase_fc_2 = cpl_table_get_data_double (vismet_table, "PHASE_FC_DRS");
2762 CPLCHECK_MSG ("Cannot get PHASE_FC_TAC from vismet_table");
2763
2764 phase_tel = gravi_table_get_data_array_double (vismet_table, "PHASE_TEL_TAC");
2765 CPLCHECK_MSG ("Cannot get PHASE_TEL_DRS from vismet_table");
2766 }
2767
2768 amplitude_tel = gravi_table_get_data_array_double (vismet_table, "AMPLITUDE_TEL");
2769 CPLCHECK_MSG ("Cannot get AMPLITUDE_TEL from vismet_table");
2770
2771 cpl_msg_info (cpl_func,"G-FAINT: create flagged amplitude");
2772 met_flag = cpl_table_get_data_int (vismet_table, "FLAG_DRS");
2773 CPLCHECK_MSG ("Cannot get FLAG_DRS from vismet_table");
2774
2775 gravi_table_init_column_array (vismet_table, "AMPLITUDE_TEL_FLAG", "V", CPL_TYPE_DOUBLE, ndiode);
2776 double ** amplitude_tel_flagged = gravi_table_get_data_array_double (vismet_table, "AMPLITUDE_TEL_FLAG");
2777 for (cpl_size row = 0; row < nrow_met*ntel; row++) {
2778 for (int diode = 0; diode < ndiode; diode++)
2779 amplitude_tel_flagged[row][diode]= amplitude_tel[row][diode] * met_flag[row];
2780 }
2781 CPLCHECK_MSG ("Cannot create flagged amplitude");
2782
2783 /* get the laser wavelength data */
2784 double lambda_met_mean = gravi_pfits_get_met_wavelength_mean(header, metrology_table);
2785 cpl_msg_info (cpl_func,"Lambda met mean :%f nm", lambda_met_mean*1e9);
2786 CPLCHECK_MSG ("Cannot get laser wavelength data");
2787
2788 /* Calculate difference in phase between TAC and DRS */
2789
2790 for (int tel = 0; tel < ntel; tel++) {
2791 double phase_diff = 0;
2792 for (cpl_size row = 0; row < nrow_met; row++)
2793 phase_diff += phase_fc[row*ntel+tel] - phase_fc_2[row*ntel+tel];
2794 phase_diff /= nrow_met;
2795 sprintf (qc_name, "ESO QC MET DRS DIFF%i", tel+1);
2796 if (phase_diff<0.02)
2797 cpl_msg_info (cpl_func, "%s = %f", qc_name, phase_diff);
2798 else
2799 cpl_msg_warning (cpl_func, "%s = %f", qc_name, phase_diff);
2800 cpl_propertylist_update_double (header, qc_name,phase_diff);
2801 cpl_propertylist_set_comment (header, qc_name, "[rad] TAC and DRS phase difference");
2802 }
2803
2804
2805 /************************************************/
2806 /* PART O: OPD_TEL and OPD_FC */
2807 /************************************************/
2808 /* Populate the OPD_TEL and OPD_FC columns */
2809 /* this is important to do it here because there should be an */
2810 /* agreement between these columns and the TAC or DRS choice above */
2811 /* also with the exact metrology value used */
2812 /************************************************/
2813
2814 gravi_table_new_column (vismet_table, "OPD_FC", "m", CPL_TYPE_DOUBLE);
2815 double * opd_fc = cpl_table_get_data_double (vismet_table, "OPD_FC");
2816
2817 gravi_table_init_column_array (vismet_table, "OPD_TEL", "m", CPL_TYPE_DOUBLE, ndiode);
2818 double ** opd_tel = gravi_table_get_data_array_double (vismet_table, "OPD_TEL");
2819
2820 for (cpl_size row = 0; row < nrow_met*ntel; row++) {
2821 opd_fc[row]= - lambda_met_mean / CPL_MATH_2PI * phase_fc[row];
2822 for (int diode = 0; diode < ndiode; diode++)
2823 opd_tel[row][diode]= - lambda_met_mean / CPL_MATH_2PI * phase_tel[row][diode];
2824 }
2825
2826
2827 /************************************************/
2828 /* PART I: OPD_FC_CORR */
2829 /************************************************/
2830 /* Correction to Fiber Coupler Metrology OPD */
2831 /* from pupil motion, fiber pickup, and defocus */
2832 /************************************************/
2833
2834
2835 cpl_msg_info (cpl_func,"Calculate OPD_FC_CORR from OPD_PUPIL and pickup/defocus.");
2836
2837 /* Create array in OI_VIS_MET table, fill with zeros, and get pointer */
2838 gravi_table_new_column (vismet_table, "OPD_FC_CORR", "m", CPL_TYPE_DOUBLE);
2839 double * opd_fc_corr = cpl_table_get_data_double (vismet_table, "OPD_FC_CORR");
2840
2841 /* Pupil motion correction already calculated before and available from OPD_PUPIL */
2842 double * opd_pupil = NULL;
2843 if ( cpl_table_has_column(vismet_table, "OPD_PUPIL") ) {
2844 opd_pupil = cpl_table_get_data_double (vismet_table, "OPD_PUPIL");
2845 }
2846 else {
2847 cpl_msg_warning(cpl_func,"Cannot get the OPD_PUPIL (not computed) so will"
2848 "not correct for pupil opd (check the --reduce-acq-cam option)");
2849 }
2850
2851 /* Retrieve object separation */
2852 double dx_in = gravi_pfits_get_sobj_x (header);
2853 double dy_in = gravi_pfits_get_sobj_y (header);
2854 double rho_in = sqrt(dx_in*dx_in + dy_in*dy_in);
2855 /* Also retreive the seperation from gvctu */
2856 double dx_gvctu = gravi_pfits_get_gvctu_x (header);
2857 double dy_gvctu = gravi_pfits_get_gvctu_y (header);
2858
2859 CPLCHECK_MSG ("Cannot get separation");
2860 cpl_msg_info(cpl_func,"X and Y OFFSET from gvctu = %.2f, %.2f mas",dx_gvctu,dy_gvctu);
2861 cpl_msg_info(cpl_func,"X and Y OFFSET from header = %.2f, %.2f mas",dx_in,dy_in);
2862
2863 /* Force separation to zero in SINGLE */
2864 if (gravi_pfits_get_mode (header) == MODE_SINGLE) {
2865 cpl_msg_info (cpl_func,"Mode SINGLE thus separation forced to 0.0");
2866 rho_in = 0.;
2867 }
2868 cpl_msg_info (cpl_func,"FE: separation in mas: %g ", rho_in );
2869
2870 /* Apply pupil, focus and pickup offsets */
2871 for (int tel = 0; tel < ntel; tel++) {
2872
2873 double opd_fc_focus = gravi_metrology_get_fc_focus (header, tel, static_param_data); // [m]
2874 double opd_fc_shift = gravi_metrology_get_fc_shift (header, tel, static_param_data) * (rho_in/1000.0); // [m/arcsec] x [arcsec] = [m]
2875 CPLCHECK_INT("Cannot get Static calib");
2876
2877 cpl_msg_info (cpl_func, "FE: Tel %d Pupil %g Focus %g Shift %g",
2878 tel,
2879 (opd_pupil?opd_pupil[0*ntel+tel]:0)*1e9,
2880 opd_fc_focus*1e9,
2881 opd_fc_shift*1e9 );
2882
2883 for (cpl_size row = 0; row < nrow_met; row++) {
2884 opd_fc_corr[row*ntel+tel] = (opd_pupil?opd_pupil[row*ntel+tel]:0)
2885 + opd_fc_focus
2886 + opd_fc_shift;
2887 }
2888 }
2889 CPLCHECK_MSG ("Cannot calculate OPD_FC_CORR");
2890
2891
2892 /*****************************************************************/
2893 /* PART II: OPD_TEL_CORR */
2894 /*****************************************************************/
2895 /* Deprojection for telescope diodes to center of telescope, */
2896 /* based on prior knowledge of object separation, */
2897 /* acqcam pointing offset, and astigmatism. */
2898 /*****************************************************************/
2899
2900
2901 cpl_msg_info (cpl_func,"Deproject telescope diodes to center of telescope in OPD_TEL_CORR.");
2902
2903 /* Create array in OI_VIS_MET table, fill with zeros, and get list of pointer.
2904 * Note that this list has to be free */
2905 gravi_table_init_column_array (vismet_table, "OPD_TEL_CORR", "m", CPL_TYPE_DOUBLE, ndiode);
2906 double ** opd_tel_corr = gravi_table_get_data_array_double (vismet_table, "OPD_TEL_CORR");
2907 CPLCHECK_MSG ("Cannot init the OPD_TEL_CORR column");
2908
2909 /* Read metrology receiver positions from MET_POS,
2910 * rec_az = [HIERARCH ESO MET AT<tel> REC<diode>X]
2911 * rec_zd = [HIERARCH ESO MET AT<tel> REC<diode>Y]
2912 * 1st index = beam/tel, 2nd index = diode */
2913 double rec_az[4][4];
2914 double rec_zd[4][4];
2915 for (int tel=0;tel<4;tel++) {
2916 for (int diode=0;diode<4;diode++) {
2917 rec_az[tel][diode] = gravi_metrology_get_posx (header, tel, diode); /* in [m] in order GV 1,2,3,4 */
2918 rec_zd[tel][diode] = gravi_metrology_get_posy (header, tel, diode); /* in [m] in order GV 1,2,3,4 */
2919 }
2920 }
2921 CPLCHECK_MSG ("Cannot read metrology reciever position from MET_POS");
2922
2923 /*----- PART II.a: Separation deprojection with acqcam correction -----*/
2924
2925 /* (rec_az E_AZ + rec_zd E_ZD) . (sobj_U E_U + sobj_V E_V) */
2926
2927 /* Declare some variables */
2928 double deproject, northangle, field_dU, field_dV;
2929 char card[100];
2930
2931 /* Load North angle array */
2932 cpl_array * northangle_array = cpl_array_new(4,CPL_TYPE_DOUBLE);
2933
2934 /* get North angle array, unwrapped */
2935 double north_angle_tmp = gravi_pfits_get_northangle_acqcam (header, 0);
2936 cpl_array_set_double (northangle_array,0,north_angle_tmp);
2937 CPLCHECK_MSG ("Cannot get north angle");
2938
2939 for (cpl_size tel=1;tel < 4; tel++)
2940 {
2941 double north_angle_tmp2 = gravi_pfits_get_northangle_acqcam (header, tel);
2942 if (north_angle_tmp2-north_angle_tmp >= 180) north_angle_tmp2 -= 360.0;
2943 if (north_angle_tmp2-north_angle_tmp < -180) north_angle_tmp2 += 360.0;
2944 cpl_array_set_double (northangle_array,tel,north_angle_tmp2);
2945 CPLCHECK_MSG ("Cannot get north angle for tel");
2946 }
2947
2948 /* Vectors used in Julien's formula */
2949 cpl_vector * rec = cpl_vector_new (3);
2950 cpl_vector * sobj = cpl_vector_new (3);
2951
2952 /* read E_U,V,AZ,ZD from OI_VIS_MET table */
2953 cpl_array ** E_U = cpl_table_get_data_array (vismet_table,"E_U");
2954 cpl_array ** E_V = cpl_table_get_data_array (vismet_table,"E_V");
2955 cpl_array ** E_AZ = cpl_table_get_data_array (vismet_table,"E_AZ");
2956 cpl_array ** E_ZD = cpl_table_get_data_array (vismet_table,"E_ZD");
2957 CPLCHECK_MSG ("Cannot get E_U, E_V, E_AZ, E_ZD");
2958
2959 /* read the field fiber offset */
2960 double * field_dX = NULL;
2961 double * field_dY = NULL;
2962 if ( cpl_table_has_column (vismet_table, "FIELD_FIBER_DX")
2963 && cpl_table_has_column (vismet_table, "FIELD_FIBER_DY") ) {
2964 field_dX = cpl_table_get_data_double (vismet_table, "FIELD_FIBER_DX");
2965 field_dY = cpl_table_get_data_double (vismet_table, "FIELD_FIBER_DY");
2966 }
2967 CPLCHECK_MSG ("Cannot get FIELD_FIBER_DX/DY");
2968
2969 /* some debug messages */
2970 cpl_msg_info (cpl_func,"FE: E_U = [%g, %g, %g].",
2971 cpl_array_get (E_U[0*ntel+0], 0, NULL),
2972 cpl_array_get (E_U[0*ntel+0], 1, NULL),
2973 cpl_array_get (E_U[0*ntel+0], 2, NULL));
2974 cpl_msg_info (cpl_func,"FE: E_V = [%g, %g, %g].",
2975 cpl_array_get (E_V[0*ntel+0], 0, NULL),
2976 cpl_array_get (E_V[0*ntel+0], 1, NULL),
2977 cpl_array_get (E_V[0*ntel+0], 2, NULL));
2978 cpl_msg_info (cpl_func,"FE: E_AZ = [%g, %g, %g].",
2979 cpl_array_get (E_AZ[0*ntel+0], 0, NULL),
2980 cpl_array_get (E_AZ[0*ntel+0], 1, NULL),
2981 cpl_array_get (E_AZ[0*ntel+0], 2, NULL));
2982 cpl_msg_info (cpl_func,"FE: E_ZD = [%g, %g, %g].",
2983 cpl_array_get (E_ZD[0*ntel+0], 0, NULL),
2984 cpl_array_get (E_ZD[0*ntel+0], 1, NULL),
2985 cpl_array_get (E_ZD[0*ntel+0], 2, NULL));
2986
2987 /* use reference value for X and Y fiber offset */
2988 /* if separation tracking is on , it is better to use X and Y position */
2989 /* if separation tracking is off, it is better to use the gvctu estimate */
2990 double dx_sep_diode = dx_gvctu;
2991 double dy_sep_diode = dy_gvctu;
2992 if (cpl_propertylist_has (header, "ESO FT KAL SEPTRK"))
2993 {
2994 if (cpl_propertylist_get_int (header, "ESO FT KAL SEPTRK") == 1)
2995 {
2996 dx_sep_diode = dx_in;
2997 dy_sep_diode = dy_in;
2998 cpl_msg_info (cpl_func,"Using X/Y offset from header");
2999 } else {
3000 cpl_msg_info (cpl_func,"Using X/Y offset from gvctu");
3001 }
3002 } else {
3003 cpl_msg_info (cpl_func,"Using X/Y offset from gvctu");
3004 }
3005 CPLCHECK_MSG ("Cannot get ESO FT KAL SEPTRK");
3006
3007 /* Verbose about option before the loop */
3008 cpl_msg_info (cpl_func,"Use fiber dxy is set to %s", use_fiber_dxy ? "TRUE" : "FALSE");
3009
3010 /* loop over all column and diodes */
3011 for (int tel = 0; tel < ntel; tel++) {
3012
3013 double RA=0.0;
3014 double RZ=0.0;
3015 for (int diode = 0; diode < ndiode; diode++) {
3016 RA+=rec_az[tel][diode];
3017 RZ+=rec_zd[tel][diode];
3018 }
3019 cpl_msg_info (cpl_func,"Tel %i; diodes RA=%g, RZ=%g, ABS=%g", tel, RA, RZ, sqrt(RA*RA+RZ*RZ));
3020
3021 /* compute the north angle on acqcam [deg] */
3022 northangle = cpl_array_get_double(northangle_array,tel,NULL);
3023 CPLCHECK_MSG ("Cannot get North Angle");
3024
3025 /* If available, get average image scale on acqcam [mas/pix] */
3026 double scale = 0.0;
3027 sprintf (card,"ESO QC ACQ FIELD%d SCALE", tel+1);
3028 if (cpl_propertylist_has (header, card))
3029 scale = cpl_propertylist_get_double (header, card);
3030 CPLCHECK_MSG ("Cannot get ACQ FIELD SCALE");
3031
3032 for (cpl_size row = 0; row < nrow_met; row++) {
3033
3034 /* If available, transform field fiber offset from (x,y) acqcam [pix] to (U,V) sky [mas].
3035 FE 20190504: add the parameter use_fiber_dxy to ignore this correction */
3036 if ((field_dX && field_dY) && use_fiber_dxy != 0) {
3037 field_dU = (field_dX[row*ntel+tel] * sin( (northangle+90.) * CPL_MATH_RAD_DEG )
3038 +field_dY[row*ntel+tel] * cos( (northangle+90.) * CPL_MATH_RAD_DEG ))*scale;
3039 field_dV = (field_dX[row*ntel+tel] * sin( (northangle ) * CPL_MATH_RAD_DEG )
3040 +field_dY[row*ntel+tel] * cos( (northangle ) * CPL_MATH_RAD_DEG ))*scale;
3041 } else {
3042 field_dU = 0.0;
3043 field_dV = 0.0;
3044 }
3045
3046 for (int diode = 0; diode < ndiode; diode++) {
3047
3048 /* Filling vectors of Julien's formula */
3049 cpl_vector_set (rec, 0, rec_az[tel][diode] * cpl_array_get (E_AZ[row*ntel+tel], 0, NULL)
3050 +rec_zd[tel][diode] * cpl_array_get (E_ZD[row*ntel+tel], 0, NULL));
3051 cpl_vector_set (rec, 1, rec_az[tel][diode] * cpl_array_get (E_AZ[row*ntel+tel], 1, NULL)
3052 +rec_zd[tel][diode] * cpl_array_get (E_ZD[row*ntel+tel], 1, NULL));
3053 cpl_vector_set (rec, 2, rec_az[tel][diode] * cpl_array_get (E_AZ[row*ntel+tel], 2, NULL)
3054 +rec_zd[tel][diode] * cpl_array_get (E_ZD[row*ntel+tel], 2, NULL));
3055 cpl_vector_set (sobj, 0, (dx_sep_diode-field_dU) * cpl_array_get (E_U[row*ntel+tel], 0, NULL)
3056 +(dy_sep_diode-field_dV) * cpl_array_get (E_V[row*ntel+tel], 0, NULL));
3057 cpl_vector_set (sobj, 1, (dx_sep_diode-field_dU) * cpl_array_get (E_U[row*ntel+tel], 1, NULL)
3058 +(dy_sep_diode-field_dV) * cpl_array_get (E_V[row*ntel+tel], 1, NULL));
3059 cpl_vector_set (sobj, 2, (dx_sep_diode-field_dU) * cpl_array_get (E_U[row*ntel+tel], 2, NULL)
3060 +(dy_sep_diode-field_dV) * cpl_array_get (E_V[row*ntel+tel], 2, NULL));
3061
3062 /* calculate deprojection */
3063 deproject = cpl_vector_product(rec, sobj); /* in m * mas */
3064 deproject = deproject / 1000. / 3600. / 360. * TWOPI; /* convert in meter */
3065
3066 /* some debug messages */
3067 if (row == 0 && tel == 0 && diode == 0) {
3068 cpl_msg_debug (cpl_func,"FE: Julien deproject diode 0 in nm: %g", deproject*1e9);
3069 }
3070 if (row == 0 && tel == 0 && diode == 1) {
3071 cpl_msg_debug (cpl_func,"FE: Julien deproject diode 1 in nm: %g", deproject*1e9);
3072 }
3073 if (row == 0 && tel == 0 && diode == 2) {
3074 cpl_msg_debug (cpl_func,"FE: Julien deproject diode 2 in nm: %g", deproject*1e9);
3075 }
3076 if (row == 0 && tel == 0 && diode == 3) {
3077 cpl_msg_debug (cpl_func,"FE: Julien deproject diode 3 in nm: %g", deproject*1e9);
3078 }
3079
3080 /* store deprojection in opd_tel_corr */
3081 opd_tel_corr[row*ntel+tel][diode] = deproject;
3082 }
3083 }
3084 }
3085
3086 /* Free memory */
3087 FREE (cpl_vector_delete, rec);
3088 FREE (cpl_vector_delete, sobj);
3089 CPLCHECK_MSG ("Cannot calculate OPD_TEL_CORR (SX and SY)");
3090
3091 /*----- Part II.b: Astigmatism -----*/
3092
3093 /* local variables */
3094 double AstigmAmplitude; // [nm]
3095 double AstigmTheta; // [rad]
3096 double AstigmRadius; // [m]
3097 double astigm;
3098 double diodeang;
3099 double astang;
3100 double astradius;
3101
3102 /* Paralactic angle is averaged from fitsheader */
3103
3104 double parang1 = cpl_propertylist_get_double(header, "ESO ISS PARANG START") * (TWOPI/360.0); // in [rad]
3105 double parang2 = cpl_propertylist_get_double(header, "ESO ISS PARANG END" ) * (TWOPI/360.0); // in [rad]
3106 if (parang2-parang1 > TWOPI/2) parang2-=TWOPI;
3107 if (parang2-parang1 < -TWOPI/2) parang2+=TWOPI;
3108
3109 double parang = (parang1 + parang2) / 2;
3110 CPLCHECK_MSG ("Cannot get paralactic angle");
3111 cpl_msg_info (cpl_func,"FE: paralactic angle in degrees: %g ", parang * (360.0/TWOPI) );
3112
3113 /* Posangle is calculated from SOBJX ansd SOBJY already read before
3114 x,y are exchanged following coordinate systems in Stefan's slide */
3115 int flag;
3116
3117 /* loop over all diodes and beams */
3118 for (int tel = 0; tel < ntel; tel++) {
3119
3120 gravi_metrology_get_astig (header, tel, &AstigmAmplitude, &AstigmTheta, &AstigmRadius, static_param_data);
3121
3122 cpl_msg_info (cpl_func,"Astigmatism params GV%i : %5g [nm], %5g [deg]", tel, AstigmAmplitude ,AstigmTheta);
3123
3124 for (cpl_size row = 0; row < nrow_met; row++) {
3125 for (int diode = 0; diode < ndiode; diode++) {
3126 diodeang = myAtan(-rec_zd[tel][diode],-rec_az[tel][diode], &flag); // [rad]
3127 double parang3 = parang1 + (parang2-parang1) * row / nrow_met; // [rad]
3128 astang = (AstigmTheta - cpl_array_get_double(northangle_array,tel,NULL)) * (TWOPI/360.) - parang3 - diodeang ; // [rad]
3129 astradius = sqrt(rec_az[tel][diode]*rec_az[tel][diode] + rec_zd[tel][diode]*rec_zd[tel][diode]) / AstigmRadius; /* normalized */
3130 astigm = (AstigmAmplitude*1e-9) * sqrt(6) * astradius * astradius * sin(2. * astang); // [m]
3131
3132 /* some debug messages */
3133 if (row == 0 && tel == 0 && diode == 0) {
3134 cpl_msg_debug (cpl_func,"FE: Frank diode angle [deg]: %g", diodeang / TWOPI * 360.);
3135 cpl_msg_debug (cpl_func,"FE: Frank astigmatism angle [deg]: %g", astang / TWOPI * 360.);
3136 cpl_msg_debug (cpl_func,"FE: Frank normalized astradius: %g", astradius);
3137 cpl_msg_debug (cpl_func,"FE: Frank astigmatism diode 0 in nm: %g", astigm*1e9);
3138 }
3139 if (row == 0 && tel == 0 && diode == 1) {
3140 cpl_msg_debug (cpl_func,"FE: Frank astigmatism diode 1 in nm: %g", astigm*1e9);
3141 }
3142 if (row == 0 && tel == 0 && diode == 2) {
3143 cpl_msg_debug (cpl_func,"FE: Frank astigmatism diode 2 in nm: %g", astigm*1e9);
3144 }
3145 if (row == 0 && tel == 0 && diode == 3) {
3146 cpl_msg_debug (cpl_func,"FE: Frank astigmatism diode 3 in nm: %g", astigm*1e9);
3147 }
3148
3149
3150 /* subtracting astigmatism */
3151 opd_tel_corr[row*ntel+tel][diode] -= astigm;
3152 }
3153 }
3154 }
3155 CPLCHECK_MSG ("Cannot calculate OPD_TEL_CORR (Astigmatism part)");
3156
3157 /*****************************************************************/
3158 /* PART III: PHASE_TELFC_CORR */
3159 /*****************************************************************/
3160 /* Calculation of the difference between the best correction of */
3161 /* the metrology receivers and the best correction of the fiber */
3162 /* coupler. This quantity should be close to zero. */
3163 /*****************************************************************/
3164
3165 cpl_msg_info (cpl_func,"Calculate difference between corrected telescope diodes and fiber coupler in OPD_TELFC_CORR.");
3166
3167 /* Create OPD_TELFC_CORR array in OI_VIS_MET table, fill with zeros, and get pointer */
3168 gravi_table_init_column_array (vismet_table, "OPD_TELFC_CORR", "m", CPL_TYPE_DOUBLE, ndiode);
3169 double ** opd_telfc_corr = gravi_table_get_data_array_double (vismet_table, "OPD_TELFC_CORR");
3170
3171 gravi_table_init_column_array (vismet_table, "PHASE_TELFC_CORR", "rad", CPL_TYPE_DOUBLE, ndiode);
3172 double ** phase_telfc_corr = gravi_table_get_data_array_double (vismet_table, "PHASE_TELFC_CORR");
3173
3174
3175 /* Compute PHASE_TELFC_CORR as (PHASE_TEL + OPD_TEL_CORR*2*pi/lambda) - (PHASE_FC + OPD_FC_CORR*2*pi/lambda) */
3176 /* going to complex phasor notation and then back to phase */
3177 double phi;
3178 double phifc;
3179 double complex phasor;
3180 double complex phasorfc;
3181 double complex dphasor;
3182
3183 for (int tel = 0; tel < ntel; tel++) {
3184 for (cpl_size row = 0; row < nrow_met; row++) {
3185 phifc = - phase_fc[row*ntel+tel]
3186 + opd_fc_corr[row*ntel+tel] / lambda_met_mean * TWOPI; /* fiber coupler */
3187 phasorfc = cos(phifc) + sin(phifc) * I;
3188 for (int diode = 0; diode < ndiode; diode++) {
3189 phi = - phase_tel[row*ntel+tel][diode] + (opd_tel_corr[row*ntel+tel][diode])
3190 / lambda_met_mean * TWOPI; /* telescope receiver */
3191 phasor = cos(phi) + sin(phi) * I;
3192 dphasor = phasor*conj(phasorfc);
3193 phase_telfc_corr[row*ntel+tel][diode] = carg(dphasor);
3194 }
3195 }
3196 }
3197
3198 CPLCHECK_MSG ("Cannot calculate PHASE_TELFC_CORR");
3199
3200 /*****************************************************************
3201 * PART IV: OPD_TELFC_CORR_XY
3202 ****************************************************************
3203 * Calculation of the actual astimatism (smoothed over 5s)
3204 * Calculation of the separation of the 2 fibers (smoothed over 5s)
3205 * Correction of both are stored in OPD_TELFC_CORR_XY
3206 * but calculations are all done in phase (hence the name phase_telfc_corr_xy)
3207 * Correction of both are subtracted to create PHASE_TELFC_CORR
3208 * PHASE_TELFC_CORR is wrapped around its mean value
3209 * PHASE_TELFC_CORR is then stored in OPD_TELFC_CORR
3210 *****************************************************************/
3211
3212 double tilt_mean[4];
3213 double tip_mean[4];
3214 double astig_mean[4];
3215
3216 gravi_table_init_column_array (vismet_table, "OPD_TELFC_CORR_XY", "rad", CPL_TYPE_DOUBLE, ndiode);
3217 double ** phase_telfc_corr_xy = gravi_table_get_data_array_double (vismet_table, "OPD_TELFC_CORR_XY");
3218 /* the variable name phase_telfc_corr_xy is used because calculated in phase */
3219 /* It is changed at the last moment to an OPD (for verification by user) */
3220
3221 int Nsmooth_astig = 2500; /* 10 seconds smoothing */
3222 int Nsmooth_sep = 1250; /* 5 seconds smoothing */
3223 int Nsmooth_tel = 50; /* 200 ms smoothing */
3224
3225 /* if faint mode activated, the smoothing length is increased to the SC DIT */
3226 if (gravi_pfits_get_met_mode(header)==MET_FAINT_HEADER)
3227 {
3228 double dit_sc = gravi_pfits_get_dit_sc (header);
3229 Nsmooth_astig = (int) 250 * 4 * (dit_sc +1.0);
3230 Nsmooth_sep = (int) 250 * 2 * (dit_sc + 1.0);
3231 if ( use_faint_met == 1) {
3232 Nsmooth_tel = 400;
3233 } else {
3234 Nsmooth_tel = (int) 250 * 2* (dit_sc + 1.0);
3235 }
3236 }
3237
3238 for (int tel = 0; tel < ntel; tel++)
3239 {
3240 cpl_msg_info (cpl_func,"G-FAINT: Use flagged amplitude for astigmatism");
3241 cpl_msg_info (cpl_func,"Smoothing astigmatism by %d metrology DITS (tel=%i)",Nsmooth_astig*2+1,tel);
3242 cpl_array * phasor_array = cpl_array_new(nrow_met,CPL_TYPE_DOUBLE_COMPLEX);
3243 for (cpl_size row = 0; row < nrow_met; row++)
3244 {
3245 double amplitude_ast = amplitude_tel_flagged[row*ntel+tel][3] *
3246 amplitude_tel_flagged[row*ntel+tel][2] *
3247 amplitude_tel_flagged[row*ntel+tel][1] *
3248 amplitude_tel_flagged[row*ntel+tel][0];
3249 double phase_ast = phase_telfc_corr[row*ntel+tel][3] -
3250 phase_telfc_corr[row*ntel+tel][2] +
3251 phase_telfc_corr[row*ntel+tel][1] -
3252 phase_telfc_corr[row*ntel+tel][0];
3253 double complex complex_ast = amplitude_ast * (cos(phase_ast)+I*sin(phase_ast));
3254
3255 cpl_array_set_double_complex(phasor_array,row,complex_ast);
3256 }
3257 CPLCHECK_MSG("Error at computation of phasors");
3258
3259 double complex mean_phasor = cpl_array_get_mean_complex (phasor_array);
3260 cpl_array_multiply_scalar_complex (phasor_array, conj(mean_phasor));
3261 cpl_array * phase_astig_corr_array = gravi_array_smooth (phasor_array, Nsmooth_astig);
3262 cpl_array_arg(phase_astig_corr_array);
3263 cpl_array_add_scalar(phase_astig_corr_array, carg(mean_phasor));
3264
3265 astig_mean[tel]=carg(mean_phasor);
3266 cpl_msg_info(cpl_func,"Mean astigmatism for tel %i : %.2f radians",tel,carg(mean_phasor));
3267
3268 CPLCHECK_MSG ("Failed when smoothing astigmatism");
3269 cpl_array_delete(phasor_array);
3270
3271
3272 cpl_array * phasor_array_X = cpl_array_new(nrow_met,CPL_TYPE_DOUBLE_COMPLEX);
3273 cpl_array * phasor_array_Y = cpl_array_new(nrow_met,CPL_TYPE_DOUBLE_COMPLEX);
3274
3275 /* get smoothed value of X separation angle */
3276 cpl_msg_info (cpl_func,"G-FAINT: Use flagged amplitude for sep. angle");
3277 cpl_msg_info (cpl_func,"Smoothing X and Y sep by %d metrology DITS (tel=%i)",Nsmooth_sep*2+1,tel);
3278 for (cpl_size row = 0; row < nrow_met; row++)
3279 {
3280 double phase_ast_2=cpl_array_get_double(phase_astig_corr_array,row,NULL)/2.;
3281
3282 double amplitude_sep = amplitude_tel_flagged[row*ntel+tel][3] *
3283 amplitude_tel_flagged[row*ntel+tel][2] *
3284 amplitude_tel_flagged[row*ntel+tel][1] *
3285 amplitude_tel_flagged[row*ntel+tel][0];
3286
3287 double phase_sepX1 = phase_telfc_corr[row*ntel+tel][3] -
3288 phase_telfc_corr[row*ntel+tel][2] - phase_ast_2;
3289 double phase_sepX2 = phase_telfc_corr[row*ntel+tel][0] -
3290 phase_telfc_corr[row*ntel+tel][1] + phase_ast_2;
3291
3292 double phase_sepY1 = phase_telfc_corr[row*ntel+tel][3] -
3293 phase_telfc_corr[row*ntel+tel][0] - phase_ast_2;
3294 double phase_sepY2 = phase_telfc_corr[row*ntel+tel][2] -
3295 phase_telfc_corr[row*ntel+tel][1] + phase_ast_2;
3296
3297 double complex phase_sepX = amplitude_sep * (cos(phase_sepX1) + cos(phase_sepX2)
3298 + I*( sin(phase_sepX1) + sin(phase_sepX2) ));
3299 double complex phase_sepY = amplitude_sep * (cos(phase_sepY1) + cos(phase_sepY2)
3300 + I*( sin(phase_sepY1) + sin(phase_sepY2) ));
3301
3302 cpl_array_set_double_complex(phasor_array_X,row,phase_sepX);
3303 cpl_array_set_double_complex(phasor_array_Y,row,phase_sepY);
3304
3305 }
3306
3307 double complex mean_phasor_X = cpl_array_get_mean_complex (phasor_array_X);
3308 cpl_array_multiply_scalar_complex (phasor_array_X, conj(mean_phasor_X));
3309 cpl_array * phase_sepX_corr_array = gravi_array_smooth (phasor_array_X, Nsmooth_sep);
3310 cpl_array_arg(phase_sepX_corr_array);
3311 cpl_array_add_scalar(phase_sepX_corr_array, carg(mean_phasor_X));
3312
3313 double complex mean_phasor_Y = cpl_array_get_mean_complex (phasor_array_Y);
3314 cpl_array_multiply_scalar_complex (phasor_array_Y, conj(mean_phasor_Y));
3315 cpl_array * phase_sepY_corr_array = gravi_array_smooth (phasor_array_Y, Nsmooth_sep);
3316 cpl_array_arg(phase_sepY_corr_array);
3317 cpl_array_add_scalar(phase_sepY_corr_array, carg(mean_phasor_Y));
3318
3319 tip_mean[tel]=carg(mean_phasor_X);
3320 tilt_mean[tel]=carg(mean_phasor_Y);
3321 cpl_msg_info(cpl_func,"Mean X separation for tel %i : %.2f radians",tel,carg(mean_phasor_X));
3322 cpl_msg_info(cpl_func,"Mean Y separation for tel %i : %.2f radians",tel,carg(mean_phasor_Y));
3323
3324
3325 CPLCHECK_MSG ("Failed when smoothing separation vectors");
3326 cpl_array_delete(phasor_array_X);
3327 cpl_array_delete(phasor_array_Y);
3328
3329 /* compute the values to be stored in phase_telfc_corr_xy */
3330 /* subtract the values to phase_telfc_corr */
3331
3332 cpl_array * phasor_telfc_corr_array[ndiode];
3333 for (int diode = 0; diode < ndiode; diode++)
3334 phasor_telfc_corr_array[diode] = cpl_array_new(nrow_met,CPL_TYPE_DOUBLE_COMPLEX);
3335
3336 cpl_msg_info(cpl_func,"Computing the PHASE_TELFC_CORR_XY column");
3337 if ( use_faint_met == 1) {
3338 cpl_msg_info(cpl_func,"G-FAINT: Use faint amplitude for PHASE_TELFC_CORR_XY");
3339 } else {
3340 cpl_msg_info(cpl_func,"G-FAINT: Use flagged amplitude for PHASE_TELFC_CORR_XY");
3341 }
3342 for (cpl_size row = 0; row < nrow_met; row++)
3343 {
3344 double amplitude_diodes;
3345 if ( use_faint_met == 1) {
3346 amplitude_diodes = amplitude_tel[row*ntel+tel][3] +
3347 amplitude_tel[row*ntel+tel][2] +
3348 amplitude_tel[row*ntel+tel][1] +
3349 amplitude_tel[row*ntel+tel][0];
3350 } else {
3351 amplitude_diodes = amplitude_tel_flagged[row*ntel+tel][3] +
3352 amplitude_tel_flagged[row*ntel+tel][2] +
3353 amplitude_tel_flagged[row*ntel+tel][1] +
3354 amplitude_tel_flagged[row*ntel+tel][0];
3355 }
3356 phase_telfc_corr_xy[row*ntel+tel][0]=
3357 -cpl_array_get_double(phase_astig_corr_array,row,NULL)/4
3358 +cpl_array_get_double(phase_sepX_corr_array,row,NULL)/2
3359 -cpl_array_get_double(phase_sepY_corr_array,row,NULL)/2;
3360 phase_telfc_corr_xy[row*ntel+tel][1]=
3361 cpl_array_get_double(phase_astig_corr_array,row,NULL)/4
3362 -cpl_array_get_double(phase_sepX_corr_array,row,NULL)/2
3363 -cpl_array_get_double(phase_sepY_corr_array,row,NULL)/2;
3364 phase_telfc_corr_xy[row*ntel+tel][2]=
3365 -cpl_array_get_double(phase_astig_corr_array,row,NULL)/4
3366 -cpl_array_get_double(phase_sepX_corr_array,row,NULL)/2
3367 +cpl_array_get_double(phase_sepY_corr_array,row,NULL)/2;
3368 phase_telfc_corr_xy[row*ntel+tel][3]=
3369 cpl_array_get_double(phase_astig_corr_array,row,NULL)/4
3370 +cpl_array_get_double(phase_sepX_corr_array,row,NULL)/2
3371 +cpl_array_get_double(phase_sepY_corr_array,row,NULL)/2;
3372
3373 for (int diode = 0; diode < ndiode; diode++) {
3374 double phase_telfc_tmp = phase_telfc_corr[row*ntel+tel][diode]-phase_telfc_corr_xy[row*ntel+tel][diode];
3375 double complex amp_telfc_corr = amplitude_diodes*(cos(phase_telfc_tmp)+I*sin(phase_telfc_tmp));
3376 cpl_array_set_double_complex(phasor_telfc_corr_array[diode],row,amp_telfc_corr);
3377 }
3378
3379 /* changing phase_telfc_corr_xy to opd value to store into fits file as OPD*/
3380 for (int diode = 0; diode < ndiode; diode++) {
3381 phase_telfc_corr_xy[row*ntel+tel][diode]*=lambda_met_mean /TWOPI ;
3382 }
3383 }
3384 CPLCHECK_MSG ("Failed when computing phasor for PHASE_TELFC_CORR_XY");
3385
3386 cpl_msg_info (cpl_func,"Smoothing phase_telfc_corr by %d metrology DITS (tel=%i)",Nsmooth_tel*2+1,tel);
3387
3388 for (int diode = 0; diode < ndiode; diode++) {
3389 double complex mean_phasor = cpl_array_get_mean_complex (phasor_telfc_corr_array[diode]);
3390 cpl_array_multiply_scalar_complex (phasor_telfc_corr_array[diode], conj(mean_phasor));
3391 cpl_array * phase_telfc_corr_array = gravi_array_smooth (phasor_telfc_corr_array[diode], Nsmooth_tel);
3392 cpl_array_arg(phase_telfc_corr_array);
3393 cpl_array_add_scalar(phase_telfc_corr_array, carg(mean_phasor));
3394
3395 for (cpl_size row = 0; row < nrow_met; row++)
3396 {
3397 phase_telfc_corr[row*ntel+tel][diode] = cpl_array_get_double(phase_telfc_corr_array,row,NULL);
3398 opd_telfc_corr[row*ntel+tel][diode] = phase_telfc_corr[row*ntel+tel][diode] / TWOPI * lambda_met_mean ;
3399 }
3400 cpl_array_delete(phase_telfc_corr_array);
3401 cpl_array_delete(phasor_telfc_corr_array[diode]);
3402 }
3403
3404 CPLCHECK_MSG ("Failed when smoothing PHASE_TELFC_CORR_XY");
3405 cpl_array_delete(phase_astig_corr_array);
3406 cpl_array_delete(phase_sepX_corr_array);
3407 cpl_array_delete(phase_sepY_corr_array);
3408
3409 }
3410
3411 /* from now on, old code from Frank. Here for historical reason */
3412 /* It is not needed for Part V */
3413
3414 cpl_vector * tmp_vector;
3415 tmp_vector = cpl_vector_new (nrow_met);
3416 double tmp_median;
3417 double mmet[4][4];
3418
3419 /* report final median for each telescope and diode */
3420 for (int tel = 0; tel < ntel; tel++) {
3421 for (int diode = 0; diode < ndiode; diode++) {
3422 /* fill tmp_vector with opd_telfc for given telescope and diode */
3423 for (cpl_size row = 0; row < nrow_met; row++) {
3424 cpl_vector_set (tmp_vector, row, opd_telfc_corr[row*ntel+tel][diode]);
3425 }
3426 /* calculate median */
3427 tmp_median = cpl_vector_get_median_const(tmp_vector);
3428 cpl_msg_info (cpl_func,"Median TEL-FC in nm for Tel %d Diode %d : %g ",
3429 tel, diode, tmp_median*1e9);
3430 /* FE: put results in a local variable for below residual tilt calculation */
3431 mmet[tel][diode] = tmp_median;
3432 }
3433 }
3434
3435 /* Put QC parameters with residual Mean and simple TT and ASTIG estimates;
3436 simple means, not taking into acount detailed receiver positions, but working
3437 on sum and differences of opposite diodes;
3438 convention such that results give "residual in nm per diode"
3439 OFFS = Mean of all diodes
3440 = (D0 + D1 + D2 + D3) / 4
3441 TIP/TILT = 1/2 * difference of opposite diodes
3442 TIP = (D0 - D2) / 2
3443 TILT = (D1 - D3) / 2
3444 ASTIG = 1/2 * difference of Mean of opposite Diodes) =
3445 = [ (D0 + D2) / 2 - (D1 + D3) / 2 ] / 2
3446 = (D0 - D1 + D2 - D3) / 4
3447 */
3448
3449 for (int tel=0; tel<ntel; tel++) {
3450 /* old way (weighted for faint mode) to calculate the tip/tilt/astigmatism residuals,mean over phase */
3451 sprintf (qc_name, "ESO QC MET OFF%i", tel+1);
3452 cpl_msg_info (cpl_func, "%s = %f [nm]", qc_name, (mmet[tel][0]+mmet[tel][1]+mmet[tel][2]+mmet[tel][3])/4.*1e9);
3453 cpl_propertylist_update_double (header, qc_name, (mmet[tel][0]+mmet[tel][1]+mmet[tel][2]+mmet[tel][3])/4.*1e9);
3454 cpl_propertylist_set_comment (header, qc_name, "[nm] residual metrology offset");
3455 sprintf (qc_name, "ESO QC MET TIP%i", tel+1);
3456 cpl_msg_info (cpl_func, "%s = %f [nm]", qc_name, (mmet[tel][0]-mmet[tel][2])/2.*1e9);
3457 sprintf (qc_name, "ESO QC MET TILT%i", tel+1);
3458 cpl_msg_info (cpl_func, "%s = %f [nm]", qc_name, (mmet[tel][1]-mmet[tel][3])/2.*1e9);
3459 sprintf (qc_name, "ESO QC MET ASTIG%i", tel+1);
3460 cpl_msg_info (cpl_func, "%s = %f [nm]", qc_name,(mmet[tel][0]-mmet[tel][1]+mmet[tel][2]-mmet[tel][3])/4.*1e9);
3461 /* new way (weighted for faint mode) to calculate the tip/tilt/astigmatism residuals,mean over phasor */
3462 sprintf (qc_name, "ESO QC MET TIP%i", tel+1);
3463 cpl_msg_info (cpl_func, "%s = %f [nm]", qc_name, tip_mean[tel]*lambda_met_mean*1e9/TWOPI);
3464 cpl_propertylist_update_double (header, qc_name, tip_mean[tel]*lambda_met_mean*1e9/TWOPI);
3465 cpl_propertylist_set_comment (header, qc_name, "[nm] residual metrology tip");
3466 sprintf (qc_name, "ESO QC MET TILT%i", tel+1);
3467 cpl_msg_info (cpl_func, "%s = %f [nm]", qc_name, tilt_mean[tel]*lambda_met_mean*1e9/TWOPI);
3468 cpl_propertylist_update_double (header, qc_name, tilt_mean[tel]*lambda_met_mean*1e9/TWOPI);
3469 cpl_propertylist_set_comment (header, qc_name, "[nm] residual metrology tilt");
3470 sprintf (qc_name, "ESO QC MET ASTIG%i", tel+1);
3471 cpl_msg_info (cpl_func, "%s = %f [nm]", qc_name,astig_mean[tel]*lambda_met_mean*1e9/TWOPI);
3472 cpl_propertylist_update_double (header, qc_name, astig_mean[tel]*lambda_met_mean*1e9/TWOPI);
3473 cpl_propertylist_set_comment (header, qc_name, "[nm] residual metrology astigmatism");
3474
3475 }
3476
3477 /* FE: calculate residual tilt of metrology signal for each telescope and convert to
3478 equivalent dx,dy on detector. In other words, output the calibration error of fc_fiber_dxy
3479 */
3480
3481 double mttx[4];
3482 double mtty[4];
3483 double mttn[4];
3484 double mtte[4];
3485 double mttdx[4];
3486 double mttdy[4];
3487 for (int tel = 0; tel < ntel; tel++) {
3488 mttx[tel] = (mmet[tel][3]*rec_az[tel][3]*pow(rec_zd[tel][0],2) - mmet[tel][0]*rec_az[tel][1]*rec_zd[tel][0]*rec_zd[tel][1] + mmet[tel][0]*rec_az[tel][0]*pow(rec_zd[tel][1],2) +
3489 mmet[tel][3]*rec_az[tel][3]*pow(rec_zd[tel][1],2) - mmet[tel][0]*rec_az[tel][2]*rec_zd[tel][0]*rec_zd[tel][2] + mmet[tel][0]*rec_az[tel][0]*pow(rec_zd[tel][2],2) +
3490 mmet[tel][3]*rec_az[tel][3]*pow(rec_zd[tel][2],2) - mmet[tel][3]*rec_az[tel][0]*rec_zd[tel][0]*rec_zd[tel][3] - mmet[tel][0]*rec_az[tel][3]*rec_zd[tel][0]*rec_zd[tel][3] -
3491 mmet[tel][3]*rec_az[tel][1]*rec_zd[tel][1]*rec_zd[tel][3] - mmet[tel][3]*rec_az[tel][2]*rec_zd[tel][2]*rec_zd[tel][3] + mmet[tel][0]*rec_az[tel][0]*pow(rec_zd[tel][3],2) -
3492 mmet[tel][2]*rec_zd[tel][2]*(rec_az[tel][0]*rec_zd[tel][0] + rec_az[tel][1]*rec_zd[tel][1] + rec_az[tel][3]*rec_zd[tel][3]) -
3493 mmet[tel][1]*rec_zd[tel][1]*(rec_az[tel][0]*rec_zd[tel][0] + rec_az[tel][2]*rec_zd[tel][2] + rec_az[tel][3]*rec_zd[tel][3]) +
3494 mmet[tel][2]*rec_az[tel][2]*(pow(rec_zd[tel][0],2) + pow(rec_zd[tel][1],2) + pow(rec_zd[tel][3],2)) +
3495 mmet[tel][1]*rec_az[tel][1]*(pow(rec_zd[tel][0],2) + pow(rec_zd[tel][2],2) + pow(rec_zd[tel][3],2))) /
3496 (pow(rec_az[tel][3],2)*pow(rec_zd[tel][0],2) + pow(rec_az[tel][0],2)*pow(rec_zd[tel][1],2) +
3497 pow(rec_az[tel][3],2)*pow(rec_zd[tel][1],2) + pow(rec_az[tel][0],2)*pow(rec_zd[tel][2],2) +
3498 pow(rec_az[tel][3],2)*pow(rec_zd[tel][2],2) - 2*rec_az[tel][0]*rec_az[tel][3]*rec_zd[tel][0]*rec_zd[tel][3] +
3499 pow(rec_az[tel][0],2)*pow(rec_zd[tel][3],2) - 2*rec_az[tel][2]*rec_zd[tel][2]*(rec_az[tel][0]*rec_zd[tel][0] + rec_az[tel][3]*rec_zd[tel][3]) -
3500 2*rec_az[tel][1]*rec_zd[tel][1]*(rec_az[tel][0]*rec_zd[tel][0] + rec_az[tel][2]*rec_zd[tel][2] + rec_az[tel][3]*rec_zd[tel][3]) +
3501 pow(rec_az[tel][2],2)*(pow(rec_zd[tel][0],2) + pow(rec_zd[tel][1],2) + pow(rec_zd[tel][3],2)) +
3502 pow(rec_az[tel][1],2)*(pow(rec_zd[tel][0],2) + pow(rec_zd[tel][2],2) + pow(rec_zd[tel][3],2)));
3503 mtty[tel] = (-(mmet[tel][2]*rec_az[tel][0]*rec_az[tel][2]*rec_zd[tel][0]) - mmet[tel][3]*rec_az[tel][0]*rec_az[tel][3]*rec_zd[tel][0] - mmet[tel][2]*rec_az[tel][1]*rec_az[tel][2]*rec_zd[tel][1] -
3504 mmet[tel][3]*rec_az[tel][1]*rec_az[tel][3]*rec_zd[tel][1] + mmet[tel][2]*pow(rec_az[tel][0],2)*rec_zd[tel][2] + mmet[tel][2]*pow(rec_az[tel][1],2)*rec_zd[tel][2] -
3505 mmet[tel][3]*rec_az[tel][2]*rec_az[tel][3]*rec_zd[tel][2] + mmet[tel][2]*pow(rec_az[tel][3],2)*rec_zd[tel][2] + mmet[tel][3]*pow(rec_az[tel][0],2)*rec_zd[tel][3] +
3506 mmet[tel][3]*pow(rec_az[tel][1],2)*rec_zd[tel][3] + mmet[tel][3]*pow(rec_az[tel][2],2)*rec_zd[tel][3] - mmet[tel][2]*rec_az[tel][2]*rec_az[tel][3]*rec_zd[tel][3] +
3507 mmet[tel][0]*(pow(rec_az[tel][1],2)*rec_zd[tel][0] + pow(rec_az[tel][2],2)*rec_zd[tel][0] + pow(rec_az[tel][3],2)*rec_zd[tel][0] -
3508 rec_az[tel][0]*rec_az[tel][1]*rec_zd[tel][1] - rec_az[tel][0]*rec_az[tel][2]*rec_zd[tel][2] - rec_az[tel][0]*rec_az[tel][3]*rec_zd[tel][3]) +
3509 mmet[tel][1]*(-(rec_az[tel][0]*rec_az[tel][1]*rec_zd[tel][0]) + pow(rec_az[tel][0],2)*rec_zd[tel][1] + pow(rec_az[tel][2],2)*rec_zd[tel][1] +
3510 pow(rec_az[tel][3],2)*rec_zd[tel][1] - rec_az[tel][1]*rec_az[tel][2]*rec_zd[tel][2] - rec_az[tel][1]*rec_az[tel][3]*rec_zd[tel][3]))/
3511 (pow(rec_az[tel][3],2)*pow(rec_zd[tel][0],2) + pow(rec_az[tel][0],2)*pow(rec_zd[tel][1],2) +
3512 pow(rec_az[tel][3],2)*pow(rec_zd[tel][1],2) + pow(rec_az[tel][0],2)*pow(rec_zd[tel][2],2) +
3513 pow(rec_az[tel][3],2)*pow(rec_zd[tel][2],2) - 2*rec_az[tel][0]*rec_az[tel][3]*rec_zd[tel][0]*rec_zd[tel][3] +
3514 pow(rec_az[tel][0],2)*pow(rec_zd[tel][3],2) - 2*rec_az[tel][2]*rec_zd[tel][2]*(rec_az[tel][0]*rec_zd[tel][0] + rec_az[tel][3]*rec_zd[tel][3]) -
3515 2*rec_az[tel][1]*rec_zd[tel][1]*(rec_az[tel][0]*rec_zd[tel][0] + rec_az[tel][2]*rec_zd[tel][2] + rec_az[tel][3]*rec_zd[tel][3]) +
3516 pow(rec_az[tel][2],2)*(pow(rec_zd[tel][0],2) + pow(rec_zd[tel][1],2) + pow(rec_zd[tel][3],2)) +
3517 pow(rec_az[tel][1],2)*(pow(rec_zd[tel][0],2) + pow(rec_zd[tel][2],2) + pow(rec_zd[tel][3],2)));
3518 // above should be in radians , convert to mas
3519 mttx[tel] = mttx[tel] / TWOPI * 360. * 3600. * 1000.;
3520 mtty[tel] = mtty[tel] / TWOPI * 360. * 3600. * 1000.;
3521 // now convert tp pixel above should be in radians, divide by pixel scale to get pixel
3522 // FE 2018-02-09 using default 18 mas/pixel from UTs, because division by 0.0 would cause NAN problems below
3523 // double scale = 0.0;
3524 double scale = 18.0;
3525 sprintf (card,"ESO QC ACQ FIELD%d SCALE", tel+1);
3526 if (cpl_propertylist_has (header, card))
3527 scale = cpl_propertylist_get_double (header, card);
3528
3529 // FE 2018-02-09 check if telescope given, otherwise
3530 // receiver pos are all zero and we get NaN
3531 /* Get telname */
3532 const char * telname = gravi_conf_get_telname (tel, header);
3533 if (telname == NULL) {
3534 mttx[tel] = 0.;
3535 mtty[tel] = 0.;
3536 }
3537
3538 mttx[tel] = mttx[tel] / scale;
3539 mtty[tel] = mtty[tel] / scale;
3540 // rotate to sky (north, east) using parang
3541 mtte[tel] = mttx[tel]*cos(parang) + mtty[tel]*sin(parang);
3542 mttn[tel] = - mttx[tel]*sin(parang) + mtty[tel]*cos(parang);
3543 // rotate to acqcam
3544 northangle = cpl_array_get_double(northangle_array,tel,NULL);
3545 mttdx[tel] = mttn[tel]*sin(northangle * CPL_MATH_RAD_DEG ) + mtte[tel]*cos(northangle * CPL_MATH_RAD_DEG );
3546 mttdy[tel] = mttn[tel]*cos(northangle * CPL_MATH_RAD_DEG ) - mtte[tel]*sin(northangle * CPL_MATH_RAD_DEG );
3547 }
3548
3549 // FE 20190509 changed uvw -> telescope, beause I think mttxy are in telescope coordinate system, not uvw (which is aligned to N/E
3550 cpl_msg_debug (cpl_func,"FE: mttx, mtty in telescope coordinates of GV1: %g %g pixel", mttx[0], mtty[0]);
3551 cpl_msg_debug (cpl_func,"FE: mttx, mtty in telescope coordinates of GV2: %g %g pixel", mttx[1], mtty[1]);
3552 cpl_msg_debug (cpl_func,"FE: mttx, mtty in telescope coordinates of GV3: %g %g pixel", mttx[2], mtty[2]);
3553 cpl_msg_debug (cpl_func,"FE: mttx, mtty in telescope coordinates of GV4: %g %g pixel", mttx[3], mtty[3]);
3554 cpl_msg_debug (cpl_func,"FE: mttx, mtty in RA DEC coordinates of GV1: %g %g pixel", mtte[0], mttn[0]);
3555 cpl_msg_debug (cpl_func,"FE: mttx, mtty in RA DEC coordinates of GV2: %g %g pixel", mtte[1], mttn[1]);
3556 cpl_msg_debug (cpl_func,"FE: mttx, mtty in RA DEC coordinates of GV3: %g %g pixel", mtte[2], mttn[2]);
3557 cpl_msg_debug (cpl_func,"FE: mttx, mtty in RA DEC coordinates of GV4: %g %g pixel", mtte[3], mttn[3]);
3558 cpl_msg_debug (cpl_func,"FE: correction for dxy of GV1: %g %g pixel", mttdx[0], mttdy[0]);
3559 cpl_msg_debug (cpl_func,"FE: correction for dxy of GV2: %g %g pixel", mttdx[1], mttdy[1]);
3560 cpl_msg_debug (cpl_func,"FE: correction for dxy of GV3: %g %g pixel", mttdx[2], mttdy[2]);
3561 cpl_msg_debug (cpl_func,"FE: correction for dxy of GV4: %g %g pixel", mttdx[3], mttdy[3]);
3562
3563 /* Add QC parameters for the corrections */
3564 // FE: declared above
3565 // char qc_name[100];
3566 for (int tel=0; tel<ntel; tel++) {
3567 sprintf (qc_name, "ESO QC MET FIBER SC%iDX", tel+1);
3568 cpl_msg_info (cpl_func, "%s = %f", qc_name, mttdx[tel]);
3569 cpl_propertylist_update_double (header, qc_name, mttdx[tel]);
3570 cpl_propertylist_set_comment (header, qc_name, "[pix] SC fiber offset measured by tel metrology");
3571 sprintf (qc_name, "ESO QC MET FIBER SC%iDY", tel+1);
3572 cpl_msg_info (cpl_func, "%s = %f", qc_name, mttdy[tel]);
3573 cpl_propertylist_update_double (header, qc_name, mttdy[tel]);
3574 cpl_propertylist_set_comment (header, qc_name, "[pix] SC fiber offset measured by tel metrology");
3575 }
3576
3577 /* FE 2018-02-10 Calculate best guess offset of SC fiber from
3578 desired SOBJXY given by SC_FIBER_DX - QC.MET.FIBER.SCiDX:
3579 sobj_dx/y in pixel on acqcam
3580 sobj_dra/dec in mas on sky
3581 sign convention is such that sobj_dxy is identical to value measured by acqcam
3582 when residual metrology tilt is zero, and - QC.MET.FIBER.SCiDX in case the acqcam
3583 doesn't see any error, but the metrology does */
3584
3585 double sobj_dx[4]; /* total offset in pixel in acqcam coordinates */
3586 double sobj_dy[4]; /* total offset in pixel in acqcam coordinates */
3587 double sobj_dra[4]; /* total RA offset in mas */
3588 double sobj_ddec[4]; /* total DEC offset in mas */
3589 double sc_fiber_dx[4]; /* pixel offset measured by AcqCam (without taking into account metrology tilt) in pixel */
3590 double sc_fiber_dy[4]; /* pixel offset measured by AcqCam (without taking into account metrology tilt) in pixel */
3591
3592 for (int tel=0; tel<ntel; tel++) {
3593
3594 /* In acquisition camera X/Y pixel */
3595 /* FE 20190504: if use_fiber_dxy is FALSE then ignore fiber_dxy */
3596 sprintf (card,"ESO QC ACQ FIELD%d SC_FIBER_DX", tel+1);
3597 if (cpl_propertylist_has (header, card) && use_fiber_dxy != 0) {
3598 sc_fiber_dx[tel] = cpl_propertylist_get_double (header, card);
3599 sobj_dx[tel] = sc_fiber_dx[tel] - mttdx[tel];
3600 } else {
3601 sobj_dx[tel] = - mttdx[tel];
3602 }
3603 sprintf (card,"ESO QC ACQ FIELD%d SC_FIBER_DY", tel+1);
3604 if (cpl_propertylist_has (header, card) && use_fiber_dxy != 0) {
3605 sc_fiber_dy[tel] = cpl_propertylist_get_double (header, card);
3606 sobj_dy[tel] = sc_fiber_dy[tel] - mttdy[tel];
3607 } else {
3608 sobj_dy[tel] = - mttdy[tel];
3609 }
3610
3611 sprintf (qc_name, "ESO QC MET SOBJ DX%i", tel+1);
3612 cpl_msg_info (cpl_func, "%s = %f", qc_name, sobj_dx[tel]);
3613 cpl_propertylist_update_double (header, qc_name, sobj_dx[tel]);
3614 cpl_propertylist_set_comment (header, qc_name, "[pixel] x offset from SOBJ");
3615 sprintf (qc_name, "ESO QC MET SOBJ DY%i", tel+1);
3616 cpl_msg_info (cpl_func, "%s = %f", qc_name, sobj_dy[tel]);
3617 cpl_propertylist_update_double (header, qc_name, sobj_dy[tel]);
3618 cpl_propertylist_set_comment (header, qc_name, "[pixel] y offset from SOBJ");
3619
3620 /* on Sky in RA/DEC mas, undoing above ttdx calculation now for
3621 offset including offset measured on acqcam */
3622 // rotate from acqcam to sky
3623 // FE 2018-02-10: below declaration should be done one in the beginning
3624 northangle = cpl_array_get_double(northangle_array,tel,NULL);
3625 sobj_dra[tel] = sobj_dx[tel]*cos(northangle * CPL_MATH_RAD_DEG ) - sobj_dy[tel]*sin(northangle * CPL_MATH_RAD_DEG );
3626 sobj_ddec[tel] = sobj_dx[tel]*sin(northangle * CPL_MATH_RAD_DEG ) + sobj_dy[tel]*cos(northangle * CPL_MATH_RAD_DEG );
3627 // convert to mas
3628 // FE 2018-02-10 using the same dafault as above, if no acqcam measurement available
3629 double scale = 18.0;
3630 sprintf (card,"ESO QC ACQ FIELD%d SCALE", tel+1);
3631 if (cpl_propertylist_has (header, card))
3632 scale = cpl_propertylist_get_double (header, card);
3633 sobj_dra[tel] = sobj_dra[tel] * scale;
3634 sobj_ddec[tel] = sobj_ddec[tel] * scale;
3635
3636 sprintf (qc_name, "ESO QC MET SOBJ DRA%i", tel+1);
3637 cpl_msg_info (cpl_func, "%s = %f", qc_name, sobj_dra[tel]);
3638 cpl_propertylist_update_double (header, qc_name, sobj_dra[tel]);
3639 cpl_propertylist_set_comment (header, qc_name, "[mas] RA offset from SOBJ");
3640 sprintf (qc_name, "ESO QC MET SOBJ DDEC%i", tel+1);
3641 cpl_msg_info (cpl_func, "%s = %f", qc_name, sobj_ddec[tel]);
3642 cpl_propertylist_update_double (header, qc_name, sobj_ddec[tel]);
3643 cpl_propertylist_set_comment (header, qc_name, "[mas] DEC offset from SOBJ");
3644 }
3645
3646
3647 /*****************************************************************/
3648 /* PART V: OPD_TELFC_MCORR */
3649 /*****************************************************************/
3650 /* Average the four telescope diodes. */
3651 /* Note: could also be other combinations, e.g. the less noisy */
3652 /* mean of opposite diodes. */
3653 /*****************************************************************/
3654 // FE 2019-05-15 remark: using less noisy mean of opposite diodes
3655 // show systematic differences in the one case I carefully checked
3656
3657 cpl_msg_info (cpl_func,"Calculate OPD_TELFC_MCORR.");
3658
3659 /* Create array in OI_VIS_MET table, fill with zeros, and get pointer */
3660 gravi_table_new_column (vismet_table, "OPD_TELFC_MCORR", "m", CPL_TYPE_DOUBLE);
3661 double * opd_telfc_mcorr = cpl_table_get_data_double (vismet_table, "OPD_TELFC_MCORR");
3662
3663 /* Calculate mean correction and store in table */
3664 for (int tel = 0; tel < ntel; tel++) {
3665 for (cpl_size row = 0; row < nrow_met; row++) {
3666 opd_telfc_mcorr[row*ntel+tel] =
3667 (opd_telfc_corr[row*ntel+tel][0]
3668 + opd_telfc_corr[row*ntel+tel][1]
3669 + opd_telfc_corr[row*ntel+tel][2]
3670 + opd_telfc_corr[row*ntel+tel][3]) / 4. ;
3671 }
3672 }
3673
3674
3675 FREE (cpl_vector_delete, tmp_vector);
3676
3677
3678 /*****************************************************************/
3679 /* PART VI: OPD_TTPUP */
3680 /*****************************************************************/
3681 /* Calculate the combined tip/tilt*pupil error */
3682 /*****************************************************************/
3683 /* FE 2019-05-07: calculate metrology tip/tilt error = fiber mispointing
3684 currently as above for the QC parameter */
3685
3686 gravi_table_new_column (vismet_table, "FIBER_DU", "rad", CPL_TYPE_DOUBLE);
3687 double * fiber_du = cpl_table_get_data_double (vismet_table, "FIBER_DU");
3688 gravi_table_new_column (vismet_table, "FIBER_DV", "rad", CPL_TYPE_DOUBLE);
3689 double * fiber_dv = cpl_table_get_data_double (vismet_table, "FIBER_DV");
3690 double met_ttx, met_tty; /* temporary variables */
3691
3692 /* to speed up calculation, pre-calculate sine and cosine of paralactic angle */
3693 float sparang = sin(parang);
3694 float cparang = cos(parang);
3695
3696 for (int tel = 0; tel < ntel; tel++) {
3697 for (cpl_size row = 0; row < nrow_met; row++) {
3698 /* calculate terology tip/tilt error telescope coordinate system in radians */
3699 met_ttx = (opd_telfc_corr[row*ntel+tel][3]*rec_az[tel][3]*pow(rec_zd[tel][0],2) - opd_telfc_corr[row*ntel+tel][0]*rec_az[tel][1]*rec_zd[tel][0]*rec_zd[tel][1] + opd_telfc_corr[row*ntel+tel][0]*rec_az[tel][0]*pow(rec_zd[tel][1],2) +
3700 opd_telfc_corr[row*ntel+tel][3]*rec_az[tel][3]*pow(rec_zd[tel][1],2) - opd_telfc_corr[row*ntel+tel][0]*rec_az[tel][2]*rec_zd[tel][0]*rec_zd[tel][2] + opd_telfc_corr[row*ntel+tel][0]*rec_az[tel][0]*pow(rec_zd[tel][2],2) +
3701 opd_telfc_corr[row*ntel+tel][3]*rec_az[tel][3]*pow(rec_zd[tel][2],2) - opd_telfc_corr[row*ntel+tel][3]*rec_az[tel][0]*rec_zd[tel][0]*rec_zd[tel][3] - opd_telfc_corr[row*ntel+tel][0]*rec_az[tel][3]*rec_zd[tel][0]*rec_zd[tel][3] -
3702 opd_telfc_corr[row*ntel+tel][3]*rec_az[tel][1]*rec_zd[tel][1]*rec_zd[tel][3] - opd_telfc_corr[row*ntel+tel][3]*rec_az[tel][2]*rec_zd[tel][2]*rec_zd[tel][3] + opd_telfc_corr[row*ntel+tel][0]*rec_az[tel][0]*pow(rec_zd[tel][3],2) -
3703 opd_telfc_corr[row*ntel+tel][2]*rec_zd[tel][2]*(rec_az[tel][0]*rec_zd[tel][0] + rec_az[tel][1]*rec_zd[tel][1] + rec_az[tel][3]*rec_zd[tel][3]) -
3704 opd_telfc_corr[row*ntel+tel][1]*rec_zd[tel][1]*(rec_az[tel][0]*rec_zd[tel][0] + rec_az[tel][2]*rec_zd[tel][2] + rec_az[tel][3]*rec_zd[tel][3]) +
3705 opd_telfc_corr[row*ntel+tel][2]*rec_az[tel][2]*(pow(rec_zd[tel][0],2) + pow(rec_zd[tel][1],2) + pow(rec_zd[tel][3],2)) +
3706 opd_telfc_corr[row*ntel+tel][1]*rec_az[tel][1]*(pow(rec_zd[tel][0],2) + pow(rec_zd[tel][2],2) + pow(rec_zd[tel][3],2))) /
3707 (pow(rec_az[tel][3],2)*pow(rec_zd[tel][0],2) + pow(rec_az[tel][0],2)*pow(rec_zd[tel][1],2) +
3708 pow(rec_az[tel][3],2)*pow(rec_zd[tel][1],2) + pow(rec_az[tel][0],2)*pow(rec_zd[tel][2],2) +
3709 pow(rec_az[tel][3],2)*pow(rec_zd[tel][2],2) - 2*rec_az[tel][0]*rec_az[tel][3]*rec_zd[tel][0]*rec_zd[tel][3] +
3710 pow(rec_az[tel][0],2)*pow(rec_zd[tel][3],2) - 2*rec_az[tel][2]*rec_zd[tel][2]*(rec_az[tel][0]*rec_zd[tel][0] + rec_az[tel][3]*rec_zd[tel][3]) -
3711 2*rec_az[tel][1]*rec_zd[tel][1]*(rec_az[tel][0]*rec_zd[tel][0] + rec_az[tel][2]*rec_zd[tel][2] + rec_az[tel][3]*rec_zd[tel][3]) +
3712 pow(rec_az[tel][2],2)*(pow(rec_zd[tel][0],2) + pow(rec_zd[tel][1],2) + pow(rec_zd[tel][3],2)) +
3713 pow(rec_az[tel][1],2)*(pow(rec_zd[tel][0],2) + pow(rec_zd[tel][2],2) + pow(rec_zd[tel][3],2)));
3714 met_tty = (-(opd_telfc_corr[row*ntel+tel][2]*rec_az[tel][0]*rec_az[tel][2]*rec_zd[tel][0]) - opd_telfc_corr[row*ntel+tel][3]*rec_az[tel][0]*rec_az[tel][3]*rec_zd[tel][0] - opd_telfc_corr[row*ntel+tel][2]*rec_az[tel][1]*rec_az[tel][2]*rec_zd[tel][1] -
3715 opd_telfc_corr[row*ntel+tel][3]*rec_az[tel][1]*rec_az[tel][3]*rec_zd[tel][1] + opd_telfc_corr[row*ntel+tel][2]*pow(rec_az[tel][0],2)*rec_zd[tel][2] + opd_telfc_corr[row*ntel+tel][2]*pow(rec_az[tel][1],2)*rec_zd[tel][2] -
3716 opd_telfc_corr[row*ntel+tel][3]*rec_az[tel][2]*rec_az[tel][3]*rec_zd[tel][2] + opd_telfc_corr[row*ntel+tel][2]*pow(rec_az[tel][3],2)*rec_zd[tel][2] + opd_telfc_corr[row*ntel+tel][3]*pow(rec_az[tel][0],2)*rec_zd[tel][3] +
3717 opd_telfc_corr[row*ntel+tel][3]*pow(rec_az[tel][1],2)*rec_zd[tel][3] + opd_telfc_corr[row*ntel+tel][3]*pow(rec_az[tel][2],2)*rec_zd[tel][3] - opd_telfc_corr[row*ntel+tel][2]*rec_az[tel][2]*rec_az[tel][3]*rec_zd[tel][3] +
3718 opd_telfc_corr[row*ntel+tel][0]*(pow(rec_az[tel][1],2)*rec_zd[tel][0] + pow(rec_az[tel][2],2)*rec_zd[tel][0] + pow(rec_az[tel][3],2)*rec_zd[tel][0] -
3719 rec_az[tel][0]*rec_az[tel][1]*rec_zd[tel][1] - rec_az[tel][0]*rec_az[tel][2]*rec_zd[tel][2] - rec_az[tel][0]*rec_az[tel][3]*rec_zd[tel][3]) +
3720 opd_telfc_corr[row*ntel+tel][1]*(-(rec_az[tel][0]*rec_az[tel][1]*rec_zd[tel][0]) + pow(rec_az[tel][0],2)*rec_zd[tel][1] + pow(rec_az[tel][2],2)*rec_zd[tel][1] +
3721 pow(rec_az[tel][3],2)*rec_zd[tel][1] - rec_az[tel][1]*rec_az[tel][2]*rec_zd[tel][2] - rec_az[tel][1]*rec_az[tel][3]*rec_zd[tel][3]))/
3722 (pow(rec_az[tel][3],2)*pow(rec_zd[tel][0],2) + pow(rec_az[tel][0],2)*pow(rec_zd[tel][1],2) +
3723 pow(rec_az[tel][3],2)*pow(rec_zd[tel][1],2) + pow(rec_az[tel][0],2)*pow(rec_zd[tel][2],2) +
3724 pow(rec_az[tel][3],2)*pow(rec_zd[tel][2],2) - 2*rec_az[tel][0]*rec_az[tel][3]*rec_zd[tel][0]*rec_zd[tel][3] +
3725 pow(rec_az[tel][0],2)*pow(rec_zd[tel][3],2) - 2*rec_az[tel][2]*rec_zd[tel][2]*(rec_az[tel][0]*rec_zd[tel][0] + rec_az[tel][3]*rec_zd[tel][3]) -
3726 2*rec_az[tel][1]*rec_zd[tel][1]*(rec_az[tel][0]*rec_zd[tel][0] + rec_az[tel][2]*rec_zd[tel][2] + rec_az[tel][3]*rec_zd[tel][3]) +
3727 pow(rec_az[tel][2],2)*(pow(rec_zd[tel][0],2) + pow(rec_zd[tel][1],2) + pow(rec_zd[tel][3],2)) +
3728 pow(rec_az[tel][1],2)*(pow(rec_zd[tel][0],2) + pow(rec_zd[tel][2],2) + pow(rec_zd[tel][3],2)));
3729 // rotate to sky (north, east) using parang (like for QC parameter)
3730 fiber_du[row*ntel+tel] = met_ttx*cparang + met_tty*sparang;
3731 fiber_dv[row*ntel+tel] = - met_ttx*sparang + met_tty*cparang;
3732 }
3733 }
3734
3735 /* calculate correction for combined tip/tilt*pupil term
3736 were we are doing the vector product in uvw coordinates (rahter than sky coordinates),
3737 because both pupil and fiber mispointing already avalable */
3738
3739 double * pupil_u = NULL;
3740 double * pupil_v = NULL;
3741 if ( cpl_table_has_column(vismet_table, "PUPIL_U") && cpl_table_has_column(vismet_table, "PUPIL_V") ) {
3742 pupil_u = cpl_table_get_data_double (vismet_table, "PUPIL_U");
3743 pupil_v = cpl_table_get_data_double (vismet_table, "PUPIL_V");
3744 }
3745 else {
3746 cpl_msg_warning(cpl_func,"Cannot get the PUPIL_U/V (not computed) so will"
3747 "not correct for pupil opd (check the --reduce-acq-cam option)");
3748 }
3749
3750
3751 gravi_table_new_column (vismet_table, "OPD_TTPUP", "m", CPL_TYPE_DOUBLE);
3752 double * opd_ttpup = cpl_table_get_data_double (vismet_table, "OPD_TTPUP");
3753
3754 for (int tel = 0; tel < ntel; tel++) {
3755 for (cpl_size row = 0; row < nrow_met; row++) {
3756 opd_ttpup[row*ntel+tel] = (pupil_u && pupil_v) ?
3757 fiber_du[row*ntel+tel] * pupil_u[row*ntel+tel] +
3758 fiber_dv[row*ntel+tel] * pupil_v[row*ntel+tel] : 0. ;
3759 }
3760 }
3761
3762
3763 /*----------------------------------------------------------------*/
3764 /* FE end */
3765 /*----------------------------------------------------------------*/
3766
3767 CPLCHECK_MSG ("Cannot fill result of TEL vs FC computation");
3768
3769 /* Free the pointer to pointer to data */
3770 cpl_array_delete(northangle_array);
3771 FREE (cpl_free, phase_tel);
3772 FREE (cpl_free, amplitude_tel);
3773 FREE (cpl_free, amplitude_tel_flagged);
3774 FREE (cpl_free, opd_tel);
3775 FREE (cpl_free, opd_tel_corr);
3776 FREE (cpl_free, opd_telfc_corr);
3777 FREE (cpl_free, phase_telfc_corr);
3778 FREE (cpl_free, phase_telfc_corr_xy);
3779
3780 gravi_msg_function_exit(1);
3781 return CPL_ERROR_NONE;
3782}
3783
3784
3785/* -------------------------------------------------------------------- */
3795/* -------------------------------------------------------------------- */
3796static int sin_lambda(const double x[], const double a[], double *result){
3797
3798 *result = a[0] + a[1] * sin(x[0] * (2 * M_PI) / a[3]) +
3799 a[2] * cos(x[0] * (2 * M_PI) / a[3]);
3800
3801 return (0);
3802}
3803
3804static int dfda_sin(const double x[], const double a[], double result[]){
3805 result[0] = 1;
3806 result[1] = sin(x[0] * (2 * M_PI) / a[3]);
3807 result[2] = cos(x[0] * (2 * M_PI) / a[3]);
3808 result[3] = a[1] * x[0] * (2 * M_PI) / (a[3] * a[3]) *
3809 sin(x[0] * (2 * M_PI) / a[3]) - a[1] * x[0] *
3810 (2 * M_PI) / (a[3] * a[3]) * cos(x[0] * (2 * M_PI) / a[3]);
3811 return (0);
3812}
3813
3814/* -------------------------------------------------------------------- */
3828/* -------------------------------------------------------------------- */
3829
3830cpl_table * gravi_metrology_compute_p2vm (cpl_table * metrology_table, double wave_met)
3831{
3832 gravi_msg_function_start(1);
3833 cpl_ensure (metrology_table, CPL_ERROR_NULL_INPUT, NULL);
3834
3835 cpl_vector * vectA, * vectB,
3836 * y_sigma, * init_val, *opl_vect;
3837 const cpl_array * volt;
3838 cpl_array * coherence_fit, * phase_fit;
3839 cpl_matrix * opd_matrix;
3840 int infos = 0;
3841 cpl_vector ** vect;
3842 int val_to_fit2[] = {1,1,1,0}, nv;
3843 double mse, red_chisq;
3844 cpl_array * trans;
3845
3846 int ntel = 4 ;
3847 int n_diode = 80;
3848
3849 /* Get data */
3850 cpl_size nrow = cpl_table_get_nrow (metrology_table);
3851 CPLCHECK_NUL ("Cannot get data");
3852
3853 /* Create p2vm_met : (REGNAME, TRANSMISSION, COHERENCE and PHASE) */
3854 cpl_table * p2vm_met = cpl_table_new (ntel * 2);
3855 cpl_table_new_column_array (p2vm_met,"TRANSMISSION", CPL_TYPE_DOUBLE, 2);
3856 cpl_table_new_column_array (p2vm_met,"COHERENCE", CPL_TYPE_DOUBLE, 2);
3857 cpl_table_new_column_array (p2vm_met,"PHASE", CPL_TYPE_DOUBLE, 2);
3858
3859 CPLCHECK_NUL ("Allocate the metrology output");
3860
3861 /* Loop on SC and FT */
3862 for (int gravi_type = 0; gravi_type < 2; gravi_type++ ){
3863 int comb = (gravi_type == GRAVI_SC ? 1 : 2);
3864
3865 /* Loop on diodes */
3866 for (int tel = 0; tel < ntel; tel++){
3867
3868 /* load vectA and vectB from metrology */
3869 vectA = cpl_vector_new (nrow);
3870 vectB = cpl_vector_new (nrow);
3871
3872 for (cpl_size row = 0; row < nrow; row ++){
3873 volt = cpl_table_get_array ( metrology_table, "VOLT", row);
3874 cpl_vector_set (vectA, row, cpl_array_get (volt , n_diode - 2*(comb*ntel - tel), &nv));
3875 cpl_vector_set (vectB, row, cpl_array_get (volt , n_diode - 2*(comb*ntel - tel) + 1, &nv));
3876 }
3877
3878 /* Compute phase from vectA and vectB*/
3879 opl_vect = gravi_ellipse_phase_create (vectA, vectB, NULL);
3880 cpl_vector_multiply_scalar (opl_vect, wave_met/(2.0*M_PI));
3881 CPLCHECK_NUL("Compute OPD");
3882
3883 /* put OPD into opd_matrix for P2VM_MET calib */
3884 opd_matrix = cpl_matrix_new(nrow, 1);
3885 for (cpl_size row = 0; row < nrow; row++){
3886 cpl_matrix_set (opd_matrix, row, 0, cpl_vector_get (opl_vect, row));
3887 }
3888
3889 cpl_vector_delete (opl_vect);
3890
3891 /* fit on a central window of 5*lambda width */
3892
3893 vect = cpl_malloc (2*sizeof(cpl_vector*));
3894 vect[0] = vectA;
3895 vect[1] = vectB;
3896 phase_fit = cpl_array_new (2, CPL_TYPE_DOUBLE);
3897 coherence_fit = cpl_array_new (2, CPL_TYPE_DOUBLE);
3898 trans = cpl_array_new (2, CPL_TYPE_DOUBLE);
3899
3900 for (int i = 0; i < 2; i++){
3901
3902 /* Get the spectrum of the vector region */
3903 y_sigma = cpl_vector_new(nrow);
3904 cpl_vector_fill(y_sigma, 1);
3905
3906 /* Define and initialize all variables to make a FIT */
3907 mse = 0;
3908 red_chisq = 0;
3909 init_val = cpl_vector_new(4);
3910 cpl_vector_set(init_val, 0, 1);
3911 cpl_vector_set(init_val, 1, 1);
3912 cpl_vector_set(init_val, 2, 1);
3913 cpl_vector_set(init_val, 3, wave_met);
3914
3915 cpl_errorstate prestate = cpl_errorstate_get();
3916 cpl_fit_lvmq(opd_matrix, NULL, vect[i],
3917 y_sigma, init_val, val_to_fit2, &sin_lambda,
3918 &dfda_sin, CPL_FIT_LVMQ_TOLERANCE, CPL_FIT_LVMQ_COUNT,
3919 CPL_FIT_LVMQ_MAXITER, &mse, &red_chisq, NULL);
3920
3921 if (cpl_error_get_code()){
3922 printf("error %f %s and %s \n", 6.0, cpl_error_get_message(), cpl_error_get_where());
3923 return NULL;
3924 }
3925
3926 if (!strcmp("The iterative process did not converge",
3927 cpl_error_get_message())){
3928 if (infos)
3929 cpl_msg_info(cpl_func, "The iterative process "
3930 "did not converge");
3931 cpl_errorstate_set (prestate);
3932 }
3933 if (infos)
3934 cpl_msg_info(cpl_func, "tel = %d : mse "
3935 "%g chi2 %g", tel, mse, red_chisq);
3936
3937
3938 /* Compute the P2VM value */
3939 cpl_array_set_double (coherence_fit, i,
3940 sqrt( pow( cpl_vector_get(init_val, 2), 2) +
3941 pow( cpl_vector_get(init_val, 1), 2)));
3942
3943 cpl_array_set_double (phase_fit, i, atan2( cpl_vector_get(init_val, 2),
3944 cpl_vector_get(init_val, 1)));
3945
3946 cpl_array_set_double (trans, i, cpl_vector_get(init_val, 0));
3947
3948 if (cpl_error_get_code()){
3949 printf("error %f %s and %s \n", 7.0, cpl_error_get_message(), cpl_error_get_where());
3950 return NULL;
3951 }
3952
3953 cpl_vector_delete(init_val);
3954 cpl_vector_delete(y_sigma);
3955 cpl_vector_delete (vect[i]);
3956 }
3957
3958 FREE (cpl_free, vect);
3959 FREE (cpl_matrix_delete, opd_matrix);
3960
3961 cpl_table_set_array (p2vm_met, "COHERENCE", ntel*(1-gravi_type) + tel, coherence_fit);
3962 cpl_array_subtract_scalar (phase_fit,
3963 cpl_array_get_double (phase_fit, 0, &nv));
3964 cpl_table_set_array (p2vm_met, "PHASE", ntel*(1-gravi_type) + tel, phase_fit);
3965 cpl_table_set_array (p2vm_met, "TRANSMISSION", ntel*(1-gravi_type) + tel, trans);
3966
3967 FREE (cpl_array_delete, trans);
3968 FREE (cpl_array_delete, coherence_fit);
3969 FREE (cpl_array_delete, phase_fit);
3970 CPLCHECK_NUL("End loop on tel");
3971 }
3972 }
3973
3974 /* Verbose */
3975 gravi_msg_function_exit(1);
3976 return p2vm_met;
3977}
3978
3979
3980/* -------------------------------------------------------------------- */
3996/* -------------------------------------------------------------------- */
3997cpl_error_code gravi_metrology_reduce (gravi_data * data,
3998 gravi_data * eop_data,
3999 gravi_data * static_param_data,
4000 gravi_data * met_pos,
4001 const cpl_parameterlist * parlist)
4002{
4003 gravi_msg_function_start(1);
4004 cpl_ensure_code (data, CPL_ERROR_NULL_INPUT);
4005
4006 /* Load data */
4007 cpl_table * metrology_table = gravi_data_get_table (data, GRAVI_METROLOGY_EXT);
4008 cpl_propertylist * header = gravi_data_get_header (data);
4009 CPLCHECK_MSG ("Cannot load met extension");
4010
4011 /* Update receiver position */
4012 if (met_pos) {
4013 cpl_table * pos_table = gravi_data_get_table (met_pos, "RECEIVER_POSITION");
4014 gravi_metrology_update_receiverpos (header, pos_table);
4015 CPLCHECK_MSG ("Cannot update receiver positions");
4016 }
4017
4018 /* Create the table */
4019 cpl_table * vismet_table = NULL;
4020 vismet_table = gravi_metrology_create (metrology_table, header);
4021 CPLCHECK_MSG ("Cannot create vismet_table");
4022
4023 /* Compute pointing directions, but do not calculate projected baseline */
4024 int save_pointing = 1;
4025 gravi_eop_pointing_uv (vismet_table, header,
4026 (eop_data ? gravi_data_get_table_x (eop_data, 0) : NULL),
4027 (eop_data ? gravi_data_get_header (eop_data) : NULL),
4028 save_pointing, NULL);
4029
4030 /* If VIS_ACQ table exist, we compute the OPD_PUPIL */
4031 if (gravi_data_has_extension (data, GRAVI_OI_VIS_ACQ_EXT)) {
4032 cpl_table * visacq_table;
4033 visacq_table = gravi_data_get_table (data, GRAVI_OI_VIS_ACQ_EXT);
4034 double delay = gravi_param_get_double_default (parlist, "gravity.metrology.acq-correction-delay",0.0);
4035 double period = gravi_pfits_get_period_acqcam (header);
4036 /* delay is increase by half of the camera period */
4037 delay += period/2;
4038 gravi_metrology_acq (visacq_table, vismet_table, delay, header);
4039 }
4040
4041 /* Reduce the metrology with the DRS algorithm, this
4042 * creates the OI_VIS_MET table */
4043 gravi_metrology_drs (metrology_table, vismet_table, header, parlist);
4044 CPLCHECK_MSG ("Cannot reduce metrology with DRS algo");
4045
4046 /* Add the columns from TAC algorithm */
4047 gravi_metrology_tac (metrology_table, vismet_table, header);
4048 CPLCHECK_MSG ("Cannot reduce metrology with TAC algo");
4049
4050 /* Compute TEL vs FC corrections */
4051 gravi_metrology_telfc (metrology_table, vismet_table, static_param_data, header, parlist);
4052 CPLCHECK_MSG ("Cannot compute TEL vs FC reference");
4053
4054 /* Add the VISMET_TABLE table to the gravi_data */
4055 gravi_data_add_table (data, NULL, GRAVI_OI_VIS_MET_EXT, vismet_table);
4056 CPLCHECK_MSG ("Cannot add OI_VIS_MET in p2vmred_data");
4057
4058 gravi_msg_function_exit(1);
4059 return CPL_ERROR_NONE;
4060}
4061
TWOPI
#define TWOPI
Definition gravi_astrometry.c:47
PI
#define PI
Definition gravi_astrometry.c:46
gravi_cpl.h
gravi_table_get_value
#define gravi_table_get_value(table, name, row, value)
Definition gravi_cpl.h:49
gravi_data.h
gravi_data
typedefCPL_BEGIN_DECLS struct _gravi_data_ gravi_data
Definition gravi_data.h:39
gravi_data_get_header
#define gravi_data_get_header(data)
Definition gravi_data.h:75
ntel
const cpl_size ntel
Definition gravi_demodulate.c:84
ndiode
const cpl_size ndiode
Definition gravi_demodulate.c:85
SC
#define SC
Definition gravi_demodulate.c:55
FT
#define FT
Definition gravi_demodulate.c:54
gravi_dfs.h
gravi_ellipse.h
gravi_eop.h
gravi_metrology.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_pfits.h
GRAVI_OI_VIS_ACQ_EXT
#define GRAVI_OI_VIS_ACQ_EXT
Definition gravi_pfits.h:89
GRAVI_OI_VIS_MET_EXT
#define GRAVI_OI_VIS_MET_EXT
Definition gravi_pfits.h:88
GRAVI_SC
#define GRAVI_SC
Definition gravi_pfits.h:165
MET_FAINT_FLAG
#define MET_FAINT_FLAG
Definition gravi_pfits.h:188
MET_FAINT_HEADER
#define MET_FAINT_HEADER
Definition gravi_pfits.h:186
MODE_SINGLE
#define MODE_SINGLE
Definition gravi_pfits.h:177
GRAVI_METROLOGY_EXT
#define GRAVI_METROLOGY_EXT
Definition gravi_pfits.h:60
MODE_ONAXIS
#define MODE_ONAXIS
Definition gravi_pfits.h:179
gravi_signal.h
gravi_utils.h
CPLCHECK_INT
#define CPLCHECK_INT(msg)
Definition gravi_utils.h:51
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
FREELOOP
#define FREELOOP(function, variable, n)
Definition gravi_utils.h:72
gravi_array_smooth
cpl_array * gravi_array_smooth(cpl_array *input_array, int DIT_smooth)
Definition gravi_cpl.c:1284
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_array_phase_unwrap
cpl_error_code gravi_array_phase_unwrap(cpl_array *input)
Definition gravi_cpl.c:1087
gravi_table_init_column_array
cpl_error_code gravi_table_init_column_array(cpl_table *table, const char *name, const char *unit, cpl_type type, cpl_size size)
Definition gravi_cpl.c:1691
gravi_table_get_data_array_double
double ** gravi_table_get_data_array_double(cpl_table *table, const char *name)
Definition gravi_cpl.c:472
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_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_get_table_x
cpl_table * gravi_data_get_table_x(gravi_data *self, int i)
Get the table of an extension by position.
Definition gravi_data.c:1901
gravi_param_get_bool
int gravi_param_get_bool(const cpl_parameterlist *parlist, const char *name)
Definition gravi_dfs.c:1534
gravi_param_get_int
int gravi_param_get_int(const cpl_parameterlist *parlist, const char *name)
Definition gravi_dfs.c:1521
gravi_param_get_double_default
double gravi_param_get_double_default(const cpl_parameterlist *parlist, const char *name, double def)
Get the parameter from the parameter list.
Definition gravi_dfs.c:1460
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_eop_pointing_uv
cpl_error_code gravi_eop_pointing_uv(cpl_table *input_table, cpl_propertylist *header, cpl_table *eop_table, cpl_propertylist *eop_header, int save_pointing, cpl_table *array_table)
Compute the pointing directions and projected baselines.
Definition gravi_eop.c:228
metrology_unwrap
int metrology_unwrap(double raw_phase, double previous_phase, double max_allowed_phase_diff, double *unwrapped_phase)
Definition gravi_metrology.c:630
gravi_metrology_compute_p2vm
cpl_table * gravi_metrology_compute_p2vm(cpl_table *metrology_table, double wave_met)
Calibrate the P2VM of the metrology.
Definition gravi_metrology.c:3830
metrology_makeDefaultTacConfiguration
structTacConfiguration * metrology_makeDefaultTacConfiguration(double lambda_met)
Definition gravi_metrology.c:326
gravi_metrology_get_posx
double gravi_metrology_get_posx(cpl_propertylist *header, int tel, int diode)
Read the receiver position from header.
Definition gravi_metrology.c:1392
FLAG_JUMP
#define FLAG_JUMP
Definition gravi_metrology.c:128
myAtan
double myAtan(double x, double y, int *flag)
Definition gravi_metrology.c:664
gravi_metrology_telfc
cpl_error_code gravi_metrology_telfc(cpl_table *metrology_table, cpl_table *vismet_table, gravi_data *static_param_data, cpl_propertylist *header, const cpl_parameterlist *parlist)
Best knowledge correction for referencing TEL to FC.
Definition gravi_metrology.c:2714
FLAG_FLUX
#define FLAG_FLUX
Definition gravi_metrology.c:124
gravi_round
long gravi_round(double number)
Definition gravi_metrology.c:107
FLAG_VOLT
#define FLAG_VOLT
Definition gravi_metrology.c:134
gravi_metrology_get_astig
cpl_error_code gravi_metrology_get_astig(cpl_propertylist *header, int gv, double *ampltiude, double *angle, double *radius, gravi_data *static_param_data)
Read the astigmatism amplitude and angle from the fits header.
Definition gravi_metrology.c:1593
gravi_metrology_get_posy
double gravi_metrology_get_posy(cpl_propertylist *header, int tel, int diode)
Definition gravi_metrology.c:1415
gravi_metrology_reduce
cpl_error_code gravi_metrology_reduce(gravi_data *data, gravi_data *eop_data, gravi_data *static_param_data, gravi_data *met_pos, const cpl_parameterlist *parlist)
Reduce the metrology.
Definition gravi_metrology.c:3997
metrology_makeDefaultTacData
structTacData * metrology_makeDefaultTacData(double lambda_met)
Definition gravi_metrology.c:387
MET_MAX_HISTORY
#define MET_MAX_HISTORY
Definition gravi_metrology.c:122
metrology_algorithm
int metrology_algorithm(structTacData *tacData)
Definition gravi_metrology.c:675
metrology_sq
double metrology_sq(double x)
Definition gravi_metrology.c:660
metrology_read_voltages
int metrology_read_voltages(structTacData *tacData, double *volt)
Definition gravi_metrology.c:526
FLAG_SNR
#define FLAG_SNR
Definition gravi_metrology.c:125
FLAG_UNWRAP
#define FLAG_UNWRAP
Definition gravi_metrology.c:126
gravi_metrology_get_fc_focus
double gravi_metrology_get_fc_focus(cpl_propertylist *header, int gv, gravi_data *static_param_data)
Read the fiber coupler focus offset from the fits header.
Definition gravi_metrology.c:1448
FLAG_SPEED
#define FLAG_SPEED
Definition gravi_metrology.c:127
gravi_metrology_update_receiverpos
cpl_error_code gravi_metrology_update_receiverpos(cpl_propertylist *header, cpl_table *receiver_table)
Update the receiver position from header from external calibration.
Definition gravi_metrology.c:1333
sin_lambda
static int sin_lambda(const double x[], const double a[], double *result)
Create the P2VM of the metrology.
Definition gravi_metrology.c:3796
SIN
#define SIN
Definition gravi_metrology.c:142
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
FLAG_PHASOR0
#define FLAG_PHASOR0
Definition gravi_metrology.c:135
FLAG_LOST_PHASE
#define FLAG_LOST_PHASE
Definition gravi_metrology.c:136
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
dfda_sin
static int dfda_sin(const double x[], const double a[], double result[])
Definition gravi_metrology.c:3804
COS
#define COS
Definition gravi_metrology.c:143
gravi_metrology_get_fc_shift
double gravi_metrology_get_fc_shift(cpl_propertylist *header, int gv, gravi_data *static_param_data)
Read the fiber coupler shift parameter, resulting from offset between fiber coupler pickup and pupil ...
Definition gravi_metrology.c:1519
gravi_metrology_acq
cpl_error_code gravi_metrology_acq(cpl_table *visacq_table, cpl_table *vismet_table, double delay, cpl_propertylist *header)
Fill the VIS_MET table with the OPD_PUPIL column.
Definition gravi_metrology.c:1787
gravi_pfits_get_met_ph
const char * gravi_pfits_get_met_ph(const cpl_propertylist *plist)
Definition gravi_pfits.c:534
gravi_pfits_get_gvctu_x
double gravi_pfits_get_gvctu_x(const cpl_propertylist *plist)
Definition gravi_pfits.c:479
gravi_pfits_get_sobj_y
double gravi_pfits_get_sobj_y(const cpl_propertylist *plist)
Definition gravi_pfits.c:461
gravi_pfits_get_mode
int gravi_pfits_get_mode(const cpl_propertylist *plist)
Definition gravi_pfits.c:193
gravi_convert_to_mjd
double gravi_convert_to_mjd(const char *start)
Definition gravi_pfits.c:1188
gravi_pfits_get_dit_sc
double gravi_pfits_get_dit_sc(const cpl_propertylist *plist)
Definition gravi_pfits.c:664
gravi_pfits_get_northangle_acqcam
double gravi_pfits_get_northangle_acqcam(const cpl_propertylist *plist, int tel)
Definition gravi_pfits.c:816
gravi_pfits_get_gvctu_y
double gravi_pfits_get_gvctu_y(const cpl_propertylist *plist)
Definition gravi_pfits.c:489
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_dit_acqcam
double gravi_pfits_get_dit_acqcam(const cpl_propertylist *plist)
Definition gravi_pfits.c:655
gravi_pfits_get_met_faint_params
cpl_vector * gravi_pfits_get_met_faint_params(const cpl_propertylist *plist)
Extract metrology faint settings.
Definition gravi_pfits.c:380
gravi_pfits_get_mjd
double gravi_pfits_get_mjd(const cpl_propertylist *plist)
Definition gravi_pfits.c:526
gravi_pfits_get_axis
int gravi_pfits_get_axis(const cpl_propertylist *plist)
Definition gravi_pfits.c:213
gravi_pfits_get_period_acqcam
double gravi_pfits_get_period_acqcam(const cpl_propertylist *plist)
Definition gravi_pfits.c:691
gravi_pfits_get_start_prcacq
const char * gravi_pfits_get_start_prcacq(const cpl_propertylist *plist)
Definition gravi_pfits.c:647
gravi_pfits_get_met_mode
int gravi_pfits_get_met_mode(const cpl_propertylist *plist)
Definition gravi_pfits.c:359
gravi_pfits_get_sobj_x
double gravi_pfits_get_sobj_x(const cpl_propertylist *plist)
Definition gravi_pfits.c:455
gravi_signal_create_sync
cpl_error_code gravi_signal_create_sync(cpl_table *vis_SC, int nbase_sc, double dit_sc, cpl_table *vis_FT, int nbase_ft, const char *name)
Compute synchronisation indices between OIFITS tables.
Definition gravi_signal.c:849
gravi_conf_get_telname
const char * gravi_conf_get_telname(int gravi_beam, cpl_propertylist *header)
Definition gravi_utils.c:1191
structTacConfiguration::min_allowed_snr_diff
double min_allowed_snr_diff[4][4][2]
Definition gravi_metrology.c:155
structTacConfiguration::min_allowed_flux_telescope
double min_allowed_flux_telescope[4][4][2]
Definition gravi_metrology.c:151
structTacConfiguration::lockin_constant_telescope
double lockin_constant_telescope[4][4][2][2]
Definition gravi_metrology.c:167
structTacConfiguration::min_allowed_flux_fiber_coupler
double min_allowed_flux_fiber_coupler[4][2]
Definition gravi_metrology.c:152
structTacConfiguration::min_allowed_snr_telescope
double min_allowed_snr_telescope[4][4][2]
Definition gravi_metrology.c:153
structTacConfiguration::offset_volt_fiber_coupler
double offset_volt_fiber_coupler[4][2][2]
Definition gravi_metrology.c:164
structTacConfiguration::rms_flux_fiber_coupler
double rms_flux_fiber_coupler[4][2]
Definition gravi_metrology.c:166
structTacConfiguration::offset_volt_telescope
double offset_volt_telescope[4][4][2][2]
Definition gravi_metrology.c:163
structTacConfiguration::min_allowed_snr_fiber_coupler
double min_allowed_snr_fiber_coupler[4][2]
Definition gravi_metrology.c:154
structTacConfiguration::rms_flux_telescope
double rms_flux_telescope[4][4][2]
Definition gravi_metrology.c:165
structTacConfiguration::lockin_constant_fiber_coupler
double lockin_constant_fiber_coupler[4][2][2]
Definition gravi_metrology.c:168
structTacData::flag_volt_telescope
int flag_volt_telescope[4][4][2]
Definition gravi_metrology.c:209
structTacData::delta_phasor_fiber_coupler
double delta_phasor_fiber_coupler[4][2][2]
Definition gravi_metrology.c:224
structTacData::sum_sq_flux_diff
double sum_sq_flux_diff[4][4][2]
Definition gravi_metrology.c:274
structTacData::buffer_delta_phasor_fiber_coupler
double buffer_delta_phasor_fiber_coupler[4][2][2][MET_MAX_HISTORY]
Definition gravi_metrology.c:223
structTacData::start_phase_fiber_coupler
double start_phase_fiber_coupler[4][2]
Definition gravi_metrology.c:290
structTacData::sum_volt_lockin_telescope
double sum_volt_lockin_telescope[4][4][2][2]
Definition gravi_metrology.c:187
structTacData::flag_flux_diff
int flag_flux_diff[4][4][2]
Definition gravi_metrology.c:268
structTacData::buffer_sq_flux_diff
double buffer_sq_flux_diff[4][4][2][MET_MAX_HISTORY]
Definition gravi_metrology.c:276
structTacData::flag_phasor0_diff
int flag_phasor0_diff[4][4][2]
Definition gravi_metrology.c:266
structTacData::raw_phase_telescope
double raw_phase_telescope[4][4][2]
Definition gravi_metrology.c:242
structTacData::raw_phase_diff
double raw_phase_diff[4][4][2]
Definition gravi_metrology.c:262
structTacData::rms_flux_telescope
double rms_flux_telescope[4][4][2]
Definition gravi_metrology.c:194
structTacData::flag_unwrap_telescope
int flag_unwrap_telescope[4][4][2]
Definition gravi_metrology.c:245
structTacData::buffer_delta_phasor_telescope
double buffer_delta_phasor_telescope[4][4][2][2][MET_MAX_HISTORY]
Definition gravi_metrology.c:196
structTacData::opl_zero_telescope_diode
double opl_zero_telescope_diode[4][4]
Definition gravi_metrology.c:287
structTacData::rms_flux_diff
double rms_flux_diff[4][4][2]
Definition gravi_metrology.c:277
structTacData::prev_phase_fiber_coupler
double prev_phase_fiber_coupler[4][2]
Definition gravi_metrology.c:253
structTacData::total_flag_telescope
int total_flag_telescope[4][4][2]
Definition gravi_metrology.c:280
structTacData::flag_snr_fiber_coupler
int flag_snr_fiber_coupler[4][2]
Definition gravi_metrology.c:240
structTacData::buffer_phasor_diff
double buffer_phasor_diff[4][4][2][2][MET_MAX_HISTORY]
Definition gravi_metrology.c:270
structTacData::buffer_sq_flux_telescope
double buffer_sq_flux_telescope[4][4][2][MET_MAX_HISTORY]
Definition gravi_metrology.c:193
structTacData::buffer_sq_delta_phase_fiber_coupler
double buffer_sq_delta_phase_fiber_coupler[4][2][MET_MAX_HISTORY]
Definition gravi_metrology.c:231
structTacData::sum_sq_delta_phase_telescope
double sum_sq_delta_phase_telescope[4][4][2]
Definition gravi_metrology.c:204
structTacData::volt_lockin_telescope
double volt_lockin_telescope[4][4][2][2]
Definition gravi_metrology.c:186
structTacData::freeze_speed_fiber_coupler
double freeze_speed_fiber_coupler[4][2]
Definition gravi_metrology.c:295
structTacData::sum_sq_flux_telescope
double sum_sq_flux_telescope[4][4][2]
Definition gravi_metrology.c:191
structTacData::sum_flux_telescope
double sum_flux_telescope[4][4][2]
Definition gravi_metrology.c:190
structTacData::flux_diff
double flux_diff[4][4][2]
Definition gravi_metrology.c:272
structTacData::unwrapped_phase_diff
double unwrapped_phase_diff[4][4][2]
Definition gravi_metrology.c:264
structTacData::snr_flux_fiber_coupler
double snr_flux_fiber_coupler[4][2]
Definition gravi_metrology.c:222
structTacData::best_estimate_phasor_fiber_coupler
double best_estimate_phasor_fiber_coupler[4][2][2]
Definition gravi_metrology.c:225
structTacData::flag_jump_fiber_coupler
int flag_jump_fiber_coupler[4][2]
Definition gravi_metrology.c:257
structTacData::opl_zero_telescope
double opl_zero_telescope[4]
Definition gravi_metrology.c:286
structTacData::buffer_used_speed_fiber_coupler
double buffer_used_speed_fiber_coupler[4][2][MET_MAX_HISTORY]
Definition gravi_metrology.c:236
structTacData::sum_speed_fiber_coupler
double sum_speed_fiber_coupler[4][2]
Definition gravi_metrology.c:234
structTacData::tacConfiguration
structTacConfiguration * tacConfiguration
Definition gravi_metrology.c:299
structTacData::flag_jump_telescope
int flag_jump_telescope[4][4][2]
Definition gravi_metrology.c:246
structTacData::flag_flux_fiber_coupler
int flag_flux_fiber_coupler[4][2]
Definition gravi_metrology.c:239
structTacData::flag_unwrap_diff
int flag_unwrap_diff[4][4][2]
Definition gravi_metrology.c:265
structTacData::delta_phasor_telescope
double delta_phasor_telescope[4][4][2][2]
Definition gravi_metrology.c:197
structTacData::flag_flux_telescope
int flag_flux_telescope[4][4][2]
Definition gravi_metrology.c:210
structTacData::opl_telescope
double opl_telescope[4]
Definition gravi_metrology.c:283
structTacData::freeze_phasor_fiber_coupler
double freeze_phasor_fiber_coupler[4][2][2]
Definition gravi_metrology.c:294
structTacData::flag_phasor0_fiber_coupler
int flag_phasor0_fiber_coupler[4][2]
Definition gravi_metrology.c:259
structTacData::sum_sq_flux_fiber_coupler
double sum_sq_flux_fiber_coupler[4][2]
Definition gravi_metrology.c:218
structTacData::raw_phase_fiber_coupler
double raw_phase_fiber_coupler[4][2]
Definition gravi_metrology.c:252
structTacData::opl_telescope_diode
double opl_telescope_diode[4][4]
Definition gravi_metrology.c:284
structTacData::buffer_flux_diff
double buffer_flux_diff[4][4][2][MET_MAX_HISTORY]
Definition gravi_metrology.c:275
structTacData::mean_phase_telescope
double mean_phase_telescope[4]
Definition gravi_metrology.c:250
structTacData::best_estimate_flux_telescope
double best_estimate_flux_telescope[4][4][2]
Definition gravi_metrology.c:199
structTacData::volt_lockin_fiber_coupler
double volt_lockin_fiber_coupler[4][2][2]
Definition gravi_metrology.c:213
structTacData::flag_volt_fiber_coupler
int flag_volt_fiber_coupler[4][2]
Definition gravi_metrology.c:238
structTacData::buffer_unwrapped_phase_fiber_coupler
double buffer_unwrapped_phase_fiber_coupler[4][2][MET_MAX_HISTORY]
Definition gravi_metrology.c:255
structTacData::sum_delta_phase_telescope
double sum_delta_phase_telescope[4][4][2]
Definition gravi_metrology.c:202
structTacData::sum_delta_phase_fiber_coupler
double sum_delta_phase_fiber_coupler[4][2]
Definition gravi_metrology.c:230
structTacData::best_estimate_phasor_telescope
double best_estimate_phasor_telescope[4][4][2][2]
Definition gravi_metrology.c:198
structTacData::unwrapped_phase_fiber_coupler
double unwrapped_phase_fiber_coupler[4][2]
Definition gravi_metrology.c:254
structTacData::rms_flux_fiber_coupler
double rms_flux_fiber_coupler[4][2]
Definition gravi_metrology.c:221
structTacData::sum_flux_fiber_coupler
double sum_flux_fiber_coupler[4][2]
Definition gravi_metrology.c:217
structTacData::buffer_speed_telescope
double buffer_speed_telescope[4][4][2][MET_MAX_HISTORY]
Definition gravi_metrology.c:207
structTacData::prev_phase_telescope
double prev_phase_telescope[4][4][2]
Definition gravi_metrology.c:243
structTacData::buffer_sq_delta_phase_telescope
double buffer_sq_delta_phase_telescope[4][4][2][MET_MAX_HISTORY]
Definition gravi_metrology.c:203
structTacData::buffer_flux_fiber_coupler
double buffer_flux_fiber_coupler[4][2][MET_MAX_HISTORY]
Definition gravi_metrology.c:219
structTacData::sum_sq_delta_phase_fiber_coupler
double sum_sq_delta_phase_fiber_coupler[4][2]
Definition gravi_metrology.c:232
structTacData::flag_snr_telescope
int flag_snr_telescope[4][4][2]
Definition gravi_metrology.c:211
structTacData::buffer_best_estimate_phasor_fiber_coupler
double buffer_best_estimate_phasor_fiber_coupler[4][2][2][MET_MAX_HISTORY]
Definition gravi_metrology.c:227
structTacData::freeze_count_fiber_coupler
long freeze_count_fiber_coupler[4][2]
Definition gravi_metrology.c:293
structTacData::opl_fiber_coupler
double opl_fiber_coupler[4]
Definition gravi_metrology.c:285
structTacData::buffer_delta_phase_telescope
double buffer_delta_phase_telescope[4][4][2][MET_MAX_HISTORY]
Definition gravi_metrology.c:201
structTacData::opl_zero_fiber_coupler
double opl_zero_fiber_coupler[4]
Definition gravi_metrology.c:288
structTacData::used_unwrapped_phase_fiber_coupler
double used_unwrapped_phase_fiber_coupler[4][2]
Definition gravi_metrology.c:297
structTacData::delta_phase_fiber_coupler
double delta_phase_fiber_coupler[4]
Definition gravi_metrology.c:260
structTacData::buffer_volt_lockin_fiber_coupler
double buffer_volt_lockin_fiber_coupler[4][2][2][MET_MAX_HISTORY]
Definition gravi_metrology.c:215
structTacData::sum_speed_telescope
double sum_speed_telescope[4][4][2]
Definition gravi_metrology.c:206
structTacData::sum_volt_lockin_fiber_coupler
double sum_volt_lockin_fiber_coupler[4][2][2]
Definition gravi_metrology.c:214
structTacData::buffer_speed_fiber_coupler
double buffer_speed_fiber_coupler[4][2][MET_MAX_HISTORY]
Definition gravi_metrology.c:235
structTacData::snr_flux_telescope
double snr_flux_telescope[4][4][2]
Definition gravi_metrology.c:195
structTacData::buffer_volt_lockin_telescope
double buffer_volt_lockin_telescope[4][4][2][2][MET_MAX_HISTORY]
Definition gravi_metrology.c:188
structTacData::best_estimate_flux_fiber_coupler
double best_estimate_flux_fiber_coupler[4][2]
Definition gravi_metrology.c:226
structTacData::sum_flux_diff
double sum_flux_diff[4][4][2]
Definition gravi_metrology.c:273
structTacData::flux_telescope
double flux_telescope[4][4][2]
Definition gravi_metrology.c:189
structTacData::buffer_flux_telescope
double buffer_flux_telescope[4][4][2][MET_MAX_HISTORY]
Definition gravi_metrology.c:192
structTacData::sum_phasor_diff
double sum_phasor_diff[4][4][2][2]
Definition gravi_metrology.c:271
structTacData::snr_flux_diff
double snr_flux_diff[4][4][2]
Definition gravi_metrology.c:278
structTacData::buffer_sq_flux_fiber_coupler
double buffer_sq_flux_fiber_coupler[4][2][MET_MAX_HISTORY]
Definition gravi_metrology.c:220
structTacData::flag_phasor0_telescope
int flag_phasor0_telescope[4][4][2]
Definition gravi_metrology.c:248
structTacData::flag_speed_telescope
int flag_speed_telescope[4][4][2]
Definition gravi_metrology.c:247
structTacData::delta_phase_telescope
double delta_phase_telescope[4][4]
Definition gravi_metrology.c:249
structTacData::buffer_delta_phase_fiber_coupler
double buffer_delta_phase_fiber_coupler[4][2][MET_MAX_HISTORY]
Definition gravi_metrology.c:229
structTacData::start_phase_telescope
double start_phase_telescope[4][4][2]
Definition gravi_metrology.c:291
structTacData::flag_speed_fiber_coupler
int flag_speed_fiber_coupler[4][2]
Definition gravi_metrology.c:258
structTacData::freeze_unwrapped_phase_fiber_coupler
double freeze_unwrapped_phase_fiber_coupler[4][2]
Definition gravi_metrology.c:296
structTacData::prev_phase_diff
double prev_phase_diff[4][4][2]
Definition gravi_metrology.c:263
structTacData::unwrapped_phase_telescope
double unwrapped_phase_telescope[4][4][2]
Definition gravi_metrology.c:244
structTacData::flux_fiber_coupler
double flux_fiber_coupler[4][2]
Definition gravi_metrology.c:216
structTacData::flag_snr_diff
int flag_snr_diff[4][4][2]
Definition gravi_metrology.c:267
structTacData::total_flag_fiber_coupler
int total_flag_fiber_coupler[4][2]
Definition gravi_metrology.c:281
structTacData::flag_unwrap_fiber_coupler
int flag_unwrap_fiber_coupler[4][2]
Definition gravi_metrology.c:256
Generated on Fri Apr 24 2026 17:14:45 for GRAVI Pipeline Reference Manual by doxygen 1.13.2