CR2RE Pipeline Reference Manual 1.7.0
cr2res_pol.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 <string.h>
29#include <math.h>
30#include <cpl.h>
31
32#include "cr2res_pol.h"
33#include "cr2res_dfs.h"
34#include "cr2res_calib.h"
35#include "cr2res_trace.h"
36#include "cr2res_pfits.h"
37
38/*-----------------------------------------------------------------------------
39 Defines
40 -----------------------------------------------------------------------------*/
41
42int cr2res_pol_resample(cpl_vector ** intens,
43 cpl_vector ** wl,
44 cpl_vector ** errors,
45 int n,
46 cpl_vector ** xmin,
47 cpl_vector ** xmax);
48
49#define CR2RES_POL_MODE_1U "1u"
50#define CR2RES_POL_MODE_1D "1d"
51#define CR2RES_POL_MODE_2U "2u"
52#define CR2RES_POL_MODE_2D "2d"
53#define CR2RES_POL_MODE_3U "3u"
54#define CR2RES_POL_MODE_3D "3d"
55#define CR2RES_POL_MODE_4U "4u"
56#define CR2RES_POL_MODE_4D "4d"
57#define CR2RES_POL_MODE_ERROR "??"
58
59#define CR2RES_POL_MODE(n) \
60 (((n) == 0) ? CR2RES_POL_MODE_1U : \
61 (((n) == 1) ? CR2RES_POL_MODE_1D : \
62 (((n) == 2) ? CR2RES_POL_MODE_2U : \
63 (((n) == 3) ? CR2RES_POL_MODE_2D : \
64 (((n) == 4) ? CR2RES_POL_MODE_3U : \
65 (((n) == 5) ? CR2RES_POL_MODE_3D : \
66 (((n) == 6) ? CR2RES_POL_MODE_4U : \
67 (((n) == 7) ? CR2RES_POL_MODE_4D : \
68 CR2RES_POL_MODE_ERROR \
69 ))))))))
70
71/*-----------------------------------------------------------------------------
72 Functions prototypes
73 -----------------------------------------------------------------------------*/
74
75/*----------------------------------------------------------------------------*/
79/*----------------------------------------------------------------------------*/
80
82
83/* TODO demod functions : rebin all 8 vectors to the same wl */
84/* --> wl[0] as reference */
85/*----------------------------------------------------------------------------*/
105/*----------------------------------------------------------------------------*/
107 cpl_vector ** intens,
108 cpl_vector ** wl,
109 cpl_vector ** errors,
110 int n)
111{
112 cpl_bivector * result;
113 cpl_vector * outspec;
114 cpl_vector * outerr;
115 cpl_size size;
116 cpl_vector * xmin, * xmax;
117
118 cpl_vector ** intens_local;
119 cpl_vector ** errors_local;
120 cpl_vector ** wl_local;
121 cpl_size i, j;
122
123 /* Check entries */
124 if (intens == NULL || wl == NULL || errors == NULL) return NULL;
125 if (n != 8) {
126 cpl_msg_error(__func__,
127 "Demodulation fail: Got %i spectra, but expected 8!", n);
128 return NULL;
129 }
130
131 if (intens[0]== NULL) {
132 cpl_msg_error(__func__,
133 "Demodulation fail: Spectrum %s is missing", CR2RES_POL_MODE(0));
134 return NULL;
135 }
136 size = cpl_vector_get_size(intens[0]);
137 for (i = 0; i < n; i++) {
138 if (intens[i] == NULL) {
139 cpl_msg_error(__func__,
140 "Demodulation fail: Spectrum %s is missing", CR2RES_POL_MODE(i));
141 return NULL;
142 }
143 if (wl[i] == NULL) {
144 cpl_msg_error(__func__,
145 "Demodulation fail: Wavelength %s is missing", CR2RES_POL_MODE(i));
146 return NULL;
147 }
148 if (errors[i] == NULL) {
149 cpl_msg_error(__func__,
150 "Demodulation fail: Errors %s are missing", CR2RES_POL_MODE(i));
151 return NULL;
152 }
153 if (cpl_vector_get_size(intens[i]) != size) {
154 cpl_msg_error(__func__, "Spectra have different sizes");
155 return NULL;
156 }
157 if (cpl_vector_get_size(wl[i]) != size) {
158 cpl_msg_error(__func__, "Wavelengths have different sizes");
159 return NULL;
160 }
161 if (cpl_vector_get_size(errors[i]) != size) {
162 cpl_msg_error(__func__, "Errors have different sizes");
163 return NULL;
164 }
165 }
166
167 // Resample to common wavelength grid
168 // This modifies intens and errors, so we should copy them first
169 intens_local = cpl_malloc(n * sizeof(cpl_vector*));
170 errors_local = cpl_malloc(n * sizeof(cpl_vector*));
171 wl_local = cpl_malloc(n * sizeof(cpl_vector*));
172 for (i = 0; i < n; i++){
173 intens_local[i] = cpl_vector_duplicate(intens[i]);
174 errors_local[i] = cpl_vector_duplicate(errors[i]);
175 wl_local[i] = cpl_vector_duplicate(wl[i]);
176 }
177 xmin = cpl_vector_new(n);
178 xmax = cpl_vector_new(n);
179 if (cr2res_pol_resample(intens_local, wl_local,
180 errors_local, n, &xmin, &xmax) == -1){
181 cpl_msg_error(__func__,
182 "Could not resample polarimetry spectra to the same wavelength scale");
183 for (i = 0; i < n; i++)
184 {
185 cpl_vector_delete(intens_local[i]);
186 cpl_vector_delete(errors_local[i]);
187 cpl_vector_delete(wl_local[i]);
188 }
189 cpl_free(intens_local);
190 cpl_free(errors_local);
191 cpl_free(wl_local);
192 cpl_vector_delete(xmin);
193 cpl_vector_delete(xmax);
194 return NULL;
195 }
196
197 result = cpl_bivector_new(size);
198 outspec = cpl_bivector_get_x(result);
199 outerr = cpl_bivector_get_y(result);
200
201 // Initialize to 0
202 for (i = 0; i < size; i++)
203 {
204 cpl_vector_set(outspec, i, 0.);
205 cpl_vector_set(outerr, i, 0.);
206 }
207
208 for (i = cpl_vector_get_max(xmin);
209 i < cpl_vector_get_min(xmax) + 1; i++){
210 double r, s, e;
211 // Calculate R
212 r = cpl_vector_get(intens_local[0], i); // 1u
213 r /= cpl_vector_get(intens_local[2], i); // 2u
214 r *= cpl_vector_get(intens_local[3], i); // 2d
215 r /= cpl_vector_get(intens_local[1], i); // 1d
216 r *= cpl_vector_get(intens_local[5], i); // 3d
217 r /= cpl_vector_get(intens_local[7], i); // 4d
218 r *= cpl_vector_get(intens_local[6], i); // 4u
219 r /= cpl_vector_get(intens_local[4], i); // 3u
220 r = pow(r, 0.25); // 0.25 = 2/n
221
222 // Calculate Spectrum
223 s = (r - 1) / (r + 1);
224
225 // Calculate Error
226 // sum((err / spec)**2)
227 e = 0;
228 for (j = 0; j < n; j++) {
229 double tmp;
230 tmp = cpl_vector_get(errors_local[j], i) /
231 cpl_vector_get(intens_local[j], i);
232 e += tmp * tmp;
233 }
234
235 // 0.5 * R / (R+1)**2 * sqrt(err)
236 e = sqrt(e);
237 e *= 0.5 * r / ((r + 1) * (r + 1));
238
239 cpl_vector_set(outspec, i, s);
240 cpl_vector_set(outerr, i, e);
241 }
242
243 // We need to remove pointless values outside the wavelength overlap
244 for (j = 0; j < cpl_vector_get_max(xmin); j++)
245 {
246 cpl_vector_set(outspec, j, NAN);
247 cpl_vector_set(outerr, j, NAN);
248 }
249 for (j = cpl_vector_get_min(xmax) + 1; j < size; j++)
250 {
251 cpl_vector_set(outspec, j, NAN);
252 cpl_vector_set(outerr, j, NAN);
253 }
254
255 for (i = 0; i < n; i++){
256 cpl_vector_delete(intens_local[i]);
257 cpl_vector_delete(errors_local[i]);
258 cpl_vector_delete(wl_local[i]);
259 }
260 cpl_free(intens_local);
261 cpl_free(errors_local);
262 cpl_free(wl_local);
263 cpl_vector_delete(xmin);
264 cpl_vector_delete(xmax);
265
266 if (cpl_error_get_code() != CPL_ERROR_NONE) {
267 cpl_msg_error(__func__, "Error message: %s", cpl_error_get_message());
268 cpl_bivector_delete(result);
269 cpl_error_reset();
270 return NULL;
271 }
272
273 return result;
274}
275
276/*----------------------------------------------------------------------------*/
288/*----------------------------------------------------------------------------*/
289int cr2res_pol_resample(cpl_vector ** intens,
290 cpl_vector ** wl,
291 cpl_vector ** errors,
292 int n,
293 cpl_vector ** xmin,
294 cpl_vector ** xmax)
295{
296 // TODO: Use some other grid as baseline?
297 double wmin, wmax;
298 cpl_vector * master_wl;
299 cpl_size i, j, m;
300
301 // Construct master_wl
302 // so that we only include points from within ALL wavelength ranges
303 master_wl = cpl_vector_duplicate(wl[0]);
304 wmin = -INFINITY;
305 wmax = INFINITY;
306 for (i = 0; i < n; i++)
307 {
308 if (cpl_vector_get(wl[i], 0) > wmin){
309 wmin = cpl_vector_get(wl[i], 0);
310 }
311 if (cpl_vector_get(wl[i], cpl_vector_get_size(wl[i]) - 1) < wmax){
312 wmax = cpl_vector_get(wl[i], cpl_vector_get_size(wl[i]) - 1);
313 }
314 }
315
316 for (i = 0; i < n; i++){
317 cpl_vector_set(*xmin, i, 0);
318 cpl_vector_set(*xmax, i, cpl_vector_get_size(wl[i]) - 1);
319
320 for (j = 0; j < cpl_vector_get_size(wl[i]); j++){
321 if (cpl_vector_get(wl[i], j) >= wmin){
322 cpl_vector_set(*xmin, i, j);
323 break;
324 }
325 }
326
327 for (j = cpl_vector_get_size(wl[i]) - 1; j >= 0; j--){
328 if (cpl_vector_get(wl[i], j) <= wmax){
329 cpl_vector_set(*xmax, i, j);
330 break;
331 }
332 }
333 }
334
335 // We need to limit master_wl to the valid range, since bivector interpolate
336 // does not extrapolate at all, but gives an error (?)
337 m = cpl_vector_get_size(wl[0]);
338 for (i = 0; i < m; i++)
339 {
340 if (cpl_vector_get(master_wl, i) < wmin){
341 cpl_vector_set(master_wl, i, wmin);
342 } else if (cpl_vector_get(master_wl, i) > wmax){
343 cpl_vector_set(master_wl, i, wmax);
344 }
345 }
346
347 // Interpolate from copies: in place, the output would feed back into the
348 // input, smoothing the spectra and spreading NaNs
349 for (i = 0; i < n; i++) {
350 cpl_bivector *tmp;
351 cpl_bivector *target;
352 cpl_vector *src;
353 src = cpl_vector_duplicate(intens[i]);
354 tmp = cpl_bivector_wrap_vectors(wl[i], src);
355 target = cpl_bivector_wrap_vectors(master_wl, intens[i]);
356 cpl_bivector_interpolate_linear(target, tmp);
357 cpl_bivector_unwrap_vectors(tmp);
358 cpl_bivector_unwrap_vectors(target);
359 cpl_vector_delete(src);
360
361 src = cpl_vector_duplicate(errors[i]);
362 tmp = cpl_bivector_wrap_vectors(wl[i], src);
363 target = cpl_bivector_wrap_vectors(master_wl, errors[i]);
364 cpl_bivector_interpolate_linear(target, tmp);
365 cpl_bivector_unwrap_vectors(tmp);
366 cpl_bivector_unwrap_vectors(target);
367 cpl_vector_delete(src);
368
369 // We copy the wavelength here, even though it is currently noy used in
370 // the rest of the code
371 cpl_vector_copy(wl[i], master_wl);
372
373 for (j = 0; j < cpl_vector_get(*xmin, i); j++) {
374 cpl_vector_set(intens[i], j, 1);
375 cpl_vector_set(errors[i], j, 1);
376 }
377 for (j = cpl_vector_get(*xmax, i) + 1; j < cpl_vector_get_size(intens[i]);
378 j++) {
379 cpl_vector_set(intens[i], j, 1);
380 cpl_vector_set(errors[i], j, 1);
381 }
382 }
383
384 cpl_vector_delete(master_wl);
385
386 return 0;
387}
388
414/*----------------------------------------------------------------------------*/
416 cpl_vector ** intens,
417 cpl_vector ** wl,
418 cpl_vector ** errors,
419 int n)
420{
421 cpl_vector ** swapintens;
422 cpl_vector ** swapwl;
423 cpl_vector ** swaperrors;
424 cpl_vector * tmpintens;
425 cpl_vector * tmpwl;
426 cpl_vector * tmperrors;
427 cpl_bivector * out ;
428 cpl_size size;
429 int i;
430
431 /* Check entries */
432 if (intens == NULL || wl == NULL || errors == NULL) return NULL;
433 if (n != 8) {
434 cpl_msg_error(__func__,
435 "Demodulation fail: Got %i spectra, but expected 8!", n);
436 return NULL;
437 }
438
439 if (intens[0]== NULL) {
440 cpl_msg_error(__func__,
441 "Demodulation fail: Spectrum %s is missing", CR2RES_POL_MODE(0));
442 return NULL;
443 }
444 size = cpl_vector_get_size(intens[0]);
445 for (i = 0; i < n; i++) {
446 if (intens[i] == NULL) {
447 cpl_msg_error(__func__,
448 "Demodulation fail: Spectrum %s is missing", CR2RES_POL_MODE(i));
449 return NULL;
450 }
451 if (wl[i] == NULL) {
452 cpl_msg_error(__func__,
453 "Demodulation fail: Wavelength %s is missing", CR2RES_POL_MODE(i));
454 return NULL;
455 }
456 if (errors[i] == NULL) {
457 cpl_msg_error(__func__,
458 "Demodulation fail: Errors %s are missing", CR2RES_POL_MODE(i));
459 return NULL;
460 }
461 if (cpl_vector_get_size(intens[i]) != size) {
462 cpl_msg_error(__func__, "Spectra have different sizes");
463 return NULL;
464 }
465 if (cpl_vector_get_size(wl[i]) != size) {
466 cpl_msg_error(__func__, "Wavelengths have different sizes");
467 return NULL;
468 }
469 if (cpl_vector_get_size(errors[i]) != size) {
470 cpl_msg_error(__func__, "Errors have different sizes");
471 return NULL;
472 }
473 }
474
475 // copy list to leave original unchanged
476 swapintens = cpl_malloc(n * sizeof(cpl_vector *));
477 swapwl = cpl_malloc(n * sizeof(cpl_vector *));
478 swaperrors = cpl_malloc(n * sizeof(cpl_vector *));
479 for (i=0 ; i<n ; i++) {
480 swapintens[i]=intens[i];
481 swapwl[i]=wl[i];
482 swaperrors[i]=errors[i];
483 }
484
485 // swap index 2 and 3, i.e. 2u and 2d
486 tmpintens = swapintens[2];
487 tmpwl = swapwl[2];
488 tmperrors = swaperrors[2];
489 swapintens[2] = swapintens[3];
490 swapwl[2] = swapwl[3];
491 swaperrors[2] = swaperrors[3];
492 swapintens[3] = tmpintens;
493 swapwl[3] = tmpwl;
494 swaperrors[3] = tmperrors;
495
496 // swap index 7 and 6, i.e. 4d and 4u
497 tmpintens = swapintens[7];
498 tmpwl = swapwl[7];
499 tmperrors = swaperrors[7];
500 swapintens[7] = swapintens[6];
501 swapwl[7] = swapwl[6];
502 swaperrors[7] = swaperrors[6];
503 swapintens[6] = tmpintens;
504 swapwl[6] = tmpwl;
505 swaperrors[6] = tmperrors;
506
507 out = cr2res_pol_demod_stokes(swapintens, swapwl, swaperrors, n);
508
509 cpl_free(swapintens);
510 cpl_free(swapwl);
511 cpl_free(swaperrors);
512
513 return out ;
514}
515
516/*----------------------------------------------------------------------------*/
532/*----------------------------------------------------------------------------*/
534 cpl_vector ** intens,
535 cpl_vector ** wl,
536 cpl_vector ** errors,
537 int n)
538{
539 cpl_bivector * result;
540 cpl_vector * outspec;
541 cpl_vector * outerr;
542 cpl_size size, i, j, imin, imax;
543 cpl_vector * xmin, * xmax;
544
545 cpl_vector ** intens_local;
546 cpl_vector ** errors_local;
547 cpl_vector ** wl_local;
548
549 /* Check entries */
550 if (intens == NULL || wl == NULL || errors == NULL) return NULL;
551 if (n != 8) {
552 cpl_msg_error(__func__,
553 "Demodulation fail: Got %i spectra, but expected 8!", n);
554 return NULL;
555 }
556
557 if (intens[0]== NULL) {
558 cpl_msg_error(__func__,
559 "Demodulation fail: Spectrum %s is missing", CR2RES_POL_MODE(0));
560 return NULL;
561 }
562 size = cpl_vector_get_size(intens[0]);
563 for (i = 0; i < n; i++) {
564 if (intens[i] == NULL) {
565 cpl_msg_error(__func__,
566 "Demodulation fail: Spectrum %s is missing", CR2RES_POL_MODE(i));
567 return NULL;
568 }
569 if (wl[i] == NULL) {
570 cpl_msg_error(__func__,
571 "Demodulation fail: Wavelength %s is missing", CR2RES_POL_MODE(i));
572 return NULL;
573 }
574 if (errors[i] == NULL) {
575 cpl_msg_error(__func__,
576 "Demodulation fail: Errors %s are missing", CR2RES_POL_MODE(i));
577 return NULL;
578 }
579 if (cpl_vector_get_size(intens[i]) != size) {
580 cpl_msg_error(__func__, "Spectra have different sizes");
581 return NULL;
582 }
583 if (cpl_vector_get_size(wl[i]) != size) {
584 cpl_msg_error(__func__, "Wavelengths have different sizes");
585 return NULL;
586 }
587 if (cpl_vector_get_size(errors[i]) != size) {
588 cpl_msg_error(__func__, "Errors have different sizes");
589 return NULL;
590 }
591 }
592
593 // Resample to common wavelength grid
594 // This modifies intens and errors, so we should copy them first
595 intens_local = cpl_malloc(n * sizeof(cpl_vector*));
596 errors_local = cpl_malloc(n * sizeof(cpl_vector*));
597 wl_local = cpl_malloc(n * sizeof(cpl_vector*));
598 for (i = 0; i < n; i++){
599 intens_local[i] = cpl_vector_duplicate(intens[i]);
600 errors_local[i] = cpl_vector_duplicate(errors[i]);
601 wl_local[i] = cpl_vector_duplicate(wl[i]);
602 }
603 xmin = cpl_vector_new(n);
604 xmax = cpl_vector_new(n);
605
606 if (cr2res_pol_resample(intens_local, wl_local,
607 errors_local, n, &xmin, &xmax) == -1){
608 cpl_msg_error(__func__,
609 "Could not resample polarimetry spectra to the same wavelength scale");
610 for (i = 0; i < n; i++)
611 {
612 cpl_vector_delete(intens_local[i]);
613 cpl_vector_delete(errors_local[i]);
614 cpl_vector_delete(wl_local[i]);
615 }
616 cpl_free(intens_local);
617 cpl_free(errors_local);
618 cpl_free(wl_local);
619 cpl_vector_delete(xmin);
620 cpl_vector_delete(xmax);
621 return NULL;
622 }
623
624 result = cpl_bivector_new(size);
625 outspec = cpl_bivector_get_x(result);
626 outerr = cpl_bivector_get_y(result);
627
628 imin = (cpl_size)cpl_vector_get_max(xmin);
629 imax = (cpl_size)cpl_vector_get_min(xmax);
630 for (i = 0; i < size; i++) {
631 double s = 0.0, v = 0.0;
632 if (i < imin || i > imax) {
633 s = v = NAN;
634 } else {
635 for (j = 0; j < n; j++) {
636 double e = cpl_vector_get(errors_local[j], i);
637 s += cpl_vector_get(intens_local[j], i);
638 v += e * e;
639 }
640 }
641 cpl_vector_set(outspec, i, s);
642 cpl_vector_set(outerr, i, sqrt(v));
643 }
644
645 for (i = 0; i < n; i++){
646 cpl_vector_delete(intens_local[i]);
647 cpl_vector_delete(errors_local[i]);
648 cpl_vector_delete(wl_local[i]);
649 }
650 cpl_free(intens_local);
651 cpl_free(errors_local);
652 cpl_free(wl_local);
653 cpl_vector_delete(xmin);
654 cpl_vector_delete(xmax);
655
656 if (cpl_error_get_code() != CPL_ERROR_NONE) {
657 cpl_msg_error(__func__, "Error message: %s", cpl_error_get_message());
658 cpl_bivector_delete(result);
659 cpl_error_reset();
660 return NULL;
661 }
662
663 return result;
664}
665
666/*----------------------------------------------------------------------------*/
677/*----------------------------------------------------------------------------*/
679 int * orders,
680 cpl_vector ** wl,
681 cpl_bivector ** stokes,
682 cpl_bivector ** null,
683 cpl_bivector ** intens,
684 int norders)
685{
686 cpl_table * out ;
687 char * col_name ;
688 const double * px ;
689 const double * py ;
690 int all_null, i, nrows ;
691
692 /* Check entries */
693 if (orders==NULL || wl==NULL || stokes==NULL || null==NULL || intens==NULL)
694 return NULL ;
695
696 /* Check if all bivectors are not null */
697 all_null = 1 ;
698 for (i=0 ; i<norders ; i++)
699 if (wl[i]!=NULL && stokes[i]!=NULL && null[i]!=NULL && intens[i]!=NULL){
700 nrows = cpl_vector_get_size(wl[i]) ;
701 if (cpl_bivector_get_size(stokes[i]) != nrows ||
702 cpl_bivector_get_size(null[i]) != nrows ||
703 cpl_bivector_get_size(intens[i]) != nrows) {
704 cpl_msg_error(__func__, "Invalid Input Sizes") ;
705 return NULL ;
706 }
707 all_null = 0 ;
708 }
709 if (all_null == 1) return NULL ;
710
711 /* Create the table */
712 out = cpl_table_new(nrows);
713 for (i=0 ; i<norders ; i++) {
714 if (wl[i]!=NULL && stokes[i]!=NULL && null[i]!=NULL && intens[i]!=NULL){
715 /* Create POL_WAVELENGTH column */
716 col_name = cr2res_dfs_POL_WAVELENGTH_colname(orders[i]) ;
717 cpl_table_new_column(out, col_name, CPL_TYPE_DOUBLE);
718 cpl_free(col_name) ;
719
720 /* Create POL_STOKES column */
721 col_name = cr2res_dfs_POL_STOKES_colname(orders[i]) ;
722 cpl_table_new_column(out, col_name, CPL_TYPE_DOUBLE);
723 cpl_free(col_name) ;
724
725 /* Create POL_STOKES_ERROR column */
726 col_name = cr2res_dfs_POL_STOKES_ERROR_colname(orders[i]) ;
727 cpl_table_new_column(out, col_name, CPL_TYPE_DOUBLE);
728 cpl_free(col_name) ;
729
730 /* Create POL_NULL column */
731 col_name = cr2res_dfs_POL_NULL_colname(orders[i]) ;
732 cpl_table_new_column(out, col_name, CPL_TYPE_DOUBLE);
733 cpl_free(col_name) ;
734
735 /* Create POL_NULL_ERROR column */
736 col_name = cr2res_dfs_POL_NULL_ERROR_colname(orders[i]) ;
737 cpl_table_new_column(out, col_name, CPL_TYPE_DOUBLE);
738 cpl_free(col_name) ;
739
740 /* Create POL_INTENS column */
741 col_name = cr2res_dfs_POL_INTENS_colname(orders[i]) ;
742 cpl_table_new_column(out, col_name, CPL_TYPE_DOUBLE);
743 cpl_free(col_name) ;
744
745 /* Create POL_INTENS_ERROR column */
746 col_name = cr2res_dfs_POL_INTENS_ERROR_colname(orders[i]) ;
747 cpl_table_new_column(out, col_name, CPL_TYPE_DOUBLE);
748 cpl_free(col_name) ;
749 }
750 }
751 /* Fill the table */
752 for (i=0 ; i<norders ; i++) {
753 if (wl[i]!=NULL && stokes[i]!=NULL && null[i]!=NULL && intens[i]!=NULL){
754 /* Fill POL_WAVELENGTH column */
755 px = cpl_vector_get_data_const(wl[i]) ;
756 col_name = cr2res_dfs_POL_WAVELENGTH_colname(orders[i]) ;
757 cpl_table_copy_data_double(out, col_name, px) ;
758 cpl_free(col_name) ;
759
760 /* Fill POL_STOKES column */
761 px = cpl_bivector_get_x_data_const(stokes[i]) ;
762 col_name = cr2res_dfs_POL_STOKES_colname(orders[i]) ;
763 cpl_table_copy_data_double(out, col_name, px) ;
764 cpl_free(col_name) ;
765 /* Fill POL_STOKES_ERROR column */
766 py = cpl_bivector_get_y_data_const(stokes[i]) ;
767 col_name = cr2res_dfs_POL_STOKES_ERROR_colname(orders[i]) ;
768 cpl_table_copy_data_double(out, col_name, py) ;
769 cpl_free(col_name) ;
770
771 /* Fill POL_NULL column */
772 px = cpl_bivector_get_x_data_const(null[i]) ;
773 col_name = cr2res_dfs_POL_NULL_colname(orders[i]) ;
774 cpl_table_copy_data_double(out, col_name, px) ;
775 cpl_free(col_name) ;
776 /* Fill POL_NULL_ERROR column */
777 py = cpl_bivector_get_y_data_const(null[i]) ;
778 col_name = cr2res_dfs_POL_NULL_ERROR_colname(orders[i]) ;
779 cpl_table_copy_data_double(out, col_name, py) ;
780 cpl_free(col_name) ;
781
782 /* Fill POL_INTENS column */
783 px = cpl_bivector_get_x_data_const(intens[i]) ;
784 col_name = cr2res_dfs_POL_INTENS_colname(orders[i]) ;
785 cpl_table_copy_data_double(out, col_name, px) ;
786 cpl_free(col_name) ;
787 /* Fill POL_INTENS_ERROR column */
788 py = cpl_bivector_get_y_data_const(intens[i]) ;
789 col_name = cr2res_dfs_POL_INTENS_ERROR_colname(orders[i]) ;
790 cpl_table_copy_data_double(out, col_name, py) ;
791 cpl_free(col_name) ;
792 }
793 }
794 return out ;
795}
796
797/*----------------------------------------------------------------------------*/
814/*----------------------------------------------------------------------------*/
816 const cpl_frame * frame1,
817 const cpl_frame * frame2,
818 const cpl_frame * frame3,
819 const cpl_frame * frame4)
820{
821 int * idx_order ;
822 const char * poltype_first ;
823 const char * fname ;
824 cpl_propertylist * pri_head ;
825 const cpl_frame * frames_holder[4] ;
826 int nframes, i ;
827
828 /* Check Inputs */
829 if (frame1 == NULL || frame2 == NULL || frame3 == NULL || frame4 == NULL)
830 return NULL;
831
832 /* Initialise */
833 nframes = 4 ;
834 frames_holder[0] = frame1 ;
835 frames_holder[1] = frame2 ;
836 frames_holder[2] = frame3 ;
837 frames_holder[3] = frame4 ;
838
839 /* Read POL.TYPE */
840 fname = cpl_frame_get_filename(frames_holder[0]);
841 pri_head = cpl_propertylist_load(fname, 0);
842 poltype_first = cr2res_pfits_get_poltype(pri_head);
843 if (cpl_error_get_code()) {
844 cpl_propertylist_delete(pri_head);
845 cpl_msg_error(__func__, "Missing POL.TYPE in header") ;
846 return NULL ;
847 }
848
849 /* Make sure the others are the same */
850 for (i=1 ; i<nframes ; i++) {
851 const char *poltype_curr;
852 cpl_propertylist *tmp_head;
853 fname = cpl_frame_get_filename(frames_holder[i]);
854 tmp_head = cpl_propertylist_load(fname, 0);
855 poltype_curr = cr2res_pfits_get_poltype(tmp_head);
856 if (cpl_error_get_code()) {
857 cpl_propertylist_delete(tmp_head);
858 cpl_propertylist_delete(pri_head);
859 cpl_msg_error(__func__, "Missing POL.TYPE in header") ;
860 return NULL ;
861 }
862 if (strcmp(poltype_first, poltype_curr)) {
863 cpl_msg_error(__func__, "Different POL.TYPE! %s != %s",
864 poltype_first, poltype_curr) ;
865 return NULL ;
866 }
867 cpl_propertylist_delete(tmp_head);
868 }
869
870 idx_order = cpl_malloc(nframes * sizeof(int)) ;
871
872 if (!strcmp(poltype_first, "V")){
873 cpl_msg_info(__func__, "Set up file order for circular pol.");
874 idx_order[0] = 0 ;
875 idx_order[1] = 1 ;
876 idx_order[2] = 2 ;
877 idx_order[3] = 3 ;
878 } else {
879 cpl_msg_info(__func__, "Set up file order for linear pol.");
880 idx_order[0] = 0 ;
881 idx_order[1] = 3 ;
882 idx_order[2] = 2 ;
883 idx_order[3] = 1 ;
884 }
885 cpl_propertylist_delete(pri_head);
886 return idx_order ;
887}
888
889/*----------------------------------------------------------------------------*/
898/*----------------------------------------------------------------------------*/
900 const cpl_table ** pol_spec_list,
901 int pol_spec_nb)
902{
903 int i, j, k;
904 int n_cols, n_rows;
905 cpl_table * merged_table;
906 double sum, avg;
907 cpl_array * column_names_array;
908
909 /* Check Inputs */
910 if (pol_spec_list == NULL || pol_spec_nb <= 0) return NULL;
911 for (i = 0; i < pol_spec_nb; i++)
912 if (pol_spec_list[i] == NULL) return NULL;
913
914 n_cols = cpl_table_get_ncol(pol_spec_list[0]);
915 n_rows = cpl_table_get_nrow(pol_spec_list[0]);
916 merged_table = cpl_table_new(n_rows);
917
918 column_names_array = cpl_table_get_column_names(pol_spec_list[0]);
919
920 for (i = 0; i < n_cols; i++) {
921 const char *col_name = cpl_array_get_string(column_names_array, i);
922 size_t len = strlen(col_name);
923 int is_err = len >= 4 && !strcmp(col_name + len - 4, "_ERR");
924 cpl_table_new_column(merged_table, col_name, CPL_TYPE_DOUBLE);
925
926 /* Loop through the rows */
927 for (j = 0; j < n_rows; j++) {
928 sum = 0.0;
929 for (k = 0; k < pol_spec_nb; k++) {
930 double val = cpl_table_get_double(pol_spec_list[k], col_name,
931 j, NULL);
932 sum += is_err ? val * val : val;
933 }
934 avg = is_err ? sqrt(sum) / pol_spec_nb : sum / pol_spec_nb;
935 cpl_table_set_double(merged_table, col_name, j, avg);
936 }
937 }
938
939 /* Clean up */
940 cpl_array_delete(column_names_array);
941
942 return merged_table;
943}
944
945/*----------------------------------------------------------------------------*/
956/*----------------------------------------------------------------------------*/
957
959 const cpl_table *tw_in,
960 cr2res_decker decker_position,
961 int up_or_down)
962{
963
964 cpl_table * tw_out = NULL;
965 cpl_polynomial * wl_poly;
966 cpl_polynomial * corr_poly;
967 double wl, y_corr, newSF;
968 cpl_array * tmp_arr;
969 cpl_size pow0=0;
970 cpl_size pow1=1;
971 char * band;
972 int nb_traces, i;
973 cpl_vector * newSFs;
974
975
976 wl_poly = cr2res_get_trace_wave_poly(tw_in, CR2RES_COL_WAVELENGTH, 5, 1);
977 wl = cpl_polynomial_eval_1d(wl_poly, 1024, NULL);
978 cpl_polynomial_delete(wl_poly);
979
980 /* Correction polynomial gets filled with hardcoded values. */
981 /* These were derived as linear difference from the order mid-line */
982 /* for the upper/lower beams and A/B nodding position, in each band.*/
983 corr_poly = cpl_polynomial_new(1);
984 if (wl < 1125) {
985 band=cpl_sprintf("Y");
986 /*
987 'YuA': array([ 0.03988915, -22.78903764]),
988 'YuB': array([ 0.04808708, 13.12469913]),
989 'YdA': array([-0.04965059, -6.51155269]),
990 'YdB': array([-0.04019722, 28.22704312]),
991 */
992 if (decker_position == CR2RES_DECKER_2_4){ // nodd A
993 if (up_or_down == 1){ // UP
994 cpl_polynomial_set_coeff(corr_poly, &pow0, -22.78903764);
995 cpl_polynomial_set_coeff(corr_poly, &pow1, 0.03988915);
996 } else { // DOWN
997 cpl_polynomial_set_coeff(corr_poly, &pow0, -6.51155269);
998 cpl_polynomial_set_coeff(corr_poly, &pow1, -0.04965059);
999 }
1000 } else if (decker_position == CR2RES_DECKER_1_3){ // nodd B
1001 if (up_or_down == 1){ // UP
1002 cpl_polynomial_set_coeff(corr_poly, &pow0, 13.12469913);
1003 cpl_polynomial_set_coeff(corr_poly, &pow1, 0.04808708);
1004 } else { // DOWN
1005 cpl_polynomial_set_coeff(corr_poly, &pow0, 28.22704312);
1006 cpl_polynomial_set_coeff(corr_poly, &pow1, -0.04019722);
1007 }
1008 }
1009 }
1010 else if (wl < 1360) {
1011 band=cpl_sprintf("J");
1012 /*
1013 'JuA': array([ 0.03903698, -26.55176271]),
1014 'JuB': array([ 0.04496791, 10.27999294]),
1015 'JdA': array([ -0.047482 , -10.87508078]),
1016 'JdB': array([-0.04101687, 25.44306868]),
1017 */
1018 if (decker_position == CR2RES_DECKER_2_4){ // nodd A
1019 if (up_or_down == 1){ // UP
1020 cpl_polynomial_set_coeff(corr_poly, &pow0, -26.55176271);
1021 cpl_polynomial_set_coeff(corr_poly, &pow1, 0.03903698);
1022 } else { // DOWN
1023 cpl_polynomial_set_coeff(corr_poly, &pow0, -10.87508078);
1024 cpl_polynomial_set_coeff(corr_poly, &pow1, -0.047482);
1025 }
1026 } else if (decker_position == CR2RES_DECKER_1_3){ // nodd B
1027 if (up_or_down == 1){ // UP
1028 cpl_polynomial_set_coeff(corr_poly, &pow0, 10.27999294);
1029 cpl_polynomial_set_coeff(corr_poly, &pow1, 0.04496791);
1030 } else { // DOWN
1031 cpl_polynomial_set_coeff(corr_poly, &pow0, 25.44306868);
1032 cpl_polynomial_set_coeff(corr_poly, &pow1, -0.04101687);
1033 }
1034 }
1035
1036 }
1037 else if (wl < 1850) {
1038 band=cpl_sprintf("H");
1039 /*
1040 'HuA': array([ 1.59902724e-02, -2.17975984e+01]),
1041 'HuB': array([ 0.0199041 , 17.44054878]),
1042 'HdA': array([ -0.02831508, -23.23946728]),
1043 'HdB': array([-0.02552662, 18.31347609]),
1044 */
1045 if (decker_position == CR2RES_DECKER_2_4){ // nodd A
1046 if (up_or_down == 1){ // UP
1047 cpl_polynomial_set_coeff(corr_poly, &pow0, -2.17975984e+01);
1048 cpl_polynomial_set_coeff(corr_poly, &pow1, 1.59902724e-02);
1049 } else { // DOWN
1050 cpl_polynomial_set_coeff(corr_poly, &pow0, -23.23946728);
1051 cpl_polynomial_set_coeff(corr_poly, &pow1, -0.02831508);
1052 }
1053 } else if (decker_position == CR2RES_DECKER_1_3){ // nodd B
1054 if (up_or_down == 1){ // UP
1055 cpl_polynomial_set_coeff(corr_poly, &pow0, 17.44054878);
1056 cpl_polynomial_set_coeff(corr_poly, &pow1, 0.0199041);
1057 } else { // DOWN
1058 cpl_polynomial_set_coeff(corr_poly, &pow0, 18.31347609);
1059 cpl_polynomial_set_coeff(corr_poly, &pow1, -0.02552662);
1060 }
1061 }
1062 }
1063 else if (wl < 2600) {
1064 band=cpl_sprintf("K");
1065 /*
1066 'KuA': array([ 2.24337404e-02, -2.30870052e+01]),
1067 'KuB': array([ 0.02467974, 16.43404529]),
1068 'KdA': array([ -0.02667301, -13.55392217]),
1069 'KdB': array([-2.42857387e-02, 2.58313952e+01])
1070 */
1071 if (decker_position == CR2RES_DECKER_2_4){ // nodd A
1072 if (up_or_down == 1){ // UP
1073 cpl_polynomial_set_coeff(corr_poly, &pow0, -2.30870052e+01);
1074 cpl_polynomial_set_coeff(corr_poly, &pow1, 2.24337404e-02);
1075 } else { // DOWN
1076 cpl_polynomial_set_coeff(corr_poly, &pow0, -13.55392217);
1077 cpl_polynomial_set_coeff(corr_poly, &pow1, -0.02667301);
1078 }
1079 } else if (decker_position == CR2RES_DECKER_1_3){ // nodd B
1080 if (up_or_down == 1){ // UP
1081 cpl_polynomial_set_coeff(corr_poly, &pow0, 16.43404529);
1082 cpl_polynomial_set_coeff(corr_poly, &pow1, 0.02467974);
1083 } else { // DOWN
1084 cpl_polynomial_set_coeff(corr_poly, &pow0, 2.58313952e+01);
1085 cpl_polynomial_set_coeff(corr_poly, &pow1, -2.42857387e-02);
1086 }
1087 }
1088 }
1089 else {
1090 cpl_msg_error(__func__, "No proper WL found to identify band.");
1091 cpl_polynomial_delete(corr_poly);
1092 return NULL;
1093 }
1094 cpl_msg_info(__func__, "Found wl=%g nm, therefore applying beam "
1095 "correction for %s-band.", wl, band);
1096
1097
1098 /* Loop through all traces in input TW */
1099 nb_traces = cpl_table_get_nrow(tw_in) ;
1100 newSFs = cpl_vector_new(nb_traces);
1101 for (i = 0; i < nb_traces; i++)
1102 {
1103 int trace_id;
1104 cpl_table_get(tw_in, CR2RES_COL_ORDER, i, NULL) ;
1105 trace_id = cpl_table_get(tw_in, CR2RES_COL_TRACENB, i, NULL);
1106 if (trace_id != 1){
1107 cpl_msg_error(__func__, "More than one input-trace per order"
1108 "is not supported.");
1109 cpl_table_delete(tw_out);
1110 cpl_free(band);
1111 cpl_polynomial_delete(corr_poly);
1112 return NULL;
1113 }
1114
1115 /* Evaluate at middle of detector to get average positions */
1116 y_corr = cpl_polynomial_eval_1d(corr_poly, wl, NULL);
1117 cpl_vector_set(newSFs, i, 0.5 + (y_corr) / CR2RES_APPROX_SLIT_HEIGHT);
1118 }
1119 newSF = cpl_vector_get_median(newSFs);
1120 cpl_vector_delete(newSFs);
1121 cpl_msg_info(__func__, "Median new slit-fraction: %g", newSF);
1122
1123 tmp_arr = cpl_array_new(3,CPL_TYPE_DOUBLE);
1124 cpl_array_set(tmp_arr, 0, newSF-0.15);
1125 cpl_array_set(tmp_arr, 1, newSF);
1126 cpl_array_set(tmp_arr, 2, newSF+0.15);
1127 tw_out = cr2res_trace_new_slit_fraction(tw_in, tmp_arr);
1128 //tw_out = cpl_table_duplicate(tw_in);
1129 cpl_array_delete(tmp_arr);
1130
1131 for (i=0 ; i<nb_traces ; i++) {
1132 double slope_corr, halfSF;
1133 int o;
1134 o = cpl_table_get(tw_in, CR2RES_COL_ORDER, i, NULL) ;
1135 //trace_id = cpl_table_get(tw_in, CR2RES_COL_TRACENB, i, NULL);
1136
1137 /* Calculate new slit-fraction, important for WL-correction */
1138 /* wl is at middle of detector from above */
1139 y_corr = cpl_polynomial_eval_1d(corr_poly, wl, NULL);
1140 tmp_arr = cpl_array_new(3,CPL_TYPE_DOUBLE);
1141 newSF = 0.5 + (y_corr / CR2RES_APPROX_SLIT_HEIGHT);
1142 halfSF = CR2RES_POLARIMETRY_DEFAULT_HEIGHT /
1143 CR2RES_APPROX_SLIT_HEIGHT / 2;
1144 cpl_array_set(tmp_arr, 0, newSF-halfSF);
1145 cpl_array_set(tmp_arr, 1, newSF);
1146 cpl_array_set(tmp_arr, 2, newSF+halfSF);
1147 cpl_table_set_array(tw_out, CR2RES_COL_SLIT_FRACTION, i, tmp_arr);
1148 cpl_array_delete(tmp_arr);
1149
1150 /* Evaluate WL at detector edges, because corr-poly is P(WL)*/
1151 wl_poly = cr2res_get_trace_wave_poly(tw_in, CR2RES_COL_WAVELENGTH, o,1);
1152 wl = cpl_polynomial_eval_1d(wl_poly, 1, NULL);
1153 y_corr = cpl_polynomial_eval_1d(corr_poly, wl, NULL);
1154 wl = cpl_polynomial_eval_1d(wl_poly, CR2RES_DETECTOR_SIZE, NULL);
1155 slope_corr = (cpl_polynomial_eval_1d(corr_poly, wl, NULL)
1156 - y_corr) / CR2RES_DETECTOR_SIZE; // SLOPE!
1157
1158 cpl_msg_debug(__func__,
1159 "newSF, halfSF, slope_corr, y_corr: "
1160 "%g, %g, %g, %g ",
1161 newSF, halfSF, slope_corr, y_corr);
1162
1163 /* Add the correction to the existing trace polys */
1164 tmp_arr = cpl_array_duplicate (
1165 cpl_table_get_array(tw_in, CR2RES_COL_ALL,i));
1166 cpl_array_set(tmp_arr, 0, cpl_array_get(tmp_arr, 0, NULL) + y_corr);
1167 cpl_array_set(tmp_arr, 1, cpl_array_get(tmp_arr, 1, NULL) + slope_corr);
1168 cpl_table_set_array(tw_out, CR2RES_COL_ALL, i, tmp_arr);
1169
1170 /* Edge polynomials are just shifted up/down by height/2 */
1171 cpl_array_set(tmp_arr, 0, cpl_array_get(tmp_arr, 0, NULL) +
1172 CR2RES_POLARIMETRY_DEFAULT_HEIGHT / 2);
1173 cpl_table_set_array(tw_out, CR2RES_COL_UPPER, i, tmp_arr);
1174 cpl_array_set(tmp_arr, 0, cpl_array_get(tmp_arr, 0, NULL) -
1175 CR2RES_POLARIMETRY_DEFAULT_HEIGHT);
1176 cpl_table_set_array(tw_out, CR2RES_COL_LOWER, i, tmp_arr);
1177 cpl_array_delete(tmp_arr);
1178
1179 cpl_polynomial_delete(wl_poly);
1180 }
1181
1182
1183 cpl_free(band);
1184 cpl_polynomial_delete(corr_poly);
1185 return tw_out;
1186}
1187
1188//
1189
char * cr2res_dfs_POL_NULL_colname(int order_idx)
Get the POL_NULL column name for a given order.
Definition cr2res_dfs.c:258
char * cr2res_dfs_POL_STOKES_ERROR_colname(int order_idx)
Get the POL_STOKES_ERROR column name for a given order.
Definition cr2res_dfs.c:241
char * cr2res_dfs_POL_INTENS_ERROR_colname(int order_idx)
Get the POL_INTENS_ERROR column name for a given order.
Definition cr2res_dfs.c:306
char * cr2res_dfs_POL_NULL_ERROR_colname(int order_idx)
Get the POL_NULL_ERROR column name for a given order.
Definition cr2res_dfs.c:274
char * cr2res_dfs_POL_STOKES_colname(int order_idx)
Get the POL_STOKES column name for a given order.
Definition cr2res_dfs.c:225
char * cr2res_dfs_POL_INTENS_colname(int order_idx)
Get the POL_INTENS column name for a given order.
Definition cr2res_dfs.c:290
char * cr2res_dfs_POL_WAVELENGTH_colname(int order_idx)
Get the POL_WAVELENGTH column name for a given order.
Definition cr2res_dfs.c:209
const char * cr2res_pfits_get_poltype(const cpl_propertylist *plist)
find out the POL TYPE value
cpl_bivector * cr2res_pol_demod_intens(cpl_vector **intens, cpl_vector **wl, cpl_vector **errors, int n)
Combine extracted spectra into Intensity spectrum.
Definition cr2res_pol.c:533
cpl_table * cr2res_pol_POL_SPEC_create(int *orders, cpl_vector **wl, cpl_bivector **stokes, cpl_bivector **null, cpl_bivector **intens, int norders)
Create the POL_SPEC table to be saved.
Definition cr2res_pol.c:678
cpl_table * cr2res_pol_spec_pol_merge(const cpl_table **pol_spec_list, int pol_spec_nb)
Merge several POL_SPEC tables together, by averaging.
Definition cr2res_pol.c:899
int * cr2res_pol_sort_frames(const cpl_frame *frame1, const cpl_frame *frame2, const cpl_frame *frame3, const cpl_frame *frame4)
Compute the positions of the passed frames.
Definition cr2res_pol.c:815
cpl_table * cr2res_pol_get_beam_trace(const cpl_table *tw_in, cr2res_decker decker_position, int up_or_down)
Compute the traces for the polarimetric beams.
Definition cr2res_pol.c:958
cpl_bivector * cr2res_pol_demod_stokes(cpl_vector **intens, cpl_vector **wl, cpl_vector **errors, int n)
Demodulate extracted spectra into Stokes parameter.
Definition cr2res_pol.c:106
int cr2res_pol_resample(cpl_vector **intens, cpl_vector **wl, cpl_vector **errors, int n, cpl_vector **xmin, cpl_vector **xmax)
Resample all spectra to the same wavelength grid.
Definition cr2res_pol.c:289
cpl_bivector * cr2res_pol_demod_null(cpl_vector **intens, cpl_vector **wl, cpl_vector **errors, int n)
Demodulate extracted spectra into Null spectrum.
Definition cr2res_pol.c:415
cpl_table * cr2res_trace_new_slit_fraction(const cpl_table *traces, const cpl_array *new_slit_fraction)
Recompute the traces at a newly specified slit fraction.
cpl_polynomial * cr2res_get_trace_wave_poly(const cpl_table *trace_wave, const char *poly_column, int order_idx, int trace_nb)
Get a polynomial from a TRACE_WAVE table.