GRAVI Pipeline Reference Manual 1.10.1
Loading...
Searching...
No Matches
gravi_old.c
Go to the documentation of this file.
1
2/* $Id: gravi_calib.c,v 1.10 2012/03/23 15:10:40 nazouaoui Exp $
3 *
4 * This file is part of the GRAVI Pipeline
5 * Copyright (C) 2002,2003 European Southern Observatory
6 *
7 * This program is free software; you can redistribute it and/or modify
8 * it under the terms of the GNU General Public License as published by
9 * the Free Software Foundation; either version 2 of the License, or
10 * (at your option) any later version.
11 *
12 * This program is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with this program; if not, write to the Free Software
19 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 */
21
22
23/*
24 * This file contains outdated function, not compiled
25 * nor included in H files.
26 */
27
28
29int gravi_get_base_from_tel (int t1, int t2)
30{
31 for ( int i=0 ; i<6; i++ ) {
32 if ( GRAVI_BASE_TEL[i][0] == t1 && GRAVI_BASE_TEL[i][1] == t2 ) return i;
33 if ( GRAVI_BASE_TEL[i][0] == t2 && GRAVI_BASE_TEL[i][1] == t1 ) return i;
34 }
35 return 99;
36}
37
38int gravi_get_basesign_from_tel (int t1, int t2)
39{
40 for ( int i=0 ; i<6; i++ ) {
41 if ( GRAVI_BASE_TEL[i][0] == t1 && GRAVI_BASE_TEL[i][1] == t2 ) return 1;
42 if ( GRAVI_BASE_TEL[i][0] == t2 && GRAVI_BASE_TEL[i][1] == t1 ) return -1;
43 }
44 return 99;
45}
46
47cpl_error_code gravi_lazer_get_wavelength (gravi_data * lazer_data){
48 gravi_msg_function_start(1);
49
50 cpl_table * spectrum_table, * detector_table, * wavelength_table;
51 char * polarisation, * data_x;
52 const char * regname;
53 cpl_propertylist * primary_hdr;
54 int type_data, pol, npol, nregion, reg, nwave, size, wave, nv;
55
56 primary_hdr = gravi_data_get_header (lazer_data);
57
58 for (type_data = 0; type_data < 2; type_data ++) {
59
60 /* Extract the IMAGING_DETECTOR table, OI_WAVELENGTH property list and the
61 * primary header */
62
63 if (type_data == 0){
64 detector_table = gravi_data_get_table (lazer_data,
65 GRAVI_IMAGING_DETECTOR_SC_EXT);
66 /* Find the spectrum data field */
67 spectrum_table = gravi_data_get_table (lazer_data,
68 GRAVI_SPECTRUM_DATA_SC_EXT);
69 }
70 else {
71 detector_table = gravi_data_get_table (lazer_data,
72 GRAVI_IMAGING_DETECTOR_FT_EXT);
73
74 /* Find the spectrum data field */
75 spectrum_table = gravi_data_get_table (lazer_data,
76 GRAVI_SPECTRUM_DATA_FT_EXT);
77 }
78
79 if ((spectrum_table == NULL) || (primary_hdr == NULL) ||
80 (wavelength_table == NULL) || (detector_table == NULL)){
81
82 return cpl_error_set_message(cpl_func, CPL_ERROR_ILLEGAL_OUTPUT,
83 "The preproc data must contrain "
84 "SPECTRUM_DATA field");
85 }
86
87
88 /* Get the number of wavelength and the region */
89 nwave = cpl_table_get_column_dimension(spectrum_table, "DATA1", 1);
90 size = cpl_table_get_column_dimension(spectrum_table, "DATA1", 0);
91
92 nregion = cpl_table_get_nrow(detector_table);
93 if (nregion > 24)
94 npol = 2;
95 else
96 npol = 1;
97
98
99 wavelength_table = gravi_data_get_oi_table(lazer_data,
100 GRAVI_OI_WAVE_EXT(type_data),
101 GRAVI_INSNAME(type_data,0,npol));
102
103 for (pol = 0; pol < npol; pol++){
104
105 if (pol == 0)
106 if (npol == 2)
107 polarisation = POLAR_2;
108 else {
109 polarisation = POLAR_3;
110 pol = 1;
111 }
112 else
113 polarisation = POLAR_1;
114
115 cpl_image * mean_pol = cpl_image_new(size, nwave, CPL_TYPE_DOUBLE);
116 cpl_image_fill_window (mean_pol, 1, 1, size, nwave, 0.0);
117
118 for (reg = 0; reg < nregion; reg ++){
119
120 regname = cpl_table_get_string (detector_table, "REGANAME", reg);
121
122 if (!strstr(regname, polarisation))
123 continue;
124
125
126 data_x = cpl_sprintf("DATA%d", reg+1);
127
128 cpl_imagelist * spec_data = gravi_table_data_to_imagelist(
129 spectrum_table, reg+1);
130
131 cpl_image * mean_spec = cpl_imagelist_collapse_create(spec_data);
132
133 cpl_image_add (mean_pol, mean_spec);
134 cpl_imagelist_delete (spec_data);
135 cpl_image_delete (mean_spec);
136
137 }
138
139 int ind = 0;
140 double val_spec = cpl_image_get (mean_pol, 1, 1, &nv);
141 for (wave = 1; wave < nwave; wave ++){
142 if (val_spec < cpl_image_get (mean_pol, 1, wave, &nv)){
143 val_spec = cpl_image_get (mean_pol, 1, wave, &nv);
144 ind = wave;
145 }
146
147 }
148
149 float wave_lazer = cpl_table_get_float (wavelength_table,
150 "EFF_WAVE", ind, &nv);
151 char * key = cpl_sprintf("LASER %s %s", GRAVI_TYPE(type_data), polarisation);
152 cpl_propertylist_append_float (primary_hdr, key, wave_lazer);
153 cpl_free (key);
154
155 cpl_image_delete (mean_pol);
156 }
157 }
158
159 /* Verbose */
160 gravi_msg_function_exit(1);
161 return CPL_ERROR_NONE;
162}
163
164
165gravi_data * gravi_compute_disp_old (gravi_data * vis_data){
166
167 cpl_table * vistable_sc, * wavelenght_sc, * fddl_met;
168 cpl_propertylist * plist, * vis_FT_p, * vis_SC_p, * wavesc_plist;
169 gravi_data * disp_data;
170 int wave, nwave_sc, base, nbase = 6, row , nfile, nv;
171 cpl_array * phase_sc_base;
172 const cpl_array * phase_sc,
173 * phase_met1, * phase_met2, * fddl_sc, * fddl_ft;
174 cpl_vector * wave_sc;
175 cpl_bivector * fout, * fref;
176 int tel_1[6] = {0,0,0,1,1,2}, pol = 0;
177 int tel_2[6] = {1,2,3,2,3,3}, size_matrix, size_opd, npol;
178 double pos_sc, pos_ft, pos_mean, ph_sc, opd_sc, opd_met,
179 ph_met1, ph_met2;
180 double time_sc, time_ft, exptime_sc;
181 gsl_matrix * phase_disp;
182 gsl_vector ** opd_disp;
183 char * pol_FT_check[2], * pol_SC_check[2];
184
185 /* Verbose */
186 gravi_msg_function_start(1);
187
188 vis_FT_p = gravi_data_get_plist (vis_data, GRAVI_OI_VIS_FT_EXT);
189 vis_SC_p = gravi_data_get_plist (vis_data, GRAVI_OI_VIS_SC_EXT);
190
191 /* Check if there is 2 polarization */
192 if (!(strcmp(gravi_pfits_get_insname (vis_FT_p), INSNAME_FT_P1)) &&
193 (strcmp(gravi_pfits_get_insname (vis_FT_p), INSNAME_FT_P2))){
194 pol_FT_check[0] = INSNAME_FT_P1;
195 pol_FT_check[1] = INSNAME_FT_P2;
196
197 npol = 2;
198// polFT_test = 1;
199 }
200 else if (!strcmp(gravi_pfits_get_insname (vis_FT_p), INSNAME_FT)){
201
202 pol_FT_check[0] = INSNAME_FT;
203 pol_FT_check[1] = INSNAME_FT;
204 npol = 1;
205 }
206
207 if (!(strcmp(gravi_pfits_get_insname (vis_SC_p), INSNAME_SC_P1)) &&
208 (strcmp(gravi_pfits_get_insname (vis_SC_p), INSNAME_SC_P2))){
209 pol_SC_check[0] = INSNAME_SC_P1;
210 pol_SC_check[1] = INSNAME_SC_P2;
211 npol = 2;
212// polSC_test = 1;
213 }
214 else if (!strcmp(gravi_pfits_get_insname (vis_SC_p), INSNAME_SC)){
215 pol_SC_check[0] = INSNAME_SC;
216 pol_SC_check[1] = INSNAME_SC;
217 npol = 1;
218 }
219
220 /* Get the FDDL_MET table and load the matrix fddl matrix*/
221 fddl_met = gravi_data_get_table(vis_data, "FDDL_MET_MEAN");
222 nfile=cpl_table_get_nrow(fddl_met);
223
224 disp_data = gravi_data_new(0);
225
226 gravi_data_append_header (disp_data, gravi_data_get_header (vis_data));
227
228
229// cpl_table_delete (p2vm_met);
230 for (pol = 0; pol < npol; pol ++){
231 vistable_sc = gravi_data_get_oi_table (vis_data,
232 GRAVI_OI_VIS_SC_EXT, pol_SC_check[pol]);
233 wavelenght_sc = gravi_data_get_oi_table(vis_data,
234 GRAVI_OI_WAVELENGTH_SC_EXT, pol_SC_check[pol]);
235 wavesc_plist = gravi_data_get_oi_plist(vis_data,
236 GRAVI_OI_WAVELENGTH_SC_EXT, pol_SC_check[pol]);
237 nwave_sc = gravi_pfits_get_nwave (wavesc_plist);
238 wave_sc = cpl_vector_new(nwave_sc);
239
240 plist = gravi_data_get_header (vis_data);
241
242 // not used ? exptime_sc = gravi_pfits_get_dit (plist, GRAVI_SC)*pow(10, 6);
243 // cpl_matrix * matrix = cpl_matrix_new (nbase, ncol);
244 // cpl_matrix_fill (matrix, 0.0);
245 // for (i = 0; i < nbase; i ++){
246 // cpl_matrix_set (matrix, i, tel_1[i], 1);
247 // cpl_matrix_set (matrix, i, tel_2[i], -1);
248 // cpl_matrix_set (matrix, i, 4 + tel_1[i], 1);
249 // cpl_matrix_set (matrix, i, 4 + tel_2[i], -1);
250 // }
251
252 nfile = cpl_table_get_nrow (vistable_sc) / 6;
253
254 // gsl_vector * bis;
255 // gsl_matrix * A;
256
257 phase_disp = gsl_matrix_calloc (nfile * 6, 8);
258
259 opd_disp = cpl_malloc(nwave_sc * sizeof(gsl_vector *));
260 for (wave = 0; wave < nwave_sc; wave ++){
261
262 cpl_vector_set (wave_sc, wave,
263 cpl_table_get_float (wavelenght_sc, "EFF_WAVE", wave, &nv));
264 opd_disp[wave] = gsl_vector_alloc (nfile * 6);
265 }
266
267 /* Get the mean of the metrology phase */
268 cpl_array * met_phase = cpl_array_new(4, CPL_TYPE_DOUBLE);
269 cpl_array_fill_window (met_phase, 0, 4, 0.0);
270 for (row = 0; row < cpl_table_get_nrow(fddl_met); row ++){
271 cpl_array_add (met_phase, cpl_table_get_array(fddl_met, "PHASE_FC", row));
272 }
273
274 cpl_array_divide_scalar(met_phase, cpl_table_get_nrow(fddl_met));
275
276 for (base = 0; base < 6; base++){
277
278 ph_met1 = cpl_array_get_double (met_phase, tel_1[base], &nv);
279 ph_met2 = cpl_array_get_double (met_phase, tel_2[base], &nv);
280 opd_met = (ph_met1 - ph_met2) * 2 * M_PI * LAMBDA_MET;
281
282 /* Construction of each mean complex visibilities and the mean
283 * squared visibilities array */
284
285
286 /* Compute the average phase of FT to evaluate the phase between
287 * the interval [n DIT, (n+1)DIT] */
288 for (row = 0; row < nfile; row ++){
289
290 // printf("time_SC = %d time_FT = %d\n", time_SC, time_FT);
291
292 phase_sc = cpl_table_get_array (vistable_sc, "VISPHI", row * 6 + base);
293
294 /* Get the visibilities product by the reduced p2vm */
295
296 /* Compute the opd dispertion */
297
298 /* Interpolate the phase ft */
299 phase_sc_base = cpl_array_duplicate (phase_sc);
300
301 if (cpl_error_get_code()){
302 cpl_bivector_unwrap_vectors(fout);
303 cpl_bivector_unwrap_vectors(fref);
304 cpl_error_set_message(cpl_func, CPL_ERROR_NULL_INPUT,
305 "Error during the interpolation of phase ft");
306 return NULL;
307 }
308
309 /*Construction of the phase matrix */
310 fddl_sc = cpl_table_get_array (fddl_met, "FT_POS", row);
311 fddl_ft = cpl_table_get_array (fddl_met, "SC_POS", row);
312 pos_sc = cpl_array_get_float (fddl_sc, tel_1[base], &nv);
313 pos_ft = cpl_array_get_float (fddl_ft, tel_1[base], &nv);
314
315 pos_mean = (pos_sc + pos_ft) / 2;
316 gsl_matrix_set (phase_disp, row * nbase + base, tel_1[base], ph_met1);
317 gsl_matrix_set (phase_disp, row * nbase + base, tel_2[base], ph_met2);
318 gsl_matrix_set (phase_disp, row * nbase + base, 4 + tel_1[base], pos_mean);
319
320 pos_sc = cpl_array_get_float (fddl_sc, tel_2[base], &nv);
321 pos_ft = cpl_array_get_float (fddl_ft, tel_2[base], &nv);
322
323 pos_mean = (pos_sc + pos_ft) / 2;
324
325 gsl_matrix_set (phase_disp, row * nbase + base, 4 + tel_2[base], pos_mean);
326
327 /* Construction of the opd vector */
328 for (wave = 0; wave < nwave_sc; wave++){
329 ph_sc = cpl_array_get_double (phase_sc_base, wave, &nv);
330 opd_sc = ph_sc * cpl_vector_get (wave_sc, wave)/ (2 * M_PI);
331 // printf("cpl_vector_get (wave_sc, wave) = %e\n", cpl_vector_get (wave_sc, wave));
332 // printf("wave = %d\n", wave);
333 // printf("ph_sc = %e opd_sc \n", ph_sc,opd_sc);
334 gsl_vector_set (opd_disp[wave], base + nbase * row,
335 opd_sc);
336
337
338
339 }
340 if (cpl_error_get_code()){
341 cpl_bivector_unwrap_vectors(fout);
342 cpl_bivector_unwrap_vectors(fref);
343 cpl_error_set_message(cpl_func, CPL_ERROR_NULL_INPUT,
344 "Error during the construction of the opd vector");
345 return NULL;
346 }
347
348 cpl_array_delete (phase_sc_base);
349 }
350
351
352 }
353 cpl_vector_unwrap (wave_sc);
354 cpl_array_delete (met_phase);
355 cpl_table * disp_table = cpl_table_new (nwave_sc);
356 cpl_table_new_column_array (disp_table, "COEFF",
357 CPL_TYPE_DOUBLE, 8);
358 cpl_array * coeff;
359
360 gsl_vector * S;
361 gsl_vector * work;
362 gsl_vector * x;
363 gsl_matrix * u = gsl_matrix_calloc (nfile * 6, 8), * V;
364 int i;
365
366 for (wave = 0; wave < nwave_sc; wave ++){
367
368 gsl_matrix_memcpy(u, phase_disp);
369 x = gsl_vector_alloc (8);
370 work = gsl_vector_alloc (8);
371 S = gsl_vector_alloc (8);
372 coeff = cpl_array_new (8, CPL_TYPE_DOUBLE);
373 V = gsl_matrix_alloc (8, 8);
374
375 gsl_linalg_SV_decomp (u, V, S, work);
376 gsl_linalg_SV_solve (u, V, S, opd_disp[wave], x);
377
378 for (i = 0; i < 8; i ++) {
379 cpl_array_set_double (coeff, i, gsl_vector_get (x, i));
380 }
381
382 cpl_table_set_array (disp_table, "COEFF", wave, coeff);
383
384 cpl_array_delete (coeff);
385 gsl_vector_free (opd_disp[wave]);
386 gsl_matrix_free (V);
387 gsl_vector_free (x);
388 gsl_vector_free (work);
389 gsl_vector_free (S);
390 }
391 gsl_matrix_free (u);
392 gsl_matrix_free (phase_disp);
393 cpl_free (opd_disp);
394
395
396 cpl_propertylist * met_plist = cpl_propertylist_duplicate (wavesc_plist);
397 char * name_tb = cpl_sprintf("DISP_%s", pol_SC_check[pol]);
398 cpl_propertylist_append_string (met_plist, "EXTNAME",
399 name_tb);
400 cpl_free (name_tb);
401 gravi_data_add (disp_data, met_plist, disp_table);
402 cpl_propertylist_delete (met_plist);
403 } // End loop on polarisation
404
405
406 gravi_msg_function_exit(1);
407 return disp_data;
408}
409
410
411
412/*----------------------------------------------------------------------------*/
428/*----------------------------------------------------------------------------*/
429
430gravi_data * gravi_fit_argon (gravi_data * argon_data, gravi_data * wave_data,
431 gravi_data * profile_map, gravi_data * dark_map,
432 gravi_data * bad_map) {
433 gravi_msg_function_start(1);
434
435 cpl_propertylist * primary_hdr, * plist;
436 cpl_table * spectrum_table, * imaging_detector, * wave_data_sc,
437 * img_output;
438 cpl_matrix * all_coord;
439 cpl_array * wavelength;
440 cpl_polynomial * fit2d;
441 cpl_vector * residuals;
442 double rechisq;
443 // Degrees of wavelength fit polynomial {2 in vertical direction, 4 in disp direction}
444 const cpl_size deg2d[2] = {2, 2};
445 const cpl_size deg1d[1] = {2};
446 cpl_vector * coord_X, * coord_Y, * all_wavelength;
447 double result;
448 int wave;
449 double slope;
450 int y_corner;
451 int nwave, n_region, region;
452 char * data_x;
453 int nv;
454 cpl_bivector * plot;
455 cpl_vector *pos;
456 double minwave = 0, maxwave = 1e10;
457 cpl_image * img_profile, * image_wave, * profile_image;
458 cpl_imagelist * imglist_wave;
459 cpl_table * profile_table;
460 cpl_table * wave_fibre;
461 int sizey, sizex, ind;
462
463 /* Check the inputs data */
464 if (argon_data == NULL){
465 cpl_error_set_message(cpl_func, CPL_ERROR_NULL_INPUT, "The data is NULL");
466 return NULL;
467 }
468
469 clock_t start, end;
470 start = clock();
471 cpl_imagelist * imglist = gravi_data_get_cube (argon_data, GRAVI_IMAGING_DATA_SC_EXT);
472 cpl_image * img_median = cpl_imagelist_collapse_median_create (imglist);
473
474 cpl_imagelist * imglist_med = cpl_imagelist_new ();
475 cpl_imagelist_set (imglist_med, img_median, 0);
476
477 /* Put the data in the output table : dark_map */
478 gravi_data_set_cube (argon_data, GRAVI_IMAGING_DATA_SC_EXT, imglist_med);
479
480 gravi_data * output_data = gravi_extract_spectrum(argon_data, profile_map,
481 dark_map, bad_map, NULL);
482
483 end = clock();
484 cpl_msg_info (cpl_func, "Execution time gravi_extract_spectrum : %f", (end - start) / (double)CLOCKS_PER_SEC);
485
486 CPLCHECK_NUL("Connot extract the spectrum");
487
488 /* Wavelengths of the argon emission lines (microns) */
489 int nlines = 10;
490 double line_wave[] = {/*1.982291,*/
491 1.997118e-6,
492 2.032256e-6,
493 2.062186e-6,
494 /*2.065277e-6,*/
495 /*2.073922e-6,*/
496 /*2.081672e-6,*/
497 2.099184e-6,
498 2.133871e-6,
499 2.154009e-6,
500 2.208321e-6,
501 2.313952e-6,
502 2.385154e-6,
503 2.397306e-6};
504
505 slope = 0.02;
506 /* Correction of the wavelength scale in glass to go to vaccuum scale */
507 imaging_detector = gravi_data_get_table(output_data, GRAVI_IMAGING_DETECTOR_SC_EXT);
508 spectrum_table = gravi_data_get_table (output_data, GRAVI_SPECTRUM_DATA_SC_EXT);
509
510 if ((spectrum_table == NULL) || (imaging_detector == NULL)) {
511 gravi_data_delete(output_data);
512 cpl_error_set_message(cpl_func,
513 CPL_ERROR_ILLEGAL_OUTPUT, "Data must contain SPECTRUM_DATA");
514 return NULL;
515 }
516
517 wave_data_sc = cpl_table_duplicate (gravi_data_get_table (wave_data, GRAVI_WAVE_DATA_SC_EXT));
518 double wavestep = cpl_array_get_double (cpl_table_get_array (wave_data_sc, "DATA1", 0), 1, &nv) -
519 cpl_array_get_double (cpl_table_get_array (wave_data_sc, "DATA1", 0), 0, &nv);
520
521 cpl_array ** wave_array;
522 double * wave_ ;
523// for (region = 0; region < n_region; region ++){
524//
525// data_x = cpl_sprintf("DATA%d", region + 1);
526// wave_array = cpl_table_get_data_array (wave_data_sc, data_x);
527// wave_ = cpl_array_get_data_double (wave_array[0]);
528// for (wave = 0; wave < nwave; wave ++){
529// wave_ [wave] -= wave * wavestep * slope;
530// }
531//
532// cpl_free (data_x);
533// }
534
535
536 /* Get the number of wavelength and the region */
537 nwave = cpl_table_get_column_dimension(spectrum_table, "DATA1", 1);
538 n_region = cpl_table_get_nrow(imaging_detector);
539
540 /* Loop region */
541 coord_X = cpl_vector_new(n_region * nlines);
542 coord_Y = cpl_vector_new(n_region * nlines);
543 cpl_vector * coord_X_fit = cpl_vector_new(n_region * nlines);
544 all_wavelength = cpl_vector_new(n_region * nlines);
545 cpl_vector * fitsigm = cpl_vector_new(n_region * nlines);
546
547 cpl_vector_fill (coord_X, -1);
548 cpl_vector_fill (coord_Y, -1);
549 cpl_vector_fill (all_wavelength, -1);
550
551 int fitwidth, fit_in, comp = 0, list;
552
553 primary_hdr = gravi_data_get_header (argon_data);
554 const char * resolution = gravi_pfits_get_resolution (primary_hdr);
555
556 if (! (strcmp(resolution, "LOW") && strcmp(resolution, "MED")) )
557 fitwidth = 3;
558 else
559 fitwidth = 10;
560
561 /* Fit of each emission line position using a quadratic model of their
562 * position as a function of channel number */
563 const cpl_array * argon;
564
565 for (region = 0; region < n_region; region ++) {
566
567 data_x = cpl_sprintf("DATA%d", region + 1);
568 wave_array = cpl_table_get_data_array (wave_data_sc, data_x);
569 wave_ = cpl_array_get_data_double (wave_array[0]);
570
571 for (wave = 0; wave < nwave; wave ++){
572 wave_ [wave] -= wave * wavestep * slope;
573 }
574
575 argon = cpl_table_get_array (spectrum_table, data_x, 0);
576
577 cpl_free (data_x);
578
579 wave = 0;
580 for (list = 0; list < nlines; list ++) {
581
582 while (wave_[wave] < line_wave[list]){
583 wave ++;
584 }
585
586// wave += 20;
587
588 if (wave >= nwave){
589 cpl_error_set_message(cpl_func,
590 CPL_ERROR_ILLEGAL_OUTPUT, "The argon wavelength does "
591 "not much with the calibration wavelength ");
592 return NULL;
593 }
594
595 double flux = 0;
596 double coord_flux = 0, coord_err = 0;
597 cpl_vector * vector_x=cpl_vector_new(fitwidth*2);
598 cpl_vector * vector_y=cpl_vector_new(fitwidth*2);
599 int i=0;
600 for (fit_in = wave - fitwidth; fit_in < wave + fitwidth; fit_in ++){
601 cpl_vector_set(vector_x, i, fit_in);
602 cpl_vector_set(vector_y, i, cpl_array_get_double (argon, fit_in, &nv));
603 i++;
604 flux += cpl_array_get_double (argon, fit_in, &nv);
605 coord_flux += fit_in * cpl_array_get_double (argon, fit_in, &nv);
606// printf("region %d; line %d; %d, flux= %g \n", region, list, fit_in, cpl_array_get_double (argon, fit_in, &nv));
607 }
608
609 cpl_errorstate prestate = cpl_errorstate_get();
610 double x0, sigma, area, offset, mse;
611 cpl_vector_fit_gaussian (vector_x, NULL, vector_y, NULL,
612 CPL_FIT_ALL, &x0, &sigma, &area,
613 &offset, &mse, NULL, NULL);
614
615 coord_flux = (coord_flux / flux);
616 coord_err = fabs (flux);
617// printf("Barycentre = %g Gaussian FIT :%g (+-%g) \n", coord_flux, x0, sigma);
618
619 if (cpl_error_get_code() == CPL_ERROR_CONTINUE){
620 cpl_errorstate_set (prestate);
621 cpl_msg_warning(cpl_func, "The gaussian fit did not converge while fitting argon line");
622 x0=coord_flux;
623 sigma=100;
624
625 }
626
627 CPLCHECK_NUL("Error during the computation");
628
629// cpl_vector_set (coord_X, list * n_region + region, coord_flux);
630 cpl_vector_set (coord_X, list * n_region + region, x0);
631 cpl_vector_set (coord_Y, list * n_region + region, region + 1);
632// cpl_vector_set (fitsigm, list * n_region + region, coord_err);
633 cpl_vector_set (fitsigm, list * n_region + region, sigma);
634 cpl_vector_set (all_wavelength, list * n_region + region, line_wave[list]);
635
636 comp++;
637 } /* End loop on list of lines */
638
639 } /* End loop on regions */
640
641 cpl_table_delete (wave_data_sc);
642 int i;
643
644 /* Fit of the barycenters of each argon wavelength of all the regions */
645 cpl_vector * coord_x_list, * fitsigm_list,
646 * coord_y_list, * coord_x_fit = cpl_vector_new (n_region);
647
648 cpl_matrix * matrix3, * matrix2;
649
650 gsl_matrix * X = gsl_matrix_alloc (n_region, 3);
651 gsl_matrix * X_bis = gsl_matrix_alloc (n_region, 3);
652 gsl_vector * y = gsl_vector_alloc (n_region);
653 gsl_vector * w = gsl_vector_alloc (n_region);
654
655 gsl_vector * c = gsl_vector_alloc (3);
656 gsl_matrix * cov = gsl_matrix_alloc (3, 3);
657 double chisq, result_gsl;
658
659 /* Loop on the list of lines */
660 for (list = 0; list < nlines; list++) {
661
662 coord_x_list = cpl_vector_extract (coord_X, list * n_region, (list + 1) * n_region - 1, 1);
663 fitsigm_list = cpl_vector_extract (fitsigm, list * n_region, (list + 1) * n_region - 1, 1);
664 coord_y_list = cpl_vector_extract (coord_Y, list * n_region, (list + 1) * n_region - 1, 1);
665
666 for (i = 0; i < n_region; i++){
667 gsl_matrix_set (X, i, 0, 1.0);
668 gsl_matrix_set (X, i, 1, cpl_vector_get (coord_y_list, i));
669 gsl_matrix_set (X, i, 2, cpl_vector_get (coord_y_list, i) *
670 cpl_vector_get (coord_y_list, i));
671
672 gsl_vector_set (y, i, cpl_vector_get (coord_x_list, i));
673 gsl_vector_set (w, i, cpl_vector_get (fitsigm_list, i));
674 }
675
676 gsl_matrix_memcpy (X_bis, X);
677 gsl_multifit_linear_workspace * work = gsl_multifit_linear_alloc (n_region, 3);
678 gsl_multifit_wlinear (X_bis, w, y, c, cov,
679 &chisq, work);
680 gsl_multifit_linear_free (work);
681
682 CPLCHECK_NUL("Cannot compute the barycneter of all regions");
683
684 for (region = 0; region < n_region; region ++){
685
686 result_gsl = gsl_vector_get (c, 0) + gsl_vector_get (c, 1) *
687 gsl_matrix_get (X, region, 1) + gsl_vector_get (c, 2) *
688 gsl_matrix_get (X, region, 2);
689 cpl_vector_set (coord_X_fit, list * n_region + region, result_gsl);
690 cpl_vector_set (coord_x_fit, region, result_gsl);
691 }
692
693 cpl_vector_delete (coord_x_list);
694 cpl_vector_delete (fitsigm_list);
695 cpl_vector_delete (coord_y_list);
696
697 } /* End loop on the list of lines */
698
699 gsl_matrix_free (cov);
700 gsl_matrix_free (X);
701 gsl_matrix_free (X_bis);
702 gsl_vector_free (y);
703 gsl_vector_free (w);
704 gsl_vector_free (c);
705
706 cpl_vector_delete (coord_x_fit);
707
708 /* Compute the polynomial model between the position and the
709 * wavelength = F(x, y) where F is polynomial function */
710
711 fit2d = cpl_polynomial_new(2);
712
713 matrix3 = cpl_matrix_wrap(1, nlines * n_region, cpl_vector_get_data(coord_X));
714 matrix2 = cpl_matrix_wrap(1, nlines * n_region, cpl_vector_get_data(coord_Y));
715
716 all_coord = cpl_matrix_duplicate (matrix3) ;
717 cpl_matrix_append(all_coord, matrix2, 1);
718
719 /* fit CPL */
720 cpl_polynomial_fit(fit2d, all_coord, NULL, all_wavelength, NULL,
721 CPL_TRUE, NULL, deg2d);
722 residuals = cpl_vector_new(nlines*n_region);
723 cpl_vector_fill_polynomial_fit_residual (residuals, all_wavelength, NULL,
724 fit2d, all_coord, &rechisq );
725
726 /* fits gsl */
727 X = gsl_matrix_alloc (n_region*nlines, 6);
728 X_bis = gsl_matrix_alloc (n_region*nlines, 6);
729 y = gsl_vector_alloc (n_region*nlines);
730 w = gsl_vector_alloc (n_region*nlines);
731
732 c = gsl_vector_alloc (6);
733 cov = gsl_matrix_alloc (6, 6);
734
735 for (list = 0; list < nlines; list++) {
736
737 coord_x_list = cpl_vector_extract (coord_X, list * n_region, (list + 1) * n_region - 1, 1);
738 fitsigm_list = cpl_vector_extract (fitsigm, list * n_region, (list + 1) * n_region - 1, 1);
739 coord_y_list = cpl_vector_extract (coord_Y, list * n_region, (list + 1) * n_region - 1, 1);
740
741 for (i = 0; i < n_region; i++){
742 gsl_matrix_set (X, i+list*n_region, 0, 1.0);
743 gsl_matrix_set (X, i+list*n_region, 1, cpl_vector_get (coord_x_list, i));
744 gsl_matrix_set (X, i+list*n_region, 2, cpl_vector_get (coord_x_list, i) *
745 cpl_vector_get (coord_x_list, i));
746 gsl_matrix_set (X, i+list*n_region, 3, cpl_vector_get (coord_y_list, i));
747 gsl_matrix_set (X, i+list*n_region, 4, cpl_vector_get (coord_y_list, i) *
748 cpl_vector_get (coord_y_list, i));
749 gsl_matrix_set (X, i+list*n_region, 5, cpl_vector_get (coord_x_list, i) *
750 cpl_vector_get (coord_y_list, i));
751
752 gsl_vector_set (y, i+list*n_region, line_wave[list]);
753 gsl_vector_set (w, i+list*n_region, 1./pow(cpl_vector_get (fitsigm_list, i),2));
754 }
755 }
756
757 gsl_multifit_linear_workspace * work = gsl_multifit_linear_alloc (n_region*nlines, 6);
758 gsl_matrix_memcpy (X_bis, X);
759 gsl_multifit_wlinear (X_bis, w, y, c, cov,
760 &chisq, work);
761 gsl_vector * r=gsl_vector_alloc (n_region*nlines);
762 gsl_multifit_linear_residuals (X_bis, y, c, r);
763 gsl_multifit_linear_free (work);
764
765 for (list = 0; list < nlines; list++) {
766 for (i = 0; i < n_region; i++){
767// printf("line %d region %d residual %g\n", list, i, gsl_vector_get(r, i+list*n_region));
768 }
769 }
770
771
772
773 CPLCHECK_NUL("Error in the fit of the wavelength and position");
774
775
776// if (PLOT_WAVELENGTH)
777// {
778// vectors=malloc(3 * sizeof(cpl_vector*));
779// vectors[0]=NULL;
780// vectors[1]=all_wavelength;
781// cpl_vector *XYpos=cpl_vector_duplicate(coord_X);
782// cpl_vector_multiply_scalar(XYpos, nlines/5);
783// cpl_vector_add(XYpos, coord_Y);
784// cpl_vector_divide_scalar(XYpos, nlines/5);
785// vectors[0]=XYpos;
786// cpl_bivector *error=cpl_bivector_wrap_vectors( XYpos, cpl_vector_duplicate(residuals));
787// cpl_vector_subtract(residuals, all_wavelength);
788// cpl_vector_multiply_scalar(residuals, -1000000);
789// cpl_vector_multiply_scalar(all_wavelength, 1000000);
790// vectors[2]=residuals;
792// if (POSTSCRIPT_PLOT) ps_string=cpl_sprintf("set term png; set output 'plot_wavelength_%s.png';", "SC_Argon");
793// else ps_string=cpl_sprintf(" ");
794// cpl_plot_vectors(cpl_sprintf("%s set xlabel \"Index\";set ylabel \"Wavelength\";", ps_string)," title 'mesured';", NULL, vectors, 3);
795// cpl_plot_bivector(cpl_sprintf("%s set xlabel \"Index\";set ylabel \"Wavelength\";", ps_string)," title 'error';", NULL, error);
796// free(vectors);
797// cpl_bivector_delete(error);
798// cpl_free(ps_string);
799// }
800
801 cpl_vector_delete(coord_X_fit);
802 cpl_vector_delete(all_wavelength);
803 cpl_vector_delete(coord_X);
804 cpl_vector_delete(coord_Y);
805 cpl_matrix_delete (all_coord);
806 cpl_matrix_unwrap(matrix2);
807 cpl_matrix_unwrap(matrix3);
808 cpl_vector_delete(residuals);
809 cpl_vector_delete(fitsigm);
810
811 /* Calcul of the new wavelength vector after calibration */
812
813 /* Save the result wavelength on the associated bases depending
814 * of the polarization */
815 cpl_array * dimension = cpl_array_new(2, CPL_TYPE_INT);
816 cpl_array_set(dimension, 0, 1);
817 cpl_array_set(dimension, 1, nwave);
818
819 profile_table = gravi_data_get_table (profile_map, GRAVI_PROFILE_DATA_EXT);
820 sizex = cpl_table_get_column_dimension (profile_table, "DATA1", 0);
821 sizey = cpl_table_get_column_dimension (profile_table, "DATA1", 1);
822 profile_image = cpl_image_new (sizex, sizey, CPL_TYPE_DOUBLE);
823 cpl_image_fill_window (profile_image, 1, 1, sizex, sizey, 0.0);
824 image_wave = cpl_image_new (sizex, sizey, CPL_TYPE_DOUBLE);
825 cpl_image_fill_window (image_wave, 1, 1, sizex, sizey, 0.0);
826
827 img_output = cpl_table_new (1);
828
829 for (region = 0 ; region < n_region; region ++) {
830
831 data_x = cpl_sprintf("DATA%d", region + 1);
832
833 imglist_wave = gravi_table_data_to_imagelist(profile_table, region + 1);
834 img_profile = cpl_imagelist_get(imglist_wave, 0);
835
836 cpl_table_new_column_array (img_output, data_x, CPL_TYPE_DOUBLE, nwave);
837 cpl_table_set_column_dimensions(img_output, data_x, dimension);
838
839 wavelength = cpl_array_new(nwave, CPL_TYPE_DOUBLE);
840
841 /* Calculate of the new wavelength vector after calibration */
842 y_corner = 1;
843
844 cpl_image_add (profile_image, img_profile);
845 pos = cpl_vector_new(2);
846 for (wave = 0; wave < nwave; wave ++) {
847
848 cpl_vector_set(pos, 0, wave);
849 cpl_vector_set(pos, 1, region + 1);
850 result=gsl_vector_get (c, 0) +
851 gsl_vector_get (c, 1) * wave +
852 gsl_vector_get (c, 2) * wave*wave +
853 gsl_vector_get (c, 3) * (region+1) +
854 gsl_vector_get (c, 4) * (region+1)*(region+1) +
855 gsl_vector_get (c, 5) * (region+1)*wave;
856
857 //result = cpl_polynomial_eval(fit2d, pos);
858
859 cpl_array_set(wavelength, wave, result);
860 y_corner ++;
861
862 for (ind = 0; ind < sizey; ind ++){
863
864 if (cpl_image_get (img_profile, wave+1, ind+1, &nv) > 0.01)
865 cpl_image_set (image_wave, wave+1, ind+1, result);
866
867 CPLCHECK_NUL("The corner and image_wave");
868 }
869 }
870 cpl_vector_delete(pos);
871
872 /* Get the miniumum and maximum wavelength */
873 minwave = CPL_MAX (minwave, cpl_array_get_min(wavelength));
874 maxwave = CPL_MIN (maxwave, cpl_array_get_max(wavelength));
875
876 cpl_table_set_array(img_output, data_x, 0, wavelength);
877
878 cpl_imagelist_delete (imglist_wave);
879 cpl_array_delete(wavelength);
880 cpl_free(data_x);
881
882 } /* End loop on regions */
883
884
885 gravi_data * wave_map;
886
887 wave_map = gravi_data_new (0);
888 plist = gravi_data_get_header (wave_data);
889 gravi_data_append_header (wave_map, plist);
890 int nb_ext = gravi_data_get_size (wave_data), ext, type_ext;
891 cpl_propertylist * img_plist;
892
893 for (ext = 0; ext < nb_ext; ext++ ){
894
895 /*
896 * Load the FT or SC
897 */
898 img_plist = gravi_data_get_plist_x (wave_data, ext);
899
900 const char * plist_name = gravi_pfits_get_extname (img_plist);
901 /* Check if the needed extentions are there */
902 type_ext = gravi_pfits_get_extension_type (img_plist);
903 if (!(strcmp (plist_name, "WAVE_FIBRE_SC") &&
904 strcmp (plist_name, GRAVI_IMAGING_DETECTOR_SC_EXT) &&
905 strcmp (plist_name, GRAVI_WAVE_DATA_SC_EXT) &&
906 strcmp (plist_name, "TEST_WAVE"))){
907 /* Load the FT table */
908 if (type_ext == 2){
909 gravi_data_add (wave_map, img_plist,
910 cpl_table_duplicate (gravi_data_get_table (wave_data, plist_name)));
911 if(! strcmp (plist_name, GRAVI_WAVE_DATA_SC_EXT))
912 cpl_propertylist_set_string (gravi_data_get_plist (wave_map,
913 GRAVI_WAVE_DATA_SC_EXT), "EXTNAME", GRAVI_WAVE_ARGON_EXT);
914 }
915
916 /* Load the SC image_list */
917 else if (type_ext == 3)
918 gravi_data_add_cube (wave_map, img_plist,
919 cpl_imagelist_duplicate (gravi_data_get_cube (wave_data, plist_name)));
920 }
921 }
922
923 /* Set the QC parameters */
924 primary_hdr = gravi_data_get_header (wave_map);
925 cpl_propertylist_append_double (primary_hdr, QC_MINWAVE_SC, minwave);
926 cpl_propertylist_append_double (primary_hdr, QC_MAXWAVE_SC, maxwave);
927 cpl_msg_info (cpl_func, "QC_MINWAVE_SC = %e QC_MAXWAVE_SC = %e", minwave, maxwave);
928
929 /* Add the img_output and wave_fibre */
930 gravi_data_set_table (wave_map, GRAVI_WAVE_ARGON_EXT, img_output);
931
932 /* Add the image_wave and profile_image */
933 imglist_wave = cpl_imagelist_new ();
934 cpl_imagelist_set(imglist_wave, image_wave, 0);
935 cpl_imagelist_set(imglist_wave, profile_image, 1);
936
937 gravi_data_set_cube (wave_map, "TEST_WAVE", imglist_wave);
938
939 cpl_polynomial_delete (fit2d);
940 cpl_array_delete (dimension);
941 gravi_data_delete (output_data);
942
943 gravi_msg_function_exit(1);
944 return wave_map;
945}
946
947gravi_data * gravi_compute_wave_offset (gravi_data * argon_wave, gravi_data * wave_data) {
948
949 cpl_propertylist * primary_hdr, * plist, * wavePlist, * detectorPlist;
950 cpl_table * detector, * waveArgon_table, * waveData_table,
951 * img_output;
952 cpl_matrix * all_coord;
953 cpl_array * wavelength;
954 cpl_polynomial * fit_slope;
955 cpl_polynomial * fit2d;
956 cpl_vector * residuals;
957 double rechisq;
958 // Degrees of wavelength fit polynomial {2 in vertical direction, 4 in disp direction}
959 const cpl_size deg2d[2] = {4, 4};
960 cpl_vector * coord_X, * coord_Y, * all_wavelength;
961 double result;
962 int wave;
963 double slope;
964 int y_corner;
965 int nwave, n_region, region;
966 char * data_x;
967 int nv;
968// const cpl_vector ** vectors;
969 cpl_bivector * plot;
970 cpl_vector *pos;
971 double minwave, maxwave;
972 cpl_image * image_wave;
973 cpl_imagelist * imglist_wave;
974 int sizey, sizex, ind, npol;
975 /* Tel for base */
976 //int tel1[6] = {0,0,0,1,1,2};
977 //int tel2[6] = {1,2,3,2,3,3};
978 char * baseString [6] = {"12", "13", "14", "23", "24", "34"};
979 char * baseString_bis [6] = {"21", "31", "41", "32", "42", "43"};
980 const cpl_size deg = 0, degM = 1;
981 cpl_vector * position;
982 cpl_matrix * matFit2;
983
984
985 if ((argon_wave == NULL) || (wave_data == NULL)){
986 cpl_error_set_message(cpl_func, CPL_ERROR_NULL_INPUT,
987 "The data is NULL");
988 return NULL;
989 }
990
991 int nbase = 6, base, type_data;
992 primary_hdr = gravi_data_get_header (argon_wave);
993 const char * resolArg = gravi_pfits_get_resolution (primary_hdr);
994 const char * resolWave = gravi_pfits_get_resolution (gravi_data_get_header (wave_data));
995
996 if (strcmp(resolArg, resolWave)){
997 cpl_error_set_message(cpl_func, CPL_ERROR_ILLEGAL_OUTPUT,
998 "The argon file doesn t have the same"
999 "resultion as the input wave");
1000 return NULL;
1001 }
1002
1003 /* Compare the wavelength for each base */
1004
1005 type_data = 0;
1006
1007 detector = gravi_data_get_table (argon_wave, GRAVI_IMAGING_DETECTOR_SC_EXT);
1008 detectorPlist = gravi_data_get_plist (wave_data, GRAVI_IMAGING_DETECTOR_SC_EXT);
1009 wavePlist = gravi_data_get_plist (wave_data, GRAVI_WAVE_DATA_SC_EXT);
1010
1011 n_region = cpl_table_get_nrow (detector) ;
1012
1013 npol = (n_region > 24)?2:1;
1014
1015 char * pola = (npol == 2)?POLAR_1:POLAR_3;
1016
1017 waveArgon_table = gravi_data_get_table (argon_wave, GRAVI_WAVE_ARGON_EXT);
1018
1019 waveData_table = gravi_data_get_table (wave_data, GRAVI_WAVE_FIBRE_SC_EXT);
1020
1021 nwave = cpl_table_get_column_depth(waveArgon_table, "DATA1");
1022 char * base_fiber = cpl_sprintf("BASE_12_%s", pola);
1023 if (nwave!= cpl_table_get_column_depth(waveData_table, base_fiber)){
1024 cpl_error_set_message(cpl_func, CPL_ERROR_ILLEGAL_OUTPUT,
1025 "The argon and the wave data do not have"
1026 " the same number of wavelength");
1027 return NULL;
1028 }
1029 cpl_array * wave_;
1030 int comp;
1031
1032 cpl_imagelist * imglist = gravi_data_get_cube (wave_data, "TEST_WAVE");
1033 cpl_image * img_profile = cpl_imagelist_get (imglist, 1);
1034 sizex = cpl_image_get_size_x (img_profile);
1035 sizey = cpl_image_get_size_y (img_profile);
1036 image_wave = cpl_image_new (sizex, sizey, CPL_TYPE_DOUBLE);
1037 cpl_image * image_real = cpl_image_new (sizex, sizey, CPL_TYPE_DOUBLE);
1038 cpl_image_fill_window (image_real, 1, 1, sizex, sizey, 0.0);
1039 matFit2 = cpl_matrix_new (2,n_region*nwave);
1040 all_wavelength = cpl_vector_new (n_region*nwave);
1041 int regionA, regionD, regionC, regionB;
1042 cpl_array * waveToCal;
1043
1044 cpl_table * cooef_table = cpl_table_new (n_region);
1045 cpl_table_new_column(cooef_table, "SLOPE", CPL_TYPE_DOUBLE);
1046 cpl_table_new_column(cooef_table, "OFFSET", CPL_TYPE_DOUBLE);
1047
1048 int startx = gravi_pfits_get_startx(wavePlist) +
1049 gravi_pfits_get_window_start(primary_hdr) -2, i;
1050 if (cpl_error_get_code ()){
1051 cpl_error_set_message(cpl_func, CPL_ERROR_NULL_INPUT, "Problem to get "
1052 "the start window detector");
1053 return NULL;
1054 }
1055 int np = 0, nq = 1;
1056 cpl_vector * vect_index = cpl_vector_new (nwave), * vect;
1057
1058 for (i = startx; i < startx + nwave; i++){
1059 if ((i/64)%2 == 0){
1060 np ++;
1061 cpl_vector_set (vect_index, i - startx, 0);
1062 }
1063 else {
1064 nq ++;
1065 cpl_vector_set (vect_index, i - startx, 1);
1066 }
1067 }
1068
1069 cpl_vector * residual_p = cpl_vector_new(np * n_region),
1070 * residual_q = cpl_vector_new(nq * n_region);
1071
1072 for (region = 0; region < n_region; region ++) {
1073 /* Look for the region with the same base none looking for
1074 * the polarization */
1075// comp = 0;
1076//
1077// wavelength = cpl_array_new (nwave, CPL_TYPE_DOUBLE);
1078// waveToCal = cpl_array_new (nwave, CPL_TYPE_DOUBLE);
1079// cpl_array_fill_window (wavelength, 0, nwave, 0.0);
1080// cpl_array_fill_window (waveToCal, 0, nwave, 0.0);
1081//
1082// regionA = gravi_table_get_regindex(detector, tel1[base] + 1,
1083// tel2[base] + 1, "A", pola);
1084//
1085// if (regionA != 0) {
1086// data_x = cpl_sprintf ("DATA%d", regionA);
1087//
1088//
1089// wave_ = cpl_array_duplicate (cpl_table_get_array (waveData_table,
1090// data_x, 0));
1091// cpl_array_add (waveToCal, wave_);
1092// cpl_array_subtract (wave_, cpl_table_get_array (waveArgon_table,
1093// data_x, 0));
1094// cpl_free (data_x);
1095//
1096// cpl_array_add (wavelength, wave_);
1097// cpl_array_delete (wave_);
1098// comp ++;
1099// }
1100//
1101// regionB = gravi_table_get_regindex(detector, tel1[base] + 1,
1102// tel2[base] + 1, "B", pola);
1103//
1104// if (regionB != 0) {
1105// data_x = cpl_sprintf ("DATA%d", regionB);
1106//
1107//
1108// wave_ = cpl_array_duplicate (cpl_table_get_array (waveData_table,
1109// data_x, 0));
1110// cpl_array_add (waveToCal, wave_);
1111// cpl_array_subtract (wave_, cpl_table_get_array (waveArgon_table,
1112// data_x, 0));
1113// cpl_free (data_x);
1114//
1115// cpl_array_add(wavelength, wave_);
1116// cpl_array_delete (wave_);
1117// comp ++;
1118// }
1119//
1120// regionC = gravi_table_get_regindex(detector, tel1[base] + 1,
1121// tel2[base] + 1, "C", pola);
1122// if (regionC != 0) {
1123// data_x = cpl_sprintf ("DATA%d", regionC);
1124//
1125//
1126// wave_ = cpl_array_duplicate (cpl_table_get_array (waveData_table,
1127// data_x, 0));
1128// cpl_array_add (waveToCal, wave_);
1129// cpl_array_subtract (wave_, cpl_table_get_array (waveArgon_table,
1130// data_x, 0));
1131// cpl_free (data_x);
1132//
1133// cpl_array_add(wavelength, wave_);
1134// cpl_array_delete (wave_);
1135// comp ++;
1136// }
1137//
1138// regionD = gravi_table_get_regindex(detector, tel1[base] + 1,
1139// tel2[base] + 1, "D", pola);
1140// if (regionD != 0) {
1141// data_x = cpl_sprintf ("DATA%d", regionD);
1142//
1143//
1144// wave_ = cpl_array_duplicate (cpl_table_get_array (waveData_table,
1145// data_x, 0));
1146// cpl_array_add (waveToCal, wave_);
1147// cpl_array_subtract (wave_, cpl_table_get_array (waveArgon_table,
1148// data_x, 0));
1149// cpl_free (data_x);
1150//
1151// cpl_array_add(wavelength, wave_);
1152// cpl_array_delete (wave_);
1153// comp ++;
1154// }
1155//
1156// if (comp != 0){
1157// cpl_array_divide_scalar (wavelength, comp);
1158// cpl_array_divide_scalar (waveToCal, comp);
1159// }
1160//
1161// else {
1162// cpl_msg_info (cpl_func, "The base %d %d does "
1163// "not exist", tel1[base]+1, tel2[base]+1);
1164// cpl_array_delete(wavelength);
1165// cpl_array_delete(waveToCal);
1166// continue;
1167// }
1168
1169 const char * regname = cpl_table_get_string (detector, "REGNAME", region);
1170 for (base = 0; base < nbase; base++){
1171// printf("baseString[%d] = %s baseString_bis[%d] = %s\n",base,
1172// baseString[base],base, baseString_bis[base]);
1173 if (strstr(regname, baseString[base]) ||
1174 strstr(regname, baseString_bis[base]))
1175 break;
1176 }
1177
1178 if (npol == 2) {
1179 if (strstr (regname, POLAR_1))
1180 pola = POLAR_1;
1181 else
1182 pola = POLAR_2;
1183 }
1184
1185 base_fiber = cpl_sprintf ("BASE_%s_%s", baseString[base], pola);
1186 data_x = cpl_sprintf ("DATA%d", region + 1);
1187 wavelength = cpl_array_duplicate (cpl_table_get_array (waveData_table,
1188 base_fiber, 0));
1189 waveToCal = cpl_array_duplicate (wavelength);
1190 cpl_array_subtract (wavelength, cpl_table_get_array (waveArgon_table,
1191 data_x, 0));
1192
1193 /* Fit to get the slope of the polynomial wavelength */
1194
1195 cpl_matrix * matrix = cpl_matrix_new (1, nwave);
1196 position = cpl_vector_wrap (nwave, cpl_array_get_data_double(wavelength));
1197
1198 for (ind = 0; ind < nwave; ind++) {
1199 cpl_matrix_set (matFit2,0, region*nwave + ind,
1200 region);
1201 cpl_matrix_set (matFit2,1, region*nwave + ind,
1202 ind);
1203 cpl_matrix_set(matrix, 0, ind, ind);
1204
1205 }
1206
1207// // fit with cpl_fit_lvmq
1208// int
1209// cpl_fit_lvmq(matrix, NULL, position,
1210// NULL, init_val, val_to_fit, &polystep,
1211// NULL, CPL_FIT_LVMQ_TOLERANCE, CPL_FIT_LVMQ_COUNT,
1212// CPL_FIT_LVMQ_MAXITER, &mse, &red_chisq, NULL);
1213
1214
1215 fit_slope = cpl_polynomial_new(1);
1216 cpl_polynomial_fit(fit_slope, matrix, NULL, position, NULL,
1217 CPL_FALSE, &deg, &degM);
1218
1219 vect = cpl_vector_new (nwave);
1220 cpl_vector_fill_polynomial_fit_residual (vect, position, NULL,
1221 fit_slope, matrix, &rechisq );
1222 if (cpl_error_get_code ()){
1223 cpl_error_set_message(cpl_func, CPL_ERROR_NULL_INPUT, "fit the residual");
1224 return NULL;
1225 }
1226
1227 int ip = 0, iq = 0;
1228
1229 for (wave = 0; wave < nwave; wave ++)
1230 if (cpl_vector_get (vect_index, wave) == 0) {
1231 cpl_vector_set (residual_p, ip + np * region ,
1232 cpl_vector_get (vect, wave));
1233 ip ++;
1234 }
1235 else {
1236 cpl_vector_set (residual_q, iq + nq * region,
1237 cpl_vector_get (vect, wave));
1238 iq ++;
1239 }
1240
1241
1242
1243 const cpl_size pow_slope = 1;
1244 const cpl_size pow_slope2 = 0;
1245 cpl_msg_info (cpl_func, "region %s : y = %e * x + %e", regname,
1246 cpl_polynomial_get_coeff(fit_slope,&pow_slope ),
1247 cpl_polynomial_get_coeff(fit_slope,
1248 &pow_slope2 ));
1249
1250 cpl_table_set_double (cooef_table, "SLOPE", region,
1251 cpl_polynomial_get_coeff(fit_slope,&pow_slope ));
1252 cpl_table_set_double (cooef_table, "OFFSET", region,
1253 cpl_polynomial_get_coeff(fit_slope,&pow_slope2 ));
1254 if (cpl_error_get_code()){
1255
1256 cpl_error_set_message(cpl_func, CPL_ERROR_NULL_INPUT, "fit slope");
1257 return NULL;
1258 }
1259
1260 pos = cpl_vector_new (1);
1261 for (wave = 0; wave < nwave; wave++){
1262 cpl_vector_set (pos, 0, wave);
1263 cpl_vector_set(all_wavelength, region*nwave + wave,
1264 cpl_array_get_double (waveToCal, wave, &nv) -
1265 cpl_polynomial_eval(fit_slope, pos));
1266 }
1267
1268 cpl_vector_delete (pos);
1269 cpl_polynomial_delete (fit_slope);
1270 cpl_matrix_delete (matrix);
1271 cpl_vector_unwrap (position);
1272 cpl_array_delete(wavelength);
1273 cpl_array_delete(waveToCal);
1274 cpl_free (data_x);
1275 cpl_free (base_fiber);
1276 }
1277
1278 /* Computation of the step */
1279 double step = cpl_vector_get_median (residual_p) - cpl_vector_get_median (residual_q);
1280
1281
1282 cpl_msg_info(cpl_func, "step = %g", step);
1283
1284 /* Fit */
1285 fit2d = cpl_polynomial_new(2);
1286 cpl_polynomial_fit(fit2d, matFit2, NULL, all_wavelength, NULL,
1287 CPL_TRUE, NULL, deg2d);
1288 cpl_matrix_delete (matFit2);
1289
1290 if (cpl_error_get_code()){
1291
1292 cpl_error_set_message(cpl_func, CPL_ERROR_NULL_INPUT, "fit 2D");
1293 return NULL;
1294 }
1295
1296 /* Calcul of the new wavelength vector after calibration */
1297
1298 /* Save the result wavelength on the associated bases depending
1299 * of the polarization */
1300 cpl_array * dimension = cpl_array_new(2, CPL_TYPE_INT);
1301 cpl_array_set(dimension, 0, 1);
1302 cpl_array_set(dimension, 1, nwave);
1303
1304 cpl_image_fill_window (image_wave, 1, 1, sizex, sizey, 0.0);
1305
1306 img_output = NULL;
1307 img_output = cpl_table_new (1);
1308
1309 for (region = 0 ; region < n_region; region ++){
1310
1311 data_x = cpl_sprintf("DATA%d", region + 1);
1312
1313 cpl_table_new_column_array (img_output, data_x, CPL_TYPE_DOUBLE, nwave);
1314 cpl_table_set_column_dimensions(img_output, data_x,
1315 dimension);
1316 const char * regname = cpl_table_get_string (detector, "REGNAME", region);
1317
1318// printf("regname = %s\n", regname);
1319 for (base = 0; base < nbase; base++){
1320// printf("baseString[%d] = %s baseString_bis[%d] = %s\n",base,
1321// baseString[base],base, baseString_bis[base]);
1322 if (strstr(regname, baseString[base]) ||
1323 strstr(regname, baseString_bis[base]))
1324 break;
1325 }
1326// printf("base = %d\n", base);
1327 if (base == nbase){
1328 cpl_error_set_message(cpl_func, CPL_ERROR_NULL_INPUT,
1329 "Check the base");
1330 return NULL;
1331 }
1332// printf("base = %d\n", base);
1333 wavelength = cpl_array_new(nwave, CPL_TYPE_DOUBLE);
1334 cpl_vector *wavelength_vect = cpl_vector_new(nwave);
1335// vectors[region + 1] = cpl_vector_new (nwave);
1336 /* Calculate of the new wavelength vector after calibration */
1337 y_corner = 1;
1338
1339 pos = cpl_vector_new(2);
1340 for (wave = 0; wave < nwave; wave ++){
1341
1342
1343 cpl_vector_set(pos, 0, region + 1);
1344 cpl_vector_set(pos, 1, wave);
1345
1346 result = cpl_polynomial_eval(fit2d, pos);
1347// result = cpl_vector_get (all_wavelength, base*nwave + wave) +
1348 result = cpl_vector_get (all_wavelength, region*nwave + wave) -
1349 cpl_vector_get (vect_index, wave) * step;
1350 cpl_array_set(wavelength, wave, result);
1351 cpl_vector_set (wavelength_vect, wave, result );
1352 y_corner ++;
1353
1354 for (ind = 0; ind < sizey; ind ++){
1355
1356 if (cpl_image_get (img_profile, wave+1, ind+1, &nv) > 0.01){
1357
1358 cpl_image_set (image_wave, wave+1, ind+1, result);
1359
1360 cpl_image_set (image_real, wave+1, ind+1,
1361// cpl_vector_get (all_wavelength, base*nwave + wave));
1362 cpl_vector_get (all_wavelength, region*nwave + wave));
1363 }
1364
1365 if (cpl_error_get_code()){
1366
1367 cpl_error_set_message(cpl_func, CPL_ERROR_NULL_INPUT, "The corner and "
1368 "image_wave");
1369 return NULL;
1370 }
1371
1372 }
1373
1374
1375 }
1376 cpl_vector_delete(pos);
1377
1378
1379 /*
1380 * Fit the wave_fibre, and add the OFFSET, SLOPE and STEP
1381 */
1382 cpl_matrix * matrix = cpl_matrix_new (1, nwave);
1383 const cpl_size deg_2=2;
1384
1385 for (ind = 0; ind < nwave; ind++) {
1386 cpl_matrix_set(matrix, 0, ind, ind);
1387 }
1388
1389 fit_slope = cpl_polynomial_new(1);
1390 cpl_polynomial_fit(fit_slope, matrix, NULL, wavelength_vect, NULL,
1391 CPL_FALSE, &deg, &deg_2);
1392
1393 pos=cpl_vector_new(1);
1394 for (wave = 0; wave < nwave; wave ++){
1395 cpl_vector_set(pos, 0, wave);
1396 result = cpl_polynomial_eval(fit_slope, pos)-
1397 cpl_vector_get (vect_index, wave) * step;
1398 cpl_array_set(wavelength, wave, result);
1399 }
1400 cpl_vector_delete(pos);
1401 cpl_polynomial_delete(fit_slope);
1402
1403 /* Get the miniumum and maximum wavelength */
1404 if (region == 0){
1405 minwave = cpl_array_get_min(wavelength);
1406 maxwave = cpl_array_get_max(wavelength);
1407 }
1408 else {
1409
1410 if (minwave < cpl_array_get_min(wavelength))
1411 minwave = cpl_array_get_min(wavelength);
1412
1413
1414 if (maxwave > cpl_array_get_max(wavelength))
1415 maxwave = cpl_array_get_max(wavelength);
1416 }
1417
1418
1419 cpl_table_set_array(img_output, data_x, 0, wavelength);
1420
1421 cpl_array_delete(wavelength);
1422 cpl_free(data_x);
1423 }
1424
1425 cpl_vector_delete (all_wavelength);
1426 gravi_data * wave_map;
1427 wave_map = gravi_data_new (0);
1428 plist = gravi_data_get_header (wave_data);
1429 gravi_data_append_header (wave_map, plist);
1430 int nb_ext = gravi_data_get_size (wave_data), ext, type_ext;
1431 cpl_propertylist * img_plist;
1432
1433 for (ext = 0; ext < nb_ext; ext++ ){
1434
1435 /*
1436 * Load the FT or SC
1437 */
1438 img_plist = gravi_data_get_plist_x (wave_data, ext);
1439
1440 const char * plist_name = gravi_pfits_get_extname (img_plist);
1441 /* Check if the needed extentions are there */
1442 type_ext = gravi_pfits_get_extension_type (img_plist);
1443 if (!(strcmp (plist_name, GRAVI_WAVE_FIBRE_SC_EXT) &&
1444 strcmp (plist_name, GRAVI_IMAGING_DETECTOR_SC_EXT) &&
1445 strcmp (plist_name, GRAVI_WAVE_DATA_SC_EXT) &&
1446 strcmp (plist_name, "TEST_WAVE") &&
1447 strcmp (plist_name, GRAVI_WAVE_DATA_FT_EXT) &&
1448 strcmp (plist_name, GRAVI_IMAGING_DETECTOR_FT_EXT))){
1449 /* Load the FT table */
1450 if (type_ext == 2)
1451 gravi_data_add (wave_map, img_plist,
1452 cpl_table_duplicate (gravi_data_get_table (wave_data, plist_name)));
1453
1454 /* Load the SC image_list */
1455 else if (type_ext == 3)
1456 gravi_data_add_cube (wave_map, img_plist,
1457 cpl_imagelist_duplicate (gravi_data_get_cube (wave_data, plist_name)));
1458 }
1459 }
1460
1461// wave_map = gravi_data_duplicate (wave_data);
1462
1463 /* Add the offset table */
1464 plist = cpl_propertylist_duplicate (gravi_data_get_plist (wave_map,
1465 GRAVI_WAVE_DATA_SC_EXT));
1466 cpl_propertylist_append_string (plist, "EXTNAME",
1467 "WAVE_OFFSET");
1468 gravi_data_add (wave_map, plist, cooef_table);
1469 cpl_propertylist_delete (plist);
1470
1471 primary_hdr = gravi_data_get_header (wave_map);
1472
1473 cpl_polynomial_delete(fit2d);
1474 cpl_propertylist_append_double (primary_hdr, QC_MINWAVE_SC, minwave);
1475 cpl_propertylist_append_double (primary_hdr, QC_MAXWAVE_SC, maxwave);
1476
1477 cpl_msg_info (cpl_func, "wave corrected : QC_MINWAVE_SC = %e QC_MAXWAVE_SC = %e",
1478 minwave, maxwave);
1479
1480 gravi_data_set_table (wave_map, GRAVI_WAVE_DATA_SC_EXT, img_output);
1481
1482 imglist_wave = cpl_imagelist_new ();
1483 cpl_imagelist_set(imglist_wave,image_wave , 0);
1484 cpl_imagelist_set(imglist_wave, cpl_image_duplicate (img_profile), 1);
1485 cpl_imagelist_set(imglist_wave, image_real, 2);
1486 gravi_data_set_cube (wave_map, "TEST_WAVE", imglist_wave);
1487
1488 cpl_array_delete (dimension);
1489
1490 return wave_map;
1491}
1492
1493cpl_error_code gravi_data_save( gravi_data * self,
1494 cpl_frameset * allframes,
1495 const char * filename,
1496 const cpl_parameterlist * parlist,
1497 cpl_frameset * usedframes,
1498 cpl_frame * frame,
1499 const char * recipe,
1500 cpl_propertylist * applist){
1501 cpl_frameset * frameset;
1502 int ext;
1503 /* Check the inputs */
1504 if ((filename == NULL) || (self == NULL)){
1505 cpl_error_set_message(cpl_func, CPL_ERROR_NULL_INPUT,
1506 "one of the inputs at least is NULL");
1507 return CPL_ERROR_NULL_INPUT;
1508 }
1509
1510 /* Create the file and save the primary header. */
1511
1512 frameset = cpl_frameset_duplicate(usedframes);
1513// cpl_frame * frame = cpl_frameset_get_first(frameset);
1514
1515 if (cpl_dfs_save_propertylist (allframes, self->primary_hdr, parlist,
1516 frameset, frame, recipe, applist,
1517 NULL, PACKAGE_STRING, filename ) != CPL_ERROR_NONE){
1518 cpl_error_set_message(cpl_func, CPL_ERROR_ILLEGAL_OUTPUT,
1519 "Cannot save the first extension primary header");
1520 return CPL_ERROR_NULL_INPUT;
1521 }
1522
1523 /* Save the extensions */
1524 for (ext = 0; ext < self->nb_ext; ext ++){
1525 if (self->exts_tbs[ext] != NULL)
1526 cpl_table_save(self->exts_tbs[ext], NULL, self->exts_hdrs[ext],
1527 filename, CPL_IO_EXTEND);
1528 else if (self->exts_imgl[ext] != NULL)
1529 cpl_imagelist_save (self->exts_imgl[ext], filename,
1530 cpl_image_get_type (cpl_imagelist_get (self->exts_imgl[ext], 0)),
1531 self->exts_hdrs[ext], CPL_IO_EXTEND);
1532 }
1533 cpl_frameset_delete (frameset);
1534
1535 return CPL_ERROR_NONE;
1536}
1537
1538/*
1539 * Normalize the OI_FLUX (FLUX and FLUXERR columns) by the mean of FLUX(lbd)
1540 * Independly for each row.
1541 */
1542
1543cpl_error_code gravi_normalize_flux (gravi_data * data)
1544{
1545 gravi_msg_function_start(0);
1546
1547 /* Get the header */
1548 cpl_propertylist * hdr_data = gravi_data_get_header (data);
1549
1550 /* Loop on FT and SC */
1551 for (int type_data = 0; type_data < 2; type_data ++) {
1552
1553 /* Loop on polarisation */
1554 int npol = gravi_pfits_get_pola_num (hdr_data, type_data);
1555 for ( int pol= 0 ; pol < npol ; pol++ ) {
1556
1557 cpl_msg_info (cpl_func, "Normalize the flux of %s (pol %i over %i)",GRAVI_TYPE(type_data),pol+1,npol);
1558
1559 cpl_table * oi_flux = gravi_data_get_oi_flux (data, type_data, pol, npol);
1560 cpl_array ** flux = cpl_table_get_data_array (oi_flux, "FLUX");
1561 cpl_array ** fluxErr = cpl_table_get_data_array (oi_flux, "FLUXERR");
1562 cpl_size nrow = cpl_table_get_nrow (oi_flux);
1563 CPLCHECK_MSG("Cannot get data");
1564
1565 /* Loop on rows to normalize */
1566 for (cpl_size row = 0; row < nrow; row++) {
1567 double flux_mean = cpl_array_get_mean (flux[row]);
1568 cpl_array_divide_scalar (flux[row], flux_mean);
1569 cpl_array_divide_scalar (fluxErr[row], flux_mean);
1570 }
1571
1572 CPLCHECK_MSG("Cannot divide by the mean flux");
1573 } /* End loop on pol */
1574 } /* End loop on SC/FT */
1575
1576 gravi_msg_function_exit(0);
1577 return CPL_ERROR_NONE;
1578}
1579
1580/*---------------------------------------------------------------------------*/
1581
1582cpl_vector * gravi_image_to_vector(cpl_image * img)
1583{
1584 cpl_ensure (img, CPL_ERROR_NULL_INPUT, NULL);
1585
1586 cpl_vector * vector;
1587 cpl_type type_img;
1588 cpl_image * _image;
1589 int x, y;
1590 double * data;
1591
1592 type_img = cpl_image_get_type (img);
1593
1594 _image = cpl_image_cast (img, CPL_TYPE_DOUBLE);
1595
1596 x = cpl_image_get_size_x (img);
1597 y = cpl_image_get_size_y (img);
1598 data = cpl_image_get_data_double(_image);
1599 vector = cpl_vector_new(x*y);
1600 for (int i = 0 ; i < x*y; i++)
1601 cpl_vector_set(vector, i, data[i]);
1602
1603 cpl_image_delete(_image);
1604 return vector;
1605}
1606
1607/* Transform tangent plane coordinates into spherical
1608 * (double precision). All in [rad]. Return RA in
1609 *
1610 * Given:
1611 * XI,ETA dp tangent plane rectangular coordinates
1612 * RAZ,DECZ dp spherical coordinates of tangent point
1613 *
1614 * Returned:
1615 * RA,DEC dp spherical coordinates (0-2pi,+/-pi/2)
1616 *
1617 * From: sla_DRANRM
1618 * P.T.Wallace Starlink 24 July 1995
1619 * Copyright (C) 1995 Rutherford Appleton Laboratory
1620 */
1621cpl_error_code gravi_dtps(double xi, double eta, double raz, double decz, double * ra, double * dec)
1622{
1623 double sdecz = sin(decz);
1624 double cdecz = cos(decz);
1625 double denom = cdecz - eta * sdecz;
1626 /* Compute new coordinates */
1627 *dec = atan2 (sdecz+eta*cdecz, sqrt(xi*xi + denom*denom));
1628 *ra = atan2 (xi,denom) + raz;
1629
1630 /* Make ra within 0-2pi */
1631 *ra = fmod (*ra, CPL_MATH_2PI);
1632 if ( *ra < 0.0 ) *ra += CPL_MATH_2PI;
1633
1634 return CPL_ERROR_NONE;
1635}
1636
1637
1638int my_gsl_matrix_complex_fprintf(FILE *stream,gsl_matrix_complex *m,char *fmt)
1639{
1640 size_t rows=m->size1;
1641 size_t cols=m->size2;
1642 size_t row,col,ml;
1643 int fill;
1644 char buf[100];
1645 gsl_vector *maxlen;
1646 gsl_complex buff;
1647
1648 maxlen=gsl_vector_alloc(cols);
1649 for (col=0;col<cols;++col) {
1650 ml=0;
1651 for (row=0;row<rows;++row) {
1652 buff=gsl_matrix_complex_get(m,row,col);
1653 sprintf(buf,"%g + i%g", GSL_REAL (buff), GSL_IMAG(buff));
1654 if (strlen(buf)>ml)
1655 ml=strlen(buf);
1656 }
1657 gsl_vector_set(maxlen,col,ml);
1658 }
1659
1660 for (row=0;row<rows;++row) {
1661 for (col=0;col<cols;++col) {
1662 buff=gsl_matrix_complex_get(m,row,col);
1663 sprintf(buf,"%g + i%g", GSL_REAL (buff), GSL_IMAG(buff));
1664 fprintf(stream,"%s",buf);
1665 fill=gsl_vector_get(maxlen,col)+1-strlen(buf);
1666 while (--fill>=0)
1667 fprintf(stream," ");
1668 }
1669 fprintf(stream,"\n");
1670 }
1671 gsl_vector_free(maxlen);
1672 return 0;
1673}
1674
1675
1676/*----------------------------------------------------------------------------*/
1688/*----------------------------------------------------------------------------*/
1689
1690cpl_error_code gravi_compute_flat_badpix(gravi_data * flat, gravi_data * dark)
1691{
1692 double threshold;
1693 cpl_ensure_code (flat, CPL_ERROR_NULL_INPUT);
1694 cpl_ensure_code (dark, CPL_ERROR_NULL_INPUT);
1695
1696 cpl_msg_info(cpl_func, "Compute SC flat bad pixels");
1697
1698 /* get the median flat */
1699 cpl_image * median_flat = cpl_imagelist_get(
1700 gravi_data_get_cube(flat, GRAVI_IMAGING_DATA_SC_EXT), 0);
1701
1702 /* Compute bad pixels from flat*/
1703 cpl_image * grad = cpl_image_duplicate(median_flat);
1704 int x, size_x = cpl_image_get_size_x (median_flat);
1705 int y, size_y = cpl_image_get_size_y (median_flat);
1706
1707 /* Mean profile */
1708 cpl_vector * vect_mean = cpl_vector_new_from_image_row (median_flat, 1);
1709 cpl_vector * vector;
1710 for (y = 2; y <= size_y; y++) {
1711 cpl_vector_add (vect_mean, cpl_vector_new_from_image_row (median_flat, y));
1712 }
1713 cpl_vector_divide_scalar (vect_mean, size_y);
1714 if (size_x < 60) { // case for LOW res
1715 vector = cpl_vector_filter_median_create(vect_mean, 2);
1716 }
1717 else{ // case for MED and HIGH
1718 vector = cpl_vector_filter_median_create (vect_mean, 5);
1719 }
1720 threshold = cpl_vector_get_max(vector)/100.;
1721
1722 /* Threshold from DARK QC */
1723 cpl_propertylist * p_dark = gravi_data_get_header (dark);
1724
1725 threshold = 20 * cpl_propertylist_get_double (p_dark, QC_DARKRMS_SC);
1726// /* mask gradian */
1727// cpl_image_shift(grad, 1, 0);
1728// cpl_image_subtract(grad, median_flat);
1729// cpl_image_abs(grad);
1730// threshold = cpl_image_get_stdev(grad)*10;
1731
1732 cpl_image_threshold(grad, threshold, threshold, 1., 0.);
1733 cpl_image * mask = cpl_image_duplicate(grad);
1734
1735 /* add mask to FLAT data */
1736 cpl_imagelist * list_mask = cpl_imagelist_new();
1737 cpl_propertylist * plist = cpl_propertylist_duplicate(
1738 gravi_data_get_plist(flat, GRAVI_IMAGING_DATA_SC_EXT));
1739 cpl_propertylist_update_string(plist, "EXTNAME", "FLAT_MASK");
1740 cpl_imagelist_set(list_mask, grad, 0);
1741 gravi_data_add_cube(flat, plist, list_mask);
1742 CPLCHECK_INT("Error mask");
1743
1744 /* Set badpix in profile */
1745 cpl_table * profiles = gravi_data_get_table(flat, "PROFILE_DATA");
1746 int region, nregion = cpl_table_get_ncol(profiles);
1747 char * regname;
1748 cpl_array ** array_data;
1749 double sum;
1750 cpl_size count;
1751 CPLCHECK_INT("Error mask");
1752
1753 printf("%g \n", cpl_image_get(mask, 1, 1, NULL));
1754 CPLCHECK_INT("Error mask");
1755
1756 for (region = 1; region <= nregion; region++) {
1757 regname = cpl_sprintf("DATA%d", region);
1758 printf("region %d \n", region);
1759 array_data=cpl_table_get_data_array(profiles, regname);
1760 count=0;
1761 for (x = 0; x < size_x; x++) {
1762 sum=0;
1763 printf("%d ", x);
1764 for (y = 0; y < size_y; y++) {
1765 printf("%d %g ", y, cpl_image_get(grad, x+1, y+1, NULL));
1766 if (cpl_image_get(grad, x+1, y+1, NULL) > 0.1){
1767 cpl_array_set(array_data[0], x+size_x*y, 0);
1768 count++;
1769 printf("1");
1770 }
1771 sum+=cpl_array_get(array_data[0], x+size_x*y, NULL);
1772 }
1773 if ( sum != 0 ){
1774 for (y = 0; y < size_y; y++) {
1775 cpl_array_set(array_data[0], x+size_x*y,
1776 cpl_array_get(array_data[0], x+size_x*y, NULL)/sum);
1777 }
1778 }
1779 }
1780 }
1781 cpl_msg_info(cpl_func, "Number of flat bad pixels : %lld", count);
1782
1783 cpl_vector_delete(vect_mean);
1784 return(CPL_ERROR_NONE);
1785}
1786
1787
1788
1789/*----------------------------------------------------------------------------*/
1800/*----------------------------------------------------------------------------*/
1801
1802int gravi_write_yorick_batch(const char* filename, const char* input_file, const char* output_file)
1803{
1804 gravi_msg_function_start(0);
1805 FILE *fp;
1806 //char buffer[MAX_LEN_BUFFER];
1807
1808 fp = fopen(filename,"w");
1809
1810 fprintf(fp,
1811"/* \n"
1812"* mira-batch.i -\n"
1813"*\n"
1814"* Mira batch generated by gravi_mira recipe \n"
1815"*/\n"
1816"include, \"mira.i\";\n"
1817"dataset = mira_new(\"%s\");\n"
1818"mira_config, dataset, dim=100, pixelsize=0.4*MIRA_MILLIARCSECOND, xform=\"exact\";\n"
1819"rgl = rgl_new(\"smoothness\");\n"
1820"dim = mira_get_dim(dataset);\n"
1821"img0 = array(double, dim, dim);\n"
1822"img0(dim/2, dim/2) = 1.0;\n"
1823"img1 = mira_solve(dataset, img0, maxeval=500, verb=0, xmin=0.0, normalization=1,\n"
1824" regul=rgl, mu=1e6);\n"
1825"fh = fits_create(\"%s\", overwrite=1, bitpix=-64, dimlist=dimsof(img0)); \n"
1826"fits_write_header, fh; \n"
1827"fits_write_array, fh, img1;\n"
1828"fits_close, fh \n"
1829
1830
1831"quit" , input_file, output_file);
1832
1833 //fputs(buffer, fp);
1834 fclose(fp);
1835
1836 gravi_msg_function_exit(0);
1837 return 0;
1838}
1839
1840
1841static void finals2000A_read_line(FILE *pFile, char *flag, double *mjd, double *pmx, double *pmy, double *dut)
1842{
1843 char line[LINE_SIZE];
1844 fread(line, LINE_SIZE, 1, pFile);
1845 *mjd = atof(line+7);
1846 *pmx = atof(line+18);
1847 *pmy = atof(line+37);
1848 *dut = atof(line+58);
1849 *flag = line[16];
1850}
1851
1852static int finals2000A_mjd_first (FILE *pFile)
1853{
1854 char flag;
1855 double mjd, pmx, pmy, dut;
1856 fseek(pFile, 0, SEEK_SET);
1857 finals2000A_read_line(pFile, &flag, &mjd, &pmx, &pmy, &dut);
1858 return mjd;
1859}
1860
1861static int finals2000A_mjd_last_type(FILE *pFile, char type)
1862{
1863 char flag;
1864 double mjd, pmx, pmy, dut;
1865 fseek(pFile, -LINE_SIZE, SEEK_END);
1866 finals2000A_read_line(pFile, &flag, &mjd, &pmx, &pmy, &dut);
1867 while( flag != type)
1868 {
1869 fseek(pFile, -2*LINE_SIZE, SEEK_CUR);
1870 finals2000A_read_line(pFile, &flag, &mjd, &pmx, &pmy, &dut);
1871 }
1872 fseek(pFile, 0, SEEK_SET);
1873 return mjd;
1874}
1875
1876gravi_data * gravi_eop_load_finals2000A(const char *eop_file)
1877{
1878 gravi_msg_function_start(1);
1879
1880 FILE *pFile;
1881 cpl_size last_mjd, n_entries;
1882 double *mjd, *pmx, *pmy, *dut;
1883 char **flag;
1884 gravi_data * eop_data = NULL;
1885
1886 cpl_ensure (eop_file, CPL_ERROR_NULL_INPUT, NULL);
1887
1888 /* Open finals2000A.data file */
1889 pFile = fopen ((char *)eop_file, "r");
1890
1891 if(pFile != NULL)
1892 {
1893 cpl_msg_info (cpl_func, "Load the file: %s", eop_file);
1894
1895 /* Find number of entries to last predicted entry */
1896 double mjd_P = finals2000A_mjd_last_type (pFile, 'P');
1897 double mjd_I = finals2000A_mjd_last_type (pFile, 'I');
1898 double mjd_S = finals2000A_mjd_first (pFile);
1899 n_entries = mjd_P - mjd_S;
1900 cpl_msg_info(cpl_func, "Reading %lli earth orientation parameters.", n_entries);
1901 cpl_msg_info(cpl_func, " First entry: MJD=%.1f", mjd_S);
1902 cpl_msg_info(cpl_func, " Last IERS entry: MJD=%.1f", mjd_I);
1903 cpl_msg_info(cpl_func, " Last predicted entry: MJD=%.1f", mjd_P);
1904
1905 /* Create tables */
1906 cpl_table * eop_table = cpl_table_new (n_entries);
1907
1908 /* Create columns (filled with zero, thus valid) */
1909 gravi_table_new_column (eop_table, "MJD", "d", CPL_TYPE_DOUBLE);
1910 gravi_table_new_column (eop_table, "PMX", "arcsec", CPL_TYPE_DOUBLE);
1911 gravi_table_new_column (eop_table, "PMY", "arcsec", CPL_TYPE_DOUBLE);
1912 gravi_table_new_column (eop_table, "DUT", "s", CPL_TYPE_DOUBLE);
1913 cpl_table_new_column (eop_table, "FLAG", CPL_TYPE_STRING);
1914 cpl_table_fill_column_window_string (eop_table, "FLAG", 0, n_entries, " ");
1915
1916 mjd = cpl_table_get_data_double (eop_table, "MJD");
1917 pmx = cpl_table_get_data_double (eop_table, "PMX");
1918 pmy = cpl_table_get_data_double (eop_table, "PMY");
1919 dut = cpl_table_get_data_double (eop_table, "DUT");
1920 flag = cpl_table_get_data_string (eop_table, "FLAG");
1921
1922 /* Read finals2000A */
1923 for(int i=0; i<n_entries; i++)
1924 {
1925 finals2000A_read_line (pFile, flag[i]+0, mjd+i, pmx+i, pmy+i, dut+i);
1926 }
1927
1928 /* close file */
1929 fclose (pFile);
1930
1931 /* Build the gravi_data */
1932 eop_data = gravi_data_new (0);
1933 gravi_data_add (eop_data, NULL, eop_table);
1934
1935 /* Create main header */
1936 cpl_propertylist * header = cpl_propertylist_new();
1937 cpl_propertylist_append_double (header, "ESO QC EOP_PARAM MJD_S", mjd_S);
1938 cpl_propertylist_append_double (header, "ESO QC EOP_PARAM MJD_I", mjd_I);
1939 cpl_propertylist_append_double (header, "ESO QC EOP_PARAM MJD_P", mjd_P);
1940 cpl_propertylist_append_double (header, "MJD-OBS", mjd_I);
1941 cpl_propertylist_append_string (header, "ESO PRO CATG", "EOP_PARAM");
1942 cpl_propertylist_append_string (header, "ESO PRO TECH", "CATALOG");
1943 cpl_propertylist_append_string (header, "ESO PRO TYPE", "IERS");
1944 gravi_data_append_header (eop_data, header);
1945 }
1946 else
1947 {
1948 cpl_msg_warning (cpl_func, "Cannot load the file: %s", eop_file);
1949 eop_data = NULL;
1950 }
1951
1952 gravi_msg_function_exit(1);
1953 return eop_data;
1954}
1955
1956// FLAT high frequency
1957
1958/* Compute the flat, mean image of the flat multiplied by the profile of this region */
1959cpl_image * flat_profiled = cpl_image_extract (profile_mean, xmin, ymin, xmax, ymax);
1960cpl_image_multiply (flat_profiled, profile_crop);
1961cpl_image * specMean = cpl_image_collapse_window_create (flat_profiled, 1,1,nxc,nyc,0);
1962cpl_image_delete (flat_profiled);
1963
1964CPLCHECK_NUL("Compute the flat");
1965
1966/* Keep only the high frequencies in the FLAT, that is the detector
1967 * pixel to pixel efficiency difference. So that we remove the detector
1968 * imprint in the data, while keeping the overall instrumental transmission. */
1969if ( nxc > 5000 )
1970{
1971 cpl_msg_debug (cpl_func, "Spectra has >50 pixels -> flat high frequencies");
1972 cpl_image * specFlat = cpl_image_duplicate (specMean);
1973 cpl_mask * kernel = cpl_mask_new (11, 1);
1974 /* Faulty line, replaced by cpl_mask_not(kernel) */
1975 // for (int imask=1;imask<12;imask++) cpl_mask_set (kernel, i, 1, CPL_BINARY_1);
1976 cpl_mask_not (kernel);
1977 cpl_image_filter_mask (specFlat, specMean, kernel, CPL_FILTER_MEDIAN,CPL_BORDER_FILTER);
1978 cpl_image_divide (specMean, specFlat);
1979 cpl_image_delete (specFlat);
1980 cpl_mask_delete (kernel);
1981}
1982else
1983{
1984 cpl_image_fill_window (specMean, 1, 1, nxc, 1, 1.0);
1985 cpl_msg_debug (cpl_func, "Spectra has <50 pixels -> don't flat");
1986}
1987
1988CPLCHECK_NUL("Creating spectrum flat");
1989
1990
1991
1992
1993 /* Fit the dispersion to ARGON */
1994 gravi_data * wave_data = gravi_compute_argon_wave (argon_data, profile_map, dark_map, badpix_map, parlist);
1995 CPLCHECK_CLEAN ("Cannot fit argon");
1996
1997 cpl_table * wave_data_sc = gravi_data_get_table (wave_data, "WAVE_DATA_SC");
1998 cpl_table * wave_map_sc = gravi_data_get_table (wave_map, "WAVE_DATA_SC");
1999
2000 cpl_propertylist * plist = gravi_data_get_plist (wave_data, "WAVE_DATA_SC");
2001 gravi_data_set_propertylist (wave_map, "WAVE_DATA_SC", plist);
2002 gravi_data_set_table (wave_map, "WAVE_DATA_SC", cpl_table_duplicate (wave_data_sc));
2003
2004 cpl_propertylist * header = gravi_data_get_header (wave_map);
2005 plist = gravi_data_get_header (wave_data);
2006 cpl_propertylist_update_double (header, "ESO QC MINWAVE SC", cpl_propertylist_get_double (plist, "ESO QC MINWAVE SC"));
2007 cpl_propertylist_update_double (header, "ESO QC MAXWAVE SC", cpl_propertylist_get_double (plist, "ESO QC MAXWAVE SC"));
2008
2009 gravi_data_save_data (wave_data, "wave_map_argon.fits", CPL_IO_CREATE);
2010 FREE (gravi_data_delete, wave_data);
2011
2012
2013
2014
2015
2016gravi_data * gravi_visdata_fromellipse (gravi_data * spectrum_data)
2017{
2018 gravi_msg_function_start(1);
2019 cpl_ensure (spectrum_data, CPL_ERROR_NULL_INPUT, NULL);
2020
2021 cpl_propertylist * plist = cpl_propertylist_new ();
2022
2023 /* Create the output */
2024 gravi_data * output = gravi_data_new (0);
2025 cpl_propertylist * header = gravi_data_get_header (spectrum_data);
2026 gravi_data_append_header (output, cpl_propertylist_duplicate (header));
2027
2028 /* Loop on SC/FT and polarisation*/
2029 for (int type_data = 0; type_data < 2; type_data ++ ) {
2030 int npol = gravi_pfits_get_pola_num (header, type_data);
2031 for (int pol = 0; pol < npol; pol ++) {
2032
2033 /* Get the data */
2034 cpl_table * spectrum_table = gravi_data_get_spectrum_data (spectrum_data, type_data);
2035 cpl_table * img_det = gravi_data_get_imaging_detector (spectrum_data, type_data);
2036
2037 cpl_size nwave = cpl_table_get_column_depth (spectrum_table,"DATA1");
2038 cpl_size nrow = cpl_table_get_nrow (spectrum_table);
2039 CPLCHECK_NUL ("Cannot get tables");
2040
2041 /* Create output table */
2042 cpl_size nbase = 6;
2043 cpl_table * vis_table = cpl_table_new (nbase*nrow);
2044 gravi_table_new_column (vis_table, "TIME", "us", CPL_TYPE_INT);
2045 gravi_table_new_column_array (vis_table, "VISDATA", "e", CPL_TYPE_FLOAT_COMPLEX, nwave);
2046
2047 /* Create memory for output */
2048 float complex** mVis = cpl_malloc (nbase*nrow * sizeof(float complex*));
2049 for (cpl_size row = 0; row < nbase*nrow; row ++)
2050 mVis[row] = cpl_malloc (nwave * sizeof(float complex));
2051
2052 /* Create vectors over row direction, for the ellipse */
2053 cpl_vector * vectX = cpl_vector_new (nrow);
2054 cpl_vector * vectY = cpl_vector_new (nrow);
2055
2056 /* Get pointer to speed up */
2057 cpl_array ** tVis = cpl_table_get_data_array (vis_table, "VISDATA");
2058 double * datX = cpl_vector_get_data (vectX);
2059 double * datY = cpl_vector_get_data (vectY);
2060
2061 /* Loop on base */
2062 for (int base = 0; base < nbase; base++) {
2063 cpl_msg_info_overwritable (cpl_func,"Ellipse to type=%s, pol=%i over %i, base=%i",
2064 GRAVI_TYPE(type_data), pol+1, npol, base);
2065
2066 /* Set the time */
2067 for (cpl_size row = 0; row < nrow; row ++) {
2068 cpl_size irow = row*nbase + base;
2069 cpl_table_set (vis_table, "TIME", irow,
2070 cpl_table_get (spectrum_table, "TIME", row, NULL));
2071 }
2072
2073 /* Get the sign of the phase for this baseline */
2074 int phi_sign = gravi_table_get_phase_sign (img_det, GRAVI_BASE_TEL[base][0]+1, GRAVI_BASE_TEL[base][1]+1);
2075
2076 /* Get regions */
2077 int regA = gravi_get_region (img_det, base, 'A', pol);
2078 int regB = gravi_get_region (img_det, base, 'B', pol);
2079 int regC = gravi_get_region (img_det, base, 'C', pol);
2080 int regD = gravi_get_region (img_det, base, 'D', pol);
2081 cpl_ensure (regA>=0 && regA>=0 && regA>=0 && regA>=0,
2082 CPL_ERROR_ILLEGAL_INPUT, NULL);
2083
2084 /* Get pointer to data to speed up */
2085 double ** datA = gravi_table_get_data_array_double (spectrum_table, GRAVI_DATA[regA]);
2086 double ** datB = gravi_table_get_data_array_double (spectrum_table, GRAVI_DATA[regB]);
2087 double ** datC = gravi_table_get_data_array_double (spectrum_table, GRAVI_DATA[regC]);
2088 double ** datD = gravi_table_get_data_array_double (spectrum_table, GRAVI_DATA[regD]);
2089 CPLCHECK_NUL ("Cannot get data");
2090
2091 /* Loop on wavelength because we need the
2092 * ellipse vectors in the row direction */
2093 for (cpl_size wave = 0 ; wave < nwave ; wave++) {
2094
2095 /* Create the ellipse vectors as:
2096 * X = C-A and Y = D-B */
2097 for (cpl_size row = 0; row < nrow; row ++) {
2098 datX[row] = datC[row][wave] - datA[row][wave];
2099 datY[row] = datD[row][wave] - datB[row][wave];
2100 }
2101
2102 /* Re-center the ellipse to ease fit */
2103 cpl_vector_subtract_scalar (vectY, cpl_vector_get_mean (vectY));
2104 cpl_vector_subtract_scalar (vectX, cpl_vector_get_mean (vectX));
2105
2106 /* Get the phase by fitting the ellipse */
2107 cpl_vector * phase = gravi_vectors_phase_create (vectX, vectY);
2108 cpl_vector_multiply_scalar (phase, phi_sign);
2109
2110 /* Set the VISDATA = sqrt (X^2 + Y^2) * exp (i.phase) */
2111 for (cpl_size row = 0; row < nrow; row ++) {
2112 cpl_size irow = row*nbase + base;
2113 mVis[irow][wave] = (float complex)
2114 (sqrt (datX[row]*datX[row] + datY[row]*datY[row]) *
2115 cexp (1.*I * cpl_vector_get (phase, row)));
2116 }
2117
2118 FREE (cpl_vector_delete, phase);
2119 } /* End loop on wave */
2120
2121 /* Wrap the array into the table */
2122 for (cpl_size row = 0; row < nrow*nbase; row ++)
2123 tVis[row] = cpl_array_wrap_float_complex (mVis[row], nwave);
2124
2125 /* Free pointer to data array */
2126 FREE (cpl_free, datA);
2127 FREE (cpl_free, datB);
2128 FREE (cpl_free, datD);
2129 FREE (cpl_free, datD);
2130
2131 } /* End loop on base */
2132
2133 /* Desalocate ellipses vectors */
2134 FREE (cpl_vector_delete, vectX);
2135 FREE (cpl_vector_delete, vectY);
2136
2137 /* Add table to data with EXTNAME */
2138 cpl_propertylist_append_string (plist, "EXTNAME", "ELLIPSE_VIS");
2139 cpl_propertylist_append_string (plist, "INSNAME", GRAVI_INSNAME(type_data,pol,npol));
2140 cpl_propertylist_append_int (plist, "EXTVER", GRAVI_EXTVER(type_data,pol,npol));
2141
2142 gravi_data_add (output, plist, vis_table);
2143
2144
2145 } /* End loop on pol */
2146 } /* End loop on SC/FT */
2147
2148 gravi_msg_function_exit(1);
2149 return output;
2150}
2151
2152double gravi_pfits_get_lambdamet(const cpl_propertylist * plist)
2153{
2154 cpl_errorstate prestate = cpl_errorstate_get();
2155 double value = cpl_propertylist_get_double(plist, "ESO INS MLC WAVELENG");
2156 cpl_ensure (cpl_errorstate_is_equal(prestate), cpl_error_get_code(), 0.0);
2157 return value;
2158}
2159
2160
2161/*----------------------------------------------------------------------------*/
2172cpl_error_code gravi_smooth_preproc (gravi_data * data,
2173 const cpl_parameterlist * parlist)
2174{
2175 cpl_table * detector_table, * spectrum_table;
2176 cpl_size row, n_row, reg, n_region;
2177 cpl_array * output = NULL;
2178
2179 /* Verbose */
2180 gravi_msg_function_start(1);
2181 cpl_ensure_code (data, CPL_ERROR_NULL_INPUT);
2182 cpl_ensure_code (parlist, CPL_ERROR_NULL_INPUT);
2183
2184 /* Get the nsmooth from parameter list */
2185 int nsmooth = cpl_parameter_get_int (cpl_parameterlist_find_const (parlist, "gravi.nsmooth_sc"));
2186
2187 /* Default from spectral resolution */
2188 if (nsmooth < 0) {
2189 const char * resolution = gravi_pfits_get_spec_res (gravi_data_get_header (data));
2190
2191 if ( !(strcmp(resolution, "HIGH")) ) {
2192 cpl_msg_info(cpl_func,"Default smoothing is 15 (HIGH)");
2193 nsmooth = 15;
2194 }
2195 else if ( !(strcmp(resolution, "MEDIUM")) ) {
2196 cpl_msg_info(cpl_func,"Default smoothing is 3 (MEDIUM)");
2197 nsmooth = 3;
2198 }
2199 else if ( !(strcmp(resolution, "LOW")) ) {
2200 cpl_msg_info(cpl_func,"Default is no smoothing of data (LOW)");
2201 nsmooth = 0;
2202 } else {
2203 cpl_msg_warning(cpl_func,"Unknown spectral resolution thus no smoothing of data");
2204 nsmooth = 0;
2205 }
2206 }
2207
2208 /* Case no smoothing */
2209 if ( nsmooth <= 0 ) {
2210 cpl_msg_info (cpl_func, "End function (no smoothing)");
2211 return CPL_ERROR_NONE;
2212 }
2213
2214 /* Search the number of regions for SC */
2215 detector_table = gravi_data_get_table (data, GRAVI_IMAGING_DETECTOR_SC_EXT);
2216 n_region = cpl_table_get_nrow (detector_table);
2217
2218 /* Get the SC table data */
2219 spectrum_table = gravi_data_get_table (data, GRAVI_SPECTRUM_DATA_SC_EXT);
2220 n_row = cpl_table_get_nrow (spectrum_table);
2221
2222 CPLCHECK_MSG ("Cannot get data");
2223
2224 /* Loop on SC regions */
2225 for (reg = 0; reg < n_region; reg++) {
2226
2227 /* Load pointer to all rows */
2228 cpl_msg_info_overwritable(cpl_func,"Smooth region %lld over %lld of SC", reg+1, n_region);
2229 cpl_array** arrays = cpl_table_get_data_array (spectrum_table, GRAVI_DATA[reg]);
2230
2231 /* Loop on rows */
2232 for (row = 0; row < n_row; row++) {
2233 /* Compute a smoothed version of the spectra of this row */
2234 output = gravi_array_smooth (arrays[row], nsmooth);
2235 /* Put it back inplace in the table */
2236 cpl_array_delete (arrays[row]);
2237 arrays[row] = output;
2238 }
2239 /* End loop on rows */
2240 }
2241 /* End loop on region */
2242
2243 /* Add the NSMOOTH parameter */
2244 cpl_propertylist * hdr_data = gravi_data_get_header (data);
2245 cpl_propertylist_append_int (hdr_data, "ESO QC P2VM NSMOOTH SC", nsmooth);
2246 cpl_propertylist_set_comment (hdr_data, "ESO QC P2VM NSMOOTH SC", "nb of smoothed channels in P2VM computation");
2247
2248 /* Verbose */
2249 gravi_msg_function_exit(1);
2250 return CPL_ERROR_NONE;
2251}
2252
2253
2254cpl_vector * gravi_construction_opd_phase (cpl_table * opl_table,
2255 cpl_table * phase_sc,
2256 cpl_table * phase_ft,
2257 double dit_sc)
2258{
2259 int ind_sc, nrow, row, nbase = 6, base;
2260 double exptime, exptime_sc, opl, time_sc, time_ft, time_metrology;
2261 int nv, comp, nrow_met, nrow_sc;
2262 int tel_1[6] = {0,0,0,1,1,2};
2263 int tel_2[6] = {1,2,3,2,3,3};
2264 const cpl_array * time_array_ft, * time_array_sc;
2265
2266 gravi_msg_function_start(0);
2267 cpl_ensure (opl_table, CPL_ERROR_NULL_INPUT, NULL);
2268 cpl_ensure (phase_sc, CPL_ERROR_NULL_INPUT, NULL);
2269 cpl_ensure (phase_ft, CPL_ERROR_NULL_INPUT, NULL);
2270
2271 /* Get the number of acquisitions */
2272 time_array_ft = cpl_table_get_array (phase_ft, "TIME", 0);
2273 time_array_sc = cpl_table_get_array (phase_sc, "TIME", 0);
2274 nrow = cpl_array_get_size (time_array_ft);
2275 nrow_sc = cpl_array_get_size (time_array_sc);
2276 nrow_met = cpl_table_get_nrow (opl_table);
2277
2278 CPLCHECK_NUL("Cannot get data");
2279
2280 /* allocate the matrices */
2281 gsl_matrix * A_data = gsl_matrix_calloc (nbase * nrow, 8), * A;
2282 gsl_matrix * U;
2283 gsl_vector * bis_data = gsl_vector_alloc (nbase * nrow), * bis;
2284 gsl_matrix * X = gsl_matrix_alloc (8, 8),
2285 * V = gsl_matrix_alloc (8, 8);
2286 gsl_vector * S = gsl_vector_alloc (8);
2287 gsl_vector * work = gsl_vector_alloc (8), * x = gsl_vector_alloc (8);
2288 CPLCHECK_NUL("Allocate matrix");
2289
2290 /* compute the periods of the signals */
2291 exptime = cpl_array_get_int (time_array_ft, 1, &nv) -
2292 cpl_array_get_int (time_array_ft, 0, &nv);
2293 exptime_sc = cpl_array_get_int (time_array_sc, 1, &nv) -
2294 cpl_array_get_int (time_array_sc, 0, &nv);
2295 //double exptime_met = cpl_table_get_int (opl_table, "TIME", 1, &nv) -
2296 // cpl_table_get_int (opl_table, "TIME", 0, &nv);
2297
2298 /*
2299 * Extract the OPl metrology and compute the mean of each exposure time
2300 */
2301
2302 for (base = 0; base < nbase; base ++){
2303 ind_sc = 0;
2304 int im = 0;
2305 const cpl_array * ft_array = cpl_table_get_array (phase_ft, "PHASE", base);
2306 const cpl_array * sc_array = cpl_table_get_array (phase_sc, "PHASE", base);
2307 char * opl1 = cpl_sprintf("OPL%d", tel_1[base]+1);
2308 char * opl2 = cpl_sprintf("OPL%d", tel_2[base]+1);
2309
2310
2311 cpl_array * temp = cpl_array_wrap_double (
2312 cpl_table_get_data_double (opl_table, opl1), nrow_met);
2313 cpl_array * opl_met1 = cpl_array_duplicate (temp);
2314 cpl_array_unwrap (temp);
2315 cpl_array * opl_met2 = cpl_array_wrap_double (
2316 cpl_table_get_data_double (opl_table, opl2), nrow_met);
2317 cpl_array_subtract (opl_met1, opl_met2);
2318 CPLCHECK_NUL ("Opd diff");
2319
2320 /*
2321 * resample the SC and Metrology at FT period
2322 */
2323 for (row = 0; row < nrow; row ++){ // loop on row FT
2324
2325 gsl_matrix_set (A_data, base * nrow + row, 6,
2326 - cpl_array_get_double (ft_array, row, &nv));
2327 CPLCHECK_NUL ("Set FT phase in matrix");
2328
2329 time_ft = cpl_array_get_int (time_array_ft, row, &nv);
2330 if (ind_sc < cpl_array_get_size (sc_array)){
2331 time_sc = cpl_array_get_int (time_array_sc, ind_sc, &nv);
2332 }
2333 CPLCHECK_NUL ("Get time sc");
2334
2335 /* increase ind ft if time_ft in the range of the nex DIT */
2336 while(time_ft > (time_sc+exptime_sc/2.)){
2337 ind_sc++;
2338 if (ind_sc < cpl_array_get_size (sc_array))
2339 time_sc = cpl_array_get_int (time_array_sc, ind_sc, &nv);
2340 else
2341 break;
2342 }
2343 CPLCHECK_NUL ("Get next sc dit");
2344
2345 /* get phi sc */
2346 double phi_sc=0;
2347 if (ind_sc > 0){
2348 if(ind_sc < nrow_sc) {
2349 phi_sc=cpl_array_get_double (sc_array, ind_sc, &nv);
2350 }
2351 else {
2352 phi_sc=cpl_array_get_double (sc_array, nrow_sc - 1, &nv);
2353 }
2354 }
2355 else {
2356 phi_sc = cpl_array_get_double (sc_array, 0, &nv);
2357 }
2358 CPLCHECK_NUL ("Get phi_sc");
2359
2360 gsl_matrix_set (A_data, base * nrow + row, 7, phi_sc);
2361 gsl_matrix_set (A_data, base * nrow + row, base, 1.0);
2362
2363
2364 /* Metrology case */
2365 opl = 0;
2366 comp = 0;
2367
2368 if (im < nrow_met)
2369 time_metrology = cpl_table_get_int (opl_table, "TIME", im, &nv);
2370
2371 while ((time_metrology < (time_ft + exptime))){ //((time_metrology + exptime_met)< time_ft){
2372
2373 if (im < nrow_met) {
2374 opl += cpl_array_get_double (opl_met1, im, &nv); //cpl_vector_get (vector_opd, im);
2375 comp ++;
2376 }
2377
2378 im++;
2379 if (im < nrow_met) {
2380 time_metrology = cpl_table_get_int (opl_table, "TIME", im, &nv);
2381 }
2382 else {
2383 break;
2384 }
2385
2386 CPLCHECK_NUL("Get time metrology");
2387 }
2388
2389 /* average the metrology over the FT DIT */
2390 if (comp != 0)
2391 gsl_vector_set (bis_data, base * nrow + row, opl/comp);
2392
2393 /* if no metrology signal within the FT DIT interpolate the metrology */
2394 else {
2395 if (im > 0){
2396 if(im < nrow_met) {
2397 opl = (cpl_array_get_double (opl_met1, im, &nv) -
2398 cpl_array_get_double (opl_met1, im - 1, &nv)) * (time_ft -
2399 cpl_table_get_int (opl_table, "TIME", im - 1, &nv))
2400 /(time_metrology -
2401 cpl_table_get_int (opl_table, "TIME", im - 1, &nv)) +
2402 cpl_array_get_double (opl_met1, im - 1, &nv);
2403
2404 gsl_vector_set (bis_data, base * nrow + row, opl);
2405 }
2406 else {
2407 gsl_vector_set (bis_data, base * nrow + row,
2408 cpl_array_get_double (opl_met1, nrow_met - 1, &nv));
2409 }
2410 }
2411 else {
2412 gsl_vector_set (bis_data, base * nrow + row,
2413 cpl_array_get_double (opl_met1, 0, &nv));
2414 }
2415 CPLCHECK_NUL("Interpolate the metrology");
2416 }
2417 } /* end loop on row FT */
2418
2419 cpl_array_delete (opl_met1);
2420 cpl_array_unwrap (opl_met2);
2421 cpl_free(opl1);
2422 cpl_free(opl2);
2423 }
2424
2425 /*
2426 * filter the data out of the SC integration time
2427 */
2428
2429 long n_row_A=nrow_sc*(int)(dit_sc/exptime+1);
2430 cpl_vector *i_A_vector = cpl_vector_new(n_row_A);
2431 int i_A=0;
2432 int time_mod;
2433 int t0_sc=cpl_array_get_int (time_array_sc, 0, &nv);
2434 int tend_sc=cpl_array_get_int (time_array_sc, nrow_sc-1, &nv);
2435
2436 /* Find the index of the frames within the SC integration */
2437 for (row=0; row<nrow; row++){
2438 time_ft = cpl_array_get_int (time_array_ft, row, &nv);
2439 if ((time_ft >= t0_sc-dit_sc/2) && (time_ft < tend_sc+dit_sc/2)){
2440 time_mod=((int)(time_ft-(t0_sc-dit_sc/2)) % (int)exptime_sc);
2441 if ( time_mod >= 0 && time_mod < dit_sc ) {
2442 cpl_vector_set (i_A_vector, i_A, row);
2443 i_A++;
2444 }
2445 }
2446 }
2447 n_row_A=i_A;
2448
2449 /* copy the frames within SC intergration in the matrix A */
2450 A = gsl_matrix_alloc (nbase * n_row_A, 8);
2451 bis = gsl_vector_alloc (nbase * n_row_A);
2452 cpl_vector *time_A = cpl_vector_new(n_row_A);
2453 int col;
2454
2455 for (base = 0; base < nbase; base ++)
2456 for (i_A = 0; i_A < n_row_A; i_A++){
2457 row = cpl_vector_get(i_A_vector, i_A);
2458 for (col = 0; col < 8; col ++)
2459 gsl_matrix_set (A, base * n_row_A + i_A, col, gsl_matrix_get (A_data, base * nrow + row, col));
2460
2461 gsl_vector_set (bis, base * n_row_A + i_A, gsl_vector_get (bis_data, base * nrow + row));
2462 if (base == 0)
2463 cpl_vector_set (time_A, i_A, cpl_array_get_int (time_array_ft, row, &nv));
2464 }
2465
2466 cpl_vector_delete(i_A_vector);
2467
2468 /*
2469 * Solve the linear equation Ax = b to get the coefficients to deduce the OPDs
2470 */
2471 nrow=n_row_A;
2472 U = gsl_matrix_alloc (nbase * nrow, 8);
2473 gsl_matrix_memcpy (U, A);
2474
2475
2476 gsl_linalg_SV_decomp (U, V, S, work);
2477 gsl_linalg_SV_solve (U, V, S, bis, x);
2478 cpl_vector * opd_coeff = cpl_vector_new (2);
2479 cpl_vector_set (opd_coeff, 0, gsl_vector_get (x, 7));
2480 cpl_vector_set (opd_coeff, 1, gsl_vector_get (x, 6));
2481
2482 cpl_msg_debug(cpl_func, "Wavelength base %d => SC = %g, FT = %g\n",
2483 base, -cpl_vector_get(opd_coeff,0)*2*M_PI, cpl_vector_get(opd_coeff,1)*2*M_PI);
2484
2485 // TODO
2486 if (PLOT_MET_PHASE_FIT)
2487
2488 {
2489 cpl_msg_info(cpl_func, "Plot fit residuals");
2490 cpl_errorstate prestate = cpl_errorstate_get();
2491 gsl_matrix_memcpy (U, A);
2492
2493 const cpl_vector ** vectors=malloc(3 * sizeof(cpl_vector*));
2494 vectors[0]=NULL;
2495 cpl_vector *vect_phase_sc_ft=cpl_vector_new(nrow*nbase);
2496 cpl_vector *vect_dopd_met=cpl_vector_new(nrow*nbase);
2497 cpl_vector *vect_diff=cpl_vector_new(nrow*nbase);
2498 for (row = 0; row < nrow*nbase; row ++){
2499 cpl_vector_set(vect_dopd_met, row, gsl_vector_get (bis, row));
2500 cpl_vector_set(vect_phase_sc_ft, row, gsl_matrix_get (A, row, 0)*gsl_vector_get (x, 0)+
2501 gsl_matrix_get (A, row, 1)*gsl_vector_get (x, 1)+
2502 gsl_matrix_get (A, row, 2)*gsl_vector_get (x, 2)+
2503 gsl_matrix_get (A, row, 3)*gsl_vector_get (x, 3)+
2504 gsl_matrix_get (A, row, 4)*gsl_vector_get (x, 4)+
2505 gsl_matrix_get (A, row, 5)*gsl_vector_get (x, 5)+
2506 gsl_matrix_get (A, row, 6)*gsl_vector_get (x, 6)+
2507 gsl_matrix_get (A, row, 7)*gsl_vector_get (x, 7));
2508 cpl_vector_set(vect_diff, row, cpl_vector_get(vect_dopd_met, row)-cpl_vector_get(vect_phase_sc_ft, row));
2509 }
2510
2511 vectors[1]=vect_phase_sc_ft;
2512 vectors[2]=vect_dopd_met;
2513
2514 cpl_plot_vectors (cpl_sprintf("set title 'Met fit case %d a=%g b=%g'; set xlabel 'time[10-6s]'; set ylabel 'Phase_ft+Phase_sc, dopd_met';", base+1, gsl_vector_get (x, 0), gsl_vector_get (x, 1)),
2515 "", "", vectors, 3);
2516 cpl_plot_vector(cpl_sprintf("set title 'Met fit case %d a=%g b=%g'; set xlabel 'time[10-6s]'; set ylabel 'Phase_ft+Phase_sc-dopd_met';", base+1, gsl_vector_get (x, 0), gsl_vector_get (x, 1)),
2517 "", "", vect_diff);
2518 cpl_vector_delete(vect_phase_sc_ft);
2519 cpl_vector_delete(vect_dopd_met);
2520 cpl_vector_delete(vect_diff);
2521 cpl_free(vectors);
2522 cpl_errorstate_set (prestate);
2523 }
2524
2525 gsl_matrix_free (A);
2526 gsl_matrix_free (A_data);
2527 gsl_matrix_free (U);
2528 gsl_matrix_free (V);
2529 gsl_vector_free (x);
2530 gsl_vector_free (bis);
2531 gsl_vector_free (bis_data);
2532 cpl_vector_delete(time_A);
2533 gsl_matrix_free (X);
2534 gsl_vector_free (S);
2535 gsl_vector_free (work);
2536
2537 gravi_msg_function_exit(0);
2538 return opd_coeff;
2539}
2540
2541
2542/*----------------------------------------------------------------------------*/
2560/*----------------------------------------------------------------------------*/
2561
2562cpl_table * gravi_metrology_calibration (cpl_table * metrology_table,
2563 cpl_table * opl_table,
2564 double mjd_obs)
2565{
2566 gravi_msg_function_start(1);
2567 cpl_ensure (metrology_table, CPL_ERROR_NULL_INPUT, NULL);
2568 cpl_ensure (opl_table, CPL_ERROR_NULL_INPUT, NULL);
2569
2570 cpl_table * p2vm_met;
2571
2572 int row, i, nb_row;
2573 cpl_vector * vectA, * vectB,
2574 * phase, * y_sigma, * init_val;
2575 cpl_array * volt;
2576 cpl_array * temp, * coherence_fit, * phase_fit;
2577 const cpl_array * volt0;
2578 cpl_matrix * opd_matrix;
2579 int tel, infos = 0, n_tel;
2580 cpl_vector ** vect;
2581 int val_to_fit2[] = {1,1,1,0}, nv;
2582 double mse, red_chisq;
2583 cpl_bivector *plot;
2584 char* ps_string;
2585 char * opl;
2586 cpl_array * trans;
2587
2588 /* Extract inputs */
2589
2590 nb_row = cpl_table_get_nrow (opl_table);
2591 for (tel = 0; tel < 4; tel ++){
2592 opl = cpl_sprintf("OPL%d", tel + 1);
2593 cpl_table_new_column (opl_table, opl, CPL_TYPE_DOUBLE);
2594 cpl_free (opl);
2595 }
2596
2597 /* Define the first ABCD k band and in witch case it starts */
2598 volt0 = cpl_table_get_array (metrology_table, "VOLT", 0);
2599
2600
2601 n_tel = cpl_array_get_size (volt0);
2602
2603 /* Creating the 4 columns (REGNAME, TRANSMISSION, COHERENCE and PHASE) */
2605
2606 int start_tel = n_tel/2 - 8;// 0
2607 p2vm_met = cpl_table_new (n_tel / 2 - start_tel);
2608 cpl_table_new_column_array (p2vm_met,"TRANSMISSION", CPL_TYPE_DOUBLE, 2);
2609 cpl_table_new_column_array (p2vm_met,"COHERENCE", CPL_TYPE_DOUBLE, 2);
2610 cpl_table_new_column_array (p2vm_met,"PHASE", CPL_TYPE_DOUBLE, 2);
2611
2612 /* Create a new table to stoke the OPL of each telescope */
2613 cpl_table_new_column (opl_table, "TIME", CPL_TYPE_INT);
2614
2615 if (cpl_error_get_code()){
2616
2617 printf("error %f %s and %s \n", 1.0, cpl_error_get_message(), cpl_error_get_where());
2618 return NULL;
2619 }
2620
2621
2622 for (tel = n_tel/2 - 8; tel < n_tel/2; tel++){
2623
2624 vectA = cpl_vector_new (nb_row);
2625 vectB = cpl_vector_new (nb_row);
2626
2627 /* Get the index of the k band AB */
2628
2629 /* Define the first AB k band and in witch case it starts */
2630 for (row = 0; row < nb_row; row ++){
2631 temp = cpl_array_duplicate (cpl_table_get_array (metrology_table,
2632 "VOLT", row));
2633 volt = cpl_array_cast (temp, CPL_TYPE_DOUBLE);
2634 cpl_array_delete (temp);
2635 cpl_vector_set (vectA, row, cpl_array_get_double (volt, 2*tel, &nv));
2636 cpl_vector_set (vectB, row, cpl_array_get_double (volt, 2*tel + 1, &nv));
2637
2638 cpl_array_delete (volt);
2639 }
2640
2641 /* Fit OPD */
2642 /* Compute the OPD */
2643
2644 phase = gravi_vectors_phase_create (vectA, vectB);
2645
2646 if (PLOT_MET_ELLIPSE)
2647 {
2648 if (POSTSCRIPT_PLOT) ps_string=cpl_sprintf("set term png; set output 'plot_met_ellipse_T%d.png';", tel);
2649 else ps_string=cpl_sprintf(" ");
2650 plot = cpl_bivector_wrap_vectors (vectA, vectB);
2651 cpl_plot_bivector (cpl_sprintf("set title 'Met ellipse Tel %d'; set xlabel 'C-A'; set ylabel 'D-B';", tel+1),
2652 NULL, NULL, plot);
2653
2654 if (POSTSCRIPT_PLOT) ps_string=cpl_sprintf("set term png; set output 'plot_met_phase_T%d.png';", tel);
2655 else ps_string=cpl_sprintf(" ");
2656 cpl_plot_vector (cpl_sprintf("set title 'Met phase Tel %d'; set xlabel 'index'; set ylabel 'Phase';", tel+1),
2657 NULL, NULL, phase);
2658
2659 cpl_free(ps_string);
2660 }
2661
2662
2663 if (cpl_error_get_code()){
2664 printf("error %f %s and %s \n", 4.0, cpl_error_get_message(), cpl_error_get_where());
2665
2666 return NULL;
2667 }
2668 /* Get the OPD vector from the phase */
2669
2670 opd_matrix = cpl_matrix_new(nb_row, 1);
2671
2672 for (i = 0; i < nb_row; i++){
2673
2674 cpl_matrix_set(opd_matrix, i, 0,
2675 cpl_vector_get(phase, i) * LAMBDA_MET/(2*M_PI));
2676 if (tel >= n_tel/2 - 8){
2677
2678 if (tel < n_tel/2 - 4){
2679 opl = cpl_sprintf("OPL%d", tel - (n_tel/2 - 8) + 1);
2680 cpl_table_set_double (opl_table, opl, i,
2681 cpl_vector_get (phase, i) * LAMBDA_MET/(2*M_PI));
2682 cpl_free (opl);
2683 }
2684 else {
2685 opl = cpl_sprintf("OPL%d", tel - (n_tel/2 - 4) + 1);
2686 double opl_ft = cpl_table_get_double (opl_table, opl, i, &nv);
2687 cpl_table_set_double (opl_table, opl, i, (cpl_vector_get (phase, i) * LAMBDA_MET/(2*M_PI))- opl_ft);
2688 cpl_free (opl);
2689 }
2690
2691 if (tel == n_tel/2 - 4)
2692 cpl_table_set_int (opl_table, "TIME", i,
2693 cpl_table_get_int (metrology_table, "TIME",
2694 i, &nv) + mjd_obs);
2695 }
2696 if (cpl_error_get_code()){
2697 printf("error %f %s and %s \n", 5.0, cpl_error_get_message(), cpl_error_get_where());
2698
2699 return NULL;
2700 }
2701 }
2702
2703 cpl_vector_delete(phase);
2704
2705 /* End fit opd */
2706
2707 /* fit on a central window of 5*lambda width */
2708
2709 vect = cpl_malloc(2*sizeof(cpl_vector*));
2710 vect[0] = vectA;
2711 vect[1] = vectB;
2712 phase_fit = cpl_array_new(2, CPL_TYPE_DOUBLE);
2713 coherence_fit = cpl_array_new(2, CPL_TYPE_DOUBLE);
2714 trans = cpl_array_new(2, CPL_TYPE_DOUBLE);
2715
2716 for (i = 0; i < 2; i++){
2717
2718 /* Fit to get the phase and the coherence */
2719
2720
2721 /* Get the spectrum of the vector region */
2722 y_sigma = cpl_vector_new(nb_row);
2723 cpl_vector_fill(y_sigma, 1);
2724
2725 /* Define and initialize all variables to make a FIT */
2726 mse = 0;
2727 red_chisq = 0;
2728 init_val = cpl_vector_new(4);
2729 cpl_vector_set(init_val, 0, 1);
2730 cpl_vector_set(init_val, 1, 1);
2731 cpl_vector_set(init_val, 2, 1);
2732 cpl_vector_set(init_val, 3, LAMBDA_MET);
2733
2734 cpl_errorstate prestate = cpl_errorstate_get();
2735 cpl_fit_lvmq(opd_matrix, NULL, vect[i],
2736 y_sigma, init_val, val_to_fit2, &sin_lambda,
2737 &dfda_sin, CPL_FIT_LVMQ_TOLERANCE, CPL_FIT_LVMQ_COUNT,
2738 CPL_FIT_LVMQ_MAXITER, &mse, &red_chisq, NULL);
2739
2740 if (cpl_error_get_code()){
2741 printf("error %f %s and %s \n", 6.0, cpl_error_get_message(), cpl_error_get_where());
2742 return NULL;
2743 }
2744
2745 if (!strcmp("The iterative process did not converge",
2746 cpl_error_get_message())){
2747 if (infos)
2748 cpl_msg_info(cpl_func, "The iterative process "
2749 "did not converge");
2750 cpl_errorstate_set (prestate);
2751 }
2752 if (infos)
2753 cpl_msg_info(cpl_func, "tel = %d : mse "
2754 "%g chi2 %g", tel, mse, red_chisq);
2755
2756
2757
2758 /* Compute the P2VM value */
2759 cpl_array_set_double (coherence_fit, i,
2760 sqrt( pow( cpl_vector_get(init_val, 2), 2) +
2761 pow( cpl_vector_get(init_val, 1), 2)));
2762
2763 cpl_array_set_double (phase_fit, i, atan2( cpl_vector_get(init_val, 2),
2764 cpl_vector_get(init_val, 1)));
2765
2766 cpl_array_set_double (trans, i, cpl_vector_get(init_val, 0));
2767
2768 if (cpl_error_get_code()){
2769 printf("error %f %s and %s \n", 7.0, cpl_error_get_message(), cpl_error_get_where());
2770 return NULL;
2771 }
2772
2773 cpl_vector_delete(init_val);
2774 cpl_vector_delete(y_sigma);
2775 cpl_vector_delete (vect[i]);
2776 }
2777
2778 cpl_free (vect);
2779 cpl_matrix_delete(opd_matrix);
2780
2781 cpl_table_set_array (p2vm_met, "COHERENCE", tel - start_tel, coherence_fit);
2782 cpl_array_subtract_scalar (phase_fit,
2783 cpl_array_get_double (phase_fit, 0, &nv));
2784 cpl_table_set_array (p2vm_met, "PHASE", tel - start_tel, phase_fit);
2785 cpl_table_set_array (p2vm_met, "TRANSMISSION", tel - start_tel, trans);
2786
2787 cpl_array_delete (trans);
2788 cpl_array_delete (coherence_fit);
2789 cpl_array_delete (phase_fit);
2790 }
2791
2792 /* Verbose */
2793 gravi_msg_function_exit(1);
2794 return p2vm_met;
2795
2796}
2797
2798
2799
2800/*----------------------------------------------------------------------------*/
2811/*----------------------------------------------------------------------------*/
2812
2813cpl_error_code gravi_data_mean_metFddl (gravi_data * oi_vis, gravi_data * p2vmred_data,
2814 gravi_data * preproc_data,
2815 int thread){
2816
2817 cpl_table * visMet_data, * fddl_data, * spectrum_sc;
2818 cpl_propertylist * primary_hdr;
2819 int nbrow_met, row, nbrow_sc;
2820 int nv;
2821 cpl_table * fddl_met_mean;
2822
2823 gravi_msg_function_start(1);
2824 cpl_ensure_code (oi_vis, CPL_ERROR_NULL_INPUT);
2825 cpl_ensure_code (p2vmred_data, CPL_ERROR_NULL_INPUT);
2826 cpl_ensure_code (preproc_data, CPL_ERROR_NULL_INPUT);
2827 cpl_ensure_code (thread>=0, CPL_ERROR_ILLEGAL_INPUT);
2828
2829 /*
2830 * Construction of the table containing columns of the mean FT_POS,
2831 * SC_POS and the phase FC of each dispersion file
2832 */
2833
2834 fddl_met_mean = cpl_table_new (1);
2835 cpl_table_new_column_array (fddl_met_mean, "FT_POS", CPL_TYPE_DOUBLE, 4);
2836 cpl_table_new_column_array (fddl_met_mean, "SC_POS", CPL_TYPE_DOUBLE, 4);
2837 cpl_table_new_column_array (fddl_met_mean, "PHASE_FC", CPL_TYPE_DOUBLE, 4);
2838 cpl_table_new_column (fddl_met_mean, "TIME", CPL_TYPE_INT);
2839
2840 /* Get the data */
2841 visMet_data = gravi_data_get_table (p2vmred_data, GRAVI_OI_VIS_MET_EXT);
2842
2843 primary_hdr = gravi_data_get_plist (p2vmred_data, GRAVI_PRIMARY_HDR_EXT);
2844 spectrum_sc = gravi_data_get_table (preproc_data, GRAVI_SPECTRUM_DATA_SC_EXT);
2845 fddl_data = gravi_data_get_table (p2vmred_data, GRAVI_FDDL_EXT);
2846
2847 CPLCHECK_MSG ("ERROR1");
2848
2849 nbrow_met = cpl_table_get_nrow (visMet_data);
2850 nbrow_sc = cpl_table_get_nrow (spectrum_sc);
2851
2852 CPLCHECK_MSG ("ERROR2");
2853
2854 int * time = cpl_table_get_data_int (visMet_data, "TIME");
2855 int exptime_sc = cpl_table_get_int (spectrum_sc, "TIME", 1, &nv) -
2856 cpl_table_get_int (spectrum_sc, "TIME", 0, &nv);
2857 cpl_array ** met_arrays = cpl_table_get_data_array (visMet_data, "PHASE_FC");
2858 cpl_array * met_mean = cpl_array_new (4, CPL_TYPE_DOUBLE);
2859 cpl_array_fill_window (met_mean, 0, 4, 0.0);
2860 CPLCHECK_MSG ("ERROR3");
2861
2862 /*
2863 * Compute the mean of the metrology included between the SC DIT
2864 */
2865
2866 int comp = 0;
2867 for (row = 0; row < nbrow_met; row ++){
2868
2869 if ((time[row] > cpl_table_get_int (spectrum_sc, "TIME", 0, &nv) - exptime_sc/2) &&
2870 (time[row] < cpl_table_get_int (spectrum_sc, "TIME", nbrow_sc-1, &nv) + exptime_sc/2)){
2871 cpl_array_add (met_mean, met_arrays[row]);
2872 comp ++;
2873 }
2874 CPLCHECK_MSG ("ERROR4");
2875 }
2876
2877 if (comp != 0)
2878 cpl_array_divide_scalar (met_mean, comp);
2879
2880 /*
2881 * Save the mean metrology
2882 */
2883
2884 cpl_table_set_array (fddl_met_mean, "PHASE_FC", 0, met_mean);
2885 cpl_table_set_int (fddl_met_mean, "TIME", 0,
2886 cpl_table_get_int (spectrum_sc, "TIME", 0, &nv) );
2887 cpl_array_delete (met_mean);
2888
2889 int nbrow_fddl = cpl_table_get_nrow (fddl_data);
2890
2891 /*
2892 * Compute the mean FT po and SC pos
2893 */
2894
2895 cpl_array ** fddl_ft = cpl_table_get_data_array (fddl_data, "FT_POS");
2896 cpl_array ** fddl_sc = cpl_table_get_data_array (fddl_data, "SC_POS");
2897 time = cpl_table_get_data_int (fddl_data, "TIME");
2898 cpl_array * ft_pos = cpl_array_new (4, CPL_TYPE_DOUBLE);
2899 cpl_array * sc_pos = cpl_array_new (4, CPL_TYPE_DOUBLE);
2900 cpl_array_fill_window (ft_pos, 0, 4, 0.0);
2901 cpl_array_fill_window (sc_pos, 0, 4, 0.0);
2902 comp = 0;
2903
2904 for (row = 0; row < nbrow_fddl; row ++){
2905
2906 if ((time[row] > cpl_table_get_int (spectrum_sc, "TIME", 0, &nv) - exptime_sc/2) &&
2907 (time[row] < cpl_table_get_int (spectrum_sc, "TIME", nbrow_sc-1, &nv) + exptime_sc/2)){
2908 cpl_array_add (ft_pos, fddl_ft[row]);
2909 cpl_array_add (sc_pos, fddl_sc[row]);
2910 comp ++;
2911 }
2912
2913 CPLCHECK_MSG ("Problem during the compute of the mean of the fddl");
2914 }
2915
2916 if (comp != 0) {
2917 cpl_array_divide_scalar (ft_pos, comp);
2918 cpl_array_divide_scalar (sc_pos, comp);
2919 }
2920
2921 /*
2922 * Save the mean FT_POS and SC_POS
2923 */
2924 cpl_table_set_array (fddl_met_mean, "FT_POS", 0, ft_pos);
2925 cpl_table_set_array (fddl_met_mean, "SC_POS", 0, sc_pos);
2926 cpl_array_delete (ft_pos);
2927 cpl_array_delete (sc_pos);
2928
2929 if (thread == 0){
2930 cpl_propertylist * plist_img = cpl_propertylist_duplicate (
2931 gravi_data_get_plist (p2vmred_data, GRAVI_FDDL_EXT));
2932 cpl_propertylist_set_string (plist_img, "EXTNAME", "FDDL_MET_MEAN");
2933 gravi_data_add (oi_vis, plist_img, fddl_met_mean);
2934 cpl_propertylist_delete (plist_img);
2935 }
2936 else {
2937
2938 cpl_table_insert (gravi_data_get_table (oi_vis, "FDDL_MET_MEAN"), fddl_met_mean, thread);
2939 cpl_table_delete (fddl_met_mean);
2940 }
2941 gravi_msg_function_exit(1);
2942 return CPL_ERROR_NONE;
2943}
2944
2945
2946 /* TEST: Make this cropped profile flux conservative if complex */
2947 //cpl_msg_info ("TEST:","profile is forced flux conservative don't work for boxcard!!");
2948 //cpl_msg_info ("TEST:","nxc = %lld nyc = %lld", nxc, nyc);
2949 //int null_val;
2950 //for (cpl_size x = 0; x < nxc; x++ ) {
2951 // double sum_y = 0.0;
2952 // double sum_yy = 0.0;
2953 // for (cpl_size y = 0; y < nyc; y++ ) {
2954 // sum_y += cpl_image_get (profile_crop, x+1, y+1, &null_val);
2955 // sum_yy += pow (cpl_image_get (profile_crop, x+1, y+1, &null_val), 2.0);
2956 // }
2957 // for (cpl_size y = 0; y < nyc; y++ ) {
2958 // double current = cpl_image_get (profile_crop, x+1, y+1, &null_val);
2959 // cpl_image_set (profile_crop, x+1, y+1, current * sum_y / sum_yy);
2960 // }
2961 //}
2962
2963cpl_error_code gravi_p2vm_mean_spectrum (gravi_data * p2vm_map, gravi_data * preproc_data)
2964{
2965 gravi_msg_function_start(1);
2966 cpl_ensure_code (p2vm_map, CPL_ERROR_NULL_INPUT);
2967 cpl_ensure_code (preproc_data, CPL_ERROR_NULL_INPUT);
2968
2969 /* Extract the spectrum of all region for SC only */
2970 cpl_table * spectrum_table = gravi_data_get_spectrum_data (preproc_data, GRAVI_SC);
2971 cpl_ensure_code (spectrum_table, CPL_ERROR_ILLEGAL_INPUT);
2972
2973 /* Get sizes */
2974 cpl_size nrow = cpl_table_get_nrow (spectrum_table);
2975 cpl_size nreg = gravi_spectrum_get_nregion (spectrum_table);
2976
2977 /* Build output table */
2978 cpl_table * output_table = cpl_table_extract (spectrum_table, 0, 1);
2979
2980 /* Loop on region and rows to integrate the spectrums */
2981 for (cpl_size reg = 0; reg < nreg ; reg++) {
2982 cpl_array ** arrays = cpl_table_get_data_array (spectrum_table, GRAVI_DATA[reg]);
2983 cpl_array * array = cpl_table_get_data_array (output_table, GRAVI_DATA[reg])[0];
2984 cpl_ensure_code (arrays, CPL_ERROR_ILLEGAL_INPUT);
2985 for (cpl_size row = 1; row < nrow ; row++) cpl_array_add (array, arrays[row]);
2986 cpl_table_erase_column (output_table, GRAVI_DATAERR[reg]);
2987 }
2988
2989 /* Add this table to the P2VM */
2990 gravi_data_add_table (p2vm_map, NULL, GRAVI_SPECTRUM_DATA_SC_EXT, output_table);
2991
2992 gravi_msg_function_exit(1);
2993 return CPL_ERROR_NONE;
2994}
2995
2996
2997/*----------------------------------------------------------------------------*/
3002/*----------------------------------------------------------------------------*/
3003
3004gravi_data * gravi_compute_disp (gravi_data * vis_data)
3005{
3006 gravi_data * disp_data;
3007 int ntel = 4, nbase = 6, npol = 2;
3008 FILE * file;
3009 gsl_matrix * U, * V;
3010 gsl_vector * S, * work;
3011
3012 gravi_msg_function_start(1);
3013 cpl_ensure (vis_data, CPL_ERROR_NULL_INPUT, NULL);
3014
3015 /*
3016 * Init useful matrix
3017 */
3018
3019 /* M_matrix is to go from OPL to OPD */
3020 double M_tab[24]={ 1., -1., 0.0, 0.0,
3021 1., 0.0, -1., 0.0,
3022 1., 0.0, 0.0, -1.,
3023 0.0, 1., -1., 0.0,
3024 0.0, 1., 0.0, -1.,
3025 0.0, 0.0, 1., -1.};
3026 gsl_matrix * M_matrix=gsl_matrix_alloc(6,4);
3027 memcpy(gsl_matrix_ptr(M_matrix, 0, 0), M_tab, 24*sizeof(double));
3028
3029 /* M_matrix2 is to go from OPL of FDDL FT and SC to OPD */
3030 gsl_matrix * M_matrix2=gsl_matrix_alloc(6,8);
3031 gsl_vector * M_vector=gsl_vector_alloc(6);
3032 for (int tel=0; tel<ntel; tel++){
3033 gsl_matrix_get_col(M_vector, M_matrix, tel);
3034 gsl_matrix_set_col(M_matrix2, tel, M_vector);
3035 gsl_vector_scale(M_vector, -1.);
3036 gsl_matrix_set_col(M_matrix2, tel+ntel, M_vector);
3037 }
3038 gsl_vector_free(M_vector);
3039
3040 if (INFO_DEBUG) {
3041 file = fopen("M_matrix2.txt","w");
3042 my_gsl_matrix_fprintf (file, M_matrix2, "%g");
3043 fclose(file);
3044 }
3045
3046 /* Compute inverse of M_matrix2 */
3047
3048 U = gsl_matrix_alloc (8, 6);
3049 V = gsl_matrix_alloc (6, 6);
3050 S = gsl_vector_alloc (6);
3051 work = gsl_vector_alloc (6);
3052
3053 gsl_matrix_transpose_memcpy(U, M_matrix2);
3054 gsl_linalg_SV_decomp (U, V, S, work);
3055
3056 /* Get inverse of the M matrix */
3057
3058 double wv_at;
3059 gsl_matrix * M_matrix2_inv = gsl_matrix_alloc(ntel*2, nbase);
3060 for (int j = 0; j < ntel*2; j++) {
3061 for (int i = 0; i < nbase; i++){
3062 wv_at = 0;
3063 for (int ii = 0; ii < nbase; ii++){
3064 if( gsl_vector_get(S, ii) > 1e-14)
3065 wv_at += gsl_matrix_get(V, i, ii) / gsl_vector_get(S, ii) *
3066 gsl_matrix_get(U, j, ii);
3067 }
3068 gsl_matrix_set(M_matrix2_inv, j, i, wv_at);
3069 }
3070 }
3071
3072 gsl_matrix_free (V);
3073 gsl_matrix_free (U);
3074 gsl_vector_free (S);
3075 gsl_vector_free (work);
3076
3077 CPLCHECK_NUL("Cannot inverse matrix M");
3078
3079 if (INFO_DEBUG) {
3080 file = fopen("M_matrix2_inv.txt","w");
3081 my_gsl_matrix_fprintf (file, M_matrix2_inv, "%g");
3082 fclose(file);
3083 }
3084
3085 /* Check if there is 2 polarization */
3086 npol = 2;
3087 if (!gravi_data_get_oi_vis_plist (vis_data, GRAVI_SC, 0, npol) ||
3088 !gravi_data_get_oi_vis_plist (vis_data, GRAVI_SC, 1, npol) ||
3089 !gravi_data_get_oi_vis_plist (vis_data, GRAVI_FT, 0, npol) ||
3090 !gravi_data_get_oi_vis_plist (vis_data, GRAVI_FT, 1, npol) ) {
3091
3092 cpl_error_set_message (cpl_func, CPL_ERROR_ILLEGAL_INPUT,
3093 "Missing OI_VIS (need split pol for SC and FT)");
3094 return NULL;
3095 }
3096
3097
3098 /* (1)
3099 *
3100 * Compute the position from FDDL and MET
3101 */
3102
3103
3104 cpl_msg_info (cpl_func, "*** 1 ) Compute the real FDDL position from their positions and Metrology ***");
3105 cpl_msg_info (cpl_func, "Load the FDDL table");
3106
3107 cpl_table * oi_flux;
3108 oi_flux = gravi_data_get_oi_flux (vis_data, GRAVI_SC, 0, npol);
3109 cpl_size nrow = cpl_table_get_nrow (oi_flux) / 4;
3110 CPLCHECK_NUL ("Cannot get data");
3111
3112 /* fddl = [SC_POS,FT_POS]
3113 * fddl2 = [SC_POS^2,FT_POS^2] */
3114 gsl_matrix * fddl = gsl_matrix_alloc (nrow, ntel*2);
3115 gsl_matrix * fddl2 = gsl_matrix_alloc (nrow, ntel*2);
3116
3117 for (cpl_size row=0; row<nrow; row++) {
3118 for (int tel = 0; tel<ntel; tel++) {
3119 double value;
3120 value = cpl_table_get (oi_flux, "SC_POS", row*ntel + tel, NULL);
3121 gsl_matrix_set (fddl, row,tel,value);
3122 gsl_matrix_set (fddl2,row,tel,value*value);
3123
3124 value = cpl_table_get (oi_flux, "FT_POS", row*ntel + tel, NULL);
3125 gsl_matrix_set (fddl, row,ntel+tel,value);
3126 gsl_matrix_set (fddl2,row,ntel+tel,value*value);
3127 }
3128 }
3129
3130 /* Rescale fddl */
3131 gsl_matrix_scale (fddl, 1e-3);
3132 gsl_matrix_scale (fddl, 1e-6);
3133
3134 if (INFO_DEBUG) {
3135 file = fopen("FDDL_matrix.txt","w");
3136 my_gsl_matrix_fprintf (file, fddl, "%g");
3137 fclose(file);
3138 }
3139
3140 /* MET : allocate and fill the met_data matrix */
3141 cpl_msg_info (cpl_func, "Load the MET table");
3142
3143 gsl_matrix * met_data = gsl_matrix_alloc (nrow, 4);
3144
3145 for (cpl_size row=0; row<nrow; row++) {
3146 for (int tel = 0; tel<ntel; tel++) {
3147 double value = cpl_table_get (oi_flux, "OPD_MET_FC", row*ntel + tel, NULL);
3148 gsl_matrix_set (met_data,row,tel,value);
3149 }
3150 }
3151
3152 /* Convert met in [um] */
3153 gsl_matrix_scale (met_data, 1e6);
3154
3155 if (INFO_DEBUG) {
3156 file = fopen("MET_data.txt","w");
3157 my_gsl_matrix_fprintf (file, met_data, "%g");
3158 fclose(file);
3159 }
3160
3161 /* Compute MET = (met_data*M_matrix.T).T */
3162
3163 gsl_matrix * MET = gsl_matrix_alloc(nrow, nbase);
3164 gsl_blas_dgemm (CblasNoTrans, CblasTrans, 1., met_data, M_matrix, 0, MET);
3165
3166 if (INFO_DEBUG) {
3167 file = fopen("MET.txt","w");
3168 my_gsl_matrix_fprintf (file, MET, "%g");
3169 fclose(file);
3170 }
3171
3172 /* Model FDDL=MET+Ki*FT_POSi+Kj*SC_POSj+Kk*FT_POSk^2+Kl*SC_POSl^2 */
3173 cpl_msg_info(cpl_func, "Compute the FDDL linearity factors");
3174 int Nmodel = 8*2+6;
3175 gsl_matrix * model = gsl_matrix_calloc( nrow*nbase, Nmodel);
3176 for (cpl_size row=0; row<nrow; row++){
3177 for (int base=0; base<nbase; base++){
3178 gsl_matrix_set(model, row*nbase+base, base, 1);
3179 for (int tel=0; tel<ntel*2; tel++) {
3180 gsl_matrix_set(model, row*nbase+base, nbase+tel,
3181 gsl_matrix_get(fddl, row, tel) * gsl_matrix_get(M_matrix2, base, tel));
3182 gsl_matrix_set(model, row*nbase+base, nbase+8+tel,
3183 gsl_matrix_get(fddl2, row, tel) * gsl_matrix_get(M_matrix2, base, tel));
3184 }
3185 }
3186 }
3187
3188 if (INFO_DEBUG) {
3189 file = fopen("model.txt","w");
3190 my_gsl_matrix_fprintf (file, model, "%g");
3191 fclose(file);
3192 }
3193
3194 /* invert model by SV decomp and solve */
3195 U = gsl_matrix_alloc (nrow*nbase, Nmodel);
3196 V = gsl_matrix_alloc (Nmodel, Nmodel);
3197 S = gsl_vector_alloc (Nmodel);
3198 work = gsl_vector_alloc (Nmodel);
3199 gsl_matrix_memcpy(U, model);
3200
3201 gsl_linalg_SV_decomp (U, V, S, work);
3202
3203 gsl_vector * K_coeff = gsl_vector_alloc (Nmodel);
3204 gsl_vector * MET_vector = gsl_vector_alloc(nrow*nbase);
3205 memcpy (MET_vector->data, MET->data, nbase*nrow*sizeof(double));
3206 gsl_linalg_SV_solve (U, V, S, MET_vector, K_coeff);
3207
3208 gsl_matrix_free(U);
3209 gsl_matrix_free(V);
3210 gsl_vector_free(S);
3211 gsl_vector_free(work);
3212 gsl_vector_free(MET_vector);
3213
3214 cpl_msg_info(cpl_func, "Compute the real FDDL positions");
3215 gsl_matrix * FDDL_matrix=gsl_matrix_alloc(nrow, 8);
3216 for (cpl_size row=0; row<nrow; row++){
3217 for (int tel=0; tel<8; tel++){
3218 gsl_matrix_set(FDDL_matrix, row, tel,
3219 gsl_vector_get(K_coeff, tel+6)*gsl_matrix_get(fddl, row, tel) +
3220 gsl_vector_get(K_coeff, tel+6+8)*gsl_matrix_get(fddl2, row, tel));
3221 }
3222 }
3223
3224 /* Correction of met offset */
3225 for (cpl_size row=0; row<nrow; row++){
3226 for (int base=0; base<nbase; base++){
3227 gsl_matrix_set(MET, row, base, gsl_matrix_get(MET, row, base)
3228 - gsl_matrix_get(model, row*nbase+base, base)
3229 * gsl_vector_get(K_coeff, base));
3230 //MET-=dot(MODEL[:,:6,:6],K[:6]);
3231 }
3232 }
3233
3234 /* Correction of FDDL
3235 * FDDL+=dot(pinv2(M_matrix2),(MET.T-dot(M_matrix2,FDDL.T))).T
3236 * or FDDL+=dot(pinv2(M_matrix2),(MET-dot(FDDL,M_matrix2.T)).T).T
3237 * or FDDL+=dot((MET-dot(FDDL,M_matrix2.T)), pinv2(M_matrix2).T)
3238 * */
3239 gsl_matrix * FDDLdotM = gsl_matrix_alloc(nrow, nbase);
3240 gsl_matrix * METsub = gsl_matrix_alloc(nrow, nbase);
3241 gsl_matrix * FDDLcorrection = gsl_matrix_alloc(nrow, ntel*2);
3242
3243 memcpy (METsub->data, MET->data, nbase*nrow*sizeof(double));
3244 gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1., FDDL_matrix, M_matrix2, 0, FDDLdotM);
3245 gsl_matrix_sub(METsub, FDDLdotM);
3246
3247 gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1., METsub, M_matrix2_inv, 0, FDDLcorrection);
3248 gsl_matrix_add(FDDL_matrix, FDDLcorrection);
3249
3250 gsl_matrix_free(FDDLdotM);
3251 gsl_matrix_free(METsub);
3252 gsl_matrix_free(FDDLcorrection);
3253
3254
3255 if (INFO_DEBUG) {
3256 file = fopen("MET.txt","w");
3257 my_gsl_matrix_fprintf (file, MET, "%g");
3258 fclose(file);
3259 file = fopen("coeff.txt","w");
3260 for (int Imodel=0; Imodel<Nmodel;Imodel++) fprintf(file, "%e \n", gsl_vector_get(K_coeff, Imodel));
3261 fclose(file);
3262 file = fopen("FDDL_matrix.txt","w");
3263 my_gsl_matrix_fprintf (file, FDDL_matrix, "%g");
3264 fclose(file);
3265 }
3266
3267
3268 /* (2)
3269 *
3270 * Compute the phase
3271 */
3272 /* Load wavelength */
3273 cpl_msg_info(cpl_func, "*** 2 ) Compute the Phases ***");
3274 cpl_msg_info(cpl_func, "Load the phase SC (VISDATA)");
3275
3276 // FIXME: consider same wavelength for both polar
3277 cpl_table * oi_wavelength = gravi_data_get_oi_wave (vis_data, GRAVI_SC, 0, npol);
3278
3279 cpl_size nwave = cpl_table_get_nrow(oi_wavelength);
3280 cpl_array * wavenumber = cpl_array_new(nwave, CPL_TYPE_DOUBLE);
3281
3282 /* Compute wavelength in fiber -- FIXME: use the QC parameters to get back */
3283 for (cpl_size wave=0; wave<nwave; wave++)
3284 cpl_array_set(wavenumber, wave, (2.*M_PI/(((double)cpl_table_get(oi_wavelength, "EFF_WAVE", wave, NULL)*1.e6-0.021*1.908)/(1.-0.021))));
3285 if (INFO_DEBUG) {
3286 file = fopen("wavenumber.txt","w");
3287 for (cpl_size i=0; i<nwave;i++) fprintf(file, "%5.10g \n", cpl_array_get(wavenumber, i, NULL));
3288 fclose(file);
3289 }
3290
3291 /* load SC phase */
3292 cpl_table * oi_vis;
3293 oi_vis = gravi_data_get_oi_vis (vis_data, GRAVI_SC, 0, npol);
3294 cpl_array ** visdata_p1_orig = cpl_table_get_data_array(oi_vis, "VISDATA");
3295 cpl_array ** visdata_p1 = cpl_malloc(cpl_table_get_nrow(oi_vis)*sizeof(cpl_array *));
3296
3297 oi_vis = gravi_data_get_oi_vis (vis_data, GRAVI_SC, 1, npol);
3298 cpl_array ** visdata_p2_orig = cpl_table_get_data_array(oi_vis, "VISDATA");
3299 cpl_array ** visdata_p2 = cpl_malloc(cpl_table_get_nrow(oi_vis)*sizeof(cpl_array *));
3300
3301 if (INFO_DEBUG) {
3302 file = fopen("vis_p1.txt","w");
3303 }
3304
3305 /* correct phase from metrology
3306 * vis_p1 = Vis_p1 * e i(MET*wavenumber)*/
3307 cpl_msg_info(cpl_func, "Add phase MET to phase SC");
3308 for (cpl_size row=0; row<nrow; row++){
3309 for (int base=0; base < nbase; base++){
3310 visdata_p1[row*nbase+base]=cpl_array_duplicate(visdata_p1_orig[row*nbase+base]);
3311 gravi_array_multiply_phasor(visdata_p1[row*nbase+base], I*(gsl_matrix_get(MET, row, base)), wavenumber);
3312 if (INFO_DEBUG) {
3313 for (cpl_size i=0; i<nwave;i++) fprintf(file, "%f + i%f \t", creal(cpl_array_get_complex(visdata_p1[row*nbase+base], i, NULL)),
3314 cimag(cpl_array_get_complex(visdata_p1[row*nbase+base], i, NULL)));
3315 }
3316
3317 visdata_p2[row*nbase+base]=cpl_array_duplicate(visdata_p2_orig[row*nbase+base]);
3318 gravi_array_multiply_phasor(visdata_p2[row*nbase+base], I*gsl_matrix_get(MET, row, base), wavenumber);
3319
3320 if (INFO_DEBUG) {
3321 fprintf(file, "\n");
3322 }
3323 }
3324 }
3325 if (INFO_DEBUG) {
3326 fclose(file);
3327 }
3328
3329 /* Compute Group Delay */
3330 cpl_msg_info(cpl_func, "Compute Group Delay");
3331 cpl_array * GD1 = cpl_array_new(nrow*nbase, CPL_TYPE_DOUBLE);
3332 cpl_array * GD2 = cpl_array_new(nrow*nbase, CPL_TYPE_DOUBLE);
3333 cpl_array * GD1_stat = cpl_array_new(nrow, CPL_TYPE_DOUBLE);
3334 cpl_array * GD2_stat = cpl_array_new(nrow, CPL_TYPE_DOUBLE);
3335 cpl_array * GD1_med = cpl_array_new(nbase, CPL_TYPE_DOUBLE);
3336 cpl_array * GD2_med = cpl_array_new(nbase, CPL_TYPE_DOUBLE);
3337 double std_gd1 = 0.;
3338 double std_gd2 = 0.;
3339 for (int base=0; base < nbase; base++){
3340 for (cpl_size row=0; row<nrow; row++){
3341 double complex tmp = 0.0 * I + 0.0;
3342 for (cpl_size wave=1; wave<nwave; wave++)
3343 tmp += cpl_array_get_complex (visdata_p1[row*nbase+base], wave, NULL) * conj(cpl_array_get_complex (visdata_p1[row*nbase+base], wave-1, NULL));
3344 cpl_array_set(GD1 , row*nbase + base, carg(tmp));
3345 cpl_array_set(GD1_stat , row, carg(tmp));
3346 tmp = 0.0 * I + 0.0;
3347 for (cpl_size wave=1; wave<nwave; wave++)
3348 tmp += cpl_array_get_complex (visdata_p2[row*nbase+base], wave, NULL) * conj(cpl_array_get_complex (visdata_p2[row*nbase+base], wave-1, NULL));
3349 cpl_array_set(GD2 , row*nbase + base, carg(tmp));
3350 printf("%g ",carg(tmp));
3351 cpl_array_set(GD2_stat , row, carg(tmp));
3352 }
3353 // FIXME: make sure this is /nrow and not /row
3354 std_gd1+=cpl_array_get_stdev(GD1_stat)/nrow;
3355 std_gd2+=cpl_array_get_stdev(GD2_stat)/nrow;
3356 cpl_array_set(GD1_med, base, cpl_array_get_median(GD1_stat));
3357 cpl_array_set(GD2_med, base, cpl_array_get_median(GD2_stat));
3358 }
3359 cpl_array_delete(GD1_stat);
3360 cpl_array_delete(GD2_stat);
3361
3362 cpl_msg_info(cpl_func, "GD rms average for pola 1 : %g[radians/element]", std_gd1);
3363 cpl_msg_info(cpl_func, "GD rms average for pola 2 : %g[radians/element]", std_gd2);
3364
3365 if (INFO_DEBUG) {
3366 file = fopen("GD1.txt","w");
3367 for (cpl_size row=0; row<nrow; row++){
3368 for (int base=0; base < nbase; base++){
3369 fprintf(file, "%g \t", cpl_array_get(GD1 , row*nbase + base, NULL));
3370 }
3371 fprintf(file, "\n");
3372 }
3373 fclose(file);
3374 file = fopen("GD2.txt","w");
3375 for (cpl_size row=0; row<nrow; row++){
3376 for (int base=0; base < nbase; base++){
3377 fprintf(file, "%g \t", cpl_array_get(GD1 , row*nbase + base, NULL));
3378 }
3379 fprintf(file, "\n");
3380 }
3381 fclose(file);
3382 }
3383
3384
3385 /* slope as function of metrology */
3386 cpl_msg_info(cpl_func, "Fit the metrology slopes");
3387 size_t * sort_met[6];
3388 for (int base=0; base<nbase; base++){
3389 sort_met[base]=cpl_malloc(nrow*sizeof(size_t));
3390 gsl_sort_index (sort_met[base], gsl_matrix_ptr(MET, 0,base), nbase, nrow);
3391 }
3392
3393 int nspace=1000;
3394 double step = 0.000001;
3395 cpl_array * amp_sum=cpl_array_new(nspace, CPL_TYPE_DOUBLE_COMPLEX);
3396 cpl_array_fill_window_complex(amp_sum, 0, nspace, 0.0 * I + 0.0);
3397 cpl_array * amp=cpl_array_new(nspace, CPL_TYPE_DOUBLE_COMPLEX);
3398 gsl_matrix * A1_met = gsl_matrix_alloc (nbase, nwave);
3399 gsl_matrix * A2_met = gsl_matrix_alloc (nbase, nwave);
3400 cpl_array * i2step=cpl_array_new(nspace, CPL_TYPE_DOUBLE);
3401 cpl_size max_pos;
3402 CPLCHECK_NUL("Error before fitting the metrology slopes");
3403 for (int i = 0; i < nspace; ++i) cpl_array_set(i2step, i, step*(double)i*2.);
3404 CPLCHECK_NUL("Error init istep");
3405
3406 /* POLA 1 */
3407 for (int base = 0; base < nbase; ++base) {
3408 for (cpl_size wave = 0; wave < nwave; ++wave) {
3409 /* for all parameter space compute amp=vis_data*e^i(param*MET) and sum on row*/
3410 amp_sum=cpl_array_new(nspace, CPL_TYPE_DOUBLE_COMPLEX);
3411 cpl_array_fill_window_complex(amp_sum, 0, nspace, 0.0 * I + 0.0);
3412 CPLCHECK_NUL("Error before loop row");
3413 for (cpl_size row = 0; row < nrow; ++row) {
3414 cpl_array_fill_window_complex(amp, 0, nspace, (_Complex double)cpl_array_get_float_complex(visdata_p1[sort_met[base][row]*nbase+base], wave, NULL));
3415 gravi_array_multiply_phasor(amp, I*gsl_matrix_get(MET, (sort_met[base])[row], base), i2step);
3416 cpl_array_add(amp_sum, amp);
3417 CPLCHECK_NUL("Error in loop row");
3418 }
3419 cpl_array_abs(amp_sum);
3420 /* get the best param */
3421 cpl_array_get_maxpos(amp_sum, &max_pos);
3422 gsl_matrix_set(A1_met, base, wave, max_pos*step);
3423 cpl_msg_info_overwritable(cpl_func, "Fit phase sc for polar 1/2 of base %d/6 (wave %lld/%lld)", base+1, wave+1, nwave);
3424 CPLCHECK_NUL("Error when fitting the metrology slopes");
3425 }
3426 }
3427
3428 /* POLA 2 */
3429 for (int base = 0; base < nbase; ++base) {
3430 for (cpl_size wave = 0; wave < nwave; ++wave) {
3431 /* for all parameter space compute amp=vis_data*e^i(param*MET) and sum on row*/
3432 amp_sum=cpl_array_new(nspace, CPL_TYPE_DOUBLE_COMPLEX);
3433 cpl_array_fill_window_complex(amp_sum, 0, nspace, 0.0 * I + 0.0);
3434 CPLCHECK_NUL("Error before loop row");
3435 for (cpl_size row = 0; row < nrow; ++row) {
3436 cpl_array_fill_window_complex(amp, 0, nspace, (_Complex double)cpl_array_get_float_complex(visdata_p2[sort_met[base][row]*nbase+base], wave, NULL));
3437 gravi_array_multiply_phasor(amp, I*gsl_matrix_get(MET, (sort_met[base])[row], base), i2step);
3438 cpl_array_add(amp_sum, amp);
3439 CPLCHECK_NUL("Error in loop row");
3440 }
3441 cpl_array_abs(amp_sum);
3442 /* get the best param */
3443 cpl_array_get_maxpos(amp_sum, &max_pos);
3444 gsl_matrix_set(A2_met, base, wave, max_pos*step);
3445 cpl_msg_info_overwritable(cpl_func, "Fit phase sc for polar 2/2 of base %d/6 (wave %lld/%lld)", base+1, wave+1, nwave);
3446 CPLCHECK_NUL("Error when fitting the metrology slopes");
3447 }
3448 }
3449
3450 cpl_array_delete(amp_sum);
3451 cpl_array_delete(amp);
3452 cpl_array_delete(i2step);
3453
3454
3455 if (INFO_DEBUG) {
3456 file = fopen("A1_met.txt","w");
3457 my_gsl_matrix_fprintf (file, A1_met, "%g");
3458 fclose(file);
3459 file = fopen("A2_met.txt","w");
3460 my_gsl_matrix_fprintf (file, A2_met, "%g");
3461 fclose(file);
3462 }
3463
3464 /* Reconstruction des phases unwrappés
3465 * Slopes=MET[:,:,None]*A1-arange(Nw)*median(GD1,axis=0)[:,None]
3466 * SC1b=SC1*exp(1j*Slopes)
3467 * K=unwrap(angle(SC1b.mean(axis=0)))
3468 *
3469 */
3470 cpl_msg_info(cpl_func, "Reconstruct unwrap phase ...");
3471 cpl_array * SC1b_unwrap;
3472 cpl_array * SC2b_unwrap;
3473 cpl_array * slope = cpl_array_new(nwave, CPL_TYPE_DOUBLE_COMPLEX);
3474 for (int base=0; base < nbase; base++){
3475 SC1b_unwrap = cpl_array_new (nwave, CPL_TYPE_DOUBLE_COMPLEX);
3476 cpl_array_fill_window_complex (SC1b_unwrap, 0, nwave, (double complex)(0.0 + 0.0 * I));
3477 SC2b_unwrap = cpl_array_new (nwave, CPL_TYPE_DOUBLE_COMPLEX);
3478 cpl_array_fill_window_complex (SC2b_unwrap, 0, nwave, (double complex)(0.0 + 0.0 * I));
3479 for (cpl_size row=0; row<nrow; row++){
3480 for (cpl_size wave=0; wave<nwave; wave++){
3481 cpl_array_set(slope, wave,cpl_array_get(visdata_p1[row*nbase+base], wave, NULL)*
3482 cexp(I*(
3483 gsl_matrix_get(MET, row, base)
3484 *gsl_matrix_get(A1_met, base, wave)
3485 -wave*cpl_array_get(GD1_med, base, NULL))
3486 ));
3487 }
3488 cpl_array_add(SC1b_unwrap, slope);
3489 for (cpl_size wave=0; wave<nwave; wave++){
3490 cpl_array_set(slope, wave,cpl_array_get(visdata_p2[row*nbase+base], wave, NULL)*
3491 cexp(I*(
3492 gsl_matrix_get(MET, row, base)
3493 *gsl_matrix_get(A2_met, base, wave)
3494 -wave*cpl_array_get(GD2_med, base, NULL))
3495 ));
3496 }
3497 cpl_array_add(SC2b_unwrap, slope);
3498 }
3499 cpl_array_divide_scalar(SC1b_unwrap, nrow);
3500 cpl_array_divide_scalar(SC2b_unwrap, nrow);
3501 // unwrap arg(SC1b)
3502 cpl_array_arg(SC1b_unwrap);
3503 cpl_array_arg(SC2b_unwrap);
3504 gravi_array_phase_unwrap (SC1b_unwrap);
3505 gravi_array_phase_unwrap (SC2b_unwrap);
3506 for (cpl_size row=0; row<nrow; row++){
3507 for (cpl_size wave=0; wave<nwave; wave++){
3508 cpl_array_set(slope, wave,cpl_array_get(visdata_p1[row*nbase+base], wave, NULL)*
3509 cexp(I*(
3510 gsl_matrix_get(MET, row, base)
3511 *gsl_matrix_get(A1_met, base, wave)
3512 -wave*cpl_array_get(GD1_med, base, NULL))
3513 ));
3514 }
3515
3516 }
3517
3518 }
3519
3520 cpl_array_delete(SC1b_unwrap);
3521 cpl_array_delete(SC2b_unwrap);
3522
3523 disp_data = gravi_data_new (0);
3524 cpl_propertylist * disp_header = gravi_data_get_header (disp_data);
3525
3526 /*
3527 * Save QC Parameters
3528 */
3529 /* Coeff linearity fddl K */
3530 char * qc_name;
3531 for (int tel = 0; tel<4; tel++ ) {
3532 cpl_msg_info (cpl_func, "FDDL linearity K FT%i = %f, K SC%i = %f ", tel+1, gsl_vector_get(K_coeff, tel+6),
3533 tel+1, gsl_vector_get(K_coeff, tel+4+6));
3534 qc_name=cpl_sprintf("ESO QC FDDL_LIN K_FT%i", tel+1);
3535 cpl_propertylist_append_double (disp_header, qc_name,gsl_vector_get(K_coeff, tel+6));
3536 cpl_propertylist_set_comment (disp_header, qc_name, "[-] K fddl linearity factor");
3537 qc_name=cpl_sprintf("ESO QC FDDL_LIN K_SC%i", tel+1);
3538 cpl_propertylist_append_double (disp_header, qc_name, gsl_vector_get(K_coeff, tel+4+6));
3539 cpl_propertylist_set_comment (disp_header, qc_name, "[-] K fddl linearity factor");
3540 }
3541 /* Coeff linearity fddl K2 */
3542 for (int tel = 0; tel<4; tel++ ) {
3543 cpl_msg_info (cpl_func, "FDDL linearity K2 FT%i = %f, K2 SC%i = %f ", tel,
3544 gsl_vector_get (K_coeff, tel+6+8),
3545 tel, gsl_vector_get(K_coeff, tel+4+6));
3546 qc_name=cpl_sprintf("ESO QC FDDL_LIN K2_FT%i", tel+1);
3547 cpl_propertylist_append_double (disp_header, qc_name,gsl_vector_get(K_coeff, tel+6+8));
3548 cpl_propertylist_set_comment (disp_header, qc_name, "[-] K2 fddl linearity factor");
3549 qc_name=cpl_sprintf("ESO QC FDDL_LIN K2_SC%i", tel+1);
3550 cpl_propertylist_append_double (disp_header, qc_name, gsl_vector_get(K_coeff, tel+4+6+8));
3551 cpl_propertylist_set_comment (disp_header, qc_name, "[-] K2 fddl linearity factor");
3552 }
3553 cpl_free(qc_name);
3554
3555 /*
3556 * Free memory
3557 */
3558 gsl_matrix_free(fddl);
3559 gsl_matrix_free(fddl2);
3560 gsl_matrix_free(met_data);
3561 gsl_matrix_free(M_matrix);
3562 gsl_matrix_free(M_matrix2);
3563 gsl_matrix_free(MET);
3564 gsl_matrix_free(model);
3565 gsl_vector_free(K_coeff);
3566 gsl_matrix_free(FDDL_matrix);
3567 cpl_array_delete(wavenumber);
3568// gsl_matrix_complex_free(vis_p1);
3569// gsl_matrix_complex_free(vis_p2);
3570 gsl_matrix_free(A1_met);
3571 gsl_matrix_free(A2_met);
3572 cpl_array_delete(GD1_med);
3573 cpl_array_delete(GD2_med);
3574
3575 gravi_msg_function_exit(1);
3576 return disp_data;
3577}
3578
3579int my_gsl_matrix_fprintf(FILE *stream, gsl_matrix *m, const char *fmt)
3580{
3581 size_t rows=m->size1;
3582 size_t cols=m->size2;
3583 size_t row,col,ml;
3584 int fill;
3585 char buf[100];
3586 gsl_vector *maxlen;
3587
3588 maxlen=gsl_vector_alloc(cols);
3589 for (col=0;col<cols;++col) {
3590 ml=0;
3591 for (row=0;row<rows;++row) {
3592 sprintf(buf,fmt,gsl_matrix_get(m,row,col));
3593 if (strlen(buf)>ml)
3594 ml=strlen(buf);
3595 }
3596 gsl_vector_set(maxlen,col,ml);
3597 }
3598
3599 for (row=0;row<rows;++row) {
3600 for (col=0;col<cols;++col) {
3601 sprintf(buf,fmt,gsl_matrix_get(m,row,col));
3602 fprintf(stream,"%s",buf);
3603 fill=gsl_vector_get(maxlen,col)+1-strlen(buf);
3604 while (--fill>=0)
3605 fprintf(stream," ");
3606 }
3607 fprintf(stream,"\n");
3608 }
3609 gsl_vector_free(maxlen);
3610 return 0;
3611}
3612
3613
3614
3615/*----------------------------------------------------------------------------*/
3627/*----------------------------------------------------------------------------*/
3628
3629cpl_table * gravi_opdmet_calibration (cpl_table * metrology_table,
3630 double mjd_obs)
3631{
3632 gravi_msg_function_start(1);
3633 cpl_ensure (metrology_table, CPL_ERROR_NULL_INPUT, NULL);
3634 cpl_vector * vectA, * vectB, *opl_vect;
3635 const cpl_array * volt;
3636 char name[90];
3637 char name2[90];
3638 int nb_tel = 4 ;
3639 int n_diode = 80;
3640
3641 /* Get data */
3642 cpl_size nb_row = cpl_table_get_nrow (metrology_table);
3643 CPLCHECK_NUL ("Cannot get data");
3644
3645 /* Create output table */
3646 cpl_table * opdmet_table = cpl_table_new (nb_row);
3647 cpl_table_new_column (opdmet_table, "TIME", CPL_TYPE_INT);
3648 for (cpl_size row = 0; row < nb_row; row++){
3649 cpl_table_set_int (opdmet_table, "TIME", row,
3650 cpl_table_get_int (metrology_table, "TIME",
3651 row, NULL) + mjd_obs);
3652 }
3653
3654 /* Loop on SC and FT */
3655 for (int gravi_type = 0; gravi_type < 2; gravi_type++ ){
3656 int comb = (gravi_type == GRAVI_SC ? 1 : 2);
3657
3658 /* Loop on beams */
3659 for (int tel = 0; tel < nb_tel; tel++){
3660
3661 /* load vectA and vectB from metrology */
3662 vectA = cpl_vector_new (nb_row);
3663 vectB = cpl_vector_new (nb_row);
3664
3665 for (cpl_size row = 0; row < nb_row; row ++){
3666 volt = cpl_table_get_array ( metrology_table, "VOLT", row);
3667 cpl_vector_set (vectA, row, cpl_array_get (volt , n_diode - 2*(comb*nb_tel - tel), NULL));
3668 cpl_vector_set (vectB, row, cpl_array_get (volt , n_diode - 2*(comb*nb_tel - tel) + 1, NULL));
3669 }
3670
3671 /* Compute phase from vectA and vectB */
3672 opl_vect = gravi_vectors_phase_create (vectA, vectB);
3673 cpl_vector_multiply_scalar (opl_vect, LAMBDA_MET/(2.0*M_PI));
3674
3675 FREE (cpl_vector_delete, vectA);
3676 FREE (cpl_vector_delete, vectB);
3677 CPLCHECK_NUL ("Compute OPD");
3678
3679 /* Compute OPD, OPL_FT and OPL_SC from the phase */
3680 if (gravi_type == GRAVI_SC) {
3681 sprintf (name,"OPL_FC_SC%d", tel + 1);
3682 cpl_table_wrap_double (opdmet_table, cpl_vector_get_data(opl_vect), name);
3683 cpl_table_set_column_unit (opdmet_table, name, "m");
3684 CPLCHECK_NUL("Wrap OPL_FC_SC");
3685
3686 /* duplicate the OPL_FC_SC into the OPD_FC */
3687 sprintf (name, "OPD_FC%d", tel + 1);
3688 sprintf (name2, "OPL_FC_SC%d", tel + 1);
3689 cpl_table_duplicate_column (opdmet_table, name, opdmet_table, name2);
3690 cpl_table_set_column_unit (opdmet_table, name, "m");
3691 CPLCHECK_NUL("Wrap OPD_FC");
3692 }
3693
3694 if (gravi_type == GRAVI_FT) {
3695 sprintf (name,"OPL_FC_FT%d", tel + 1);
3696 cpl_table_wrap_double (opdmet_table, cpl_vector_get_data(opl_vect), name);
3697 cpl_table_set_column_unit (opdmet_table, name, "m");
3698
3699 /* compute the OPD as OPL_SC - OPL_FT */
3700 sprintf (name2, "OPD_FC%d", tel + 1);
3701 cpl_table_subtract_columns (opdmet_table, name2, name);
3702 cpl_table_set_column_unit (opdmet_table, name, "m");
3703 CPLCHECK_NUL("Wrap OPL_FC_FT");
3704 }
3705 cpl_vector_unwrap (opl_vect);
3706
3707 } /* End loop beam */
3708 } /* End loop SC, FT */
3709
3710 /* Verbose */
3711 gravi_msg_function_exit(1);
3712 return opdmet_table;
3713}
3714
3715
3716cpl_error_code gravi_phase_correct_closures_new (cpl_table * phase_table)
3717{
3718 gravi_msg_function_start(1);
3719 cpl_ensure_code (phase_table, CPL_ERROR_NULL_INPUT);
3720
3721 cpl_table * input_table = cpl_table_duplicate (phase_table);
3722 cpl_array ** ref = cpl_table_get_data_array (input_table, "PHASE");
3723 cpl_array ** out = cpl_table_get_data_array (phase_table, "PHASE");
3724
3725 for (int base = 0; base < nbase; base++) {
3726 for (int clo = 0; clo < 2; clo ++) {
3727 for (int b = 0; b < 2; b ++) {
3728 cpl_array * tmp = ref[ GRAVI_TRI_BASE[base][clo][b] ];
3729 double sign = GRAVI_TRI_SIGN[base][clo][b];
3730 cpl_array_multiply_scalar (tmp, sign);
3731 cpl_array_add (out[base], tmp);
3732 cpl_array_multiply_scalar (tmp, sign);
3733 CPLCHECK_MSG ("Cannot correct CP");
3734 }
3735 }
3736 cpl_array_divide_scalar (out[base], 3.0);
3737 }
3738 FREE (cpl_table_delete, input_table);
3739
3740 gravi_msg_function_exit(1);
3741 return CPL_ERROR_NONE;
3742}
3743
3744
3745/*----------------------------------------------------------------------------*/
3749/*----------------------------------------------------------------------------*/
3750
3751cpl_table * gravi_metrology_reform (cpl_table * vis_met)
3752{
3753 gravi_msg_function_start(1);
3754 cpl_ensure (vis_met, CPL_ERROR_NULL_INPUT, NULL);
3755 char name[90];
3756 int ntel = 4;
3757
3758 /* Get data */
3759 cpl_size nrow = cpl_table_get_nrow (vis_met)/ntel;
3760 CPLCHECK_NUL ("Cannot get data");
3761
3762 /* Create output table */
3763 /* Create the output table for OPD_MET */
3764 cpl_table * opd_met = cpl_table_new (nrow);
3765 cpl_table_new_column (opd_met, "TIME", CPL_TYPE_INT);
3766 cpl_table_set_column_unit (opd_met, "TIME", "us");
3767
3768 /* Create and get the pointer of the table columns */
3769 double **p_opd_met;
3770 p_opd_met = cpl_malloc(ntel*sizeof(double*));
3771 for (int tel=0; tel<ntel; tel++) {
3772 sprintf (name, "OPD_FC%d", tel + 1);
3773 cpl_table_new_column (opd_met, name, CPL_TYPE_DOUBLE);
3774 cpl_table_set_column_unit (opd_met, name, "m");
3775 p_opd_met[tel]=cpl_table_get_data_double(opd_met, name);
3776 }
3777 CPLCHECK_NUL ("Cannot create data");
3778
3779 /*
3780 * Loop on the row to reform the table
3781 */
3782 double * p_vis_met=cpl_table_get_data_double(vis_met, "PHASE_FC");
3783 for (cpl_size row = 0; row < nrow; row++){
3784
3785 /* Copy the TIME column */
3786 cpl_table_set_int (opd_met, "TIME", row,
3787 cpl_table_get_int (vis_met, "TIME", row*ntel, NULL));
3788
3789 /* Reform the VIS_MET column */
3790 for (int tel=0; tel<ntel; tel++) {
3791 p_opd_met[tel][row] = p_vis_met[row*ntel+tel]*LAMBDA_MET/(2.0*M_PI);
3792 CPLCHECK_NUL ("Cannot set data");
3793 }
3794 }
3795 CPLCHECK_NUL ("Cannot reform data");
3796
3797 cpl_free(p_opd_met);
3798
3799 /* Verbose */
3800 gravi_msg_function_exit(1);
3801 return opd_met;
3802}
3803
3804
3805 if (cpl_table_get (phase_ft, "PHASE_SC", row*nbase+base, NULL) !=0
3806 && (row % (nrow/5) == 0) ) {
3807 double phase_sc = cpl_table_get (phase_ft, "PHASE_SC", row*nbase+base, NULL);
3808 double phasec_sc = gsl_matrix_get (A_data, base * nrow + row, 7);
3809 cpl_msg_info ("TEST", "FT base %i row %lld : phi_sc = %g phase_sc = %g diff = %g [rad]",
3810 base, row, phasec_sc, phase_sc, phasec_sc - phase_sc);
3811 double phase_met = cpl_table_get (phase_ft, "PHASE_MET_FC", row*nbase+base, NULL);
3812 double phasec_met = gsl_vector_get (bis_data, base * nrow + row) / LAMBDA_MET * CPL_MATH_2PI;
3813 cpl_msg_info ("TEST", "FT base %i row %lld : phi_met = %g phase_met = %g diff = %g [rad]",
3814 base, row, phasec_met, phase_met, phasec_met - phase_met);
3815 }
3816
3817
3818
3819 /* Resample MET and FT phases at the SC period */
3820 gravi_vis_create_met_ft (phase_ft, vis_met);
3821 gravi_vis_create_phasesc_ft (phase_ft, phase_sc, dit_sc);
3822 cpl_table_save (vis_met, NULL, NULL, "vismet.fits", CPL_IO_CREATE);
3823 cpl_table_save (phase_sc, NULL, NULL, "phase_sc.fits", CPL_IO_CREATE);
3824 cpl_table_save (phase_ft, NULL, NULL, "phase_ft.fits", CPL_IO_CREATE);
3825
3826
3827
3828cpl_vector * gravi_phase_fit_opdmet (cpl_table * vis_met,
3829 cpl_table * phase_sc,
3830 cpl_table * phase_ft,
3831 double dit_sc)
3832{
3833 int nbase = 6, ntel = 4;
3834 double opl, time_sc, time_ft, time_met;
3835 int nv, comp;
3836
3837 gravi_msg_function_start(1);
3838 cpl_ensure (vis_met, CPL_ERROR_NULL_INPUT, NULL);
3839 cpl_ensure (phase_sc, CPL_ERROR_NULL_INPUT, NULL);
3840 cpl_ensure (phase_ft, CPL_ERROR_NULL_INPUT, NULL);
3841
3842 /* Get the number of acquisitions */
3843 cpl_size nrow = cpl_table_get_nrow (phase_ft) / nbase;
3844 cpl_size nrow_sc = cpl_table_get_nrow (phase_sc) / nbase;
3845 cpl_size nrow_met = cpl_table_get_nrow (vis_met) / ntel;
3846 CPLCHECK_NUL("Cannot get data");
3847
3848 /* Get time */
3849 int * time_FT = cpl_table_get_data_int (phase_ft, "TIME");
3850 int * time_SC = cpl_table_get_data_int (phase_sc, "TIME");
3851 CPLCHECK_NUL ("Cannot get time");
3852
3853 /* Compute the periods of the signals */
3854 double exptime_ft = time_FT[nbase] - time_FT[0];
3855 double exptime_sc = time_SC[nbase] - time_SC[0];
3856
3857 /* Get pointer to data */
3858 double * ft_phase = cpl_table_get_data_double (phase_ft, "PHASE");
3859 double * sc_phase = cpl_table_get_data_double (phase_sc, "PHASE");
3860 double * phase_met = cpl_table_get_data_double (vis_met, "PHASE_FC");
3861 CPLCHECK_NUL ("Cannot get data");
3862
3863 /* Allocate the matrices */
3864 gsl_matrix * A_data = gsl_matrix_calloc (nbase * nrow, 8);
3865 gsl_matrix * U, * A;
3866 gsl_vector * bis_data = gsl_vector_alloc (nbase * nrow), * bis;
3867 gsl_matrix * X = gsl_matrix_alloc (8, 8);
3868 gsl_matrix * V = gsl_matrix_alloc (8, 8);
3869 gsl_vector * S = gsl_vector_alloc (8);
3870 gsl_vector * work = gsl_vector_alloc (8);
3871 gsl_vector * x = gsl_vector_alloc (8);
3872 CPLCHECK_NUL ("Allocate matrix");
3873
3874 /*
3875 * Extract the OPl metrology and compute the mean of each exposure time
3876 */
3877 int ind_sc;
3878 for (int base = 0; base < nbase; base ++) {
3879 int tel0 = GRAVI_BASE_TEL[base][0];
3880 int tel1 = GRAVI_BASE_TEL[base][1];
3881
3882 ind_sc = 0;
3883 int im = 0;
3884
3885 /*
3886 * resample the SC and Metrology at FT period
3887 */
3888 for (cpl_size row = 0; row < nrow; row ++){ // loop on row FT
3889
3890 gsl_matrix_set (A_data, base * nrow + row, 6, - ft_phase[row*nbase+base]);
3891 CPLCHECK_NUL ("Set FT phase in matrix");
3892
3893 time_ft = time_FT[row*nbase+base];
3894 if (ind_sc < nrow_sc) time_sc = time_SC[ind_sc*nbase+base];
3895 CPLCHECK_NUL ("Get time sc");
3896
3897 /* increase ind ft if time_ft in the range of the nex DIT */
3898 while(time_ft > (time_sc+exptime_sc/2.)){
3899 ind_sc++;
3900 if (ind_sc < nrow_sc) time_sc = time_SC[ind_sc*nbase+base];
3901 else break;
3902 }
3903 CPLCHECK_NUL ("Get next sc dit");
3904
3905 /* get phi sc */
3906 double phi_sc = 0;
3907 if (ind_sc > 0){
3908 if(ind_sc < nrow_sc) phi_sc = sc_phase[ind_sc*nbase+base];
3909 else phi_sc = sc_phase[(nrow_sc - 1)*nbase+base];
3910 }
3911 else phi_sc = sc_phase[0*nbase+base];
3912 CPLCHECK_NUL ("Get phi_sc");
3913
3914
3915 gsl_matrix_set (A_data, base * nrow + row, 7, phi_sc);
3916 gsl_matrix_set (A_data, base * nrow + row, base, 1.0);
3917
3918 /* Metrology case */
3919 opl = 0;
3920 comp = 0;
3921
3922 if (im < nrow_met)
3923 time_met = cpl_table_get_int (vis_met, "TIME", im*ntel, &nv);
3924
3925 while ((time_met < (time_ft + exptime_ft))){
3926
3927 if (im < nrow_met) {
3928 opl += (phase_met[im*ntel+tel0] - phase_met[im*ntel+tel1]) * LAMBDA_MET / CPL_MATH_2PI;
3929 comp ++;
3930 }
3931
3932 im++;
3933 if (im < nrow_met) {
3934 time_met = cpl_table_get_int (vis_met, "TIME", im*ntel, &nv);
3935 }
3936 else {
3937 break;
3938 }
3939
3940 CPLCHECK_NUL("Get time metrology");
3941 }
3942
3943 /* average the metrology over the FT DIT */
3944 if (comp != 0)
3945 gsl_vector_set (bis_data, base * nrow + row, opl/comp);
3946
3947 /* if no metrology signal within the FT DIT interpolate the metrology */
3948 else {
3949 if (im > 0){
3950 if(im < nrow_met) {
3951 double opd1 = (phase_met[(im-1)*ntel+tel0] - phase_met[(im-1)*ntel+tel1]) * LAMBDA_MET / CPL_MATH_2PI;
3952 double opd2 = (phase_met[im*ntel+tel0] - phase_met[im*ntel+tel1]) * LAMBDA_MET / CPL_MATH_2PI;
3953 double opd = (opd2 - opd1) *
3954 (time_ft - cpl_table_get_int (vis_met, "TIME", (im-1)*ntel, &nv)) /
3955 (time_met - cpl_table_get_int (vis_met, "TIME", (im-1)*ntel, &nv)) +
3956 opd1;
3957 gsl_vector_set (bis_data, base * nrow + row, opd);
3958 }
3959 else {
3960 double opd = (phase_met[(nrow_met-1)*ntel+tel0] - phase_met[(nrow_met-1)*ntel+tel1]) * LAMBDA_MET / CPL_MATH_2PI;
3961 gsl_vector_set (bis_data, base * nrow + row, opd);
3962 }
3963 }
3964 else {
3965 double opd = (phase_met[0*ntel+tel0] - phase_met[0*ntel+tel1]) * LAMBDA_MET / CPL_MATH_2PI;
3966 gsl_vector_set (bis_data, base * nrow + row, opd);
3967 }
3968 CPLCHECK_NUL("Interpolate the metrology");
3969 }
3970
3971 } /* end loop on row FT */
3972 } /* End loop on base */
3973
3974 /*
3975 * filter the data out of the SC integration time
3976 */
3977
3978 long n_row_A=nrow_sc*(int)(dit_sc/exptime_ft+1);
3979 cpl_vector *i_A_vector = cpl_vector_new(n_row_A);
3980 int i_A=0;
3981 int time_mod;
3982 int t0_sc = time_SC[0*nbase];
3983 int tend_sc = time_SC[(nrow_sc-1)*nbase];
3984
3985 /* Find the index of the frames within the SC integration
3986 *Assume all baseline have the same timing */
3987 for (cpl_size row=0; row<nrow; row++){
3988 time_ft = time_FT[row*nbase];
3989 if ((time_ft >= t0_sc-dit_sc/2) && (time_ft < tend_sc+dit_sc/2)){
3990 time_mod=((int)(time_ft-(t0_sc-dit_sc/2)) % (int)exptime_sc);
3991 if ( time_mod >= 0 && time_mod < dit_sc ) {
3992 cpl_vector_set (i_A_vector, i_A, row);
3993 i_A++;
3994 }
3995 }
3996 }
3997 n_row_A=i_A;
3998
3999 /* copy the frames within SC integration in the matrix A */
4000 A = gsl_matrix_alloc (nbase * n_row_A, 8);
4001 bis = gsl_vector_alloc (nbase * n_row_A);
4002 cpl_vector *time_A = cpl_vector_new(n_row_A);
4003
4004 for (int base = 0; base < nbase; base ++)
4005 for (cpl_size i_A = 0; i_A < n_row_A; i_A++){
4006 cpl_size row = cpl_vector_get(i_A_vector, i_A);
4007 for (int col = 0; col < 8; col ++)
4008 gsl_matrix_set (A, base * n_row_A + i_A, col, gsl_matrix_get (A_data, base * nrow + row, col));
4009
4010 gsl_vector_set (bis, base * n_row_A + i_A, gsl_vector_get (bis_data, base * nrow + row));
4011 if (base == 0)
4012 cpl_vector_set (time_A, i_A, time_FT[row*nbase+base]);
4013 }
4014
4015 cpl_vector_delete(i_A_vector);
4016
4017 /*
4018 * Solve the linear equation Ax = b to get the coefficients to deduce the OPDs
4019 */
4020 nrow=n_row_A;
4021 U = gsl_matrix_alloc (nbase * nrow, 8);
4022 gsl_matrix_memcpy (U, A);
4023
4024
4025 gsl_linalg_SV_decomp (U, V, S, work);
4026 gsl_linalg_SV_solve (U, V, S, bis, x);
4027 cpl_vector * opd_coeff = cpl_vector_new (3);
4028 cpl_vector_set (opd_coeff, 0, gsl_vector_get (x, 7));
4029 cpl_vector_set (opd_coeff, 1, gsl_vector_get (x, 6));
4030
4031 cpl_msg_debug (cpl_func, "Wavelength => SC = %g, FT = %g\n",
4032 cpl_vector_get(opd_coeff,0)*2*M_PI, cpl_vector_get(opd_coeff,1)*2*M_PI);
4033
4034
4035 /*
4036 * Save the A and B matrix for debug
4037 */
4038 if (INFO_DEBUG) {
4039 cpl_table * table_output = cpl_table_new (nrow * nbase);
4040 const char * table_name[9] = {"A0","A1","A2","A3","A4","A5","A6","A7","B"};
4041
4042 /* Init memory */
4043 double ** table_value = cpl_malloc (9 * sizeof(double*));
4044 for (int col = 0; col < 9; col++)
4045 table_value[col] = cpl_malloc (nrow * nbase * sizeof(double));
4046
4047 /* Fill values */
4048 cpl_msg_info (cpl_func, "fill tmp table");
4049 for (int base = 0; base < nbase; base ++)
4050 for (cpl_size row = 0; row < nrow; row++) {
4051 for (int col = 0; col < 8; col ++)
4052 table_value[col][base * nrow + row] = gsl_matrix_get (A, base * nrow + row, col);
4053
4054 table_value[8][base * nrow + row] = gsl_vector_get (bis, base * nrow + row);
4055 }
4056
4057 /* Wrap into table */
4058 cpl_msg_info (cpl_func, "wrap tmp table");
4059 for (int col = 0; col < 9; col++) {
4060 cpl_table_wrap_double (table_output, table_value[col], table_name[col]);
4061 }
4062
4063 /* Add best fit in header */
4064 cpl_propertylist * table_header = cpl_propertylist_new ();
4065 for (int col = 0; col < 8; col ++)
4066 cpl_propertylist_update_double (table_header, table_name[col], gsl_vector_get (x, col));
4067
4068 /* Save tmp table */
4069 cpl_msg_info (cpl_func, "save tmp table");
4070 cpl_table_save (table_output, NULL, table_header, "matrix_AB.fits", CPL_IO_CREATE);
4071
4072 cpl_msg_info (cpl_func, "done");
4073 }
4074 /* --- */
4075
4076
4077
4078 cpl_vector *vect_phase_sc_ft=cpl_vector_new(nrow*nbase);
4079 cpl_vector *vect_dopd_met=cpl_vector_new(nrow*nbase);
4080 cpl_vector *vect_diff=cpl_vector_new(nrow*nbase);
4081 for (cpl_size row = 0; row < nrow*nbase; row ++){
4082 cpl_vector_set(vect_dopd_met, row, gsl_vector_get (bis, row));
4083 cpl_vector_set(vect_phase_sc_ft, row,
4084 gsl_matrix_get (A, row, 0)*gsl_vector_get (x, 0)+
4085 gsl_matrix_get (A, row, 1)*gsl_vector_get (x, 1)+
4086 gsl_matrix_get (A, row, 2)*gsl_vector_get (x, 2)+
4087 gsl_matrix_get (A, row, 3)*gsl_vector_get (x, 3)+
4088 gsl_matrix_get (A, row, 4)*gsl_vector_get (x, 4)+
4089 gsl_matrix_get (A, row, 5)*gsl_vector_get (x, 5)+
4090 gsl_matrix_get (A, row, 6)*gsl_vector_get (x, 6)+
4091 gsl_matrix_get (A, row, 7)*gsl_vector_get (x, 7));
4092 cpl_vector_set(vect_diff, row, cpl_vector_get(vect_dopd_met, row)-cpl_vector_get(vect_phase_sc_ft, row));
4093 }
4094
4095 double rms_fit_phase_met = cpl_vector_get_stdev(vect_diff);
4096 cpl_msg_info (cpl_func, "RMS of residuals on fit of a.phi_ft+b.phi_sc+c = met : %g", rms_fit_phase_met);
4097 cpl_vector_set (opd_coeff, 2, rms_fit_phase_met);
4098
4099 cpl_msg_info ("TEST", "coeff SC = %.20g [um]", cpl_vector_get (opd_coeff, 0) * CPL_MATH_2PI * 1e6);
4100 cpl_msg_info ("TEST", "coeff FT = %.20g [um]", cpl_vector_get (opd_coeff, 1) * CPL_MATH_2PI * 1e6);
4101 cpl_msg_info ("TEST", "residual = %.20g [um]", rms_fit_phase_met * CPL_MATH_2PI * 1e6);
4102
4103 /*
4104 * If the residuals are too high print a warning : SC and RMN data may be desynchronized
4105 * */
4106 if (rms_fit_phase_met > 1e-7 ){
4107 cpl_msg_info (cpl_func, "*************************************************");
4108 cpl_msg_warning (cpl_func, "**** !!! residuals of the fit too high !!! ****");
4109 cpl_msg_warning (cpl_func, "**** SC and RMN may be desynchronized ****");
4110 cpl_msg_warning (cpl_func, "**** (or out of the envelope in LOW) ****");
4111 cpl_msg_info (cpl_func, "*************************************************");
4112 }
4113
4114 if (PLOT_MET_PHASE_FIT) {
4115 cpl_msg_info(cpl_func, "Plot fit residuals");
4116
4117 cpl_errorstate prestate = cpl_errorstate_get();
4118 const cpl_vector ** vectors=malloc(3 * sizeof(cpl_vector*));
4119 vectors[0]=NULL;
4120 vectors[1]=vect_phase_sc_ft;
4121 vectors[2]=vect_dopd_met;
4122
4123 cpl_plot_vectors (cpl_sprintf("set title 'Met fit case a=%g b=%g'; set xlabel 'time[10-6s]'; set ylabel 'Phase_ft+Phase_sc, dopd_met';", gsl_vector_get (x, 0), gsl_vector_get (x, 1)),
4124 "", "", vectors, 3);
4125 cpl_plot_vector(cpl_sprintf("set title 'Met fit case a=%g b=%g'; set xlabel 'time[10-6s]'; set ylabel 'Phase_ft+Phase_sc-dopd_met';", gsl_vector_get (x, 0), gsl_vector_get (x, 1)),
4126 "", "", vect_diff);
4127 cpl_free(vectors);
4128 cpl_errorstate_set (prestate);
4129 }
4130
4131 FREE(cpl_vector_delete, vect_phase_sc_ft);
4132 FREE(cpl_vector_delete, vect_dopd_met);
4133 FREE(cpl_vector_delete, vect_diff);
4134 gsl_matrix_free (A);
4135 gsl_matrix_free (A_data);
4136 gsl_matrix_free (U);
4137 gsl_matrix_free (V);
4138 gsl_vector_free (x);
4139 gsl_vector_free (bis);
4140 gsl_vector_free (bis_data);
4141 cpl_vector_delete(time_A);
4142 gsl_matrix_free (X);
4143 gsl_vector_free (S);
4144 gsl_vector_free (work);
4145
4146 gravi_msg_function_exit(1);
4147 return opd_coeff;
4148}
4149
4150
4151
4152
4153cpl_table * gravi_vis_create_opdguess_sc (cpl_table * spectrumsc_table,
4154 cpl_table * vis_FT,
4155 cpl_table * vis_MET,
4156 double dit_sc)
4157{
4158 gravi_msg_function_start(1);
4159 cpl_ensure (spectrumsc_table, CPL_ERROR_NULL_INPUT, NULL);
4160 cpl_ensure (vis_FT, CPL_ERROR_NULL_INPUT, NULL);
4161 cpl_ensure (vis_MET, CPL_ERROR_NULL_INPUT, NULL);
4162 cpl_ensure (dit_sc>0, CPL_ERROR_ILLEGAL_INPUT, NULL);
4163
4164 cpl_size nbase = 6, ntel = 4;
4165 cpl_size nrow_sc = cpl_table_get_nrow (spectrumsc_table);
4166
4167 /* Create table */
4168 cpl_table * vis_SC = cpl_table_new (nrow_sc * nbase);
4169
4170 /* Create the time colum */
4171 gravi_table_new_column (vis_SC, "TIME", "us", CPL_TYPE_INT);
4172 for (cpl_size row_sc = 0; row_sc < nrow_sc; row_sc ++) {
4173 double value = cpl_table_get (spectrumsc_table, "TIME", row_sc, NULL);
4174 for (cpl_size base = 0; base < nbase; base++)
4175 cpl_table_set (vis_SC, "TIME", row_sc*nbase+base, value);
4176 }
4177
4178 /* Create synch columns if not yet existing */
4179 gravi_signal_create_sync (vis_SC, 6, dit_sc, vis_MET, 4, "MET");
4180 gravi_signal_create_sync (vis_SC, 6, dit_sc, vis_FT, 6, "FT");
4181
4182 CPLCHECK_NUL ("Cannot create synch");
4183
4184 /* Get SC data */
4185 int * first_met = cpl_table_get_data_int (vis_SC, "FIRST_MET");
4186 int * last_met = cpl_table_get_data_int (vis_SC, "LAST_MET");
4187 int * first_ft = cpl_table_get_data_int (vis_SC, "FIRST_FT");
4188 int * last_ft = cpl_table_get_data_int (vis_SC, "LAST_FT");
4189
4190 CPLCHECK_NUL ("Cannot get data");
4191
4192 /* Get MET and data FT data */
4193 double * phase_met = cpl_table_get_data_double (vis_MET, "PHASE_FC");
4194 double * phase_ft = cpl_table_get_data_double (vis_FT, "PHASE");
4195
4196 CPLCHECK_NUL ("Cannot get direct pointer to data");
4197
4198 /* New columns */
4199 gravi_table_new_column (vis_SC, "OPD", "m", CPL_TYPE_DOUBLE);
4200 double * opd_sc = cpl_table_get_data_double (vis_SC, "OPD");
4201
4202 /* Loop on base and rows */
4203 for (cpl_size base = 0; base < nbase; base++) {
4204 for (cpl_size row_sc = 0; row_sc < nrow_sc; row_sc ++) {
4205 cpl_size nsc = row_sc * nbase + base;
4206
4207 /* Sum over synch MET frames */
4208 double opd_met = 0.0;
4209 for (cpl_size row_met = first_met[nsc] ; row_met < last_met[nsc]; row_met++) {
4210 cpl_size nmet0 = row_met * ntel + GRAVI_BASE_TEL[base][0];
4211 cpl_size nmet1 = row_met * ntel + GRAVI_BASE_TEL[base][1];
4212 opd_met += (phase_met[nmet0] - phase_met[nmet1]) * LAMBDA_MET / CPL_MATH_2PI;
4213 }
4214 cpl_size nframe_met = last_met[nsc] - first_met[nsc];
4215 if (nframe_met != 0 ) opd_met /= nframe_met;
4216
4217 /* Sum over synch FT frames */
4218 double opd_ft = 0.0;
4219 for (cpl_size row_ft = first_ft[nsc] ; row_ft < last_ft[nsc]; row_ft++) {
4220 cpl_size nft = row_ft * nbase + base;
4221 opd_ft += phase_ft[nft] * 2.2e-6 / CPL_MATH_2PI;
4222 }
4223 cpl_size nframe_ft = last_ft[nsc] - first_ft[nsc];
4224 if (nframe_ft != 0 ) opd_ft /= nframe_ft;
4225
4226 /* Set the guess */
4227 opd_sc[nsc] = opd_met + opd_ft;
4228
4229 CPLCHECK_NUL ("Cannot build opdguess");
4230 } /* End loop on SC rows */
4231 } /* End loop on base */
4232
4233 gravi_msg_function_exit(1);
4234 return vis_SC;
4235}
4236
4237
4238cpl_error_code gravi_phase_correct_closures (cpl_table * phase_table)
4239{
4240 gravi_msg_function_start(1);
4241 cpl_ensure_code (phase_table, CPL_ERROR_NULL_INPUT);
4242
4243 int ntel = 4, nbase = 6;
4244 int tel_1[6] = {0,0,0,1,1,2};
4245 int tel_2[6] = {1,2,3,2,3,3};
4246 double ai;
4247 gsl_matrix * gsl_SC_kernel = NULL;
4248 gsl_vector * phase_coeff = NULL;
4249 gsl_matrix * M = NULL;
4250 gsl_matrix * U = NULL;
4251 gsl_matrix * V = NULL;
4252 gsl_vector * S = NULL;
4253 gsl_vector * work = NULL;
4254 cpl_matrix * M_inv = NULL;
4255 cpl_matrix * SC_kernel_mat_final = NULL;
4256 gsl_vector * Ks = NULL;
4257 cpl_matrix * kernel = NULL;
4258
4259 /* Compute of the matrix M
4260 * [1 -1 0 0
4261 * M = 1 0 -1 0
4262 * 1 0 0 -1
4263 * 0 1 -1 0
4264 * 0 1 0 -1
4265 * 0 0 1 -1] */
4266 M = gsl_matrix_alloc (nbase,ntel);
4267 gsl_matrix_set_zero (M);
4268
4269 for (int j = 0; j < nbase; j++){
4270 gsl_matrix_set (M, j, tel_1[j], 1);
4271 gsl_matrix_set (M, j, tel_2[j], -1);
4272 }
4273
4274 CPLCHECK_CLEAN ("Cannot set matrix M");
4275
4276 /*
4277 * Compute of the SV dec of M
4278 */
4279
4280 U = gsl_matrix_alloc (6, 4);
4281 V = gsl_matrix_alloc (4, 4);
4282 S = gsl_vector_alloc (4);
4283 work = gsl_vector_alloc (4);
4284
4285 gsl_matrix_memcpy (U, M);
4286 gsl_linalg_SV_decomp (U, V, S, work);
4287
4288 /*
4289 * Get inverse of the M matrix
4290 */
4291 double wv_at;
4292 double * a_inv_data = cpl_malloc(ntel * nbase * sizeof(double));
4293 for (int j = 0; j < nbase; j++) {
4294 for (int i = 0; i < ntel; i++){
4295 wv_at = 0;
4296 for (int ii = 0; ii < ntel; ii++){
4297 if( gsl_vector_get(S, ii) > 1e-14)
4298 wv_at += gsl_matrix_get(V, i, ii) / gsl_vector_get(S, ii) *
4299 gsl_matrix_get(U, j, ii);
4300 }
4301 a_inv_data[j + i * nbase] = wv_at;
4302 }
4303 }
4304
4305 gsl_matrix_free (V);
4306 gsl_matrix_free (U);
4307 gsl_vector_free (S);
4308 gsl_vector_free (work);
4309
4310 M_inv = cpl_matrix_wrap (ntel, nbase, a_inv_data);
4311
4312 CPLCHECK_CLEAN ("Cannot invers matrix M");
4313
4314 /*
4315 * Compute of the matrix Kernel of M
4316 */
4317
4318 kernel = cpl_matrix_new (nbase, nbase);
4319
4320 for (int j = 0; j < nbase; j++) {
4321 for (int k = 0; k < nbase; k++) {
4322 ai = 0;
4323 for (int ii = 0; ii < ntel; ii++){
4324 ai += gsl_matrix_get (M, j, ii)*cpl_matrix_get (M_inv, ii, k);
4325 }
4326 if (j == k)
4327 cpl_matrix_set (kernel, j, k, ai - 1);
4328 else
4329 cpl_matrix_set (kernel, j, k, ai);
4330 }
4331 }
4332
4333 cpl_matrix_unwrap (M_inv);
4334 cpl_free (a_inv_data);
4335
4336 CPLCHECK_CLEAN ("Cannot compute the matrix kernel of M");
4337
4338
4339 /*
4340 * Introducing the phase of each base and compute the Ks vector
4341 */
4342
4343 int nrow = cpl_table_get_nrow (phase_table) / nbase;
4344 double * ph_data;
4345 double * pphase = cpl_table_get_data_double (phase_table, "PHASE");
4346 cpl_ensure_code (pphase, CPL_ERROR_ILLEGAL_INPUT);
4347
4348
4349 Ks = gsl_vector_alloc (nbase* nrow);
4350 for (int k = 0; k < nbase; k++) {
4351 for (cpl_size row = 0;row < nrow; row++) {
4352 ai = 0;
4353 for (int ii = 0; ii < nbase; ii++){
4354 ai += cpl_matrix_get (kernel, k, ii) * pphase[row*nbase+ii];
4355 }
4356 gsl_vector_set (Ks, row + k * nrow, ai);
4357
4358 CPLCHECK_CLEAN ("Cannot multply phase kernel to compute Ks");
4359 }
4360 }
4361
4362 /*
4363 * Construction of the big matrix T
4364 */
4365
4366 double * ones_t;
4367 cpl_array * phase_temp;
4368
4369 double * SC_kernel = cpl_malloc ( nbase*nbase * nrow * sizeof (double));
4370 double * SC_kernel2 = cpl_malloc ( nbase*nbase * nrow * sizeof (double));
4371 cpl_array * ones = cpl_array_new (nrow, CPL_TYPE_DOUBLE);
4372
4373 for (int i = 0; i < nbase; i++) {
4374 for (int k = 0; k < nbase; k ++){
4375
4376 // phase_temp = cpl_array_duplicate (phase [i]);
4377 phase_temp = cpl_array_new (nrow, CPL_TYPE_DOUBLE);
4378 for (cpl_size row = 0; row < nrow; row++)
4379 cpl_array_set (phase_temp, row, pphase[row*nbase+i]);
4380
4381 cpl_array_fill_window_double (ones, 0, nrow, 1);
4382 cpl_array_multiply_scalar (ones, cpl_matrix_get (kernel, i, k));
4383 cpl_array_multiply_scalar (phase_temp, cpl_matrix_get (kernel, i, k));
4384 ph_data = cpl_array_get_data_double (phase_temp);
4385 ones_t = cpl_array_get_data_double (ones);
4386
4387 memcpy (SC_kernel2 + (i*nbase + k) * nrow , ph_data, nrow * sizeof(double));
4388 memcpy (SC_kernel + (i*nbase + k) * nrow, ones_t, nrow * sizeof(double));
4389
4390 CPLCHECK_CLEAN ("Error in loop SC_kernel");
4391 cpl_array_delete (phase_temp);
4392 }
4393
4394 }
4395
4396 FREE (cpl_array_delete, ones);
4397 CPLCHECK_CLEAN ("Cannot compute SC_kernel");
4398
4399 cpl_matrix * temp = cpl_matrix_wrap (nbase, nbase* nrow, SC_kernel);
4400 cpl_matrix * SC_kernel_mat = cpl_matrix_duplicate (temp);
4401 cpl_matrix * SC_kernel_mat2 = cpl_matrix_wrap (nbase, nbase* nrow, SC_kernel2);
4402
4403 cpl_matrix_append (SC_kernel_mat, SC_kernel_mat2, 1);
4404
4405 SC_kernel_mat_final = cpl_matrix_extract (SC_kernel_mat, 0, 0, 1, 1, 11, nbase*nrow);
4406
4407 cpl_matrix_unwrap (temp);
4408 cpl_matrix_unwrap (SC_kernel_mat2);
4409 cpl_free (SC_kernel2);
4410 cpl_free (SC_kernel);
4411 cpl_matrix_delete (SC_kernel_mat);
4412
4413 gsl_SC_kernel = gsl_matrix_alloc (nbase*nrow, 11);
4414
4415 for (int j = 0; j < cpl_matrix_get_nrow(SC_kernel_mat_final); j++){
4416 for (int k = 0; k < cpl_matrix_get_ncol(SC_kernel_mat_final); k++)
4417 gsl_matrix_set (gsl_SC_kernel, k, j, cpl_matrix_get (SC_kernel_mat_final, j, k));
4418 }
4419
4420 /*
4421 * Resolve the linear equation T * x = Ks
4422 */
4423
4424 U = gsl_matrix_alloc (nbase*nrow, 11);
4425 V = gsl_matrix_alloc (11, 11);
4426 S = gsl_vector_alloc (11);
4427 phase_coeff = gsl_vector_alloc (11);
4428 work = gsl_vector_alloc (11);
4429 gsl_matrix_memcpy (U, gsl_SC_kernel);
4430
4431 gsl_linalg_SV_decomp (U, V, S, work);
4432 for (int i = 0; i < 11; i++)
4433 if (gsl_vector_get (S, i) < 1e-10)
4434 gsl_vector_set (S, i , 0);
4435 gsl_linalg_SV_solve (U, V, S, Ks, phase_coeff);
4436
4437 /*
4438 * Apply the correction on 5 baselines
4439 */
4440
4441 for (int base = 0; base < nbase - 1; base ++){
4442 double f = gsl_vector_get (phase_coeff, nbase + base);
4443 cpl_msg_info (cpl_func,"correction factor = 1 %+.20f", f);
4444 for (cpl_size row = 0; row < nrow; row++) {
4445 pphase[row*nbase+base] *= 1 - f;
4446 }
4447 }
4448
4449 /*
4450 * Cleanup
4451 */
4452
4453 cleanup :
4454 FREE (gsl_matrix_free, M);
4455 FREE (gsl_vector_free, Ks);
4456 FREE (cpl_matrix_delete, SC_kernel_mat_final);
4457 FREE (cpl_matrix_delete, kernel);
4458 FREE (gsl_matrix_free, gsl_SC_kernel);
4459 FREE (gsl_matrix_free, U);
4460 FREE (gsl_matrix_free, V);
4461 FREE (gsl_vector_free, S);
4462 FREE (gsl_vector_free, phase_coeff);
4463 FREE (gsl_vector_free, work);
4464
4465 gravi_msg_function_exit(1);
4466 return CPL_ERROR_NONE;
4467}
4468
4469
4470double gravi_envelope (double lambda, double delta_lambda, double opd)
4471{
4472 /* Compute coherent length [m] */
4473 double coh_len= lambda * lambda / delta_lambda / 3;
4474
4475 /* Gaussien enveloppe */
4476 double value = 1.0;
4477 if (opd != 0.0)
4478 value = exp (-1 * opd * opd / (coh_len*coh_len/2.) );
4479
4480 return value;
4481}
4482
4483
4484
4485 /* TEST */
4486 if (wave == nwave/2 && base == 0 & pol == 0) {
4487 int step = nwave == 5 ? 100 : 1;
4488
4489 cpl_vector ** array = cpl_malloc (3 * sizeof(cpl_vector*));
4490 array[0] = gravi_vector_extract (opd_vector, 0, step);
4491 cpl_vector_multiply_scalar (array[0], 1e6);
4492
4493 int iA = gravi_get_region (detector_table, base, 'A', 0);
4494 cpl_vector * tmp = gravi_table_get_vector (spectrum_table, wave, GRAVI_DATA[iA]);
4495 cpl_vector_subtract_scalar (tmp, cpl_vector_get_mean (tmp));
4496 array[1] = gravi_vector_extract (tmp, 0, step);
4497 cpl_vector_divide_scalar (array[1], cpl_vector_get_max (array[1]));
4498 cpl_vector_delete (tmp);
4499
4500 array[2] = gravi_vector_extract (envelope_vector, 0, step);
4501
4502 cpl_plot_vectors (NULL, NULL, NULL, (const cpl_vector **)array, 3);
4503 FREELOOP (cpl_vector_delete, array, 3);
4504 }
4505
4506
4507
4508// if (!(reg%4) && wave == nwave/2 && iA == 50) {
4509// int step = 100;
4510// cpl_vector ** array = cpl_calloc (4,sizeof(cpl_vector*));
4511// array[0] = gravi_vector_extract (opd_vector, 0, step);
4512// array[1] = gravi_vector_extract (X_vector, 0, step);
4513// array[2] = gravi_vector_extract (R_vector, 0, step);
4514// array[3] = gravi_vector_extract (I_vector, 0, step);
4515// cpl_plot_vectors (NULL, NULL, NULL, (const cpl_vector **)array, 4);
4516// cpl_free (array);
4517// CPLCHECK_NUL ("Cannot plot");
4518// }
4519
4520
4521 // FIXME:
4522 cpl_msg_info (cpl_func, "Compute WAVE_SCAN for %s", GRAVI_TYPE(type_data));
4523 cpl_table * wavescan_table;
4524 wavescan_table = gravi_wave_scan (spectrum_table, detector_table, opd_table);
4525
4526 gravi_data_add_table (wave_map, NULL, type_data == GRAVI_SC ? "WAVE_SCAN_SC" : "WAVE_SCAN_FT",
4527 wavescan_table);
4528
4529
4530
4531cpl_table * gravi_wave_scan (cpl_table * spectrum_table,
4532 cpl_table * detector_table,
4533 cpl_table * opd_table)
4534{
4535 gravi_msg_function_start(1);
4536 cpl_ensure (spectrum_table, CPL_ERROR_NULL_INPUT, NULL);
4537 cpl_ensure (detector_table, CPL_ERROR_NULL_INPUT, NULL);
4538 cpl_ensure (opd_table, CPL_ERROR_NULL_INPUT, NULL);
4539
4540 int nbase = 6;
4541 // char name[100];
4542
4543 /* Create the output table */
4544 cpl_table * wave_table = cpl_table_new (1);
4545
4546 /* Get the number of wavelength, region, polarisation... */
4547 cpl_size nwave = gravi_spectrum_get_nwave (spectrum_table);
4548 cpl_size nregion = cpl_table_get_nrow (detector_table);
4549 cpl_size nrow = cpl_table_get_nrow (spectrum_table);
4550
4551 /* To save results */
4552 cpl_array * wave_array = cpl_array_new (nwave, CPL_TYPE_DOUBLE);
4553
4554 /*
4555 * Calibration of each polarization and base
4556 */
4557
4558 for (int reg = 0; reg < nregion; reg++) {
4559 cpl_msg_info (cpl_func, "Compute wave region %i over %lli", reg+1, nregion);
4560
4561 /* Get the base of this region */
4562 int base = gravi_region_get_base (detector_table, reg);
4563
4564 /* Sign of this baseline */
4565 double phi_sign = gravi_region_get_base_sign (detector_table, base);
4566
4567 /* Get OPD of this region */
4568 cpl_vector * opd_vector = cpl_vector_new (nrow);
4569 for (cpl_size row = 0; row < nrow; row ++ ) {
4570 double value = cpl_table_get (opd_table, "OPD", row*nbase+base, NULL);
4571 cpl_vector_set (opd_vector, row, value);
4572 }
4573 CPLCHECK_NUL ("Cannot get opd");
4574
4575 cpl_msg_info (cpl_func,"opd min = %g max = %g [um]",
4576 cpl_vector_get_min (opd_vector)*1e6,
4577 cpl_vector_get_max (opd_vector)*1e6);
4578
4579 /* Loop on wave */
4580 for (cpl_size wave = 0; wave < nwave; wave++) {
4581
4582 cpl_vector * X_vector;
4583 X_vector = gravi_table_get_vector (spectrum_table, wave, GRAVI_DATA[reg]);
4584 cpl_vector_subtract_scalar (X_vector, cpl_vector_get_mean (X_vector));
4585 cpl_vector_divide_scalar (X_vector, 1.4 * cpl_vector_get_stdev (X_vector));
4586 CPLCHECK_NUL ("Cannot get data");
4587
4588 /* Guess wave */
4589 double lbd0 = 1.97e-6 + (2.48e-6 - 1.97e-6) / nwave * wave;
4590
4591 /* Compute envelope from OPD for this channel */
4592 cpl_vector * env_vector = gravi_compute_envelope (opd_vector, wave, nwave);
4593
4594 /* Search for a better match */
4595 cpl_size nA = (nwave == 5 ? 20 : 100);
4596 double searchA = 0.05;
4597 cpl_vector * V_vector = cpl_vector_new (nA);
4598 cpl_vector * R_vector = cpl_vector_new (nrow);
4599 cpl_vector * I_vector = cpl_vector_new (nrow);
4600
4601 for (cpl_size iA = 0; iA < nA; iA++) {
4602 double A = phi_sign * (1 - searchA + (2*searchA * iA) / nA ) * CPL_MATH_2PI / lbd0;
4603
4604 for (cpl_size row = 0; row < nrow; row++) {
4605 cpl_vector_set (R_vector, row, cos (A * cpl_vector_get (opd_vector, row)) *
4606 cpl_vector_get (env_vector, row) );
4607 cpl_vector_set (I_vector, row, sin (A * cpl_vector_get (opd_vector, row)) *
4608 cpl_vector_get (env_vector, row));
4609 }
4610
4611 cpl_vector_multiply (R_vector, X_vector);
4612 cpl_vector_multiply (I_vector, X_vector);
4613 double Rvalue = cpl_vector_get_mean (R_vector);
4614 double Ivalue = cpl_vector_get_mean (I_vector);
4615
4616 /* Fill chi2 */
4617 cpl_vector_set (V_vector, iA, Rvalue*Rvalue + Ivalue*Ivalue);
4618 CPLCHECK_NUL ("Cannot fill V_vector");
4619 }
4620
4621 /* Plot chi2 versus iA */
4622 if (!(reg%8) && (wave == 0 || wave == nwave-1)) {
4623 cpl_plot_vector (NULL, NULL, NULL, V_vector);
4624 CPLCHECK_NUL ("Cannot plot");
4625 }
4626
4627 /* Compute wavelength */
4628 double lbd = ( (1-searchA) + 2*searchA * gravi_vector_get_maxpos (V_vector) / nA ) * lbd0;
4629 cpl_array_set (wave_array, wave, lbd);
4630
4631 FREE (cpl_vector_delete, env_vector);
4632 FREE (cpl_vector_delete, V_vector);
4633 FREE (cpl_vector_delete, R_vector);
4634 FREE (cpl_vector_delete, I_vector);
4635 FREE (cpl_vector_delete, X_vector);
4636 CPLCHECK_NUL ("Cannot delete");
4637 }
4638
4639 /* Add column */
4640 gravi_table_new_column_array (wave_table, GRAVI_DATA[reg],
4641 "m", CPL_TYPE_DOUBLE, nwave);
4642 cpl_table_set_array (wave_table, GRAVI_DATA[reg], 0, wave_array);
4643
4644 cpl_vector_delete (opd_vector);
4645 }
4646
4647 cpl_array_delete (wave_array);
4648
4649 gravi_msg_function_exit(1);
4650 return wave_table;
4651}
4652
4653
4654/*
4655 * Compute the group-delay in [m] as the maximum of |Env(x)| where
4656 * Env(x) = < visdata(lbd) * exp(2i.pi * x / lbd) > with <> sum over lbd
4657 *
4658 * Not working in HIGH because the
4659 */
4660cpl_error_code gravi_array_get_group_delay_loop_new (cpl_array ** input, cpl_array * sigma,
4661 double * gd, cpl_size nrow)
4662{
4663 gravi_msg_function_start(1);
4664 cpl_ensure_code (input, CPL_ERROR_NULL_INPUT);
4665 cpl_ensure_code (sigma, CPL_ERROR_NULL_INPUT);
4666 cpl_ensure_code (gd, CPL_ERROR_ILLEGAL_OUTPUT);
4667
4668 cpl_size nsigma = cpl_array_get_size (sigma);
4669 double lbd = 1.0 / cpl_array_get (sigma,nsigma/2,NULL);
4670
4671 /* Copy data as double to secure their type */
4672 double * sigdata = cpl_malloc (sizeof(double complex) * nsigma);
4673 for (cpl_size w=0; w<nsigma; w++) sigdata[w] = cpl_array_get (sigma, w, NULL);
4674
4675 /* Coherence lenght in [m] */
4676 int size = 50;
4677 double coherence = 0.5 * nsigma / fabs (sigdata[0] - sigdata[nsigma-1]);
4678 double coherence_bb = 1 / fabs (sigdata[0] - sigdata[nsigma-1]);
4679 double width = CPL_MIN (size * coherence_bb, coherence);
4680
4681 cpl_msg_info ("TEST", "coh = %g and %g [um]", coherence*1e6, coherence_bb*1e6);
4682
4683 /* 10 maximum waveform */
4684 double complex ** waveform = cpl_calloc (10, sizeof(double complex*));
4685 double * step = cpl_calloc (10, sizeof(double));
4686
4687 /* Allocate waveforms to explore */
4688 int nloop = 0;
4689 while (width > lbd/50) {
4690 step[nloop] = width / size;
4691 cpl_msg_info (cpl_func, "Build waveform of %g [um] with step [%gum]",
4692 width *1e6, step[nloop] *1e6);
4693
4694 waveform[nloop] = cpl_malloc (sizeof(double complex) * size * nsigma);
4695 for (cpl_size s=0; s<size; s++) {
4696 double x = (s-size/2) * step[nloop];
4697 for (cpl_size w=0; w<nsigma; w++)
4698 waveform[nloop][s*nsigma+w] = cexp (-2.*I*CPL_MATH_PI * x * sigdata[w]);
4699 }
4700
4701 width /= 10.;
4702 nloop ++;
4703 }
4704
4705 /* Allocate memory for the grid search */
4706 cpl_msg_info (cpl_func, "Loop on %lli rows to compute gdelay", nrow);
4707
4708 /* Memory for data */
4709 double complex * visdata = cpl_malloc (sizeof(double complex) * nsigma);
4710
4711 /* Loop on rows */
4712 for (cpl_size row = 0; row<nrow; row++) {
4713 gd[row] = 0.0;
4714
4715 /* Copy data as double complex to secure their type */
4716 for (cpl_size w=0; w<nsigma; w++) visdata[w] = cpl_array_get_complex (input[row], w, NULL);
4717
4718 /* Loop on resolution */
4719 for (int loop = 0; loop < nloop; loop++) {
4720 double current_max = -1, current_gd = 0.0;
4721
4722 /* Loop on x to find the maximum of P(x) = |FT(input(sigma))| */
4723 for (cpl_size s=0; s<size; s++) {
4724 double complex tmp = 0.0 * I + 0.0;
4725 for (cpl_size w=0; w<nsigma; w++) {tmp += visdata[w] * waveform[loop][s*nsigma+w];}
4726 double P = cabs (tmp);
4727 if ( P > current_max) { current_max = P; current_gd = (s-size/2) * step[loop]; }
4728 }
4729
4730 /* Remove gd */
4731 for (cpl_size w=0; w<nsigma; w++) visdata[w] *= cexp (-2.*I*CPL_MATH_PI*current_gd*sigdata[w]);
4732 gd[row] += current_gd;
4733
4734 } /* End loop on resolution */
4735
4736 CPLCHECK_MSG("Cannot compute GD");
4737 } /* End loop on rows */
4738
4739 /* Clean memory */
4740 FREE (cpl_free, visdata);
4741 FREE (cpl_free, sigdata);
4742 FREE (cpl_free, step);
4743 FREELOOP (cpl_free, waveform, 10);
4744
4745 gravi_msg_function_exit(1);
4746 return CPL_ERROR_NONE;
4747}
4748
4749
4750
4751/*
4752 * Compute the group-delay in [m] as the maximum of |Env(x)| where
4753 * Env(x) = < visdata(lbd) * exp(2i.pi * x / lbd) > with <> sum over lbd
4754 *
4755 * Perform a double pass:
4756 * - crude search over the entire coherence length, resolution 1lbd
4757 * - fine search within 2lbd, resolution lbd/100
4758 */
4759double gravi_array_get_group_delay (cpl_array * input, cpl_array * sigma)
4760{
4761 cpl_ensure (input, CPL_ERROR_NULL_INPUT, 0);
4762 cpl_ensure (sigma, CPL_ERROR_NULL_INPUT, 0);
4763
4764 int nv = 0;
4765 double width, step;
4766 cpl_size nsigma = cpl_array_get_size (sigma);
4767 cpl_size ninput = cpl_array_get_size (input);
4768
4769 cpl_ensure (nsigma == ninput, CPL_ERROR_ILLEGAL_INPUT, 0);
4770
4771 /* Width of a single spectral channel in [m]
4772 * Step is 1 lambda */
4773 width = 1.0 * nsigma / fabs (cpl_array_get (sigma,0,&nv) - cpl_array_get (sigma,nsigma-1,&nv));
4774 step = 1. / cpl_array_get (sigma,nsigma/2,&nv);
4775
4776 /* Init the search for maximum */
4777 double gd0 = 0.0;
4778 double current_gd = -1e10;
4779 double current_max = -1.0;
4780
4781 /* Loop on x to find the maximum of P(x) = |FT(input(sigma))|
4782 * FIXME: the get_complex in the middle of the loop is very slow */
4783 for (double x = gd0-width/2.0 ; x < gd0+width/2.0 ; x += step) {
4784 double complex tmp = 0.0 * I + 0.0;
4785 for (cpl_size w=0; w<nsigma; w++)
4786 tmp += cpl_array_get_complex (input, w, &nv) * cexp (-2.*I*CPL_MATH_PI * x * cpl_array_get (sigma,w,&nv));
4787 if ( cabs (tmp) > current_max ) {
4788 current_max = cabs (tmp);
4789 current_gd = x;
4790 }
4791 }
4792
4793 /* Width is 4 lambda
4794 * Step is lambda/100 */
4795 width = 4.0 * step;
4796 step = step / 100.0;
4797
4798 /* Re-init the search for maximum */
4799 gd0 = current_gd;
4800 current_gd = -1e10;
4801 current_max = -1.0;
4802
4803 /* Loop on x to find the maximum of P(x) = |FT(input(sigma))| */
4804 for (double x = gd0-width/2 ; x < gd0+width/2 ; x += step) {
4805 double complex tmp = 0.0 * I + 0.0;
4806 for (cpl_size w=0; w<nsigma; w++)
4807 tmp += cpl_array_get_complex (input, w, &nv) * cexp (-2.*I*CPL_MATH_PI * x * cpl_array_get (sigma,w,&nv));
4808 if ( cabs (tmp) > current_max ) {
4809 current_max = cabs (tmp);
4810 current_gd = x;
4811 }
4812 }
4813
4814 CPLCHECK("Cannot compute GD");
4815
4816 return current_gd;
4817}
4818
4819/*
4820 * Compute the group-delay in [m] as the maximum of |Env(x)| where
4821 * Env(x) = < visdata(lbd) * exp(2i.pi * x / lbd) > with <> sum over lbd
4822 *
4823 * Perform a single pass:
4824 * - crude search over the entire coherence length, resolution 1lbd
4825 */
4826double gravi_array_get_group_delay_coarse (cpl_array * input, cpl_array * sigma)
4827{
4828 cpl_ensure (input, CPL_ERROR_NULL_INPUT, 0);
4829 cpl_ensure (sigma, CPL_ERROR_NULL_INPUT, 0);
4830
4831 int nv = 0;
4832 double width, step;
4833 cpl_size nsigma = cpl_array_get_size (sigma);
4834 cpl_size ninput = cpl_array_get_size (input);
4835
4836 cpl_ensure (nsigma == ninput, CPL_ERROR_ILLEGAL_INPUT, 0);
4837
4838 /* Width of a single spectral channel in [m]
4839 * Step is 1 lambda */
4840 width = 1.0 * nsigma / fabs (cpl_array_get (sigma,0,&nv) - cpl_array_get (sigma,nsigma-1,&nv));
4841 step = 1. / cpl_array_get (sigma,nsigma/2,&nv);
4842
4843 /* Init the search for maximum */
4844 double gd0 = 0.0;
4845 double current_gd = -1e10;
4846 double current_max = -1.0;
4847
4848 /* Loop on x to find the maximum of P(x) = |FT(input(sigma))|
4849 * FIXME: the get_complex in the middle of the loop is very slow */
4850 for (double x = gd0-width/2.0 ; x < gd0+width/2.0 ; x += step) {
4851 double complex tmp = 0.0 * I + 0.0;
4852 for (cpl_size w=0; w<nsigma; w++)
4853 tmp += cpl_array_get_complex (input, w, &nv) * cexp (-2.*I*CPL_MATH_PI * x * cpl_array_get (sigma,w,&nv));
4854 if ( cabs (tmp) > current_max ) {
4855 current_max = cabs (tmp);
4856 current_gd = x;
4857 }
4858 }
4859
4860 CPLCHECK("Cannot compute GD");
4861
4862 return current_gd;
4863}
4864
4865
4866/*
4867 * Compute the group-delay in [m] with the Pedretii algorithm
4868 * http://cdsads.u-strasbg.fr/abs/2005ApOpt..44.5173P
4869 */
4870double gravi_array_get_group_delay_iota (cpl_array * input, cpl_array * wavenumber)
4871{
4872 cpl_ensure (input, CPL_ERROR_NULL_INPUT, 0);
4873 cpl_ensure (wavenumber, CPL_ERROR_NULL_INPUT, 0);
4874
4875 double complex interspectra = 0.0 + I * 0.0;
4876 double sum_wavenumber_diff = 0.0;
4877
4878 cpl_size size = cpl_array_get_size (input);
4879 cpl_size sizew = cpl_array_get_size (wavenumber);
4880
4881 cpl_ensure (size == sizew, CPL_ERROR_ILLEGAL_INPUT, 0);
4882
4883 for (cpl_size n = 1; n < size; n++) {
4884
4885 /* InterSpectre = < visData_sc_exp(i) * visData_sc_exp(i+1) > over lbd */
4886 interspectra += cpl_array_get_complex (input, n, NULL) *
4887 conj (cpl_array_get_complex (input, n - 1, NULL));
4888
4889 /* Sum of the wavenumber differences */
4890 sum_wavenumber_diff += cpl_array_get (wavenumber, n, NULL) - cpl_array_get (wavenumber, n-1, NULL);
4891 }
4892
4893 /* Return the group delay in [m] */
4894 return carg ( interspectra ) / sum_wavenumber_diff * (size-1) / (2*M_PI);
4895}
4896
4897
4898
4899//cpl_error_code gravi_data_erase_obs (gravi_data * data,
4900// const cpl_parameterlist params)
4901//{
4902// cpl_ensure (data, CPL_ERROR_NULL_INPUT, NULL);
4903// cpl_ensure (params, CPL_ERROR_NULL_INPUT, NULL);
4904//
4905// if (params)
4906//
4907// /* Rejection flag (1 = rejected) */
4908// cpl_size max_obs = 10000;
4909// cpl_array * flag_array = cpl_array_new (max_obs, CPL_TYPE_INT);
4910// cpl_array_fill_window (flag_array, 0, max_obs, 0);
4911//
4912// /* Loop on extension */
4913// for (int ext = 0; ext < data->nb_ext ; ext++) {
4914//
4915// const char * extname = gravi_pfits_get_extname (data->exts_hdrs[ext]);
4916// if (strstr (extname, GRAVI_OI_VIS_EXT) ||
4917// strstr (extname, GRAVI_OI_VIS2_EXT) ) {
4918//
4919// }
4920// if (strstr (extname, GRAVI_OI_FLUX_EXT) ||
4921// strstr (extname, GRAVI_OI_T3_EXT) ) {
4922//
4923// }
4924//
4925// } /* End loop on extensions */
4926//
4927// return CPL_ERROR_NONE;
4928//}
4929
4930
4931
4932
4933/* TODO, FIXME
4934 * Test function to find a way to compute the bias of the column without
4935 * introducing noise between columns by fitting the bias with polynome */
4936double gravi_bivector_get_med_poly (cpl_bivector * bivector_in)
4937{
4938 cpl_ensure_code (bivector_in, CPL_ERROR_NULL_INPUT);
4939
4940 cpl_bivector_sort(bivector_in, bivector_in, CPL_SORT_ASCENDING, CPL_SORT_BY_Y);
4941
4942 int size = cpl_bivector_get_size(bivector_in);
4943 int sizeout = size*(0.5);
4944 int start = (size-sizeout)/2;
4945 cpl_matrix * coord = cpl_matrix_new(1,sizeout);
4946 cpl_vector * vector = cpl_vector_new(sizeout);
4947
4948 for (cpl_size i = 0 ; i < sizeout ; i++){
4949 cpl_vector_set(vector, i, cpl_vector_get(cpl_bivector_get_y(bivector_in), i+start));
4950 cpl_matrix_set(coord, 0, i, cpl_vector_get(cpl_bivector_get_x(bivector_in), i+start));
4951 }
4952 CPLCHECK_MSG ("Cannot clip the data");
4953
4954 cpl_size power = 2;
4955 cpl_polynomial * poly = cpl_polynomial_new(1);
4956 cpl_polynomial_fit (poly, coord, NULL, vector, NULL,
4957 CPL_FALSE, NULL, &power);
4958
4959 sizeout=cpl_vector_get (cpl_bivector_get_y (bivector_in), size-1)-cpl_vector_get (cpl_bivector_get_y (bivector_in), 0);
4960 cpl_vector * vector_mean = cpl_vector_new (sizeout);
4961 CPLCHECK_MSG ("Cannot get median");
4962 cpl_vector_fill_polynomial (vector_mean, poly, 0, sizeout);
4963 double mean = cpl_vector_get_mean(vector_mean);
4964
4965 cpl_vector_delete(vector_mean);
4966 cpl_vector_delete(vector);
4967 cpl_matrix_delete(coord);
4968
4969 return mean;
4970}
4971
4972
4973/* TODO
4974 * Test function to find a way to compute the bias of the column without
4975 * introducing noise between columns doing a median on the value inside +-n*rms
4976 * and after remooving a percent of the extrem
4977 */
4978double gravi_vector_get_med_percent(cpl_vector * vector_in, float percent)
4979{
4980 cpl_vector_sort(vector_in, CPL_SORT_ASCENDING);
4981 int size = cpl_vector_get_size(vector_in);
4982 int sizeout = size*(1-percent*2);
4983 int start = (size-sizeout)/2;
4984 cpl_vector * vector = cpl_vector_new(sizeout);
4985
4986 for (cpl_size i = 0 ; i < sizeout ; i++)
4987 cpl_vector_set(vector, i, cpl_vector_get(vector_in, i+start));
4988
4989
4990 size = cpl_vector_get_size(vector);
4991 sizeout = 0;
4992 double med = cpl_vector_get_median(vector);
4993 double rms = cpl_vector_get_stdev(vector);
4994 int nsig = 3;
4995
4996 for (cpl_size i = 0 ; i < size ; i++)
4997 if ( (cpl_vector_get(vector, i) > med-nsig*rms) && (cpl_vector_get(vector, i) < med+nsig*rms) ) sizeout++;
4998
4999 int i_out=0;
5000 cpl_vector * vector_med = cpl_vector_new(sizeout);
5001 for (cpl_size i = 0 ; i < size ; i++)
5002 if ( (cpl_vector_get(vector, i) > med-nsig*rms) && (cpl_vector_get(vector, i) < med+nsig*rms) )
5003 cpl_vector_set(vector_med, i_out++, cpl_vector_get(vector, i));
5004
5005 return cpl_vector_get_mean(vector_med);
5006}
5007
5008
5009
5010cpl_error_code gravi_flux_create_average (cpl_table * i_table, cpl_table * o_table,
5011 const char * colname, const char * sname)
5012{
5013 gravi_msg_function_start(1);
5014 cpl_ensure_code (i_table);
5015 cpl_ensure_code (o_table);
5016 cpl_ensure_code (colname);
5017
5018 /* Check input */
5019 if ( !cpl_table_has_column (i_table, colname)) {
5020 cpl_msg_info (cpl_func, "Cannot average column %s (not in data)", colname);
5021 return CPL_ERROR_NONE;
5022 }
5023
5024
5025 /* Get synch data */
5026 int ntel = 4;
5027 cpl_size nrow_i = cpl_table_get_nrow (output_o) / ntel;
5028 int * first_i = cpl_table_get_data_int (output_o, sname);
5029 int * last_i = cpl_table_get_data_int (output_o, sname);
5030 cpl_ensure_code (first_i);
5031 cpl_ensure_code (last_i);
5032
5033 const char * unit = cpl_table_get_column_unit (i_table, colname);
5034
5035 /* Get data */
5036 double * i_data = cpl_table_get_data_double (i_table, colname);
5037 cpl_ensure_code (i_data);
5038
5039 /* Create output */
5040 gravi_table_new_column (o_table, o_table, unit, CPL_TYPE_DOUBLE);
5041 double * o_data = cpl_table_get_data_double (o_table, colname);
5042
5043 /* Loop on tel and rows */
5044 for (cpl_size tel = 0; tel < ntel; tel++) {
5045 for (cpl_size row_o = 0; row_o < nrow_o; row_o ++) {
5046 cpl_size no = row_o * ntel + tel;
5047
5048 /* Sum over synch input frames */
5049 for (cpl_size row_i = first_i[no] ; row_i < last_i[no]; row_i++) {
5050 cpl_size ni = row_i * ntel + tel;
5051 o_data [no] += i_data[ni];
5052 }
5053 CPLCHECK_MSG ("Fail to compute metrology per base from metrology per tel");
5054
5055 /* Normalize the means (if nframe == 0, values are zero) */
5056 cpl_size nframe = last_met[nsc] - first_met[nsc];
5057 if (nframe != 0 ) o_data[no] /= nframe;
5058 } /* End loop on output frames */
5059 }/* End loop on tel */
5060
5061 gravi_msg_function_exit(1);
5062 return CPL_ERROR_NONE;
5063}
5064
GRAVI_PRIMARY_HDR_EXT
#define GRAVI_PRIMARY_HDR_EXT
Definition gravi-test.h:212
gravi_data
typedefCPL_BEGIN_DECLS struct _gravi_data_ gravi_data
Definition gravi_data.h:39
gravi_data_get_spectrum_data
#define gravi_data_get_spectrum_data(data, type)
Definition gravi_data.h:63
gravi_data_get_oi_flux
#define gravi_data_get_oi_flux(data, type, pol, npol)
Definition gravi_data.h:49
gravi_data_get_header
#define gravi_data_get_header(data)
Definition gravi_data.h:75
gravi_data_get_oi_wave
#define gravi_data_get_oi_wave(data, type, pol, npol)
Definition gravi_data.h:45
gravi_data_get_oi_vis_plist
#define gravi_data_get_oi_vis_plist(data, type, pol, npol)
Definition gravi_data.h:70
gravi_data_get_imaging_detector
#define gravi_data_get_imaging_detector(data, type)
Definition gravi_data.h:60
gravi_data_get_oi_vis
#define gravi_data_get_oi_vis(data, type, pol, npol)
Definition gravi_data.h:46
ntel
const cpl_size ntel
Definition gravi_demodulate.c:84
gravi_vis_create_opdguess_sc
cpl_table * gravi_vis_create_opdguess_sc(cpl_table *spectrumsc_table, cpl_table *vis_FT, cpl_table *vis_MET, double dit_sc)
Definition gravi_old.c:4153
gravi_lazer_get_wavelength
cpl_error_code gravi_lazer_get_wavelength(gravi_data *lazer_data)
Definition gravi_old.c:47
my_gsl_matrix_fprintf
int my_gsl_matrix_fprintf(FILE *stream, gsl_matrix *m, const char *fmt)
Definition gravi_old.c:3579
gravi_opdmet_calibration
cpl_table * gravi_opdmet_calibration(cpl_table *metrology_table, double mjd_obs)
Compute the opd from metrology for the FC.
Definition gravi_old.c:3629
gravi_data_set_propertylist
gravi_data_set_propertylist(wave_map, "WAVE_DATA_SC", plist)
gravi_wave_scan
cpl_table * gravi_wave_scan(cpl_table *spectrum_table, cpl_table *detector_table, cpl_table *opd_table)
Definition gravi_old.c:4531
gravi_flux_create_average
cpl_error_code gravi_flux_create_average(cpl_table *i_table, cpl_table *o_table, const char *colname, const char *sname)
Definition gravi_old.c:5010
gravi_get_basesign_from_tel
int gravi_get_basesign_from_tel(int t1, int t2)
Definition gravi_old.c:38
gravi_compute_disp_old
gravi_data * gravi_compute_disp_old(gravi_data *vis_data)
Definition gravi_old.c:165
gravi_data_mean_metFddl
cpl_error_code gravi_data_mean_metFddl(gravi_data *oi_vis, gravi_data *p2vmred_data, gravi_data *preproc_data, int thread)
Compute mean position of MET and FDDL in each DITs.
Definition gravi_old.c:2813
wavescan_table
cpl_table * wavescan_table
Definition gravi_old.c:4523
finals2000A_mjd_last_type
static int finals2000A_mjd_last_type(FILE *pFile, char type)
Definition gravi_old.c:1861
cpl_table_save
cpl_table_save(vis_met, NULL, NULL, "vismet.fits", CPL_IO_CREATE)
gravi_visdata_fromellipse
gravi_data * gravi_visdata_fromellipse(gravi_data *spectrum_data)
Definition gravi_old.c:2016
gravi_dtps
cpl_error_code gravi_dtps(double xi, double eta, double raz, double decz, double *ra, double *dec)
Definition gravi_old.c:1621
gravi_metrology_reform
cpl_table * gravi_metrology_reform(cpl_table *vis_met)
reform the PHASE_FC column of the metrology for the FC into OPD_FC
Definition gravi_old.c:3751
gravi_p2vm_mean_spectrum
cpl_error_code gravi_p2vm_mean_spectrum(gravi_data *p2vm_map, gravi_data *preproc_data)
Definition gravi_old.c:2963
wave_data
gravi_data * wave_data
Definition gravi_old.c:1994
gravi_phase_fit_opdmet
cpl_vector * gravi_phase_fit_opdmet(cpl_table *vis_met, cpl_table *phase_sc, cpl_table *phase_ft, double dit_sc)
Definition gravi_old.c:3828
gravi_metrology_calibration
cpl_table * gravi_metrology_calibration(cpl_table *metrology_table, cpl_table *opl_table, double mjd_obs)
The given output FITS file contain a p2vm table of the metrology with the values of the transmission,...
Definition gravi_old.c:2562
wave_map_sc
cpl_table * wave_map_sc
Definition gravi_old.c:1998
gravi_bivector_get_med_poly
double gravi_bivector_get_med_poly(cpl_bivector *bivector_in)
Definition gravi_old.c:4936
specMean
cpl_image * specMean
Definition gravi_old.c:1961
gravi_array_get_group_delay
double gravi_array_get_group_delay(cpl_array *input, cpl_array *sigma)
Definition gravi_old.c:4759
gravi_compute_flat_badpix
cpl_error_code gravi_compute_flat_badpix(gravi_data *flat, gravi_data *dark)
The aim of this function is to identify the bad pixels and to update profile accordingly.
Definition gravi_old.c:1690
finals2000A_mjd_first
static int finals2000A_mjd_first(FILE *pFile)
Definition gravi_old.c:1852
gravi_data_save
cpl_error_code gravi_data_save(gravi_data *self, cpl_frameset *allframes, const char *filename, const cpl_parameterlist *parlist, cpl_frameset *usedframes, cpl_frame *frame, const char *recipe, cpl_propertylist *applist)
Definition gravi_old.c:1493
finals2000A_read_line
static void finals2000A_read_line(FILE *pFile, char *flag, double *mjd, double *pmx, double *pmy, double *dut)
Definition gravi_old.c:1841
gravi_data_set_table
gravi_data_set_table(wave_map, "WAVE_DATA_SC", cpl_table_duplicate(wave_data_sc))
cpl_image_multiply
cpl_image_multiply(flat_profiled, profile_crop)
gravi_data_add_table
gravi_data_add_table(wave_map, NULL, type_data==GRAVI_SC ? "WAVE_SCAN_SC" :"WAVE_SCAN_FT", wavescan_table)
gravi_data_save_data
gravi_data_save_data(wave_data, "wave_map_argon.fits", CPL_IO_CREATE)
cpl_msg_debug
cpl_msg_debug(cpl_func, "Spectra has <50 pixels -> don't flat")
gravi_envelope
double gravi_envelope(double lambda, double delta_lambda, double opd)
Definition gravi_old.c:4470
gravi_eop_load_finals2000A
gravi_data * gravi_eop_load_finals2000A(const char *eop_file)
Definition gravi_old.c:1876
flat_profiled
cpl_image * flat_profiled
Definition gravi_old.c:1959
header
cpl_propertylist * header
Definition gravi_old.c:2004
plist
cpl_propertylist * plist
Definition gravi_old.c:2000
gravi_get_base_from_tel
int gravi_get_base_from_tel(int t1, int t2)
Definition gravi_old.c:29
gravi_phase_correct_closures
cpl_error_code gravi_phase_correct_closures(cpl_table *phase_table)
Definition gravi_old.c:4238
gravi_array_get_group_delay_iota
double gravi_array_get_group_delay_iota(cpl_array *input, cpl_array *wavenumber)
Definition gravi_old.c:4870
gravi_normalize_flux
cpl_error_code gravi_normalize_flux(gravi_data *data)
Definition gravi_old.c:1543
gravi_fit_argon
gravi_data * gravi_fit_argon(gravi_data *argon_data, gravi_data *wave_data, gravi_data *profile_map, gravi_data *dark_map, gravi_data *bad_map)
This function calibrate the wavelength using the argon.
Definition gravi_old.c:430
cpl_msg_info
cpl_msg_info(cpl_func, "Compute WAVE_SCAN for %s", GRAVI_TYPE(type_data))
gravi_construction_opd_phase
cpl_vector * gravi_construction_opd_phase(cpl_table *opl_table, cpl_table *phase_sc, cpl_table *phase_ft, double dit_sc)
Definition gravi_old.c:2254
gravi_vector_get_med_percent
double gravi_vector_get_med_percent(cpl_vector *vector_in, float percent)
Definition gravi_old.c:4978
gravi_write_yorick_batch
int gravi_write_yorick_batch(const char *filename, const char *input_file, const char *output_file)
Create a yorick batch to execute mira on a OIFITS file.
Definition gravi_old.c:1802
gravi_phase_correct_closures_new
cpl_error_code gravi_phase_correct_closures_new(cpl_table *phase_table)
Definition gravi_old.c:3716
gravi_pfits_get_lambdamet
double gravi_pfits_get_lambdamet(const cpl_propertylist *plist)
Definition gravi_old.c:2152
cpl_propertylist_update_double
cpl_propertylist_update_double(header, "ESO QC MINWAVE SC", cpl_propertylist_get_double(plist, "ESO QC MINWAVE SC"))
gravi_compute_wave_offset
gravi_data * gravi_compute_wave_offset(gravi_data *argon_wave, gravi_data *wave_data)
Definition gravi_old.c:947
gravi_image_to_vector
cpl_vector * gravi_image_to_vector(cpl_image *img)
Definition gravi_old.c:1582
my_gsl_matrix_complex_fprintf
int my_gsl_matrix_complex_fprintf(FILE *stream, gsl_matrix_complex *m, char *fmt)
Definition gravi_old.c:1638
gravi_array_get_group_delay_coarse
double gravi_array_get_group_delay_coarse(cpl_array *input, cpl_array *sigma)
Definition gravi_old.c:4826
gravi_vis_create_phasesc_ft
gravi_vis_create_phasesc_ft(phase_ft, phase_sc, dit_sc)
wave_data_sc
cpl_table * wave_data_sc
Definition gravi_old.c:1997
cpl_image_delete
cpl_image_delete(flat_profiled)
gravi_smooth_preproc
cpl_error_code gravi_smooth_preproc(gravi_data *data, const cpl_parameterlist *parlist)
This function smooth the SC preproc data in the spectral direction in place, only for preproc (per re...
Definition gravi_old.c:2172
gravi_array_get_group_delay_loop_new
cpl_error_code gravi_array_get_group_delay_loop_new(cpl_array **input, cpl_array *sigma, double *gd, cpl_size nrow)
Definition gravi_old.c:4660
gravi_vis_create_met_ft
gravi_vis_create_met_ft(phase_ft, vis_met)
GRAVI_INSNAME
#define GRAVI_INSNAME(type, pol, npol)
Definition gravi_pfits.h:198
GRAVI_OI_VIS_MET_EXT
#define GRAVI_OI_VIS_MET_EXT
Definition gravi_pfits.h:88
GRAVI_IMAGING_DETECTOR_FT_EXT
#define GRAVI_IMAGING_DETECTOR_FT_EXT
Definition gravi_pfits.h:50
INSNAME_SC_P2
#define INSNAME_SC_P2
Definition gravi_pfits.h:195
GRAVI_SC
#define GRAVI_SC
Definition gravi_pfits.h:165
GRAVI_WAVE_FIBRE_SC_EXT
#define GRAVI_WAVE_FIBRE_SC_EXT
Definition gravi_pfits.h:69
QC_DARKRMS_SC
#define QC_DARKRMS_SC
Definition gravi_pfits.h:119
GRAVI_PROFILE_DATA_EXT
#define GRAVI_PROFILE_DATA_EXT
Definition gravi_pfits.h:79
GRAVI_WAVE_DATA_FT_EXT
#define GRAVI_WAVE_DATA_FT_EXT
Definition gravi_pfits.h:64
GRAVI_FDDL_EXT
#define GRAVI_FDDL_EXT
Definition gravi_pfits.h:75
GRAVI_WAVE_DATA_SC_EXT
#define GRAVI_WAVE_DATA_SC_EXT
Definition gravi_pfits.h:65
INSNAME_SC
#define INSNAME_SC
Definition gravi_pfits.h:197
INSNAME_FT_P1
#define INSNAME_FT_P1
Definition gravi_pfits.h:192
GRAVI_EXTVER
#define GRAVI_EXTVER(type, pol, npol)
Definition gravi_pfits.h:204
GRAVI_WAVE_ARGON_EXT
#define GRAVI_WAVE_ARGON_EXT
Definition gravi_pfits.h:63
GRAVI_IMAGING_DATA_SC_EXT
#define GRAVI_IMAGING_DATA_SC_EXT
Definition gravi_pfits.h:44
GRAVI_SPECTRUM_DATA_FT_EXT
#define GRAVI_SPECTRUM_DATA_FT_EXT
Definition gravi_pfits.h:54
GRAVI_TYPE
#define GRAVI_TYPE(type)
Definition gravi_pfits.h:167
INSNAME_SC_P1
#define INSNAME_SC_P1
Definition gravi_pfits.h:194
INSNAME_FT_P2
#define INSNAME_FT_P2
Definition gravi_pfits.h:193
INSNAME_FT
#define INSNAME_FT
Definition gravi_pfits.h:196
LAMBDA_MET
#define LAMBDA_MET
Definition gravi_pfits.h:103
GRAVI_IMAGING_DETECTOR_SC_EXT
#define GRAVI_IMAGING_DETECTOR_SC_EXT
Definition gravi_pfits.h:49
GRAVI_FT
#define GRAVI_FT
Definition gravi_pfits.h:166
GRAVI_SPECTRUM_DATA_SC_EXT
#define GRAVI_SPECTRUM_DATA_SC_EXT
Definition gravi_pfits.h:53
CPLCHECK_INT
#define CPLCHECK_INT(msg)
Definition gravi_utils.h:51
CPLCHECK_CLEAN
#define CPLCHECK_CLEAN(msg)
Definition gravi_utils.h:54
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
POLAR_1
#define POLAR_1
Definition gravi_utils.h:141
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
POLAR_2
#define POLAR_2
Definition gravi_utils.h:142
POLAR_3
#define POLAR_3
Definition gravi_utils.h:143
FREELOOP
#define FREELOOP(function, variable, n)
Definition gravi_utils.h:72
CPLCHECK
#define CPLCHECK(msg)
Definition gravi_utils.h:42
INFO_DEBUG
#define INFO_DEBUG
Definition gravi_calib.c:47
gravi_array_smooth
cpl_array * gravi_array_smooth(cpl_array *input_array, int DIT_smooth)
Definition gravi_cpl.c:1284
gravi_vector_extract
cpl_vector * gravi_vector_extract(const cpl_vector *vector, int start, int step)
Extract part of a vector.
Definition gravi_cpl.c:2695
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_get_data_array_double
double ** gravi_table_get_data_array_double(cpl_table *table, const char *name)
Definition gravi_cpl.c:472
gravi_table_get_vector
cpl_vector * gravi_table_get_vector(cpl_table *spectrum_data, cpl_size index, const char *regname)
Create a vector from the row index of the column regname.
Definition gravi_cpl.c:2234
gravi_array_multiply_phasor
cpl_error_code gravi_array_multiply_phasor(cpl_array *input, double complex factor, cpl_array *phase)
Multiply a REAL phase to a COMPLEX array, in-place: input = input * cexp (factor * phase)
Definition gravi_cpl.c:1158
gravi_vector_get_maxpos
cpl_size gravi_vector_get_maxpos(cpl_vector *vector)
Return the index of maximum in a vector. If several indexes exists with the maximum value,...
Definition gravi_cpl.c:2613
gravi_table_new_column_array
cpl_error_code gravi_table_new_column_array(cpl_table *table, const char *name, const char *unit, cpl_type type, cpl_size size)
Definition gravi_cpl.c:1677
gravi_data_get_oi_plist
cpl_propertylist * gravi_data_get_oi_plist(gravi_data *self, const char *extname, const char *insname)
Get the propertylist from EXTNAME and INSNAME.
Definition gravi_data.c:2001
gravi_data_add_cube
cpl_error_code gravi_data_add_cube(gravi_data *self, cpl_propertylist *plist, const char *extname, cpl_imagelist *imglist)
Add an IMAGE (imagelist) extension in gravi_data.
Definition gravi_data.c:2353
gravi_data_get_plist
cpl_propertylist * gravi_data_get_plist(gravi_data *self, const char *extname)
Get the propertylist from EXTNAME.
Definition gravi_data.c:2049
gravi_data_new
gravi_data * gravi_data_new(int nb_ext)
Create an empty gravi_data.
Definition gravi_data.c:110
gravi_data_get_cube
cpl_imagelist * gravi_data_get_cube(gravi_data *self, const char *extname)
Return a pointer on an IMAGE extension by its EXTNAME.
Definition gravi_data.c:2131
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_oi_table
cpl_table * gravi_data_get_oi_table(gravi_data *self, const char *extname, const char *insname)
Get an OI_FITS table from EXTNAME and INSNAME.
Definition gravi_data.c:1952
gravi_data_get_plist_x
cpl_propertylist * gravi_data_get_plist_x(gravi_data *self, int i)
Get the propertylist of an extension by position.
Definition gravi_data.c:1876
gravi_data_get_size
int gravi_data_get_size(const gravi_data *self)
Get the number of extension in a gravi_data.
Definition gravi_data.c:828
gravi_data_delete
void gravi_data_delete(gravi_data *self)
Delete a gravi data.
Definition gravi_data.c:146
gravi_compute_disp
gravi_data * gravi_compute_disp(gravi_data *vis_data)
Compute the dispersion calibration map (DISP) FIXME: to be done.
Definition gravi_old.c:3004
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
dfda_sin
static int dfda_sin(const double x[], const double a[], double result[])
Definition gravi_metrology.c:3804
gravi_pfits_get_startx
int gravi_pfits_get_startx(const cpl_propertylist *plist)
find out the name of the propertylist
Definition gravi_pfits.c:68
gravi_pfits_get_insname
const char * gravi_pfits_get_insname(const cpl_propertylist *plist)
Definition gravi_pfits.c:937
gravi_pfits_get_pola_num
int gravi_pfits_get_pola_num(const cpl_propertylist *plist, int type_data)
Definition gravi_pfits.c:263
gravi_pfits_get_extension_type
int gravi_pfits_get_extension_type(const cpl_propertylist *plist)
Definition gravi_pfits.c:271
gravi_pfits_get_spec_res
const char * gravi_pfits_get_spec_res(const cpl_propertylist *plist)
Definition gravi_pfits.c:162
gravi_pfits_get_resolution
const char * gravi_pfits_get_resolution(const cpl_propertylist *plist)
Definition gravi_pfits.c:155
gravi_pfits_get_extname
const char * gravi_pfits_get_extname(const cpl_propertylist *plist)
Definition gravi_pfits.c:140
gravi_pfits_get_window_start
int gravi_pfits_get_window_start(const cpl_propertylist *plist)
Definition gravi_pfits.c:92
gravi_extract_spectrum
gravi_data * gravi_extract_spectrum(gravi_data *raw_data, gravi_data *profile_map, gravi_data *dark_map, gravi_data *bad_map, gravi_data *sky_map, const cpl_parameterlist *parlist, enum gravi_detector_type det_type)
Create the SPECTRUM gravi_data with extracted spectrum per region.
Definition gravi_preproc.c:612
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_TRI_SIGN
int GRAVI_TRI_SIGN[GRAVI_NBASE][2][2]
Definition gravi_utils.c:61
gravi_region_get_base
int gravi_region_get_base(cpl_table *imaging_detector, int region)
Return the base of a region.
Definition gravi_utils.c:361
gravi_spectrum_get_nwave
cpl_size gravi_spectrum_get_nwave(const cpl_table *table)
Definition gravi_utils.c:1013
GRAVI_DATAERR
char GRAVI_DATAERR[50][10]
Definition gravi_utils.c:79
GRAVI_DATA
char GRAVI_DATA[50][7]
Definition gravi_utils.c:71
gravi_spectrum_get_nregion
cpl_size gravi_spectrum_get_nregion(const cpl_table *table)
Definition gravi_utils.c:1018
gravi_compute_envelope
cpl_vector * gravi_compute_envelope(const cpl_vector *opd, int wave, int nwave)
Compute the envelope value.
Definition gravi_utils.c:1081
GRAVI_BASE_TEL
int GRAVI_BASE_TEL[GRAVI_NBASE][2]
Definition gravi_utils.c:56
gravi_get_region
int gravi_get_region(cpl_table *img_det, int base, char phase, int pol)
Find the region matching base, phase and pol.
Definition gravi_utils.c:584
gravi_region_get_base_sign
int gravi_region_get_base_sign(cpl_table *imaging_detector, int base)
Return the sign of a base by looking at the PORT order.
Definition gravi_utils.c:399
GRAVI_TRI_BASE
int GRAVI_TRI_BASE[GRAVI_NBASE][2][2]
Definition gravi_utils.c:60
gravi_compute_argon_wave
cpl_table * gravi_compute_argon_wave(cpl_table *spectrum_table)
Compute a WAVE calibration from the ARGON data (SC only)
Definition gravi_wave.c:148
Generated on Fri Apr 24 2026 17:14:45 for GRAVI Pipeline Reference Manual by doxygen 1.13.2