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;
352 tmp = cpl_bivector_wrap_vectors(wl[i], intens[i]);
353 target = cpl_bivector_wrap_vectors(master_wl, intens[i]);
354 cpl_bivector_interpolate_linear(target, tmp);
355 cpl_bivector_unwrap_vectors(tmp);
356 cpl_bivector_unwrap_vectors(target);
358 tmp = cpl_bivector_wrap_vectors(wl[i], errors[i]);
359 target = cpl_bivector_wrap_vectors(master_wl, errors[i]);
360 cpl_bivector_interpolate_linear(target, tmp);
361 cpl_bivector_unwrap_vectors(tmp);
362 cpl_bivector_unwrap_vectors(target);
366 cpl_vector_copy(wl[i], master_wl);
368 for (j = 0; j < cpl_vector_get(*xmin, i); j++) {
369 cpl_vector_set(intens[i], j, 1);
370 cpl_vector_set(errors[i], j, 1);
372 for (j = cpl_vector_get(*xmax, i) + 1; j < cpl_vector_get_size(intens[i]);
374 cpl_vector_set(intens[i], j, 1);
375 cpl_vector_set(errors[i], j, 1);
379 cpl_vector_delete(master_wl);
529 cpl_vector ** intens,
531 cpl_vector ** errors,
534 cpl_bivector * result;
535 cpl_vector * outspec;
541 cpl_vector * xmin, * xmax;
543 cpl_vector ** intens_local;
544 cpl_vector ** errors_local;
545 cpl_vector ** wl_local;
548 if (intens == NULL || wl == NULL || errors == NULL)
return NULL;
550 cpl_msg_error(__func__,
551 "Demodulation fail: Got %i spectra, but expected 8!", n);
555 if (intens[0]== NULL) {
556 cpl_msg_error(__func__,
557 "Demodulation fail: Spectrum %s is missing", CR2RES_POL_MODE(0));
560 size = cpl_vector_get_size(intens[0]);
561 for (i = 0; i < n; i++) {
562 if (intens[i] == NULL) {
563 cpl_msg_error(__func__,
564 "Demodulation fail: Spectrum %s is missing", CR2RES_POL_MODE(i));
568 cpl_msg_error(__func__,
569 "Demodulation fail: Wavelength %s is missing", CR2RES_POL_MODE(i));
572 if (errors[i] == NULL) {
573 cpl_msg_error(__func__,
574 "Demodulation fail: Errors %s are missing", CR2RES_POL_MODE(i));
577 if (cpl_vector_get_size(intens[i]) != size) {
578 cpl_msg_error(__func__,
"Spectra have different sizes");
581 if (cpl_vector_get_size(wl[i]) != size) {
582 cpl_msg_error(__func__,
"Wavelengths have different sizes");
585 if (cpl_vector_get_size(errors[i]) != size) {
586 cpl_msg_error(__func__,
"Errors have different sizes");
593 intens_local = cpl_malloc(n *
sizeof(cpl_vector*));
594 errors_local = cpl_malloc(n *
sizeof(cpl_vector*));
595 wl_local = cpl_malloc(n *
sizeof(cpl_vector*));
596 for (i = 0; i < n; i++){
597 intens_local[i] = cpl_vector_duplicate(intens[i]);
598 errors_local[i] = cpl_vector_duplicate(errors[i]);
599 wl_local[i] = cpl_vector_duplicate(wl[i]);
601 xmin = cpl_vector_new(n);
602 xmax = cpl_vector_new(n);
605 errors_local, n, &xmin, &xmax) == -1){
606 cpl_msg_error(__func__,
607 "Could not resample polarimetry spectra to the same wavelength scale");
608 for (i = 0; i < n; i++)
610 cpl_vector_delete(intens_local[i]);
611 cpl_vector_delete(errors_local[i]);
612 cpl_vector_delete(wl_local[i]);
614 cpl_free(intens_local);
615 cpl_free(errors_local);
617 cpl_vector_delete(xmin);
618 cpl_vector_delete(xmax);
622 result = cpl_bivector_new(size);
623 outspec = cpl_bivector_get_x(result);
624 outerr = cpl_bivector_get_y(result);
627 for (i = 0; i < size; i++) {
628 cpl_vector_set(outspec, i, 0.);
629 cpl_vector_set(outerr, i, 0.);
633 for (i=0 ; i<n ; i++) {
635 cpl_vector_copy(outspec, intens[i]);
636 cpl_vector_copy(outerr, errors[i]);
637 cpl_vector_power(outerr, 2.0);
639 cpl_vector_add(outspec, intens[i]);
640 tmp = cpl_vector_duplicate(errors[i]);
641 cpl_vector_power(tmp, 2.0);
642 cpl_vector_add(outerr, tmp);
643 cpl_vector_delete(tmp);
649 pouterr = cpl_vector_get_data(outerr) ;
650 for (i=0 ; i<size ; i++) {
651 if (isnan(pouterr[i]) || pouterr[i] < 0.0) {
656 if (ncorrections > 20)
657 cpl_msg_warning(__func__,
658 "The Errors vector contained %d negative values",
661 cpl_vector_power(outerr, 0.5);
663 for (i = 0; i < n; i++){
664 cpl_vector_delete(intens_local[i]);
665 cpl_vector_delete(errors_local[i]);
666 cpl_vector_delete(wl_local[i]);
668 cpl_free(intens_local);
669 cpl_free(errors_local);
671 cpl_vector_delete(xmin);
672 cpl_vector_delete(xmax);
674 if (cpl_error_get_code() != CPL_ERROR_NONE) {
675 cpl_msg_error(__func__,
"Error message: %s", cpl_error_get_message());
676 cpl_bivector_delete(result);
699 cpl_bivector ** stokes,
700 cpl_bivector ** null,
701 cpl_bivector ** intens,
708 int all_null, i, nrows ;
711 if (orders==NULL || wl==NULL || stokes==NULL || null==NULL || intens==NULL)
716 for (i=0 ; i<norders ; i++)
717 if (wl[i]!=NULL && stokes[i]!=NULL && null[i]!=NULL && intens[i]!=NULL){
718 nrows = cpl_vector_get_size(wl[i]) ;
719 if (cpl_bivector_get_size(stokes[i]) != nrows ||
720 cpl_bivector_get_size(null[i]) != nrows ||
721 cpl_bivector_get_size(intens[i]) != nrows) {
722 cpl_msg_error(__func__,
"Invalid Input Sizes") ;
727 if (all_null == 1)
return NULL ;
730 out = cpl_table_new(nrows);
731 for (i=0 ; i<norders ; i++) {
732 if (wl[i]!=NULL && stokes[i]!=NULL && null[i]!=NULL && intens[i]!=NULL){
735 cpl_table_new_column(out, col_name, CPL_TYPE_DOUBLE);
740 cpl_table_new_column(out, col_name, CPL_TYPE_DOUBLE);
745 cpl_table_new_column(out, col_name, CPL_TYPE_DOUBLE);
750 cpl_table_new_column(out, col_name, CPL_TYPE_DOUBLE);
755 cpl_table_new_column(out, col_name, CPL_TYPE_DOUBLE);
760 cpl_table_new_column(out, col_name, CPL_TYPE_DOUBLE);
765 cpl_table_new_column(out, col_name, CPL_TYPE_DOUBLE);
770 for (i=0 ; i<norders ; i++) {
771 if (wl[i]!=NULL && stokes[i]!=NULL && null[i]!=NULL && intens[i]!=NULL){
773 px = cpl_vector_get_data_const(wl[i]) ;
775 cpl_table_copy_data_double(out, col_name, px) ;
779 px = cpl_bivector_get_x_data_const(stokes[i]) ;
781 cpl_table_copy_data_double(out, col_name, px) ;
784 py = cpl_bivector_get_y_data_const(stokes[i]) ;
786 cpl_table_copy_data_double(out, col_name, py) ;
790 px = cpl_bivector_get_x_data_const(null[i]) ;
792 cpl_table_copy_data_double(out, col_name, px) ;
795 py = cpl_bivector_get_y_data_const(null[i]) ;
797 cpl_table_copy_data_double(out, col_name, py) ;
801 px = cpl_bivector_get_x_data_const(intens[i]) ;
803 cpl_table_copy_data_double(out, col_name, px) ;
806 py = cpl_bivector_get_y_data_const(intens[i]) ;
808 cpl_table_copy_data_double(out, col_name, py) ;
971 const cpl_table *tw_in,
972 cr2res_decker decker_position,
976 cpl_table * tw_out = NULL;
977 cpl_polynomial * wl_poly;
978 cpl_polynomial * corr_poly;
979 double wl, y_corr, newSF;
989 wl = cpl_polynomial_eval_1d(wl_poly, 1024, NULL);
990 cpl_polynomial_delete(wl_poly);
995 corr_poly = cpl_polynomial_new(1);
997 band=cpl_sprintf(
"Y");
1004 if (decker_position == CR2RES_DECKER_2_4){
1005 if (up_or_down == 1){
1006 cpl_polynomial_set_coeff(corr_poly, &pow0, -22.78903764);
1007 cpl_polynomial_set_coeff(corr_poly, &pow1, 0.03988915);
1009 cpl_polynomial_set_coeff(corr_poly, &pow0, -6.51155269);
1010 cpl_polynomial_set_coeff(corr_poly, &pow1, -0.04965059);
1012 }
else if (decker_position == CR2RES_DECKER_1_3){
1013 if (up_or_down == 1){
1014 cpl_polynomial_set_coeff(corr_poly, &pow0, 13.12469913);
1015 cpl_polynomial_set_coeff(corr_poly, &pow1, 0.04808708);
1017 cpl_polynomial_set_coeff(corr_poly, &pow0, 28.22704312);
1018 cpl_polynomial_set_coeff(corr_poly, &pow1, -0.04019722);
1022 else if (wl < 1360) {
1023 band=cpl_sprintf(
"J");
1030 if (decker_position == CR2RES_DECKER_2_4){
1031 if (up_or_down == 1){
1032 cpl_polynomial_set_coeff(corr_poly, &pow0, -26.55176271);
1033 cpl_polynomial_set_coeff(corr_poly, &pow1, 0.03903698);
1035 cpl_polynomial_set_coeff(corr_poly, &pow0, -10.87508078);
1036 cpl_polynomial_set_coeff(corr_poly, &pow1, -0.047482);
1038 }
else if (decker_position == CR2RES_DECKER_1_3){
1039 if (up_or_down == 1){
1040 cpl_polynomial_set_coeff(corr_poly, &pow0, 10.27999294);
1041 cpl_polynomial_set_coeff(corr_poly, &pow1, 0.04496791);
1043 cpl_polynomial_set_coeff(corr_poly, &pow0, 25.44306868);
1044 cpl_polynomial_set_coeff(corr_poly, &pow1, -0.04101687);
1049 else if (wl < 1850) {
1050 band=cpl_sprintf(
"H");
1057 if (decker_position == CR2RES_DECKER_2_4){
1058 if (up_or_down == 1){
1059 cpl_polynomial_set_coeff(corr_poly, &pow0, -2.17975984e+01);
1060 cpl_polynomial_set_coeff(corr_poly, &pow1, 1.59902724e-02);
1062 cpl_polynomial_set_coeff(corr_poly, &pow0, -23.23946728);
1063 cpl_polynomial_set_coeff(corr_poly, &pow1, -0.02831508);
1065 }
else if (decker_position == CR2RES_DECKER_1_3){
1066 if (up_or_down == 1){
1067 cpl_polynomial_set_coeff(corr_poly, &pow0, 17.44054878);
1068 cpl_polynomial_set_coeff(corr_poly, &pow1, 0.0199041);
1070 cpl_polynomial_set_coeff(corr_poly, &pow0, 18.31347609);
1071 cpl_polynomial_set_coeff(corr_poly, &pow1, -0.02552662);
1075 else if (wl < 2600) {
1076 band=cpl_sprintf(
"K");
1083 if (decker_position == CR2RES_DECKER_2_4){
1084 if (up_or_down == 1){
1085 cpl_polynomial_set_coeff(corr_poly, &pow0, -2.30870052e+01);
1086 cpl_polynomial_set_coeff(corr_poly, &pow1, 2.24337404e-02);
1088 cpl_polynomial_set_coeff(corr_poly, &pow0, -13.55392217);
1089 cpl_polynomial_set_coeff(corr_poly, &pow1, -0.02667301);
1091 }
else if (decker_position == CR2RES_DECKER_1_3){
1092 if (up_or_down == 1){
1093 cpl_polynomial_set_coeff(corr_poly, &pow0, 16.43404529);
1094 cpl_polynomial_set_coeff(corr_poly, &pow1, 0.02467974);
1096 cpl_polynomial_set_coeff(corr_poly, &pow0, 2.58313952e+01);
1097 cpl_polynomial_set_coeff(corr_poly, &pow1, -2.42857387e-02);
1102 cpl_msg_error(__func__,
"No proper WL found to identify band.");
1103 cpl_polynomial_delete(corr_poly);
1106 cpl_msg_info(__func__,
"Found wl=%g nm, therefore applying beam "
1107 "correction for %s-band.", wl, band);
1111 nb_traces = cpl_table_get_nrow(tw_in) ;
1112 newSFs = cpl_vector_new(nb_traces);
1113 for (i = 0; i < nb_traces; i++)
1116 cpl_table_get(tw_in, CR2RES_COL_ORDER, i, NULL) ;
1117 trace_id = cpl_table_get(tw_in, CR2RES_COL_TRACENB, i, NULL);
1119 cpl_msg_error(__func__,
"More than one input-trace per order"
1120 "is not supported.");
1121 cpl_table_delete(tw_out);
1123 cpl_polynomial_delete(corr_poly);
1128 y_corr = cpl_polynomial_eval_1d(corr_poly, wl, NULL);
1129 cpl_vector_set(newSFs, i, 0.5 + (y_corr) / CR2RES_APPROX_SLIT_HEIGHT);
1131 newSF = cpl_vector_get_median(newSFs);
1132 cpl_vector_delete(newSFs);
1133 cpl_msg_info(__func__,
"Median new slit-fraction: %g", newSF);
1135 tmp_arr = cpl_array_new(3,CPL_TYPE_DOUBLE);
1136 cpl_array_set(tmp_arr, 0, newSF-0.15);
1137 cpl_array_set(tmp_arr, 1, newSF);
1138 cpl_array_set(tmp_arr, 2, newSF+0.15);
1141 cpl_array_delete(tmp_arr);
1143 for (i=0 ; i<nb_traces ; i++) {
1144 double slope_corr, halfSF;
1146 o = cpl_table_get(tw_in, CR2RES_COL_ORDER, i, NULL) ;
1151 y_corr = cpl_polynomial_eval_1d(corr_poly, wl, NULL);
1152 tmp_arr = cpl_array_new(3,CPL_TYPE_DOUBLE);
1153 newSF = 0.5 + (y_corr / CR2RES_APPROX_SLIT_HEIGHT);
1154 halfSF = CR2RES_POLARIMETRY_DEFAULT_HEIGHT /
1155 CR2RES_APPROX_SLIT_HEIGHT / 2;
1156 cpl_array_set(tmp_arr, 0, newSF-halfSF);
1157 cpl_array_set(tmp_arr, 1, newSF);
1158 cpl_array_set(tmp_arr, 2, newSF+halfSF);
1159 cpl_table_set_array(tw_out, CR2RES_COL_SLIT_FRACTION, i, tmp_arr);
1160 cpl_array_delete(tmp_arr);
1164 wl = cpl_polynomial_eval_1d(wl_poly, 1, NULL);
1165 y_corr = cpl_polynomial_eval_1d(corr_poly, wl, NULL);
1166 wl = cpl_polynomial_eval_1d(wl_poly, CR2RES_DETECTOR_SIZE, NULL);
1167 slope_corr = (cpl_polynomial_eval_1d(corr_poly, wl, NULL)
1168 - y_corr) / CR2RES_DETECTOR_SIZE;
1170 cpl_msg_debug(__func__,
1171 "newSF, halfSF, slope_corr, y_corr: "
1173 newSF, halfSF, slope_corr, y_corr);
1176 tmp_arr = cpl_array_duplicate (
1177 cpl_table_get_array(tw_in, CR2RES_COL_ALL,i));
1178 cpl_array_set(tmp_arr, 0, cpl_array_get(tmp_arr, 0, NULL) + y_corr);
1179 cpl_array_set(tmp_arr, 1, cpl_array_get(tmp_arr, 1, NULL) + slope_corr);
1180 cpl_table_set_array(tw_out, CR2RES_COL_ALL, i, tmp_arr);
1183 cpl_array_set(tmp_arr, 0, cpl_array_get(tmp_arr, 0, NULL) +
1184 CR2RES_POLARIMETRY_DEFAULT_HEIGHT / 2);
1185 cpl_table_set_array(tw_out, CR2RES_COL_UPPER, i, tmp_arr);
1186 cpl_array_set(tmp_arr, 0, cpl_array_get(tmp_arr, 0, NULL) -
1187 CR2RES_POLARIMETRY_DEFAULT_HEIGHT);
1188 cpl_table_set_array(tw_out, CR2RES_COL_LOWER, i, tmp_arr);
1189 cpl_array_delete(tmp_arr);
1191 cpl_polynomial_delete(wl_poly);
1196 cpl_polynomial_delete(corr_poly);