107 cpl_vector ** intens,
109 cpl_vector ** errors,
112 cpl_bivector * result;
113 cpl_vector * outspec;
116 cpl_vector * xmin, * xmax;
118 cpl_vector ** intens_local;
119 cpl_vector ** errors_local;
120 cpl_vector ** wl_local;
124 if (intens == NULL || wl == NULL || errors == NULL)
return NULL;
126 cpl_msg_error(__func__,
127 "Demodulation fail: Got %i spectra, but expected 8!", n);
131 if (intens[0]== NULL) {
132 cpl_msg_error(__func__,
133 "Demodulation fail: Spectrum %s is missing", CR2RES_POL_MODE(0));
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));
144 cpl_msg_error(__func__,
145 "Demodulation fail: Wavelength %s is missing", CR2RES_POL_MODE(i));
148 if (errors[i] == NULL) {
149 cpl_msg_error(__func__,
150 "Demodulation fail: Errors %s are missing", CR2RES_POL_MODE(i));
153 if (cpl_vector_get_size(intens[i]) != size) {
154 cpl_msg_error(__func__,
"Spectra have different sizes");
157 if (cpl_vector_get_size(wl[i]) != size) {
158 cpl_msg_error(__func__,
"Wavelengths have different sizes");
161 if (cpl_vector_get_size(errors[i]) != size) {
162 cpl_msg_error(__func__,
"Errors have different sizes");
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]);
177 xmin = cpl_vector_new(n);
178 xmax = cpl_vector_new(n);
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++)
185 cpl_vector_delete(intens_local[i]);
186 cpl_vector_delete(errors_local[i]);
187 cpl_vector_delete(wl_local[i]);
189 cpl_free(intens_local);
190 cpl_free(errors_local);
192 cpl_vector_delete(xmin);
193 cpl_vector_delete(xmax);
197 result = cpl_bivector_new(size);
198 outspec = cpl_bivector_get_x(result);
199 outerr = cpl_bivector_get_y(result);
202 for (i = 0; i < size; i++)
204 cpl_vector_set(outspec, i, 0.);
205 cpl_vector_set(outerr, i, 0.);
208 for (i = cpl_vector_get_max(xmin);
209 i < cpl_vector_get_min(xmax) + 1; i++){
212 r = cpl_vector_get(intens_local[0], i);
213 r /= cpl_vector_get(intens_local[2], i);
214 r *= cpl_vector_get(intens_local[3], i);
215 r /= cpl_vector_get(intens_local[1], i);
216 r *= cpl_vector_get(intens_local[5], i);
217 r /= cpl_vector_get(intens_local[7], i);
218 r *= cpl_vector_get(intens_local[6], i);
219 r /= cpl_vector_get(intens_local[4], i);
223 s = (r - 1) / (r + 1);
228 for (j = 0; j < n; j++) {
230 tmp = cpl_vector_get(errors_local[j], i) /
231 cpl_vector_get(intens_local[j], i);
237 e *= 0.5 * r / ((r + 1) * (r + 1));
239 cpl_vector_set(outspec, i, s);
240 cpl_vector_set(outerr, i, e);
244 for (j = 0; j < cpl_vector_get_max(xmin); j++)
246 cpl_vector_set(outspec, j, NAN);
247 cpl_vector_set(outerr, j, NAN);
249 for (j = cpl_vector_get_min(xmax) + 1; j < size; j++)
251 cpl_vector_set(outspec, j, NAN);
252 cpl_vector_set(outerr, j, NAN);
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]);
260 cpl_free(intens_local);
261 cpl_free(errors_local);
263 cpl_vector_delete(xmin);
264 cpl_vector_delete(xmax);
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);
291 cpl_vector ** errors,
298 cpl_vector * master_wl;
303 master_wl = cpl_vector_duplicate(wl[0]);
306 for (i = 0; i < n; i++)
308 if (cpl_vector_get(wl[i], 0) > wmin){
309 wmin = cpl_vector_get(wl[i], 0);
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);
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);
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);
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);
337 m = cpl_vector_get_size(wl[0]);
338 for (i = 0; i < m; i++)
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);
349 for (i = 0; i < n; i++) {
351 cpl_bivector *target;
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);
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);
371 cpl_vector_copy(wl[i], master_wl);
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);
377 for (j = cpl_vector_get(*xmax, i) + 1; j < cpl_vector_get_size(intens[i]);
379 cpl_vector_set(intens[i], j, 1);
380 cpl_vector_set(errors[i], j, 1);
384 cpl_vector_delete(master_wl);
534 cpl_vector ** intens,
536 cpl_vector ** errors,
539 cpl_bivector * result;
540 cpl_vector * outspec;
542 cpl_size size, i, j, imin, imax;
543 cpl_vector * xmin, * xmax;
545 cpl_vector ** intens_local;
546 cpl_vector ** errors_local;
547 cpl_vector ** wl_local;
550 if (intens == NULL || wl == NULL || errors == NULL)
return NULL;
552 cpl_msg_error(__func__,
553 "Demodulation fail: Got %i spectra, but expected 8!", n);
557 if (intens[0]== NULL) {
558 cpl_msg_error(__func__,
559 "Demodulation fail: Spectrum %s is missing", CR2RES_POL_MODE(0));
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));
570 cpl_msg_error(__func__,
571 "Demodulation fail: Wavelength %s is missing", CR2RES_POL_MODE(i));
574 if (errors[i] == NULL) {
575 cpl_msg_error(__func__,
576 "Demodulation fail: Errors %s are missing", CR2RES_POL_MODE(i));
579 if (cpl_vector_get_size(intens[i]) != size) {
580 cpl_msg_error(__func__,
"Spectra have different sizes");
583 if (cpl_vector_get_size(wl[i]) != size) {
584 cpl_msg_error(__func__,
"Wavelengths have different sizes");
587 if (cpl_vector_get_size(errors[i]) != size) {
588 cpl_msg_error(__func__,
"Errors have different sizes");
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]);
603 xmin = cpl_vector_new(n);
604 xmax = cpl_vector_new(n);
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++)
612 cpl_vector_delete(intens_local[i]);
613 cpl_vector_delete(errors_local[i]);
614 cpl_vector_delete(wl_local[i]);
616 cpl_free(intens_local);
617 cpl_free(errors_local);
619 cpl_vector_delete(xmin);
620 cpl_vector_delete(xmax);
624 result = cpl_bivector_new(size);
625 outspec = cpl_bivector_get_x(result);
626 outerr = cpl_bivector_get_y(result);
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) {
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);
641 cpl_vector_set(outspec, i, s);
642 cpl_vector_set(outerr, i, sqrt(v));
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]);
650 cpl_free(intens_local);
651 cpl_free(errors_local);
653 cpl_vector_delete(xmin);
654 cpl_vector_delete(xmax);
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);
681 cpl_bivector ** stokes,
682 cpl_bivector ** null,
683 cpl_bivector ** intens,
690 int all_null, i, nrows ;
693 if (orders==NULL || wl==NULL || stokes==NULL || null==NULL || intens==NULL)
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") ;
709 if (all_null == 1)
return NULL ;
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){
717 cpl_table_new_column(out, col_name, CPL_TYPE_DOUBLE);
722 cpl_table_new_column(out, col_name, CPL_TYPE_DOUBLE);
727 cpl_table_new_column(out, col_name, CPL_TYPE_DOUBLE);
732 cpl_table_new_column(out, col_name, CPL_TYPE_DOUBLE);
737 cpl_table_new_column(out, col_name, CPL_TYPE_DOUBLE);
742 cpl_table_new_column(out, col_name, CPL_TYPE_DOUBLE);
747 cpl_table_new_column(out, col_name, CPL_TYPE_DOUBLE);
752 for (i=0 ; i<norders ; i++) {
753 if (wl[i]!=NULL && stokes[i]!=NULL && null[i]!=NULL && intens[i]!=NULL){
755 px = cpl_vector_get_data_const(wl[i]) ;
757 cpl_table_copy_data_double(out, col_name, px) ;
761 px = cpl_bivector_get_x_data_const(stokes[i]) ;
763 cpl_table_copy_data_double(out, col_name, px) ;
766 py = cpl_bivector_get_y_data_const(stokes[i]) ;
768 cpl_table_copy_data_double(out, col_name, py) ;
772 px = cpl_bivector_get_x_data_const(null[i]) ;
774 cpl_table_copy_data_double(out, col_name, px) ;
777 py = cpl_bivector_get_y_data_const(null[i]) ;
779 cpl_table_copy_data_double(out, col_name, py) ;
783 px = cpl_bivector_get_x_data_const(intens[i]) ;
785 cpl_table_copy_data_double(out, col_name, px) ;
788 py = cpl_bivector_get_y_data_const(intens[i]) ;
790 cpl_table_copy_data_double(out, col_name, py) ;
959 const cpl_table *tw_in,
960 cr2res_decker decker_position,
964 cpl_table * tw_out = NULL;
965 cpl_polynomial * wl_poly;
966 cpl_polynomial * corr_poly;
967 double wl, y_corr, newSF;
977 wl = cpl_polynomial_eval_1d(wl_poly, 1024, NULL);
978 cpl_polynomial_delete(wl_poly);
983 corr_poly = cpl_polynomial_new(1);
985 band=cpl_sprintf(
"Y");
992 if (decker_position == CR2RES_DECKER_2_4){
993 if (up_or_down == 1){
994 cpl_polynomial_set_coeff(corr_poly, &pow0, -22.78903764);
995 cpl_polynomial_set_coeff(corr_poly, &pow1, 0.03988915);
997 cpl_polynomial_set_coeff(corr_poly, &pow0, -6.51155269);
998 cpl_polynomial_set_coeff(corr_poly, &pow1, -0.04965059);
1000 }
else if (decker_position == CR2RES_DECKER_1_3){
1001 if (up_or_down == 1){
1002 cpl_polynomial_set_coeff(corr_poly, &pow0, 13.12469913);
1003 cpl_polynomial_set_coeff(corr_poly, &pow1, 0.04808708);
1005 cpl_polynomial_set_coeff(corr_poly, &pow0, 28.22704312);
1006 cpl_polynomial_set_coeff(corr_poly, &pow1, -0.04019722);
1010 else if (wl < 1360) {
1011 band=cpl_sprintf(
"J");
1018 if (decker_position == CR2RES_DECKER_2_4){
1019 if (up_or_down == 1){
1020 cpl_polynomial_set_coeff(corr_poly, &pow0, -26.55176271);
1021 cpl_polynomial_set_coeff(corr_poly, &pow1, 0.03903698);
1023 cpl_polynomial_set_coeff(corr_poly, &pow0, -10.87508078);
1024 cpl_polynomial_set_coeff(corr_poly, &pow1, -0.047482);
1026 }
else if (decker_position == CR2RES_DECKER_1_3){
1027 if (up_or_down == 1){
1028 cpl_polynomial_set_coeff(corr_poly, &pow0, 10.27999294);
1029 cpl_polynomial_set_coeff(corr_poly, &pow1, 0.04496791);
1031 cpl_polynomial_set_coeff(corr_poly, &pow0, 25.44306868);
1032 cpl_polynomial_set_coeff(corr_poly, &pow1, -0.04101687);
1037 else if (wl < 1850) {
1038 band=cpl_sprintf(
"H");
1045 if (decker_position == CR2RES_DECKER_2_4){
1046 if (up_or_down == 1){
1047 cpl_polynomial_set_coeff(corr_poly, &pow0, -2.17975984e+01);
1048 cpl_polynomial_set_coeff(corr_poly, &pow1, 1.59902724e-02);
1050 cpl_polynomial_set_coeff(corr_poly, &pow0, -23.23946728);
1051 cpl_polynomial_set_coeff(corr_poly, &pow1, -0.02831508);
1053 }
else if (decker_position == CR2RES_DECKER_1_3){
1054 if (up_or_down == 1){
1055 cpl_polynomial_set_coeff(corr_poly, &pow0, 17.44054878);
1056 cpl_polynomial_set_coeff(corr_poly, &pow1, 0.0199041);
1058 cpl_polynomial_set_coeff(corr_poly, &pow0, 18.31347609);
1059 cpl_polynomial_set_coeff(corr_poly, &pow1, -0.02552662);
1063 else if (wl < 2600) {
1064 band=cpl_sprintf(
"K");
1071 if (decker_position == CR2RES_DECKER_2_4){
1072 if (up_or_down == 1){
1073 cpl_polynomial_set_coeff(corr_poly, &pow0, -2.30870052e+01);
1074 cpl_polynomial_set_coeff(corr_poly, &pow1, 2.24337404e-02);
1076 cpl_polynomial_set_coeff(corr_poly, &pow0, -13.55392217);
1077 cpl_polynomial_set_coeff(corr_poly, &pow1, -0.02667301);
1079 }
else if (decker_position == CR2RES_DECKER_1_3){
1080 if (up_or_down == 1){
1081 cpl_polynomial_set_coeff(corr_poly, &pow0, 16.43404529);
1082 cpl_polynomial_set_coeff(corr_poly, &pow1, 0.02467974);
1084 cpl_polynomial_set_coeff(corr_poly, &pow0, 2.58313952e+01);
1085 cpl_polynomial_set_coeff(corr_poly, &pow1, -2.42857387e-02);
1090 cpl_msg_error(__func__,
"No proper WL found to identify band.");
1091 cpl_polynomial_delete(corr_poly);
1094 cpl_msg_info(__func__,
"Found wl=%g nm, therefore applying beam "
1095 "correction for %s-band.", wl, band);
1099 nb_traces = cpl_table_get_nrow(tw_in) ;
1100 newSFs = cpl_vector_new(nb_traces);
1101 for (i = 0; i < nb_traces; i++)
1104 cpl_table_get(tw_in, CR2RES_COL_ORDER, i, NULL) ;
1105 trace_id = cpl_table_get(tw_in, CR2RES_COL_TRACENB, i, NULL);
1107 cpl_msg_error(__func__,
"More than one input-trace per order"
1108 "is not supported.");
1109 cpl_table_delete(tw_out);
1111 cpl_polynomial_delete(corr_poly);
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);
1119 newSF = cpl_vector_get_median(newSFs);
1120 cpl_vector_delete(newSFs);
1121 cpl_msg_info(__func__,
"Median new slit-fraction: %g", newSF);
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);
1129 cpl_array_delete(tmp_arr);
1131 for (i=0 ; i<nb_traces ; i++) {
1132 double slope_corr, halfSF;
1134 o = cpl_table_get(tw_in, CR2RES_COL_ORDER, i, NULL) ;
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);
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;
1158 cpl_msg_debug(__func__,
1159 "newSF, halfSF, slope_corr, y_corr: "
1161 newSF, halfSF, slope_corr, y_corr);
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);
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);
1179 cpl_polynomial_delete(wl_poly);
1184 cpl_polynomial_delete(corr_poly);