CR2RE Pipeline Reference Manual 1.6.11
cr2res_etalon.c
1/*
2 * This file is part of the CR2RES Pipeline
3 * Copyright (C) 2002,2003 European Southern Observatory
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 2 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02111-1307 USA
18 */
19
20#ifdef HAVE_CONFIG_H
21#include <config.h>
22#endif
23
24/*-----------------------------------------------------------------------------
25 Includes
26 -----------------------------------------------------------------------------*/
27
28#include <cpl.h>
29#include <math.h>
30#include "cr2res_dfs.h"
31#include "cr2res_etalon.h"
32#include "cr2res_wave.h"
33#include "cr2res_pfits.h"
34
35/*-----------------------------------------------------------------------------
36 Defines
37 -----------------------------------------------------------------------------*/
38
39#define SPEED_OF_LIGHT 299792.458
40#define modulo(x, y) ((x) - ((y) * trunc((x)/(y))))
41#define min(a,b) (((a)<(b))?(a):(b))
42#define max(a,b) (((a)>(b))?(a):(b))
43
44#define CR2RES_ETALON_GAP_H 9.9950e6
45#define CR2RES_ETALON_GAP_J 9.9900e6
46#define CR2RES_ETALON_GAP_K 9.9999e6
47#define CR2RES_ETALON_GAP_Y 9.9913e6
48
49#define CR2RES_ETALON_DEGREE_H 2
50#define CR2RES_ETALON_DEGREE_J 1
51#define CR2RES_ETALON_DEGREE_K 1
52#define CR2RES_ETALON_DEGREE_Y 2
53
54/*-----------------------------------------------------------------------------
55 Functions prototypes
56 -----------------------------------------------------------------------------*/
57
58#ifdef CR2RES_UNUSED
59static cpl_mask * cr2res_etalon_binary_image(const cpl_image *) ;
60#endif
61static cpl_vector * cr2res_etalon_get_peaks_gaussian(
62 cpl_bivector * spectra,
63 cpl_bivector * spectra_err,
64 cpl_polynomial * wavesol_init,
65 cpl_array * wavesol_init_err,
66 cpl_vector * peaks,
67 int display,
68 cpl_vector ** sigma,
69 cpl_vector ** heights,
70 cpl_vector ** fit_error);
71
72/*----------------------------------------------------------------------------*/
76/*----------------------------------------------------------------------------*/
77
79#ifdef CR2RES_UNUSED
80/*----------------------------------------------------------------------------*/
86/*----------------------------------------------------------------------------*/
87cpl_image * cr2res_etalon_computation(const cpl_image * in)
88{
89 cpl_mask * mask ;
90 cpl_size nlabels ;
91 cpl_image * labels ;
92 cpl_apertures * aperts ;
93
94 /* Compute Binary image */
95 if ((mask = cr2res_etalon_binary_image(in)) == NULL) {
96 cpl_msg_error(__func__, "Cannot Binarise the image") ;
97 return NULL ;
98 }
99
100 /* Debugging */
101 cpl_mask_save(mask, "mask.fits", NULL, CPL_IO_CREATE) ;
102
103 /* Labelise the different detected apertures */
104 if ((labels = cpl_image_labelise_mask_create(mask, &nlabels))==NULL) {
105 cpl_msg_error(cpl_func, "Cannot Labelise") ;
106 cpl_mask_delete(mask) ;
107 return NULL ;
108 }
109 cpl_mask_delete(mask) ;
110
111 cpl_msg_debug(__func__, "Number of Apertures: %"CPL_SIZE_FORMAT, nlabels) ;
112
113 /* Create the detected apertures list */
114 if ((aperts = cpl_apertures_new_from_image(in, labels)) == NULL) {
115 cpl_msg_error(cpl_func, "Cannot Compute the apertures") ;
116 cpl_image_delete(labels) ;
117 return NULL ;
118 }
119
120 /* Display Apertures */
121 cpl_apertures_dump(aperts, stdout) ;
122
123 /* Free and return */
124 cpl_apertures_delete(aperts) ;
125 return labels ;
126}
127
129
130/*----------------------------------------------------------------------------*/
136/*----------------------------------------------------------------------------*/
137static cpl_mask * cr2res_etalon_binary_image(const cpl_image * in)
138{
139 cpl_mask * mask ;
140 cpl_vector * collapsed_vec ;
141 cpl_image * collapsed_im ;
142 cpl_vector * local_maxima ;
143 double * pmax ;
144 int i;
145
146 /* Check Entries */
147 if (in == NULL) return NULL;
148
149 /*************************/
150 /* Get Local Maxima list */
151 /*************************/
152 /* Collapse in y */
153 collapsed_im = cpl_image_collapse_create(in, 0) ;
154
155 /* Smooth */
156 /* TODO */
157
158 /* Convert as vector */
159 collapsed_vec = cpl_vector_new_from_image_row(collapsed_im, 1) ;
160 cpl_image_delete(collapsed_im) ;
161
162 /* Debugging */
163 cpl_plot_vector(NULL, "w lines", NULL, collapsed_vec) ;
164
165 /* Store local maxima positions of the plot */
166 local_maxima = cr2res_etalon_get_maxpos(collapsed_vec) ;
167 pmax = cpl_vector_get_data(local_maxima) ;
168 cpl_vector_delete(collapsed_vec) ;
169
170 /* Create the output mask */
171 mask = cpl_mask_new(cpl_image_get_size_x(in), cpl_image_get_size_y(in)) ;
172
173 /* Loop on the Maxima positions and isolate the blob */
174 for (i = 0; i < cpl_vector_get_size(local_maxima); i++) {
175 cpl_image *blob_image;
176 cpl_mask *blob_mask;
177 double threshold;
178 int j, start_x, end_x;
179
180 /* Get start_x end_x */
181 start_x = 1;
182 end_x = cpl_image_get_size_x(in) ;
183 if (i>0)
184 start_x = (int)(pmax[i-1] + (pmax[i]-pmax[i-1])/2.) ;
185 if (i<cpl_vector_get_size(local_maxima)-1)
186 end_x = (int)(pmax[i] + (pmax[i+1]-pmax[i])/2.) ;
187
188 cpl_msg_debug(__func__, "Extract %d -> %d", start_x, end_x) ;
189
190 /* Extract the current blob */
191 blob_image = cpl_image_extract(in, start_x, 1, end_x,
192 cpl_image_get_size_y(in)) ;
193
194 /* Compute the threshold */
195 /* TODO */
196 threshold = cpl_image_get_mean(blob_image) ;
197
198 cpl_msg_debug(__func__, "Threshold: %g", threshold) ;
199
200 /* Binarise the image */
201 blob_mask = cpl_mask_threshold_image_create(blob_image, threshold,
202 DBL_MAX) ;
203 cpl_image_delete(blob_image);
204
205 /* Set the left column to 0 to separate from the neighbor */
206 for (j=0 ; j<cpl_mask_get_size_y(blob_mask) ; j++)
207 cpl_mask_set(blob_mask, 1, j+1, CPL_BINARY_0) ;
208
209 /* Fill the Binary with the current blob */
210 cpl_mask_copy(mask, blob_mask, start_x, 1);
211
212 cpl_mask_delete(blob_mask) ;
213 }
214 /* Free and return */
215 cpl_vector_delete(local_maxima) ;
216
217 return mask ;
218}
219
221/*----------------------------------------------------------------------------*/
227/*----------------------------------------------------------------------------*/
228cpl_vector * cr2res_etalon_get_maxpos(const cpl_vector * in)
229{
230 const double * pin ;
231 cpl_vector * maxima_pos ;
232 double * pmax ;
233 int i, nb_max ;
234
235 /* Check Entries */
236 if (in==NULL) return NULL ;
237
238 /* Initialise */
239 pin = cpl_vector_get_data_const(in) ;
240
241 /* Count the number of Max positions */
242 nb_max = 0 ;
243 for (i=1 ; i<cpl_vector_get_size(in)-1 ; i++) {
244 if (pin[i] > pin[i+1] && pin[i] > pin[i-1]) nb_max++ ;
245 }
246
247 /* Create the output vector */
248 maxima_pos = cpl_vector_new(nb_max) ;
249 pmax = cpl_vector_get_data(maxima_pos) ;
250 nb_max = 0 ;
251 for (i=1 ; i<cpl_vector_get_size(in)-1 ; i++) {
252 if (pin[i] > pin[i+1] && pin[i] > pin[i-1]) {
253 /* Store the Max X position */
254 pmax[nb_max] = i+1 ;
255 nb_max++ ;
256 }
257 }
258 return maxima_pos ;
259}
260
261#endif
262/*----------------------------------------------------------------------------*/
277/*----------------------------------------------------------------------------*/
278cpl_vector * cr2res_etalon_get_local_maxima(
279 const cpl_vector * in,
280 cpl_vector ** left_edges,
281 cpl_vector ** right_edges){
282 /*
283 Find local maxima in a 1D array.
284 This function finds all local maxima in a 1D array and returns the indices
285 for their edges and midpoints (rounded down for even plateau sizes).
286 Parameters
287 ----------
288 x : ndarray
289 The array to search for local maxima.
290 Returns
291 -------
292 midpoints : ndarray
293 Indices of midpoints of local maxima in `x`.
294 left_edges : ndarray
295 Indices of edges to the left of local maxima in `x`.
296 right_edges : ndarray
297 Indices of edges to the right of local maxima in `x`.
298 Notes
299 -----
300 - Compared to `argrelmax` this function is significantly faster and can
301 detect maxima that are more than one sample wide. However this comes at
302 the cost of being only applicable to 1D arrays.
303 - A maxima is defined as one or more samples of equal value that are
304 surrounded on both sides by at least one smaller sample.
305 .. versionadded:: 1.1.0
306 */
307
308 cpl_vector * midpoints;
309 int m, i, i_ahead, i_max;
310
311 int width = cpl_vector_get_size(in);
312
313 // Preallocate, there can't be more maxima than half the size of `x`
314 midpoints = cpl_vector_new(width / 2);
315 *left_edges = cpl_vector_new(width / 2);
316 *right_edges = cpl_vector_new(width / 2);
317 m = 0; // Pointer to the end of valid area in allocated arrays
318
319 i = 1; // Pointer to current sample, first one can't be maxima
320 i_max = width - 1; // Last sample can't be maxima
321 while (i < i_max){
322 // Test if previous sample is smaller
323 if (cpl_vector_get(in, i - 1) < cpl_vector_get(in, i)){
324 i_ahead = i + 1; // Index to look ahead of current sample
325
326 // Find next sample that is unequal to x[i]
327 while (i_ahead < i_max && cpl_vector_get(in, i_ahead) == cpl_vector_get(in, i)){
328 i_ahead += 1;
329 }
330
331 // Maxima is found if next unequal sample is smaller than x[i]
332 if (cpl_vector_get(in, i_ahead) < cpl_vector_get(in, i)){
333 cpl_vector_set(*left_edges, m, i);
334 cpl_vector_set(*right_edges, m, i_ahead - 1);
335 cpl_vector_set(midpoints, m, (i + i_ahead - 1) / 2);
336 m += 1;
337 // Skip samples that can't be maximum
338 i = i_ahead;
339 }
340 }
341 i++;
342 }
343
344 if (m <1) {
345 cpl_msg_warning(__func__, "Cannot find maxima");
346 cpl_vector_delete(midpoints);
347 cpl_vector_delete(*left_edges);
348 cpl_vector_delete(*right_edges);
349 return NULL;
350 }
351 // Keep only valid part of array memory.
352 cpl_vector_set_size(midpoints, m);
353
354 cpl_vector_set_size(*left_edges, m);
355 cpl_vector_set_size(*right_edges, m);
356
357 return midpoints;
358}
359
360/*----------------------------------------------------------------------------*/
368/*----------------------------------------------------------------------------*/
369cpl_vector * cr2res_etalon_select_by_peak_distance(const cpl_vector * peaks,
370 const cpl_vector * peak_heights,
371 float distance)
372{
373 /*
374 Evaluate which peaks fulfill the distance condition.
375 Parameters
376 ----------
377 peaks : ndarray
378 Indices of peaks in `vector`.
379 priority : ndarray
380 An array matching `peaks` used to determine priority of each peak. A
381 peak with a higher priority value is kept over one with a lower one.
382 distance : np.float64
383 Minimal distance that peaks must be spaced.
384 Returns
385 -------
386 keep : ndarray[bool]
387 A boolean mask evaluating to true where `peaks` fulfill the distance
388 condition.
389 Notes
390 -----
391 Declaring the input arrays as C-contiguous doesn't seem to have performance
392 advantages.
393 .. versionadded:: 1.1.0
394 */
395
396
397 cpl_vector * keep;
398 cpl_vector * priority_to_position;
399 cpl_bivector * bivector_tmp;
400 int i, k, peaks_size, distance_;
401
402 peaks_size = cpl_vector_get_size(peaks);
403 // Round up because actual peak distance can only be natural number
404 distance_ = distance;
405 // Prepare array of flags
406 keep = cpl_vector_new(peaks_size);
407 cpl_vector_fill(keep, 1);
408
409 // Create map from `i` (index for `peaks` sorted by `priority`) to `j` (index
410 // for `peaks` sorted by position). This allows to iterate `peaks` and `keep`
411 // with `j` by order of `priority` while still maintaining the ability to
412 // step to neighbouring peaks with (`j` + 1) or (`j` - 1).
413
414 bivector_tmp = cpl_bivector_new(peaks_size);
415 priority_to_position = cpl_bivector_get_x(bivector_tmp);
416
417 for (i = 0; i < peaks_size; i++){
418 cpl_vector_set(priority_to_position, i, i);
419 }
420 cpl_vector_copy(cpl_bivector_get_y(bivector_tmp), peak_heights);
421
422 cpl_bivector_sort(bivector_tmp, bivector_tmp, CPL_SORT_ASCENDING, CPL_SORT_BY_Y);
423
424
425 // Highest priority first -> iterate in reverse order (decreasing)
426 for (i = peaks_size - 1; i> -1; i--){
427 int j;
428 // "Translate" `i` to `j` which points to current peak whose
429 // neighbours are to be evaluated
430 j = cpl_vector_get(priority_to_position, i);
431 if (cpl_vector_get(keep, j) == 0){
432 // Skip evaluation for peak already marked as "don't keep"
433 continue;
434 }
435
436 k = j - 1;
437 // Flag "earlier" peaks for removal until minimal distance is exceeded
438 while (0 <= k && cpl_vector_get(peaks, j) - cpl_vector_get(peaks,k) < distance_){
439 cpl_vector_set(keep, k, 0);
440 k--;
441 }
442
443 k = j + 1;
444 // Flag "later" peaks for removal until minimal distance is exceeded
445 while (k < peaks_size && cpl_vector_get(peaks, k) - cpl_vector_get(peaks, j) < distance_){
446 cpl_vector_set(keep, k, 0);
447 k++;
448 }
449 }
450 cpl_bivector_delete(bivector_tmp);
451
452 return keep;
453}
454
455/*----------------------------------------------------------------------------*/
465/*----------------------------------------------------------------------------*/
466cpl_vector * cr2res_etalon_find_peaks(
467 const cpl_vector * in,
468 double height,
469 double distance){
470
471 cpl_vector * right_edges;
472 cpl_vector * left_edges;
473 cpl_vector * peaks;
474 cpl_vector * peak_heights;
475 cpl_vector * peaks_out;
476 cpl_vector * peak_distance;
477
478 double peak, peak_height ;
479 int npeaks;
480 int k;
481
482 peaks = cr2res_etalon_get_local_maxima(in, &left_edges, &right_edges);
483 if (peaks == NULL) return NULL;
484 npeaks = cpl_vector_get_size(peaks);
485
486 peak_heights = cpl_vector_new(npeaks);
487 for (cpl_size i = 0; i < npeaks; i++)
488 {
489 peak = cpl_vector_get(peaks, i);
490 peak_height = cpl_vector_get(in, peak);
491 cpl_vector_set(peak_heights, i, peak_height);
492 }
493
494 peak_distance = cr2res_etalon_select_by_peak_distance(peaks,
495 peak_heights, distance);
496
497 // Evaluate height condition
498 peaks_out = cpl_vector_new(npeaks);
499 k = 0;
500 for (cpl_size i = 0; i < npeaks; i++)
501 {
502 peak = cpl_vector_get(peaks, i);
503 peak_height = cpl_vector_get(in, peak);
504
505 if ((peak_height > height) && (cpl_vector_get(peak_distance, i))){
506 cpl_vector_set(peaks_out, k, peak);
507 k++;
508 }
509 }
510 cpl_vector_set_size(peaks_out, k);
511 cpl_vector_delete(peaks);
512 cpl_vector_delete(left_edges);
513 cpl_vector_delete(right_edges);
514 cpl_vector_delete(peak_heights);
515 cpl_vector_delete(peak_distance);
516
517 return peaks_out;
518}
519
520static cpl_vector * cr2res_etalon_get_peaks_gaussian(
521 cpl_bivector * spectra,
522 cpl_bivector * spectra_err,
523 cpl_polynomial * wavesol_init,
524 cpl_array * wavesol_init_err,
525 cpl_vector * peaks,
526 int display,
527 cpl_vector ** sigma,
528 cpl_vector ** heights,
529 cpl_vector ** fit_error)
530{
531 cpl_bivector * linelist;
532 cpl_matrix * px;
533 cpl_vector * py;
534 cpl_vector * sigma_loc;
535 cpl_vector * heights_loc;
536 cpl_vector * fit_error_loc;
537 cpl_vector * wavelength;
538 cpl_vector * peak_height;
539 cpl_vector * new_peaks;
540 const cpl_vector * in;
541 cpl_size j, npeaks;
542
543 in = cpl_bivector_get_y_const(spectra);
544 npeaks = cpl_vector_get_size(peaks);
545 linelist = cpl_bivector_new(npeaks);
546
547 wavelength = cpl_bivector_get_x(linelist);
548 peak_height = cpl_bivector_get_y(linelist);
549 for (j = 0; j < npeaks; j++)
550 {
551 double wave, height;
552 wave = cpl_polynomial_eval_1d(wavesol_init,
553 cpl_vector_get(peaks, j), NULL);
554 height = cpl_vector_get(in, (cpl_size) cpl_vector_get(peaks, j));
555 cpl_vector_set(wavelength, j, wave);
556 cpl_vector_set(peak_height, j, height);
557 }
558
559 if (cr2res_wave_extract_lines(spectra, spectra_err, wavesol_init,
560 wavesol_init_err, linelist, 40, 4, display, &px, &py, &sigma_loc,
561 &heights_loc, &fit_error_loc) == -1)
562 {
563 // Abort
564 cpl_msg_warning(__func__, "No lines could be found");
565 cpl_bivector_delete(linelist);
566 return NULL;
567 };
568
569 // Make px the new peaks
570 npeaks = cpl_matrix_get_nrow(px);
571 new_peaks = cpl_vector_wrap(npeaks, cpl_matrix_get_data(px));
572 cpl_matrix_unwrap(px);
573 // Delete all the other stuff we don't need
574 cpl_bivector_delete(linelist);
575 cpl_vector_delete(py);
576
577 if (sigma != NULL) *sigma = sigma_loc;
578 else cpl_vector_delete(sigma_loc);
579 if (heights != NULL) *heights = heights_loc;
580 else cpl_vector_delete(heights_loc);
581 if (fit_error_loc != NULL) *fit_error = fit_error_loc;
582 else cpl_vector_delete(fit_error_loc);
583
584 return new_peaks;
585}
586
587#define signum(x) ((x > 0) ? 1 : ((x < 0) ? -1 : 0))
588
589double rofunc(
590 double * x,
591 double * y,
592 int ndata,
593 double * a,
594 double * abdev,
595 const double b)
596{
597 // Evaluates the right-hand side of equation (15.7.16) for a given value of b.
598 int j;
599 double sum=0.0;
600 double arr[ndata];
601 cpl_vector * tmp_vec;
602
603 for (j=0;j<ndata;j++) arr[j]=y[j]-b*x[j];
604 tmp_vec = cpl_vector_wrap(ndata, arr);
605 *a = cpl_vector_get_median(tmp_vec);
606 cpl_vector_unwrap(tmp_vec);
607
608 *abdev=0.0;
609 for (j=0;j<ndata;j++)
610 {
611 double d;
612 d=y[j]-(b*x[j]+*a);
613 *abdev += *abdev + fabs(d);
614 if (y[j] != 0.0)
615 d /= fabs(y[j]);
616 if (fabs(d) > DBL_EPSILON)
617 sum += (d >= 0.0 ? x[j] : -x[j]);
618 }
619 return sum;
620}
621
622// This is taken from Numerical Recipes Chapter 15.7
623// with some minor changes to make it work in C instead of C++
624cpl_polynomial * cr2res_robust_polynomial_fit(cpl_matrix *px, cpl_vector *py)
625{
626 // Object for fitting a straight line y = a+bx to a set of points .xi ; yi /, by the criterion of least
627 // absolute deviations. Call the constructor to calculate the fit. The answers are then available as
628 // the variables a, b, and abdev (the mean absolute deviation of the points from the line).
629 // Constructor. Given a set of data points xx[0..ndata-1], yy[0..ndata-1], sets a, b, and abdev.
630
631 double a, b, abdev;
632 int j;
633 double b1, del, f1, sigb;
634 double sx=0.0, sy=0.0, sxy=0.0, sxx=0.0, chisq=0.0;
635 double *x, *y;
636 int ndata;
637 cpl_polynomial * poly;
638 cpl_size deg;
639
640 x = cpl_matrix_get_data(px);
641 y = cpl_vector_get_data(py);
642 ndata = cpl_vector_get_size(py);
643
644 if(ndata < 1)
645 return NULL;
646
647 // As a first guess for a and b, we will find the
648 // least-squares fitting line.
649 for (j=0;j<ndata;j++)
650 {
651 sx += x[j];
652 sy += y[j];
653 sxy += x[j] * y[j];
654 sxx += x[j] * x[j];
655 }
656 del = ndata * sxx - sx * sx;
657 a = (sxx * sy - sx * sxy) / del; // Least-squares solutions.
658 b = (ndata * sxy - sx * sy) / del;
659
660 for (j = 0; j < ndata; j++){
661 double temp;
662 temp = y[j] - (a + b * x[j]);
663 chisq += temp * temp;
664 }
665 sigb = sqrt(chisq / del);
666 // The standard deviation will give some idea of
667 // how big an iteration step to take.
668 b1 = b;
669 f1 = rofunc(x, y, ndata, &a, &abdev, b1);
670 if (sigb > 0.0)
671 {
672 double f2, b2;
673 //Guess bracket as 3-sigma away, in the downhill direction known from f1.
674 b2 = b + 3.0 * sigb * signum(f1);
675 f2 = rofunc(x, y, ndata, &a, &abdev, b2);
676 if (b2 == b1)
677 {
678 abdev /= ndata;
679 poly = cpl_polynomial_new(1);
680 deg = 0;
681 cpl_polynomial_set_coeff(poly, &deg, a);
682 deg = 1;
683 cpl_polynomial_set_coeff(poly, &deg, b);
684 return poly;
685 }
686 //Bracketing
687 while (f1 * f2 > 0.0)
688 {
689 b = b2 + 1.6 * (b2-b1);
690 b1 = b2;
691 f1 = f2;
692 b2 = b;
693 f2 = rofunc(x, y, ndata, &a, &abdev, b2);
694 }
695 sigb = 0.01 * sigb;
696 while (fabs(b2 - b1) > sigb)
697 {
698 double f;
699 b = b1 + 0.5 * (b2 - b1); //Bisection.
700 if (b == b1 || b == b2)
701 break;
702 f=rofunc(x, y, ndata, &a, &abdev, b);
703 if (f * f1 >= 0.0)
704 {
705 f1 = f;
706 b1 = b;
707 } else {
708 //f2 = f; Never used
709 b2 = b;
710 }
711 }
712 }
713 abdev /= ndata;
714
715 poly = cpl_polynomial_new(1);
716 deg = 0;
717 cpl_polynomial_set_coeff(poly, &deg, a);
718 deg = 1;
719 cpl_polynomial_set_coeff(poly, &deg, b);
720
721 return poly;
722}
723
724
725/*----------------------------------------------------------------------------*/
740/*----------------------------------------------------------------------------*/
741cpl_polynomial * cr2res_etalon_wave_2d(
742 cpl_bivector ** spectra,
743 cpl_bivector ** spectra_err,
744 cpl_polynomial ** wavesol_init,
745 cpl_array ** wavesol_init_err,
746 int * orders,
747 int * traces_nb,
748 int ninputs,
749 cpl_size degree_x,
750 cpl_size degree_y,
751 int zp_order,
752 int display,
753 cpl_array ** wavelength_error,
754 cpl_table ** line_diagnostics)
755{
756 const cpl_vector * in;
757 cpl_vector * peaks;
758 cpl_vector * peaks_new;
759 cpl_vector * freq;
760 cpl_vector * wcen;
761 cpl_vector * mpos;
762 cpl_polynomial * poly;
763 cpl_matrix * px;
764 cpl_vector * pxa;
765 cpl_vector * pxb;
766 cpl_vector * py;
767 cpl_vector ** heights;
768 cpl_vector ** sigmas;
769 cpl_vector ** fit_errors;
770 cpl_vector ** fpe_xobs;
771 cpl_vector ** fpe_wobs;
772 cpl_vector ** fpe_freq;
773 cpl_vector ** fpe_mord;
774 cpl_vector ** fpe_cord;
775 cpl_polynomial * result;
776 cpl_vector * pos;
777 cpl_vector * diff;
778
779 cpl_vector * tmp_vec;
780 cpl_vector * corr;
781 cpl_size degree_2d[2];
782 cpl_size i, j, k, deg, npeaks, npeaks_total, npoints, setting_deg;
783 double wave, gap, tmp;
784 double f0, fr, m;
785 char * path;
786 cpl_table * lines_diagnostics_loc;
787 double offset;
788 int pad;
789
790 /* Check Inputs */
791 if (spectra==NULL || spectra_err==NULL || wavesol_init==NULL ||
792 orders==NULL)
793 return NULL ;
794 for (i=0 ; i<ninputs ; i++) {
795 if (spectra[i]==NULL || spectra_err[i]==NULL || wavesol_init[i]==NULL)
796 return NULL ;
797 }
798 *wavelength_error = NULL;
799
800 // Determine the setting based on the initial wavelength guess
801 deg = 0;
802 tmp = cpl_polynomial_get_coeff(wavesol_init[0], &deg);
803 switch ((int)(tmp/100))
804 {
805 case 9:
806 case 10:
807 case 11:
808 gap = CR2RES_ETALON_GAP_Y;
809 setting_deg = CR2RES_ETALON_DEGREE_Y;
810 break;
811 case 12:
812 case 13:
813 case 14:
814 gap = CR2RES_ETALON_GAP_J;
815 setting_deg = CR2RES_ETALON_DEGREE_J;
816 break;
817 case 16:
818 case 17:
819 case 18:
820 gap = CR2RES_ETALON_GAP_H;
821 setting_deg = CR2RES_ETALON_DEGREE_H;
822 break;
823 case 23:
824 case 24:
825 case 25:
826 gap = CR2RES_ETALON_GAP_K;
827 setting_deg = CR2RES_ETALON_DEGREE_K;
828 break;
829 default:
830 gap = -1;
831 setting_deg = 1;
832 break;
833 }
834
835
836 fpe_xobs = cpl_malloc(ninputs * sizeof(cpl_vector *));
837 fpe_wobs = cpl_malloc(ninputs * sizeof(cpl_vector *));
838 fpe_freq = cpl_malloc(ninputs * sizeof(cpl_vector *));
839 fpe_mord = cpl_malloc(ninputs * sizeof(cpl_vector *));
840 fpe_cord = cpl_malloc(ninputs * sizeof(cpl_vector *));
841
842 heights = cpl_malloc(ninputs * sizeof(cpl_vector *));
843 sigmas = cpl_malloc(ninputs * sizeof(cpl_vector *));
844 fit_errors = cpl_malloc(ninputs * sizeof(cpl_vector *));
845
846
847 npeaks_total = 0;
848
849 for (i = 0; i < ninputs; i++){
850 if (spectra[i] == NULL){
851 // Skip this spectrum
852 fpe_xobs[i] = NULL;
853 fpe_wobs[i] = NULL;
854 fpe_freq[i] = NULL;
855 fpe_mord[i] = NULL;
856 fpe_cord[i] = NULL;
857 sigmas[i] = NULL;
858 heights[i] = NULL;
859 fit_errors[i] = NULL;
860 continue;
861 }
862 // Find peaks in the etalon spectra
863 in = cpl_bivector_get_y_const(spectra[i]);
864 if ((peaks = cr2res_etalon_find_peaks(in, cpl_vector_get_mean(in), 3))
865 == NULL){
866 fpe_xobs[i] = NULL;
867 fpe_wobs[i] = NULL;
868 fpe_freq[i] = NULL;
869 fpe_mord[i] = NULL;
870 fpe_cord[i] = NULL;
871 sigmas[i] = NULL;
872 heights[i] = NULL;
873 fit_errors[i] = NULL;
874 continue;
875 }
876
877 // get the peak using the gaussian fit
878 if ((peaks_new = cr2res_etalon_get_peaks_gaussian(spectra[i], spectra_err[i],
879 wavesol_init[i], wavesol_init_err[i], peaks, display,
880 &sigmas[i], &heights[i], &fit_errors[i])) == NULL){
881 cpl_vector_delete(peaks);
882 fpe_xobs[i] = NULL;
883 fpe_wobs[i] = NULL;
884 fpe_freq[i] = NULL;
885 fpe_mord[i] = NULL;
886 fpe_cord[i] = NULL;
887 sigmas[i] = NULL;
888 heights[i] = NULL;
889 fit_errors[i] = NULL;
890 continue;
891 }
892 cpl_vector_delete(peaks);
893 npeaks = cpl_vector_get_size(peaks_new);
894 npeaks_total += npeaks;
895
896 /*
897 Adjust frequency sequence so that freq_i = i * const
898 Then determine FPE order number m for each peak in a given CRIRES+ order i.
899 This is done with the following smart algebra based on diffraction equation:
900
901 m * lambda_m = const (1)
902 (m+1) * lambda_(m+1) = const
903
904 Subtracting:
905
906 m * [lambda_m - lambda_(m+1)] = lambda_(m+1)
907 m = lambda_(m+1) / [lambda_m - lambda_(m+1)] = lambda_(m+1) / delta lambda_m
908
909 To appreciate how smart this is note that const is not really constant but a slow
910 function of the wavelength. For the two adjacent lines it is very much the same
911 helping us avoiding 9th order polynomial fit to the const as e.g. in Cersullo et al.,
912 2019, A&A 624, 122.
913
914 We can reformulate the formula to frequencies by using
915 freq = c_light / lambda
916 freq_i = f0 + n_i * fr
917 where we can fit f0 and fr to the detected peaks, to:
918 m = abs(f0 / fr + n_i + 1)
919 */
920
921 // Create the frequencies
922 freq = cpl_vector_new(npeaks);
923 wcen = cpl_vector_new(npeaks);
924 for ( j = 0; j < npeaks; j++)
925 {
926 wave = cpl_polynomial_eval_1d(wavesol_init[i],
927 cpl_vector_get(peaks_new, j), NULL);
928 cpl_vector_set(freq, j, SPEED_OF_LIGHT / wave);
929 cpl_vector_set(wcen, j, wave);
930 }
931
932 // get first guesses for the polynomial parameters fr
933 tmp_vec = cpl_vector_new(npeaks - 1);
934 for ( j = 0; j < npeaks - 1; j++)
935 {
936 cpl_vector_set(tmp_vec, j,
937 cpl_vector_get(freq, j + 1)
938 - cpl_vector_get(freq, j));
939 }
940 fr = fabs(cpl_vector_get_median(tmp_vec));
941 cpl_vector_delete(tmp_vec);
942
943 // and a first guess for fd
944 tmp_vec = cpl_vector_new(npeaks);
945 for (j = 0; j < npeaks; j++)
946 {
947 cpl_vector_set(tmp_vec, j,
948 modulo(cpl_vector_get(freq, j), fr));
949 }
950 f0 = cpl_vector_get_median(tmp_vec);
951 cpl_vector_delete(tmp_vec);
952
953 // Fit a 1d polynomial to the frequencies
954 // determine the peak numbers
955 px = cpl_matrix_new(1, npeaks);
956 for (j = 0; j < npeaks; j++)
957 {
958 m = (cpl_vector_get(freq, j) - f0) / fr;
959 cpl_matrix_set(px, 0, j, round(m));
960 }
961
962 deg = 1;
963 poly = cpl_polynomial_new(1);
964 cpl_polynomial_fit(poly, px, NULL, freq, NULL, CPL_FALSE, NULL, &deg);
965 // and evaluate it at each peak
966 for ( j = 0; j < npeaks; j++)
967 {
968 m = cpl_matrix_get(px, 0, j);
969 // actually the frequency, just reusing the parameter name
970 wave = cpl_polynomial_eval_1d(poly, m, NULL);
971 cpl_vector_set(freq, j, wave);
972 cpl_vector_set(wcen, j, SPEED_OF_LIGHT / wave);
973 }
974 deg = 0;
975 //f0 = cpl_polynomial_get_coeff(poly, &deg); Never used
976 deg = 1;
977 //fr = cpl_polynomial_get_coeff(poly, &deg); Never used
978 cpl_polynomial_delete(poly);
979 // Determine M
980 mpos = cpl_vector_new(npeaks);
981 for ( j = 0; j < npeaks; j++)
982 {
983 // This is really sensitive to the wavelength solution...
984 // TODO: m is awkwardly close to 0.5 before rounding...
985 m = cpl_matrix_get(px, 0, j);
986 cpl_vector_set(mpos, j, round(m));
987 // m = cpl_vector_get(wcen, j-1)/
988 // fabs(cpl_vector_get(wcen, j) - cpl_vector_get(wcen, j-1));
989 // cpl_vector_set(mpos, j, m);
990 }
991 // cpl_vector_set(mpos, 0, 1 + cpl_vector_get(mpos, 1));
992 cpl_matrix_delete(px);
993
994 //
995 tmp_vec = cpl_vector_new(npeaks);
996 for (j = 1; j < npeaks; j++){
997 m = cpl_vector_get(wcen, j-1)/
998 fabs(cpl_vector_get(wcen, j) - cpl_vector_get(wcen, j-1));
999 cpl_vector_set(tmp_vec, j, m - cpl_vector_get(mpos, j));
1000 }
1001 offset = cpl_vector_get_median(tmp_vec);
1002 cpl_vector_add_scalar(mpos, round(offset));
1003 cpl_vector_delete(tmp_vec);
1004
1005 // Store vectors for later
1006 fpe_xobs[i] = peaks_new;
1007 fpe_wobs[i] = wcen;
1008 fpe_freq[i] = freq;
1009 fpe_mord[i] = mpos;
1010 // This is only used for the debug output
1011 fpe_cord[i] = cpl_vector_new(npeaks);
1012 cpl_vector_fill(fpe_cord[i], orders[i] + zp_order);
1013 }
1014
1015 if (npeaks_total == 0){
1016 cpl_msg_error(__func__, "No peaks found for Etalon wavecal");
1017 for (i = 0; i < ninputs; i++)
1018 {
1019 if (fpe_mord[i] == NULL) continue;
1020 cpl_vector_delete(fpe_xobs[i]);
1021 cpl_vector_delete(fpe_wobs[i]);
1022 cpl_vector_delete(fpe_freq[i]);
1023 cpl_vector_delete(fpe_mord[i]);
1024 cpl_vector_delete(fpe_cord[i]);
1025 cpl_vector_delete(heights[i]);
1026 cpl_vector_delete(sigmas[i]);
1027 cpl_vector_delete(fit_errors[i]);
1028 }
1029 cpl_free(fpe_xobs);
1030 cpl_free(fpe_wobs);
1031 cpl_free(fpe_freq);
1032 cpl_free(fpe_mord);
1033 cpl_free(fpe_cord);
1034 cpl_free(heights);
1035 cpl_free(sigmas);
1036 cpl_free(fit_errors);
1037 return NULL;
1038 }
1039
1040 if (cpl_msg_get_level() == CPL_MSG_DEBUG){
1041 path = cpl_sprintf("debug_etalon_mord.fits");
1042 cpl_vector_save(fpe_mord[0], path, CPL_TYPE_DOUBLE, NULL, CPL_IO_CREATE);
1043 for (i = 1; i < ninputs; i++)
1044 {
1045 cpl_vector_save(fpe_mord[i], path, CPL_TYPE_DOUBLE, NULL, CPL_IO_EXTEND);
1046 }
1047 cpl_free(path);
1048 path = cpl_sprintf("debug_etalon_wobs.fits");
1049 cpl_vector_save(fpe_wobs[0], path, CPL_TYPE_DOUBLE, NULL, CPL_IO_CREATE);
1050 for (i = 1; i < ninputs; i++)
1051 {
1052 cpl_vector_save(fpe_wobs[i], path, CPL_TYPE_DOUBLE, NULL, CPL_IO_EXTEND);
1053 }
1054 cpl_free(path);
1055 }
1056
1057 /*
1058 If the initial wavelength solution has an offset larger than the FPE line
1059 spacing the m's defined above will not match across the CRIRES+ orders.
1060 Here with find the offsets using eq. 1 and apply them:
1061 */
1062 // gap here is the constant m * w
1063 if (gap < 0)
1064 {
1065 cpl_vector * fpe_gap;
1066 k = 0;
1067 fpe_gap = cpl_vector_new(npeaks_total);
1068 for (i = 0; i < ninputs; i++)
1069 {
1070 if (fpe_mord[i] == NULL) continue;
1071 for (j = 0; j < cpl_vector_get_size(fpe_mord[i]); j++)
1072 {
1073 cpl_vector_set(fpe_gap, k,
1074 cpl_vector_get(fpe_mord[i], j)
1075 * cpl_vector_get(fpe_wobs[i], j));
1076 k++;
1077 }
1078 }
1079 gap = cpl_vector_get_median(fpe_gap);
1080 cpl_vector_delete(fpe_gap);
1081 cpl_msg_debug(__func__, "Median gap is: %g", gap);
1082 }
1083
1084 // Calculate and apply the correction
1085 corr = cpl_vector_new(ninputs);
1086 for (i = 0; i < ninputs; i++)
1087 {
1088 if (fpe_mord[i] == NULL) continue;
1089 tmp_vec = cpl_vector_new(cpl_vector_get_size(fpe_mord[i]));
1090 for (j = 0; j < cpl_vector_get_size(fpe_mord[i]); j++)
1091 {
1092 m = cpl_vector_get(fpe_mord[i], j);
1093 wave = cpl_vector_get(fpe_wobs[i], j);
1094 tmp = gap / wave - m;
1095 cpl_vector_set(tmp_vec, j, tmp);
1096 }
1097 cpl_vector_set(corr, i, round(cpl_vector_get_median(tmp_vec)));
1098 cpl_vector_delete(tmp_vec);
1099 // Apply the correction
1100 cpl_vector_add_scalar(fpe_mord[i], cpl_vector_get(corr, i));
1101 }
1102 cpl_vector_delete(corr);
1103
1104 // In a second step we assume that m*wave varies linearly across all orders
1105 // and then correct for that
1106 npoints = 0;
1107 pad = min(2, ninputs-2);
1108 for (i = pad; i < ninputs - pad; i++)
1109 {
1110 if (fpe_mord[i] == NULL) continue;
1111 npoints += cpl_vector_get_size(fpe_mord[i]);
1112 }
1113
1114 py = cpl_vector_new(npoints);
1115 px = cpl_matrix_new(1, npoints);
1116 k = 0;
1117 for (i = pad; i < ninputs - pad; i++)
1118 {
1119 if (fpe_mord[i] == NULL) continue;
1120 for (j = 0; j < cpl_vector_get_size(fpe_mord[i]); j++)
1121 {
1122 wave = cpl_vector_get(fpe_wobs[i], j);
1123 m = cpl_vector_get(fpe_mord[i], j);
1124 cpl_matrix_set(px, 0, k, m);
1125 cpl_vector_set(py, k, m * wave);
1126 k++;
1127 }
1128 }
1129 poly = cpl_polynomial_new(1);
1130 deg = 1;
1131 cpl_polynomial_fit(poly, px, NULL, py, NULL, CPL_FALSE, NULL, &deg);
1132 cpl_matrix_delete(px);
1133 cpl_vector_delete(py);
1134
1135 // Calculate and apply the correction
1136 corr = cpl_vector_new(ninputs);
1137 for (i = 0; i < ninputs; i++)
1138 {
1139 if (fpe_mord[i] == NULL) continue;
1140 npeaks = cpl_vector_get_size(fpe_mord[i]);
1141 tmp_vec = cpl_vector_new(npeaks);
1142 for (j = 0; j < npeaks; j++)
1143 {
1144 m = cpl_vector_get(fpe_mord[i], j);
1145 wave = cpl_vector_get(fpe_wobs[i], j);
1146 gap = cpl_polynomial_eval_1d(poly, m, NULL);
1147 tmp = gap / wave - m;
1148 cpl_vector_set(tmp_vec, j, tmp);
1149 }
1150 cpl_vector_set(corr, i, round(cpl_vector_get_median(tmp_vec)));
1151 cpl_vector_delete(tmp_vec);
1152 // Apply the correction
1153 cpl_vector_add_scalar(fpe_mord[i], cpl_vector_get(corr, i));
1154 }
1155 cpl_vector_delete(corr);
1156 cpl_polynomial_delete(poly);
1157
1158 /*
1159 Fit m * wave of all orders with a linear fit, i.e. assuming it only varies
1160 slowly between orders. Then determine new wavelengths based on that.
1161 */
1162 py = cpl_vector_new(npeaks_total);
1163 px = cpl_matrix_new(1, npeaks_total);
1164 k = 0;
1165 for (i = 0; i < ninputs; i++)
1166 {
1167 if (fpe_mord[i] == NULL) continue;
1168 for (j = 0; j < cpl_vector_get_size(fpe_mord[i]); j++)
1169 {
1170 wave = cpl_vector_get(fpe_wobs[i], j);
1171 m = cpl_vector_get(fpe_mord[i], j);
1172 cpl_matrix_set(px, 0, k, m);
1173 cpl_vector_set(py, k, m * wave);
1174 k++;
1175 }
1176 }
1177
1178 if (setting_deg == 1){
1179
1180 cpl_matrix * px_tmp;
1181 cpl_vector * py_tmp;
1182
1183 double m0, ymax, ymin;
1184 // Do a robust (L1 norm) fit to the data
1185 // use robust to avoid outliers
1186 // However the method is not as numerically stable
1187 // so we normalize the x and y data first
1188 px_tmp = cpl_matrix_duplicate(px);
1189 py_tmp = cpl_vector_duplicate(py);
1190
1191 m0 = cpl_matrix_get_min(px_tmp);
1192 cpl_matrix_subtract_scalar(px_tmp, m0);
1193 ymin = cpl_vector_get_min(py_tmp);
1194 cpl_vector_subtract_scalar(py_tmp, ymin);
1195 ymax = cpl_vector_get_max(py_tmp);
1196 cpl_vector_divide_scalar(py_tmp, ymax);
1197 poly = cr2res_robust_polynomial_fit(px_tmp, py_tmp);
1198
1199 cpl_matrix_delete(px_tmp);
1200 cpl_vector_delete(py_tmp);
1201
1202 // Fix orders to be on the same linear fit
1203 k = 0;
1204 for (i = 0; i < ninputs; i++)
1205 {
1206 double tmp1, tmp2;
1207 if (fpe_mord[i] == NULL) continue;
1208 // mean(mord * wobs)
1209 tmp1 = 0;
1210 for (j = 0; j < cpl_vector_get_size(fpe_mord[i]); j++)
1211 {
1212 m = cpl_vector_get(fpe_mord[i], j);
1213 wave = cpl_vector_get(fpe_wobs[i], j);
1214 tmp1 += m * wave;
1215 }
1216 tmp1 /= cpl_vector_get_size(fpe_mord[i]);
1217
1218 // mean(poly(mord))
1219 tmp2 = 0;
1220 for (j = 0; j < cpl_vector_get_size(fpe_mord[i]); j++)
1221 {
1222 m = cpl_vector_get(fpe_mord[i], j);
1223 // Remember to use the normalization factors created above
1224 tmp2 += cpl_polynomial_eval_1d(poly, m - m0, NULL) * ymax + ymin;
1225 }
1226 tmp2 /= cpl_vector_get_size(fpe_mord[i]);
1227
1228 cpl_msg_debug(__func__,
1229 "Difference between measured and expected in order %lli: %f",
1230 i, fabs(tmp1 - tmp2));
1231
1232 for (j = 0; j < cpl_vector_get_size(fpe_mord[i]); j++)
1233 {
1234 wave = cpl_vector_get(fpe_wobs[i], j);
1235 m = cpl_vector_get(fpe_mord[i], j);
1236 wave += (tmp2 - tmp1) / m;
1237 cpl_vector_set(py, k, m * wave);
1238 cpl_vector_set(fpe_wobs[i], j, wave);
1239 cpl_vector_set(fpe_freq[i], j, SPEED_OF_LIGHT / wave);
1240 k++;
1241 }
1242 }
1243 cpl_polynomial_delete(poly);
1244 }
1245
1246 // Do yet another fit to the points
1247 // and recalculate the wavelength on now (hopefully) consistent data
1248 deg = setting_deg;
1249 poly = cpl_polynomial_new(1);
1250 cpl_polynomial_fit(poly, px, NULL, py, NULL, CPL_FALSE, NULL, &deg);
1251
1252 cpl_matrix_delete(px);
1253 cpl_vector_delete(py);
1254
1255 for (i = 0; i < ninputs; i++)
1256 {
1257 if (fpe_mord[i] == NULL) continue;
1258 for (j = 0; j < cpl_vector_get_size(fpe_mord[i]); j++)
1259 {
1260 m = cpl_vector_get(fpe_mord[i], j);
1261 wave = cpl_polynomial_eval_1d(poly, m, NULL) / m;
1262 cpl_vector_set(fpe_freq[i], j, SPEED_OF_LIGHT / wave);
1263 cpl_vector_set(fpe_wobs[i], j, wave);
1264 }
1265 }
1266 cpl_polynomial_delete(poly);
1267
1268
1269 // Do the 2d fit
1270 pxa = cpl_vector_new(npeaks_total);
1271 pxb = cpl_vector_new(npeaks_total);
1272 py = cpl_vector_new(npeaks_total);
1273 k = 0;
1274 for (i = 0; i < ninputs; i++)
1275 {
1276 if (fpe_mord[i] == NULL) continue;
1277 npeaks = cpl_vector_get_size(fpe_mord[i]);
1278 for (j = 0; j < npeaks; j++){
1279 cpl_vector_set(pxa, k, cpl_vector_get(fpe_xobs[i], j));
1280 cpl_vector_set(pxb, k, orders[i] + zp_order);
1281 cpl_vector_set(py, k, cpl_vector_get(fpe_wobs[i], j));
1282 k++;
1283 }
1284 }
1285
1286 degree_2d[0] = degree_x ;
1287 degree_2d[1] = degree_y ;
1288 result = cr2res_polyfit_2d(pxa, pxb, py, degree_2d);
1289
1290
1291 if (result == NULL){
1292 cpl_msg_error(__func__, "Error in Etalon polynomial fit");
1293 // cpl_polynomial_delete(result);
1294 cpl_vector_delete(pxa);
1295 cpl_vector_delete(pxb);
1296 cpl_vector_delete(py);
1297 for (i = 0; i < ninputs; i++)
1298 {
1299 if (fpe_mord[i] == NULL) continue;
1300 cpl_vector_delete(fpe_xobs[i]);
1301 cpl_vector_delete(fpe_wobs[i]);
1302 cpl_vector_delete(fpe_freq[i]);
1303 cpl_vector_delete(fpe_mord[i]);
1304 cpl_vector_delete(fpe_cord[i]);
1305 cpl_vector_delete(heights[i]);
1306 cpl_vector_delete(sigmas[i]);
1307 cpl_vector_delete(fit_errors[i]);
1308 }
1309 cpl_free(fpe_xobs);
1310 cpl_free(fpe_wobs);
1311 cpl_free(fpe_freq);
1312 cpl_free(fpe_mord);
1313 cpl_free(fpe_cord);
1314 cpl_free(heights);
1315 cpl_free(sigmas);
1316 cpl_free(fit_errors);
1317 return NULL;
1318 }
1319
1320 if (cpl_msg_get_level() == CPL_MSG_DEBUG){
1321 cpl_table * tmp_table;
1322 path = cpl_sprintf("debug_etalon_final_mord.fits");
1323 cpl_vector_save(fpe_mord[0], path, CPL_TYPE_DOUBLE, NULL, CPL_IO_CREATE);
1324 for (i = 1; i < ninputs; i++)
1325 {
1326 cpl_vector_save(fpe_mord[i], path, CPL_TYPE_DOUBLE, NULL, CPL_IO_EXTEND);
1327 }
1328 cpl_free(path);
1329 path = cpl_sprintf("debug_etalon_final_wobs.fits");
1330 cpl_vector_save(fpe_wobs[0], path, CPL_TYPE_DOUBLE, NULL, CPL_IO_CREATE);
1331 for (i = 1; i < ninputs; i++)
1332 {
1333 cpl_vector_save(fpe_wobs[i], path, CPL_TYPE_DOUBLE, NULL, CPL_IO_EXTEND);
1334 }
1335 cpl_free(path);
1336 path = cpl_sprintf("debug_etalon_final_xobs.fits");
1337 cpl_vector_save(fpe_xobs[0], path, CPL_TYPE_DOUBLE, NULL, CPL_IO_CREATE);
1338 for (i = 1; i < ninputs; i++)
1339 {
1340 cpl_vector_save(fpe_xobs[i], path, CPL_TYPE_DOUBLE, NULL, CPL_IO_EXTEND);
1341 }
1342 cpl_free(path);
1343 path = cpl_sprintf("debug_etalon_final_cord.fits");
1344 cpl_vector_save(fpe_cord[0], path, CPL_TYPE_DOUBLE, NULL, CPL_IO_CREATE);
1345 for (i = 1; i < ninputs; i++)
1346 {
1347 cpl_vector_save(fpe_cord[i], path, CPL_TYPE_DOUBLE, NULL, CPL_IO_EXTEND);
1348 }
1349 cpl_free(path);
1350
1351 tmp_table = cpl_table_new((degree_x + 1) * (degree_y + 1));
1352 cpl_table_new_column(tmp_table, "DEGREE_X", CPL_TYPE_INT);
1353 cpl_table_new_column(tmp_table, "DEGREE_Y", CPL_TYPE_INT);
1354 cpl_table_new_column(tmp_table, "COEFFICIENT", CPL_TYPE_DOUBLE);
1355 k = 0;
1356 for (i = 0; i <= degree_x; i++)
1357 {
1358 for (j = 0; j <= degree_y; j++)
1359 {
1360 degree_2d[0] = i;
1361 degree_2d[1] = j;
1362 tmp = cpl_polynomial_get_coeff(result, &degree_2d[0]);
1363 cpl_table_set_int(tmp_table, "DEGREE_X", k, i);
1364 cpl_table_set_int(tmp_table, "DEGREE_Y", k, j);
1365 cpl_table_set_double(tmp_table, "COEFFICIENT", k, tmp);
1366 k++;
1367 }
1368 }
1369 path = cpl_sprintf("debug_etalon_final_poly.fits");
1370 cpl_table_save(tmp_table, NULL, NULL, path, CPL_IO_CREATE);
1371 cpl_table_delete(tmp_table);
1372 cpl_free(path);
1373 }
1374
1375 /* Create / Fill / Merge the lines diagnostics table */
1376 if (line_diagnostics != NULL) {
1377 *line_diagnostics = NULL;
1378 cpl_msg_debug(__func__, "Number of lines: %"CPL_SIZE_FORMAT, npeaks_total);
1379 for (i = 0; i < ninputs; i++)
1380 {
1381 cpl_polynomial * wavesol_loc;
1382 if (fpe_mord[i] == NULL) continue;
1383 npeaks = cpl_vector_get_size(fpe_mord[i]);
1384 /* Create */
1385 lines_diagnostics_loc =
1386 cr2res_dfs_create_lines_diagnostics_table(npeaks) ;
1387 wavesol_loc = cr2res_wave_poly_2d_to_1d(result,
1388 orders[i] + zp_order);
1389 /* Fill */
1390 for (j=0 ; j < npeaks ; j++) {
1391 double pix_pos, lambda_cat, lambda_meas, line_width, line_intens, fit_error;
1392 pix_pos = cpl_vector_get(fpe_xobs[i], j);
1393 lambda_cat = cpl_vector_get(fpe_wobs[i], j) ;
1394 lambda_meas = cpl_polynomial_eval_1d(wavesol_loc, pix_pos,
1395 NULL) ;
1396 line_width = cpl_vector_get(sigmas[i], j) ;
1397 line_intens = cpl_vector_get(heights[i], j) ;
1398 fit_error = cpl_vector_get(fit_errors[i], j) ;
1399 m = cpl_vector_get(fpe_mord[i], j);
1400 cpl_table_set_int(lines_diagnostics_loc,
1401 CR2RES_COL_ORDER, j, orders[i]) ;
1402 cpl_table_set_int(lines_diagnostics_loc,
1403 CR2RES_COL_TRACENB, j, traces_nb[i]) ;
1404 cpl_table_set_double(lines_diagnostics_loc,
1405 CR2RES_COL_MEASURED_LAMBDA, j, lambda_meas) ;
1406 cpl_table_set_double(lines_diagnostics_loc,
1407 CR2RES_COL_CATALOG_LAMBDA, j, lambda_cat);
1408 cpl_table_set_double(lines_diagnostics_loc,
1409 CR2RES_COL_DELTA_LAMBDA, j, lambda_cat-lambda_meas);
1410 cpl_table_set_double(lines_diagnostics_loc,
1411 CR2RES_COL_MEASURED_PIXEL, j, pix_pos);
1412 cpl_table_set_double(lines_diagnostics_loc,
1413 CR2RES_COL_LINE_WIDTH, j, line_width) ;
1414 cpl_table_set_double(lines_diagnostics_loc,
1415 CR2RES_COL_FIT_QUALITY, j, fit_error) ;
1416 cpl_table_set_double(lines_diagnostics_loc,
1417 CR2RES_COL_INTENSITY, j, line_intens) ;
1418 cpl_table_set_double(lines_diagnostics_loc,
1419 CR2RES_COL_FPET_M, j, m) ;
1420 }
1421 cpl_polynomial_delete(wavesol_loc);
1422 /* Merge */
1423 if (*line_diagnostics == NULL) {
1424 *line_diagnostics = lines_diagnostics_loc ;
1425 lines_diagnostics_loc = NULL ;
1426 } else if (lines_diagnostics_loc != NULL) {
1427 /* Merge with previous */
1428 cpl_table_insert(*line_diagnostics, lines_diagnostics_loc,
1429 cpl_table_get_nrow(*line_diagnostics)) ;
1430 cpl_table_delete(lines_diagnostics_loc) ;
1431 }
1432 }
1433 }
1434 cpl_msg_debug(__func__,"%s",cpl_error_get_where());
1435 npeaks = cpl_vector_get_size(py);
1436 // Calculate absolute difference between polynomial and
1437 // catalog value for each line
1438 // use px and py, so that only good lines are used
1439 diff = cpl_vector_new(npeaks);
1440 pos = cpl_vector_new(2);
1441 for (i = 0; i < npeaks; i++){
1442 cpl_vector_set(pos, 0, cpl_vector_get(pxa, i));
1443 cpl_vector_set(pos, 1, cpl_vector_get(pxb, i));
1444 cpl_vector_set(diff, i, fabs(
1445 cpl_polynomial_eval(result, pos)
1446 - cpl_vector_get(py, i)));
1447 }
1448
1449 if (*wavelength_error == NULL)
1450 *wavelength_error = cpl_array_new(2, CPL_TYPE_DOUBLE);
1451 cpl_array_set_double(*wavelength_error, 0,
1452 cpl_vector_get_mean(diff));
1453 cpl_array_set_double(*wavelength_error, 1,
1454 cpl_vector_get_max(diff));
1455
1456 cpl_vector_delete(diff);
1457 cpl_vector_delete(pos);
1458
1459 cpl_vector_delete(pxa);
1460 cpl_vector_delete(pxb);
1461 cpl_vector_delete(py);
1462
1463 for (i = 0; i < ninputs; i++)
1464 {
1465 if (fpe_mord[i] == NULL) continue;
1466 cpl_vector_delete(fpe_xobs[i]);
1467 cpl_vector_delete(fpe_wobs[i]);
1468 cpl_vector_delete(fpe_freq[i]);
1469 cpl_vector_delete(fpe_mord[i]);
1470 cpl_vector_delete(fpe_cord[i]);
1471 cpl_vector_delete(heights[i]);
1472 cpl_vector_delete(sigmas[i]);
1473 cpl_vector_delete(fit_errors[i]);
1474 }
1475 cpl_free(fpe_xobs);
1476 cpl_free(fpe_wobs);
1477 cpl_free(fpe_freq);
1478 cpl_free(fpe_mord);
1479 cpl_free(fpe_cord);
1480 cpl_free(heights);
1481 cpl_free(sigmas);
1482 cpl_free(fit_errors);
1483
1484 return result;
1485
1486
1487}
cr2res_dfs_create_lines_diagnostics_table
cpl_table * cr2res_dfs_create_lines_diagnostics_table(int nrows)
Create an empty LINES DIAGNOSTICS table.
Definition cr2res_dfs.c:553
cr2res_etalon_wave_2d
cpl_polynomial * cr2res_etalon_wave_2d(cpl_bivector **spectra, cpl_bivector **spectra_err, cpl_polynomial **wavesol_init, cpl_array **wavesol_init_err, int *orders, int *traces_nb, int ninputs, cpl_size degree_x, cpl_size degree_y, int zp_order, int display, cpl_array **wavelength_error, cpl_table **line_diagnostics)
Create the 2d wavecal fit using etalon peaks.
Definition cr2res_etalon.c:741
cr2res_etalon_find_peaks
cpl_vector * cr2res_etalon_find_peaks(const cpl_vector *in, double height, double distance)
Detect peaks from a 1d periodic signal and store their positions.
Definition cr2res_etalon.c:466
cr2res_etalon_get_local_maxima
cpl_vector * cr2res_etalon_get_local_maxima(const cpl_vector *in, cpl_vector **left_edges, cpl_vector **right_edges)
Find local maxima in a 1D array. This function finds all local maxima in a 1D array and returns the i...
Definition cr2res_etalon.c:278
cr2res_etalon_select_by_peak_distance
cpl_vector * cr2res_etalon_select_by_peak_distance(const cpl_vector *peaks, const cpl_vector *peak_heights, float distance)
Selects the highest peaks that are atleast distance pixel apart.
Definition cr2res_etalon.c:369
cr2res_wave_poly_2d_to_1d
cpl_polynomial * cr2res_wave_poly_2d_to_1d(cpl_polynomial *poly_2d, int order)
Create a 1D Wavelength Polynomial out of a 2D one.
Definition cr2res_wave.c:1776
Generated on Tue Jun 16 2026 00:11:34 for CR2RE Pipeline Reference Manual by doxygen 1.13.2