CR2RE Pipeline Reference Manual 1.7.0
irplib_wlxcorr.c
1/*
2 * This file is part of the IRPLIB package
3 * Copyright (C) 2002,2003,2014 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 "irplib_wavecal_impl.h"
29
30#include "irplib_wlxcorr.h"
31
32#include <cpl.h>
33
34#include <math.h>
35#include <string.h>
36
37/*----------------------------------------------------------------------------*/
47/*----------------------------------------------------------------------------*/
48
49/*-----------------------------------------------------------------------------
50 Defines
51 -----------------------------------------------------------------------------*/
52
53#ifndef inline
54#define inline /* inline */
55#endif
56
57#define IRPLIB_MAX(A, B) ((A) > (B) ? (A) : (B))
58#define IRPLIB_MIN(A, B) ((A) < (B) ? (A) : (B))
59
60#define IRPLIB_PTR_SWAP(a, b) \
61 do { \
62 void *irplib_ptr_swap = (a); \
63 (a) = (b); \
64 (b) = irplib_ptr_swap; \
65 } while (0)
66
67/*-----------------------------------------------------------------------------
68 Private functions
69 -----------------------------------------------------------------------------*/
70
71static void irplib_wlxcorr_estimate(cpl_vector *,
72 cpl_vector *,
73 const cpl_vector *,
74 const cpl_bivector *,
75 const cpl_vector *,
76 const cpl_polynomial *,
77 double,
78 double);
79
80static int irplib_wlxcorr_signal_resample(cpl_vector *,
81 const cpl_vector *,
82 const cpl_bivector *);
83static cpl_error_code
84cpl_vector_fill_lss_profile_symmetric(cpl_vector *, double, double);
85static cpl_error_code irplib_wlcalib_fill_spectrum(cpl_vector *,
86 const cpl_bivector *,
87 const cpl_vector *,
88 const cpl_polynomial *,
89 int);
90
91static cpl_boolean irplib_wlcalib_is_lines(const cpl_vector *,
92 const cpl_polynomial *,
93 int,
94 double);
95
97
98
99/*----------------------------------------------------------------------------*/
135/*----------------------------------------------------------------------------*/
136cpl_polynomial *
137irplib_wlxcorr_best_poly(const cpl_vector *spectrum,
138 const cpl_bivector *lines_catalog,
139 int degree,
140 const cpl_polynomial *guess_poly,
141 const cpl_vector *wl_error,
142 int nsamples,
143 double slitw,
144 double fwhm,
145 double *xc,
146 cpl_table **wlres,
147 cpl_vector **xcorrs)
148{
149 const int spec_sz = cpl_vector_get_size(spectrum);
150 const int nfree = cpl_vector_get_size(wl_error);
151 int ntests = 1;
152 cpl_vector *model;
153 cpl_vector *vxc;
154 cpl_vector *init_pts_wl;
155 cpl_matrix *init_pts_x;
156 cpl_vector *pts_wl;
157 cpl_vector *vxcorrs;
158 cpl_vector *conv_kernel = NULL;
159 cpl_polynomial *poly_sol;
160 cpl_polynomial *poly_candi;
161 const double *pwl_error = cpl_vector_get_data_const(wl_error);
162 const double *dxc;
163 cpl_size degree_loc;
164 const cpl_boolean symsamp = CPL_TRUE; /* init_pts_x is symmetric */
165 const cpl_boolean is_lines =
166 irplib_wlcalib_is_lines(cpl_bivector_get_x_const(lines_catalog),
167 guess_poly, spec_sz, 1.0);
168 int i;
169
170 /* FIXME: Need mode parameter for catalogue type (lines <=> profile) */
171
172 /* In case of failure */
173 if (wlres != NULL)
174 *wlres = NULL;
175 if (xcorrs != NULL)
176 *xcorrs = NULL;
177
178 /* Useful for knowing if resampling is used */
179 cpl_msg_debug(cpl_func,
180 "Checking %d^%d dispersion polynomials (slitw=%g, "
181 "fwhm=%g) against %d-point observed spectrum with%s "
182 "catalog resampling",
183 nsamples, nfree, slitw, fwhm, spec_sz, is_lines ? "out" : "");
184
185 cpl_ensure(xc != NULL, CPL_ERROR_NULL_INPUT, NULL);
186 *xc = -1.0;
187 cpl_ensure(spectrum != NULL, CPL_ERROR_NULL_INPUT, NULL);
188 cpl_ensure(lines_catalog != NULL, CPL_ERROR_NULL_INPUT, NULL);
189 cpl_ensure(guess_poly != NULL, CPL_ERROR_NULL_INPUT, NULL);
190 cpl_ensure(wl_error != NULL, CPL_ERROR_NULL_INPUT, NULL);
191 cpl_ensure(nfree >= 2, CPL_ERROR_ILLEGAL_INPUT, NULL);
192 cpl_ensure(nsamples > 0, CPL_ERROR_ILLEGAL_INPUT, NULL);
193 /* FIXME: degree is redundant */
194 cpl_ensure(1 + degree == nfree, CPL_ERROR_ILLEGAL_INPUT, NULL);
195
196 cpl_ensure(cpl_polynomial_get_dimension(guess_poly) == 1,
197 CPL_ERROR_ILLEGAL_INPUT, NULL);
198
199 if (nsamples > 1) {
200 /* Search place must consist of more than one point */
201 /* FIXME: The bounds should probably not be negative */
202 for (i = 0; i < nfree; i++) {
203 if (pwl_error[i] != 0.0)
204 break;
205 }
206 cpl_ensure(i < nfree, CPL_ERROR_ILLEGAL_INPUT, NULL);
207 }
208
209 if (!is_lines) {
210 /* Create the convolution kernel */
211 conv_kernel = irplib_wlxcorr_convolve_create_kernel(slitw, fwhm);
212 cpl_ensure(conv_kernel != NULL, CPL_ERROR_ILLEGAL_INPUT, NULL);
213 }
214
215 /* Create initial test points */
216 init_pts_x = cpl_matrix_new(1, nfree);
217 init_pts_wl = cpl_vector_new(nfree);
218 pts_wl = cpl_vector_new(nfree);
219 for (i = 0; i < nfree; i++) {
220 const double xpos = spec_sz * i / (double)degree;
221 const double wlpos =
222 cpl_polynomial_eval_1d(guess_poly, xpos, NULL) - 0.5 * pwl_error[i];
223
224 cpl_matrix_set(init_pts_x, 0, i, xpos);
225 cpl_vector_set(init_pts_wl, i, wlpos);
226
227 ntests *= nsamples; /* Count number of tests */
228 }
229
230 vxcorrs = xcorrs != NULL ? cpl_vector_new(ntests) : NULL;
231
232 poly_sol = cpl_polynomial_new(1);
233 poly_candi = cpl_polynomial_new(1);
234 model = cpl_vector_new(spec_sz);
235 vxc = cpl_vector_new(1);
236 dxc = cpl_vector_get_data_const(vxc);
237
238 /* Create the polynomial candidates and estimate them */
239 for (i = 0; i < ntests; i++) {
240 int idiv = i;
241 int deg;
242
243 /* Update wavelength at one anchor point - and reset wavelengths
244 to their default for any anchor point(s) at higher wavelengths */
245 for (deg = degree; deg >= 0; deg--, idiv /= nsamples) {
246 const int imod = idiv % nsamples;
247 const double wlpos = cpl_vector_get(init_pts_wl, deg) +
248 imod * pwl_error[deg] / nsamples;
249
250 /* FIXME: If wlpos causes pts_wl to be non-increasing, the
251 solution will be non-physical with no need for evaluation.
252 (*xc could be set to -1 in this case). */
253 cpl_vector_set(pts_wl, deg, wlpos);
254
255 if (imod > 0)
256 break;
257 }
258
259 /* Generate */
260 degree_loc = (cpl_size)degree;
261 cpl_polynomial_fit(poly_candi, init_pts_x, &symsamp, pts_wl, NULL,
262 CPL_FALSE, NULL, &degree_loc);
263 /* *** Estimate *** */
264 irplib_wlxcorr_estimate(vxc, model, spectrum, lines_catalog,
265 conv_kernel, poly_candi, slitw, fwhm);
266 if (vxcorrs != NULL)
267 cpl_vector_set(vxcorrs, i, *dxc);
268 if (*dxc > *xc) {
269 /* Found a better solution */
270 *xc = *dxc;
271 IRPLIB_PTR_SWAP(poly_sol, poly_candi);
272 }
273 }
274
275 cpl_vector_delete(model);
276 cpl_vector_delete(vxc);
277 cpl_vector_delete(conv_kernel);
278 cpl_vector_delete(pts_wl);
279 cpl_matrix_delete(init_pts_x);
280 cpl_vector_delete(init_pts_wl);
281 cpl_polynomial_delete(poly_candi);
282
283#ifdef CPL_WLCALIB_FAIL_ON_CONSTANT
284 /* FIXME: */
285 if (cpl_polynomial_get_degree(poly_sol) == 0) {
286 cpl_polynomial_delete(poly_sol);
287 cpl_vector_delete(vxcorrs);
288 *xc = 0.0;
289 cpl_error_set_message_macro(cpl_func, CPL_ERROR_ILLEGAL_OUTPUT,
290 __FILE__, __LINE__,
291 "Found a constant "
292 "dispersion");
293 cpl_errorstate_dump(prestate, CPL_FALSE, NULL);
294 return NULL;
295 }
296#endif
297
298 if (wlres != NULL) {
299 /* FIXME: A failure in the table creation is not considered a failure
300 of the whole function call (although all outputs may be useless) */
301
302 cpl_errorstate prestate = cpl_errorstate_get();
303 /* Create the spc_table */
304 *wlres = irplib_wlxcorr_gen_spc_table(spectrum, lines_catalog, slitw,
305 fwhm, guess_poly, poly_sol);
306 if (*wlres == NULL) {
307 cpl_polynomial_delete(poly_sol);
308 cpl_vector_delete(vxcorrs);
309 *xc = -1.0;
310 cpl_error_set_message_macro(cpl_func, CPL_ERROR_ILLEGAL_OUTPUT,
311 __FILE__, __LINE__,
312 "Cannot generate "
313 "infos table");
314 /* cpl_errorstate_dump(prestate, CPL_FALSE, NULL); */
315 cpl_errorstate_set(prestate);
316 return NULL;
317 }
318 }
319
320 if (xcorrs != NULL) {
321 *xcorrs = vxcorrs;
322 }
323 else {
324 /* assert(vxcorrs == NULL); */
325 }
326
327 return poly_sol;
328}
329
330/*
331TODO :
332 Merge irplib_wlxcorr_best_poly_prop() with irplib_wlxcorr_best_poly() by
333 adding a new parameter.
334 Need to coordinate with all pipelines.
335*/
336
337/*----------------------------------------------------------------------------*/
363/*----------------------------------------------------------------------------*/
364cpl_polynomial *
365irplib_wlxcorr_best_poly_prop(const cpl_vector *spectrum,
366 const cpl_bivector *lines_catalog,
367 int degree,
368 const cpl_polynomial *guess_poly,
369 const cpl_vector *wl_error,
370 int nsamples,
371 double slitw,
372 double fwhm,
373 double *xc,
374 cpl_table **wlres,
375 cpl_vector **xcorrs)
376{
377 const int spec_sz = cpl_vector_get_size(spectrum);
378 const int nfree = cpl_vector_get_size(wl_error);
379 int ntests = 1;
380 cpl_vector *model;
381 cpl_vector *vxc;
382 cpl_vector *init_pts_wl;
383 cpl_matrix *init_pts_x;
384 cpl_vector *pts_wl;
385 cpl_vector *vxcorrs;
386 cpl_vector *conv_kernel = NULL;
387 cpl_polynomial *poly_sol;
388 cpl_polynomial *poly_candi;
389 const double *pwl_error = cpl_vector_get_data_const(wl_error);
390 const double *dxc;
391 cpl_size degree_loc;
392 const cpl_boolean symsamp = CPL_TRUE; /* init_pts_x is symmetric */
393 const cpl_boolean is_lines =
394 irplib_wlcalib_is_lines(cpl_bivector_get_x_const(lines_catalog),
395 guess_poly, spec_sz, 1.0);
396 int i;
397
398 /* FIXME: Need mode parameter for catalogue type (lines <=> profile) */
399
400 /* In case of failure */
401 if (wlres != NULL)
402 *wlres = NULL;
403 if (xcorrs != NULL)
404 *xcorrs = NULL;
405
406 /* Useful for knowing if resampling is used */
407 cpl_msg_debug(cpl_func,
408 "Checking %d^%d dispersion polynomials (slitw=%g, "
409 "fwhm=%g) against %d-point observed spectrum with%s "
410 "catalog resampling",
411 nsamples, nfree, slitw, fwhm, spec_sz, is_lines ? "out" : "");
412
413 cpl_ensure(xc != NULL, CPL_ERROR_NULL_INPUT, NULL);
414 *xc = -1.0;
415 cpl_ensure(spectrum != NULL, CPL_ERROR_NULL_INPUT, NULL);
416 cpl_ensure(lines_catalog != NULL, CPL_ERROR_NULL_INPUT, NULL);
417 cpl_ensure(guess_poly != NULL, CPL_ERROR_NULL_INPUT, NULL);
418 cpl_ensure(wl_error != NULL, CPL_ERROR_NULL_INPUT, NULL);
419 cpl_ensure(nfree >= 1, CPL_ERROR_ILLEGAL_INPUT, NULL);
420 cpl_ensure(nsamples > 0, CPL_ERROR_ILLEGAL_INPUT, NULL);
421 /* FIXME: degree is redundant */
422 cpl_ensure(1 + degree == nfree, CPL_ERROR_ILLEGAL_INPUT, NULL);
423
424 cpl_ensure(cpl_polynomial_get_dimension(guess_poly) == 1,
425 CPL_ERROR_ILLEGAL_INPUT, NULL);
426
427 if (nsamples > 1) {
428 /* Search place must consist of more than one point */
429 /* FIXME: The bounds should probably not be negative */
430 for (i = 0; i < nfree; i++) {
431 if (pwl_error[i] != 0.0)
432 break;
433 }
434 cpl_ensure(i < nfree, CPL_ERROR_ILLEGAL_INPUT, NULL);
435 }
436
437 if (!is_lines) {
438 /* Create the convolution kernel */
439 conv_kernel = irplib_wlxcorr_convolve_create_kernel(slitw, fwhm);
440 cpl_ensure(conv_kernel != NULL, CPL_ERROR_ILLEGAL_INPUT, NULL);
441 }
442
443 /* Create initial test points */
444 init_pts_x = cpl_matrix_new(1, nfree);
445 init_pts_wl = cpl_vector_new(nfree);
446 pts_wl = cpl_vector_new(nfree);
447 const double degree_loc2 = degree == 0 ? 1 : (double)degree;
448 for (i = 0; i < nfree; i++) {
449 const double xpos = spec_sz * i / degree_loc2;
450 const double wlpos =
451 cpl_polynomial_eval_1d(guess_poly, xpos, NULL) - 0.5 * pwl_error[i];
452
453 cpl_matrix_set(init_pts_x, 0, i, xpos);
454 cpl_vector_set(init_pts_wl, i, wlpos);
455
456 ntests *= nsamples; /* Count number of tests */
457 }
458
459 vxcorrs = xcorrs != NULL ? cpl_vector_new(ntests) : NULL;
460
461 poly_sol = cpl_polynomial_new(1);
462 poly_candi = cpl_polynomial_new(1);
463 model = cpl_vector_new(spec_sz);
464 vxc = cpl_vector_new(1);
465 dxc = cpl_vector_get_data_const(vxc);
466
467 /* Create the polynomial candidates and estimate them */
468 for (i = 0; i < ntests; i++) {
469 int idiv = i;
470 int deg;
471 cpl_size power;
472
473 /* Update wavelength at one anchor point - and reset wavelengths
474 to their default for any anchor point(s) at higher wavelengths */
475 for (deg = degree; deg >= 0; deg--, idiv /= nsamples) {
476 const int imod = idiv % nsamples;
477 const double wlpos = cpl_vector_get(init_pts_wl, deg) +
478 imod * pwl_error[deg] / nsamples;
479
480 /* FIXME: If wlpos causes pts_wl to be non-increasing, the
481 solution will be non-physical with no need for evaluation.
482 (*xc could be set to -1 in this case). */
483 cpl_vector_set(pts_wl, deg, wlpos);
484
485 if (imod > 0)
486 break;
487 }
488
489 /* Generate */
490 degree_loc = (cpl_size)degree;
491 cpl_polynomial_fit(poly_candi, init_pts_x, &symsamp, pts_wl, NULL,
492 CPL_FALSE, NULL, &degree_loc);
493
494 /* Use the degrees of the input guess */
495 degree_loc = cpl_polynomial_get_degree(guess_poly);
496 for (power = degree + 1; power < degree_loc + 1; power++) {
497 cpl_polynomial_set_coeff(poly_candi, &power,
498 cpl_polynomial_get_coeff(guess_poly,
499 &power));
500 }
501
502 /* *** Estimate *** */
503 irplib_wlxcorr_estimate(vxc, model, spectrum, lines_catalog,
504 conv_kernel, poly_candi, slitw, fwhm);
505 if (vxcorrs != NULL)
506 cpl_vector_set(vxcorrs, i, *dxc);
507 if (*dxc > *xc) {
508 /* Found a better solution */
509 *xc = *dxc;
510 IRPLIB_PTR_SWAP(poly_sol, poly_candi);
511 }
512 }
513
514 /* This trunctates poly_sol to degree, not wanted for propagation */
515 /* degree_loc = cpl_polynomial_get_degree(guess_poly);
516 for (cpl_size power = degree_loc; power > degree; power--){
517 cpl_polynomial_set_coeff(poly_sol, &power, 0);
518 }
519 */
520
521 cpl_vector_delete(model);
522 cpl_vector_delete(vxc);
523 cpl_vector_delete(conv_kernel);
524 cpl_vector_delete(pts_wl);
525 cpl_matrix_delete(init_pts_x);
526 cpl_vector_delete(init_pts_wl);
527 cpl_polynomial_delete(poly_candi);
528
529#ifdef CPL_WLCALIB_FAIL_ON_CONSTANT
530 /* FIXME: */
531 if (cpl_polynomial_get_degree(poly_sol) == 0) {
532 cpl_polynomial_delete(poly_sol);
533 cpl_vector_delete(vxcorrs);
534 *xc = 0.0;
535 cpl_error_set_message_macro(cpl_func, CPL_ERROR_ILLEGAL_OUTPUT,
536 __FILE__, __LINE__,
537 "Found a constant "
538 "dispersion");
539 cpl_errorstate_dump(prestate, CPL_FALSE, NULL);
540 return NULL;
541 }
542#endif
543
544 if (wlres != NULL) {
545 /* FIXME: A failure in the table creation is not considered a failure
546 of the whole function call (although all outputs may be useless) */
547
548 cpl_errorstate prestate = cpl_errorstate_get();
549 /* Create the spc_table */
550 *wlres = irplib_wlxcorr_gen_spc_table(spectrum, lines_catalog, slitw,
551 fwhm, guess_poly, poly_sol);
552 if (*wlres == NULL) {
553 cpl_polynomial_delete(poly_sol);
554 cpl_vector_delete(vxcorrs);
555 *xc = -1.0;
556 cpl_error_set_message_macro(cpl_func, CPL_ERROR_ILLEGAL_OUTPUT,
557 __FILE__, __LINE__,
558 "Cannot generate "
559 "infos table");
560 /* cpl_errorstate_dump(prestate, CPL_FALSE, NULL); */
561 cpl_errorstate_set(prestate);
562 return NULL;
563 }
564 }
565
566 if (xcorrs != NULL) {
567 *xcorrs = vxcorrs;
568 }
569 else {
570 /* assert(vxcorrs == NULL); */
571 }
572
573 return poly_sol;
574}
575
576/*----------------------------------------------------------------------------*/
594/*----------------------------------------------------------------------------*/
595cpl_table *
596irplib_wlxcorr_gen_spc_table(const cpl_vector *spectrum,
597 const cpl_bivector *lines_catalog,
598 double slitw,
599 double fwhm,
600 const cpl_polynomial *guess_poly,
601 const cpl_polynomial *corr_poly)
602{
603 cpl_vector *conv_kernel = NULL;
604 cpl_bivector *gen_init;
605 cpl_bivector *gen_corr;
606 cpl_table *spc_table;
607 const double *pgen;
608 const double xtrunc = 0.5 * slitw + 5.0 * fwhm * CPL_MATH_SIG_FWHM;
609 const int spec_sz = cpl_vector_get_size(spectrum);
610 const cpl_boolean guess_resamp =
611 !irplib_wlcalib_is_lines(cpl_bivector_get_x_const(lines_catalog),
612 guess_poly, spec_sz, 1.0);
613 const cpl_boolean corr_resamp =
614 !irplib_wlcalib_is_lines(cpl_bivector_get_x_const(lines_catalog),
615 corr_poly, spec_sz, 1.0);
616 cpl_error_code error;
617
618 cpl_msg_debug(cpl_func,
619 "Table for guess dispersion polynomial (slitw=%g, "
620 "fwhm=%g) with %d-point observed spectrum with%s catalog re"
621 "sampling",
622 slitw, fwhm, spec_sz, guess_resamp ? "out" : "");
623 cpl_msg_debug(cpl_func,
624 "Table for corr. dispersion polynomial (slitw=%g, "
625 "fwhm=%g) with %d-point observed spectrum with%s catalog re"
626 "sampling",
627 slitw, fwhm, spec_sz, corr_resamp ? "out" : "");
628
629 /* Test inputs */
630 cpl_ensure(spectrum, CPL_ERROR_NULL_INPUT, NULL);
631 cpl_ensure(lines_catalog, CPL_ERROR_NULL_INPUT, NULL);
632 cpl_ensure(guess_poly, CPL_ERROR_NULL_INPUT, NULL);
633 cpl_ensure(corr_poly, CPL_ERROR_NULL_INPUT, NULL);
634
635 /* Create the convolution kernel */
636 if (guess_resamp || corr_resamp) {
637 conv_kernel = irplib_wlxcorr_convolve_create_kernel(slitw, fwhm);
638
639 if (conv_kernel == NULL) {
640 cpl_error_set_message_macro(cpl_func, CPL_ERROR_ILLEGAL_INPUT,
641 __FILE__, __LINE__,
642 "Cannot create "
643 "convolution kernel");
644 return NULL;
645 }
646 }
647
648 /* Get the emission at initial wavelengths */
649 gen_init = cpl_bivector_new(spec_sz);
650 if (guess_resamp) {
651 error = irplib_wlcalib_fill_spectrum(cpl_bivector_get_y(gen_init),
652 lines_catalog, conv_kernel,
653 guess_poly, 0);
654 }
655 else {
656 error =
657 irplib_vector_fill_line_spectrum_model(cpl_bivector_get_y(gen_init),
658 NULL, NULL, guess_poly,
659 lines_catalog, slitw, fwhm,
660 xtrunc, 0, CPL_FALSE,
661 CPL_FALSE, NULL);
662 }
663
664 if (error || cpl_vector_fill_polynomial(cpl_bivector_get_x(gen_init),
665 guess_poly, 1, 1)) {
666 cpl_vector_delete(conv_kernel);
667 cpl_bivector_delete(gen_init);
668 cpl_error_set_message_macro(cpl_func, CPL_ERROR_ILLEGAL_INPUT, __FILE__,
669 __LINE__,
670 "Cannot get the "
671 "emission spectrum");
672 return NULL;
673 }
674
675 /* Get the emission at corrected wavelengths */
676 gen_corr = cpl_bivector_new(spec_sz);
677 if (corr_resamp) {
678 error = irplib_wlcalib_fill_spectrum(cpl_bivector_get_y(gen_corr),
679 lines_catalog, conv_kernel,
680 corr_poly, 0);
681 }
682 else {
683 error =
684 irplib_vector_fill_line_spectrum_model(cpl_bivector_get_y(gen_corr),
685 NULL, NULL, corr_poly,
686 lines_catalog, slitw, fwhm,
687 xtrunc, 0, CPL_FALSE,
688 CPL_FALSE, NULL);
689 }
690
691 if (error || cpl_vector_fill_polynomial(cpl_bivector_get_x(gen_corr),
692 corr_poly, 1, 1)) {
693 cpl_vector_delete(conv_kernel);
694 cpl_bivector_delete(gen_init);
695 cpl_bivector_delete(gen_corr);
696 cpl_error_set_message_macro(cpl_func, CPL_ERROR_ILLEGAL_INPUT, __FILE__,
697 __LINE__,
698 "Cannot get the "
699 "emission spectrum");
700 return NULL;
701 }
702 cpl_vector_delete(conv_kernel);
703
704 /* Create the ouput table */
705 spc_table = cpl_table_new(spec_sz);
706 cpl_table_new_column(spc_table, IRPLIB_WLXCORR_COL_WAVELENGTH,
707 CPL_TYPE_DOUBLE);
708 cpl_table_new_column(spc_table, IRPLIB_WLXCORR_COL_CAT_INIT,
709 CPL_TYPE_DOUBLE);
710 cpl_table_new_column(spc_table, IRPLIB_WLXCORR_COL_CAT_FINAL,
711 CPL_TYPE_DOUBLE);
712 cpl_table_new_column(spc_table, IRPLIB_WLXCORR_COL_OBS, CPL_TYPE_DOUBLE);
713
714 /* Update table */
715 pgen = cpl_bivector_get_x_data_const(gen_corr);
716 cpl_table_copy_data_double(spc_table, IRPLIB_WLXCORR_COL_WAVELENGTH, pgen);
717 pgen = cpl_bivector_get_y_data_const(gen_corr);
718 cpl_table_copy_data_double(spc_table, IRPLIB_WLXCORR_COL_CAT_FINAL, pgen);
719 pgen = cpl_vector_get_data_const(spectrum);
720 cpl_table_copy_data_double(spc_table, IRPLIB_WLXCORR_COL_OBS, pgen);
721 pgen = cpl_bivector_get_y_data_const(gen_init);
722 cpl_table_copy_data_double(spc_table, IRPLIB_WLXCORR_COL_CAT_INIT, pgen);
723 cpl_bivector_delete(gen_init);
724 cpl_bivector_delete(gen_corr);
725
726 return spc_table;
727}
728
729/*----------------------------------------------------------------------------*/
741/*----------------------------------------------------------------------------*/
742cpl_bivector *
743irplib_wlxcorr_cat_extract(const cpl_bivector *lines_catalog,
744 double wave_min,
745 double wave_max)
746{
747 const int nlines = cpl_bivector_get_size(lines_catalog);
748 int wave_min_id, wave_max_id;
749 cpl_vector *sub_cat_wl;
750 cpl_vector *sub_cat_int;
751 const cpl_vector *xlines = cpl_bivector_get_x_const(lines_catalog);
752 const double *dxlines = cpl_vector_get_data_const(xlines);
753
754 cpl_ensure(lines_catalog != NULL, CPL_ERROR_NULL_INPUT, NULL);
755
756 /* Find the 1st line */
757 wave_min_id = (int)cpl_vector_find(xlines, wave_min);
758 if (wave_min_id < 0) {
759 cpl_error_set_message_macro(cpl_func, CPL_ERROR_ILLEGAL_INPUT, __FILE__,
760 __LINE__,
761 "The starting wavelength cannot be found");
762 return NULL;
763 }
764
765 /* The first line must be greater than (at least?) wave_min */
766 if (dxlines[wave_min_id] <= wave_min)
767 wave_min_id++;
768
769 /* Find the last line */
770 wave_max_id = (int)cpl_vector_find(xlines, wave_max);
771 if (wave_max_id < 0) {
772 cpl_error_set_message_macro(cpl_func, CPL_ERROR_ILLEGAL_INPUT, __FILE__,
773 __LINE__,
774 "The ending wavelength cannot be found");
775 return NULL;
776 }
777 /* The last line must be less than wave_max */
778 if (dxlines[wave_max_id] >= wave_max)
779 wave_max_id--;
780
781 /* Checking the wavelength range at this point via the indices also
782 verifies that they were not found using non-increasing wavelengths */
783 cpl_ensure(wave_min_id <= wave_max_id, CPL_ERROR_ILLEGAL_INPUT, NULL);
784
785 if (wave_min_id < 0 || wave_max_id == nlines) {
786 cpl_error_set_message_macro(cpl_func, CPL_ERROR_ILLEGAL_INPUT, __FILE__,
787 __LINE__,
788 "The %d-line catalogue "
789 "has no lines in the range %g -> %g",
790 nlines, wave_min, wave_max);
791 return NULL;
792 }
793
794 sub_cat_wl = cpl_vector_extract(xlines, wave_min_id, wave_max_id, 1);
795 sub_cat_int = cpl_vector_extract(cpl_bivector_get_y_const(lines_catalog),
796 wave_min_id, wave_max_id, 1);
797
798 return cpl_bivector_wrap_vectors(sub_cat_wl, sub_cat_int);
799}
800
801/*----------------------------------------------------------------------------*/
818/*----------------------------------------------------------------------------*/
819cpl_vector *
820irplib_wlxcorr_convolve_create_kernel(double slitw, double fwhm)
821{
822 const double sigma = fwhm * CPL_MATH_SIG_FWHM;
823 const int size = 1 + (int)(5.0 * sigma + 0.5 * slitw);
824 cpl_vector *kernel = cpl_vector_new(size);
825
826
827 if (cpl_vector_fill_lss_profile_symmetric(kernel, slitw, fwhm)) {
828 cpl_vector_delete(kernel);
829 kernel = NULL;
830 (void)cpl_error_set_where(cpl_func);
831 }
832
833 return kernel;
834}
835
836/*----------------------------------------------------------------------------*/
849/*----------------------------------------------------------------------------*/
850int
851irplib_wlxcorr_convolve(cpl_vector *smoothed, const cpl_vector *conv_kernel)
852{
853 int nsamples;
854 int ihwidth;
855 cpl_vector *raw;
856 double *psmoothe;
857 double *praw;
858 const double *psymm;
859 int i, j;
860
861 /* Test entries */
862 cpl_ensure(smoothed, CPL_ERROR_NULL_INPUT, -1);
863 cpl_ensure(conv_kernel, CPL_ERROR_NULL_INPUT, -1);
864
865 /* Initialise */
866 nsamples = cpl_vector_get_size(smoothed);
867 ihwidth = cpl_vector_get_size(conv_kernel) - 1;
868 cpl_ensure(ihwidth < nsamples, CPL_ERROR_ILLEGAL_INPUT, -1);
869 psymm = cpl_vector_get_data_const(conv_kernel);
870 psmoothe = cpl_vector_get_data(smoothed);
871
872 /* Create raw vector */
873 raw = cpl_vector_duplicate(smoothed);
874 praw = cpl_vector_get_data(raw);
875
876 /* Convolve with the symmetric function */
877 for (i = 0; i < ihwidth; i++) {
878 psmoothe[i] = praw[i] * psymm[0];
879 for (j = 1; j <= ihwidth; j++) {
880 const int k = i - j < 0 ? 0 : i - j;
881 psmoothe[i] += (praw[k] + praw[i + j]) * psymm[j];
882 }
883 }
884
885 for (i = ihwidth; i < nsamples - ihwidth; i++) {
886 psmoothe[i] = praw[i] * psymm[0];
887 for (j = 1; j <= ihwidth; j++)
888 psmoothe[i] += (praw[i - j] + praw[i + j]) * psymm[j];
889 }
890 for (i = nsamples - ihwidth; i < nsamples; i++) {
891 psmoothe[i] = praw[i] * psymm[0];
892 for (j = 1; j <= ihwidth; j++) {
893 const int k = i + j > nsamples - 1 ? nsamples - 1 : i + j;
894 psmoothe[i] += (praw[k] + praw[i - j]) * psymm[j];
895 }
896 }
897 cpl_vector_delete(raw);
898 return 0;
899}
900
901/*----------------------------------------------------------------------------*/
911/*----------------------------------------------------------------------------*/
912int
913irplib_wlxcorr_plot_solution(const cpl_polynomial *init,
914 const cpl_polynomial *comp,
915 const cpl_polynomial *sol,
916 int pix_start,
917 int pix_stop)
918{
919 int nsamples, nplots;
920 cpl_vector **vectors;
921 int i;
922
923 /* Test entries */
924 if (init == NULL || comp == NULL)
925 return -1;
926
927 /* Initialise */
928 nsamples = pix_stop - pix_start + 1;
929 if (sol != NULL)
930 nplots = 3;
931 else
932 nplots = 2;
933
934 /* Create vectors */
935 vectors = cpl_malloc((nplots + 1) * sizeof(cpl_vector *));
936 for (i = 0; i < nplots + 1; i++)
937 vectors[i] = cpl_vector_new(nsamples);
938
939 /* First plot with the lambda/pixel relation */
940 /* Fill vectors */
941 for (i = 0; i < nsamples; i++) {
942 cpl_vector_set(vectors[0], i, pix_start + i);
943 cpl_vector_set(vectors[1], i,
944 cpl_polynomial_eval_1d(init, (double)(pix_start + i),
945 NULL));
946 cpl_vector_set(vectors[2], i,
947 cpl_polynomial_eval_1d(comp, (double)(pix_start + i),
948 NULL));
949 if (sol != NULL)
950 cpl_vector_set(vectors[3], i,
951 cpl_polynomial_eval_1d(sol, (double)(pix_start + i),
952 NULL));
953 }
954
955 /* Plot */
956 cpl_plot_vectors("set grid;set xlabel 'Position (pixels)';",
957 "t '1-Initial / 2-Computed / 3-Solution' w lines", "",
958 (const cpl_vector **)vectors, nplots + 1);
959
960 /* Free vectors */
961 for (i = 0; i < nplots + 1; i++)
962 cpl_vector_delete(vectors[i]);
963 cpl_free(vectors);
964
965 /* Allocate vectors */
966 nplots--;
967 vectors = cpl_malloc((nplots + 1) * sizeof(cpl_vector *));
968 for (i = 0; i < nplots + 1; i++)
969 vectors[i] = cpl_vector_new(nsamples);
970
971 /* Second plot with the delta-lambda/pixel relation */
972 /* Fill vectors */
973 for (i = 0; i < nsamples; i++) {
974 double diff;
975 cpl_vector_set(vectors[0], i, pix_start + i);
976 diff = cpl_polynomial_eval_1d(comp, (double)(pix_start + i), NULL) -
977 cpl_polynomial_eval_1d(init, (double)(pix_start + i), NULL);
978 cpl_vector_set(vectors[1], i, diff);
979 if (sol != NULL) {
980 diff = cpl_polynomial_eval_1d(sol, (double)(pix_start + i), NULL) -
981 cpl_polynomial_eval_1d(init, (double)(pix_start + i), NULL);
982 cpl_vector_set(vectors[2], i, diff);
983 }
984 }
985
986 /* Plot */
987 if (sol == NULL) {
988 cpl_bivector *bivector;
989 bivector = cpl_bivector_wrap_vectors(vectors[0], vectors[1]);
990 cpl_plot_bivector("set grid;set xlabel 'Position (pixels)';set ylabel "
991 "'Wavelength difference';",
992 "t 'Computed-Initial wavelenth' w lines", "",
993 bivector);
994 cpl_bivector_unwrap_vectors(bivector);
995 }
996 else {
997 cpl_plot_vectors(
998 "set grid;set xlabel 'Position (pixels)';",
999 "t '1-Computed - Initial / 2--Solution - Initial' w lines", "",
1000 (const cpl_vector **)vectors, nplots + 1);
1001 }
1002
1003 /* Free vectors */
1004 for (i = 0; i < nplots + 1; i++)
1005 cpl_vector_delete(vectors[i]);
1006 cpl_free(vectors);
1007
1008 /* Return */
1009 return 0;
1010}
1011
1012/*----------------------------------------------------------------------------*/
1023/*----------------------------------------------------------------------------*/
1024int
1025irplib_wlxcorr_plot_spc_table(const cpl_table *spc_table,
1026 const char *title,
1027 int first_plotted_line,
1028 int last_plotted_line)
1029{
1030 char title_loc[1024];
1031 cpl_vector **vectors;
1032 cpl_vector **sub_vectors;
1033 cpl_vector *tmp_vec;
1034 int nsamples;
1035 double mean1, mean3;
1036 int start_ind, stop_ind, hsize_pix;
1037 int i, j;
1038
1039 /* Test entries */
1040 if (first_plotted_line > last_plotted_line)
1041 return -1;
1042 if (spc_table == NULL)
1043 return -1;
1044
1045 /* Initialise */
1046 nsamples = cpl_table_get_nrow(spc_table);
1047 hsize_pix = 10;
1048
1049 sprintf(title_loc,
1050 "t '%s - 1-Initial catalog/2-Corrected catalog/3-Observed' w lines",
1051 title);
1052 title_loc[1023] = (char)0;
1053
1054 vectors = cpl_malloc(4 * sizeof(cpl_vector *));
1055 vectors[0] = cpl_vector_wrap(
1056 nsamples, cpl_table_get_data_double((cpl_table *)spc_table,
1057 IRPLIB_WLXCORR_COL_WAVELENGTH));
1058 vectors[1] =
1059 cpl_vector_wrap(nsamples,
1060 cpl_table_get_data_double((cpl_table *)spc_table,
1061 IRPLIB_WLXCORR_COL_CAT_INIT));
1062 vectors[2] = cpl_vector_wrap(
1063 nsamples, cpl_table_get_data_double((cpl_table *)spc_table,
1064 IRPLIB_WLXCORR_COL_CAT_FINAL));
1065 vectors[3] =
1066 cpl_vector_wrap(nsamples,
1067 cpl_table_get_data_double((cpl_table *)spc_table,
1068 IRPLIB_WLXCORR_COL_OBS));
1069
1070 /* Scale the signal for a bettre display */
1071 mean1 = cpl_vector_get_mean(vectors[1]);
1072 mean3 = cpl_vector_get_mean(vectors[3]);
1073 if (fabs(mean3) > 1)
1074 cpl_vector_multiply_scalar(vectors[3], fabs(mean1 / mean3));
1075
1076 cpl_plot_vectors("set grid;set xlabel 'Wavelength (nm)';", title_loc, "",
1077 (const cpl_vector **)vectors, 4);
1078
1079 /* Unscale the signal */
1080 if (fabs(mean3) > 1)
1081 cpl_vector_multiply_scalar(vectors[3], mean3 / mean1);
1082
1083 /* Loop on the brightest lines and zoom on them */
1084 sprintf(title_loc,
1085 "t '%s - 1-Initial catalog/2-Corrected catalog/3-Observed "
1086 "(ZOOMED)' w lines",
1087 title);
1088 title_loc[1023] = (char)0;
1089 tmp_vec = cpl_vector_duplicate(vectors[2]);
1090 for (i = 0; i < last_plotted_line; i++) {
1091 double max;
1092 /* Find the brightest line */
1093 if ((max = cpl_vector_get_max(tmp_vec)) <= 0.0)
1094 break;
1095 for (j = 0; j < nsamples; j++) {
1096 if (cpl_vector_get(tmp_vec, j) == max)
1097 break;
1098 }
1099 if (j - hsize_pix < 0)
1100 start_ind = 0;
1101 else
1102 start_ind = j - hsize_pix;
1103 if (j + hsize_pix > nsamples - 1)
1104 stop_ind = nsamples - 1;
1105 else
1106 stop_ind = j + hsize_pix;
1107 for (j = start_ind; j <= stop_ind; j++)
1108 cpl_vector_set(tmp_vec, j, 0.0);
1109
1110 if (i + 1 >= first_plotted_line) {
1111 sub_vectors = cpl_malloc(4 * sizeof(cpl_vector *));
1112 sub_vectors[0] =
1113 cpl_vector_extract(vectors[0], start_ind, stop_ind, 1);
1114 sub_vectors[1] =
1115 cpl_vector_extract(vectors[1], start_ind, stop_ind, 1);
1116 sub_vectors[2] =
1117 cpl_vector_extract(vectors[2], start_ind, stop_ind, 1);
1118 sub_vectors[3] =
1119 cpl_vector_extract(vectors[3], start_ind, stop_ind, 1);
1120
1121 cpl_plot_vectors("set grid;set xlabel 'Wavelength (nm)';",
1122 title_loc, "", (const cpl_vector **)sub_vectors,
1123 4);
1124
1125 cpl_vector_delete(sub_vectors[0]);
1126 cpl_vector_delete(sub_vectors[1]);
1127 cpl_vector_delete(sub_vectors[2]);
1128 cpl_vector_delete(sub_vectors[3]);
1129 cpl_free(sub_vectors);
1130 }
1131 }
1132 cpl_vector_delete(tmp_vec);
1133
1134 cpl_vector_unwrap(vectors[0]);
1135 cpl_vector_unwrap(vectors[1]);
1136 cpl_vector_unwrap(vectors[2]);
1137 cpl_vector_unwrap(vectors[3]);
1138 cpl_free(vectors);
1139
1140 return 0;
1141}
1142
1143/*----------------------------------------------------------------------------*/
1151/*----------------------------------------------------------------------------*/
1152int
1153irplib_wlxcorr_catalog_plot(const cpl_bivector *cat, double wmin, double wmax)
1154{
1155 int start, stop;
1156 cpl_bivector *subcat;
1157 cpl_vector *subcat_x;
1158 cpl_vector *subcat_y;
1159 const double *pwave;
1160 int nvals, nvals_tot;
1161 int i;
1162
1163 /* Test entries */
1164 if (cat == NULL)
1165 return -1;
1166 if (wmax <= wmin)
1167 return -1;
1168
1169 /* Initialise */
1170 nvals_tot = cpl_bivector_get_size(cat);
1171
1172 /* Count the nb of values */
1173 pwave = cpl_bivector_get_x_data_const(cat);
1174 if (pwave[0] >= wmin)
1175 start = 0;
1176 else
1177 start = -1;
1178 if (pwave[nvals_tot - 1] <= wmax)
1179 stop = nvals_tot - 1;
1180 else
1181 stop = -1;
1182 i = 0;
1183 while ((i < nvals_tot - 1) && (pwave[i] < wmin))
1184 i++;
1185 start = i;
1186 i = nvals_tot - 1;
1187 while ((i > 0) && (pwave[i] > wmax))
1188 i--;
1189 stop = i;
1190
1191 if (start >= stop) {
1192 cpl_msg_error(cpl_func, "Cannot plot the catalog");
1193 return -1;
1194 }
1195 nvals = stop - start + 1;
1196
1197 /* Create the bivector to plot */
1198 subcat_x =
1199 cpl_vector_extract(cpl_bivector_get_x_const(cat), start, stop, 1);
1200 subcat_y =
1201 cpl_vector_extract(cpl_bivector_get_y_const(cat), start, stop, 1);
1202 subcat = cpl_bivector_wrap_vectors(subcat_x, subcat_y);
1203
1204 /* Plot */
1205 if (nvals > 500) {
1206 cpl_plot_bivector(
1207 "set grid;set xlabel 'Wavelength (nm)';set ylabel 'Emission';",
1208 "t 'Catalog Spectrum' w lines", "", subcat);
1209 }
1210 else {
1211 cpl_plot_bivector(
1212 "set grid;set xlabel 'Wavelength (nm)';set ylabel 'Emission';",
1213 "t 'Catalog Spectrum' w impulses", "", subcat);
1214 }
1215 cpl_bivector_unwrap_vectors(subcat);
1216 cpl_vector_delete(subcat_x);
1217 cpl_vector_delete(subcat_y);
1218
1219 return 0;
1220}
1221
1223
1224/*----------------------------------------------------------------------------*/
1239/*----------------------------------------------------------------------------*/
1240static void
1241irplib_wlxcorr_estimate(cpl_vector *vxc,
1242 cpl_vector *model,
1243 const cpl_vector *spectrum,
1244 const cpl_bivector *lines_catalog,
1245 const cpl_vector *conv_kernel,
1246 const cpl_polynomial *poly_candi,
1247 double slitw,
1248 double fwhm)
1249{
1250 cpl_errorstate prestate = cpl_errorstate_get();
1251 const int hsize = cpl_vector_get_size(vxc) / 2;
1252
1253 if (conv_kernel != NULL) {
1254 irplib_wlcalib_fill_spectrum(model, lines_catalog, conv_kernel,
1255 poly_candi, hsize);
1256 }
1257 else {
1258 const double xtrunc = 0.5 * slitw + 5.0 * fwhm * CPL_MATH_SIG_FWHM;
1259
1260 irplib_vector_fill_line_spectrum_model(model, NULL, NULL, poly_candi,
1261 lines_catalog, slitw, fwhm,
1262 xtrunc, 0, CPL_FALSE, CPL_FALSE,
1263 NULL);
1264 }
1265
1266 if (cpl_errorstate_is_equal(prestate))
1267 cpl_vector_correlate(vxc, model, spectrum);
1268
1269 if (!cpl_errorstate_is_equal(prestate)) {
1270 cpl_vector_fill(vxc, 0.0);
1271
1272 /* cpl_errorstate_dump(prestate, CPL_FALSE, NULL); */
1273 cpl_errorstate_set(prestate);
1274 }
1275
1276 return;
1277}
1278
1279
1280/*----------------------------------------------------------------------------*/
1290/*----------------------------------------------------------------------------*/
1291static cpl_boolean
1292irplib_wlcalib_is_lines(const cpl_vector *wavelengths,
1293 const cpl_polynomial *disp1d,
1294 int spec_sz,
1295 double tol)
1296{
1297 const int nlines = cpl_vector_get_size(wavelengths);
1298 /* The dispersion on the detector center */
1299 const double dispersion =
1300 cpl_polynomial_eval_1d_diff(disp1d, 0.5 * spec_sz + 1.0, 0.5 * spec_sz,
1301 NULL);
1302 const double range = cpl_vector_get(wavelengths, nlines - 1) -
1303 cpl_vector_get(wavelengths, 0);
1304
1305 cpl_ensure(wavelengths != NULL, CPL_ERROR_NULL_INPUT, CPL_FALSE);
1306 cpl_ensure(disp1d != NULL, CPL_ERROR_NULL_INPUT, CPL_FALSE);
1307 cpl_ensure(cpl_polynomial_get_dimension(disp1d) == 1,
1308 CPL_ERROR_ILLEGAL_INPUT, CPL_FALSE);
1309 cpl_ensure(range > 0.0, CPL_ERROR_ILLEGAL_INPUT, CPL_FALSE);
1310
1311 return nlines * fabs(dispersion) <= tol * fabs(range) ? CPL_TRUE
1312 : CPL_FALSE;
1313}
1314
1315/*----------------------------------------------------------------------------*/
1330/*----------------------------------------------------------------------------*/
1331static cpl_error_code
1332irplib_wlcalib_fill_spectrum(cpl_vector *self,
1333 const cpl_bivector *lines_catalog,
1334 const cpl_vector *conv_kernel,
1335 const cpl_polynomial *poly,
1336 int search_hs)
1337{
1338 const int size = cpl_vector_get_size(self);
1339 const int nlines = cpl_bivector_get_size(lines_catalog);
1340 const cpl_vector *xlines = cpl_bivector_get_x_const(lines_catalog);
1341 const double *dxlines = cpl_vector_get_data_const(xlines);
1342 cpl_bivector *sub_cat;
1343 cpl_vector *sub_cat_x;
1344 cpl_vector *sub_cat_y;
1345 cpl_vector *wl_limits;
1346 double wave_min, wave_max;
1347 int wave_min_id, wave_max_id;
1348 int nsub;
1349 int error;
1350
1351 cpl_ensure_code(self != NULL, CPL_ERROR_NULL_INPUT);
1352 cpl_ensure_code(lines_catalog != NULL, CPL_ERROR_NULL_INPUT);
1353 cpl_ensure_code(conv_kernel != NULL, CPL_ERROR_NULL_INPUT);
1354 cpl_ensure_code(poly != NULL, CPL_ERROR_NULL_INPUT);
1355 cpl_ensure_code(size > 0, CPL_ERROR_ILLEGAL_INPUT);
1356
1357
1358 /* Resample the spectrum */
1359 wl_limits = cpl_vector_new(size + 1);
1360 cpl_vector_fill_polynomial(wl_limits, poly, 0.5 - search_hs, 1);
1361
1362 /* The spectrum wavelength bounds */
1363 wave_min = cpl_vector_get(wl_limits, 0);
1364 wave_max = cpl_vector_get(wl_limits, size);
1365
1366 /* Find the 1st line */
1367 wave_min_id = cpl_vector_find(xlines, wave_min);
1368 /* The first line must be less than or equal to wave_min */
1369 if (dxlines[wave_min_id] > wave_min)
1370 wave_min_id--;
1371
1372 if (wave_min_id < 0) {
1373 cpl_vector_delete(wl_limits);
1374 return cpl_error_set_message_macro(cpl_func, CPL_ERROR_ILLEGAL_INPUT,
1375 __FILE__, __LINE__,
1376 "The %d-line "
1377 "catalogue only has lines above %g",
1378 nlines, wave_min);
1379 }
1380
1381 /* Find the last line */
1382 wave_max_id = cpl_vector_find(xlines, wave_max);
1383 /* The last line must be greater than or equal to wave_max */
1384 if (dxlines[wave_max_id] < wave_max)
1385 wave_max_id++;
1386
1387 if (wave_max_id == nlines) {
1388 cpl_vector_delete(wl_limits);
1389 return cpl_error_set_message_macro(cpl_func, CPL_ERROR_ILLEGAL_INPUT,
1390 __FILE__, __LINE__,
1391 "The %d-line "
1392 "catalogue only has lines below %g",
1393 nlines, wave_max);
1394 }
1395
1396 /* Checking the wavelength range at this point via the indices also
1397 verifies that they were not found using non-increasing wavelengths */
1398 nsub = 1 + wave_max_id - wave_min_id;
1399 cpl_ensure_code(nsub > 1, CPL_ERROR_ILLEGAL_INPUT);
1400
1401 /* Wrap a new bivector around the relevant part of the catalog */
1402 /* The data is _not_ modified */
1403 sub_cat_x = cpl_vector_wrap(nsub, wave_min_id + (double *)dxlines);
1404 sub_cat_y =
1405 cpl_vector_wrap(nsub,
1406 wave_min_id + (double *)cpl_bivector_get_y_data_const(
1407 lines_catalog));
1408 sub_cat = cpl_bivector_wrap_vectors(sub_cat_x, sub_cat_y);
1409
1410 /* High resolution catalog */
1411 error = irplib_wlxcorr_signal_resample(self, wl_limits, sub_cat);
1412
1413 cpl_vector_delete(wl_limits);
1414 cpl_bivector_unwrap_vectors(sub_cat);
1415 (void)cpl_vector_unwrap(sub_cat_x);
1416 (void)cpl_vector_unwrap(sub_cat_y);
1417
1418 cpl_ensure_code(!error, CPL_ERROR_ILLEGAL_INPUT);
1419
1420 /* Smooth the instrument resolution */
1421 cpl_ensure_code(!irplib_wlxcorr_convolve(self, conv_kernel),
1422 cpl_error_get_code());
1423
1424 return CPL_ERROR_NONE;
1425}
1426
1427
1428/*----------------------------------------------------------------------------*/
1438/*----------------------------------------------------------------------------*/
1439static int
1440irplib_wlxcorr_signal_resample(cpl_vector *resampled,
1441 const cpl_vector *xbounds,
1442 const cpl_bivector *hires)
1443{
1444 const int hrsize = cpl_bivector_get_size(hires);
1445 const cpl_vector *xhires;
1446 const cpl_vector *yhires;
1447 const double *pxhires;
1448 const double *pyhires;
1449 const double *pxbounds;
1450 cpl_vector *ybounds;
1451 cpl_bivector *boundary;
1452 double *pybounds;
1453 double *presampled;
1454 int nsamples;
1455 int i, itt;
1456
1457 /* Test entries */
1458 if ((!resampled) || (!xbounds) || (!hires))
1459 return -1;
1460
1461 /* Initialise */
1462 nsamples = cpl_vector_get_size(resampled);
1463
1464 /* Initialise */
1465 presampled = cpl_vector_get_data(resampled);
1466 pxbounds = cpl_vector_get_data_const(xbounds);
1467 xhires = cpl_bivector_get_x_const(hires);
1468 yhires = cpl_bivector_get_y_const(hires);
1469 pxhires = cpl_vector_get_data_const(xhires);
1470 pyhires = cpl_vector_get_data_const(yhires);
1471
1472 /* Create a new vector */
1473 ybounds = cpl_vector_new(cpl_vector_get_size(xbounds));
1474 boundary = cpl_bivector_wrap_vectors((cpl_vector *)xbounds, ybounds);
1475 pybounds = cpl_vector_get_data(ybounds);
1476
1477 /* Test entries */
1478 if (cpl_bivector_get_size(boundary) != nsamples + 1) {
1479 cpl_bivector_unwrap_vectors(boundary);
1480 cpl_vector_delete(ybounds);
1481 return -1;
1482 }
1483
1484 /* Get the ind */
1485 itt = cpl_vector_find(xhires, pxbounds[0]);
1486
1487 /* Interpolate the signal */
1488 if (cpl_bivector_interpolate_linear(boundary, hires)) {
1489 cpl_bivector_unwrap_vectors(boundary);
1490 cpl_vector_delete(ybounds);
1491 return -1;
1492 }
1493
1494 /* At this point itt most likely points to element just below
1495 pxbounds[0] */
1496 while (pxhires[itt] < pxbounds[0])
1497 itt++;
1498
1499 for (i = 0; i < nsamples; i++) {
1500 /* The i'th signal is the weighted average of the two interpolated
1501 signals at the pixel boundaries and those table signals in
1502 between */
1503
1504 double xlow = pxbounds[i];
1505 double x = pxhires[itt];
1506
1507 if (x > pxbounds[i + 1])
1508 x = pxbounds[i + 1];
1509 /* Contribution from interpolated value at wavelength at lower pixel
1510 boundary */
1511 presampled[i] = pybounds[i] * (x - xlow);
1512
1513 /* Contribution from table values in between pixel boundaries */
1514 while ((pxhires[itt] < pxbounds[i + 1]) && (itt < hrsize)) {
1515 const double xprev = x;
1516 x = pxhires[itt + 1];
1517 if (x > pxbounds[i + 1])
1518 x = pxbounds[i + 1];
1519 presampled[i] += pyhires[itt] * (x - xlow);
1520 xlow = xprev;
1521 itt++;
1522 }
1523
1524 /* Contribution from interpolated value at wavelength at upper pixel
1525 boundary */
1526 presampled[i] += pybounds[i + 1] * (pxbounds[i + 1] - xlow);
1527
1528 /* Compute average by dividing integral by length of pixel range
1529 (the factor 2 comes from the contributions) */
1530 presampled[i] /= 2 * (pxbounds[i + 1] - pxbounds[i]);
1531 }
1532 cpl_bivector_unwrap_vectors(boundary);
1533 cpl_vector_delete(ybounds);
1534 return 0;
1535}
1536
1537
1538/*----------------------------------------------------------------------------*/
1559/*----------------------------------------------------------------------------*/
1560static cpl_error_code
1561cpl_vector_fill_lss_profile_symmetric(cpl_vector *self,
1562 double slitw,
1563 double fwhm)
1564{
1565 const double sigma = fwhm * CPL_MATH_SIG_FWHM;
1566 const int n = cpl_vector_get_size(self);
1567 int i;
1568
1569
1570 cpl_ensure_code(self != NULL, CPL_ERROR_NULL_INPUT);
1571 cpl_ensure_code(slitw > 0.0, CPL_ERROR_ILLEGAL_INPUT);
1572 cpl_ensure_code(fwhm > 0.0, CPL_ERROR_ILLEGAL_INPUT);
1573
1574 /* Cannot fail now */
1575
1576 /* Special case for i = 0 */
1577 (void)cpl_vector_set(self, 0,
1578 (irplib_erf_antideriv(0.5 * slitw + 0.5, sigma) -
1579 irplib_erf_antideriv(0.5 * slitw - 0.5, sigma)) /
1580 slitw);
1581
1582 for (i = 1; i < n; i++) {
1583 /* FIXME: Reuse two irplib_erf_antideriv() calls from previous value */
1584 const double x1p = i + 0.5 * slitw + 0.5;
1585 const double x1n = i - 0.5 * slitw + 0.5;
1586 const double x0p = i + 0.5 * slitw - 0.5;
1587 const double x0n = i - 0.5 * slitw - 0.5;
1588 const double val = 0.5 / slitw *
1589 (irplib_erf_antideriv(x1p, sigma) -
1590 irplib_erf_antideriv(x1n, sigma) -
1591 irplib_erf_antideriv(x0p, sigma) +
1592 irplib_erf_antideriv(x0n, sigma));
1593 (void)cpl_vector_set(self, i, val);
1594 }
1595
1596 return CPL_ERROR_NONE;
1597}
irplib_erf_antideriv
double irplib_erf_antideriv(double x, double sigma)
The antiderivative of erx(x/sigma/sqrt(2)) with respect to x.
Definition irplib_wavecal.c:1369
irplib_vector_fill_line_spectrum_model
cpl_error_code irplib_vector_fill_line_spectrum_model(cpl_vector *self, cpl_vector *linepix, cpl_vector *erftmp, const cpl_polynomial *disp, const cpl_bivector *lines, double wslit, double wfwhm, double xtrunc, int hsize, cpl_boolean dofast, cpl_boolean dolog, cpl_size *pulines)
Generate a 1D spectrum from (arc) lines and a dispersion relation.
Definition irplib_wavecal.c:1030
Generated on Tue Oct 13 2026 12:27:26 for CR2RE Pipeline Reference Manual by doxygen 1.13.2