00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028 #ifdef HAVE_CONFIG_H
00029 #include <config.h>
00030 #endif
00031
00032
00033
00034
00035
00036 #include <math.h>
00037 #include <cpl.h>
00038
00039 #include "irplib_wlxcorr.h"
00040 #include "irplib_plot.h"
00041
00042
00047
00048
00049 static int irplib_wlxcorr_signal_resample(cpl_vector *, const cpl_vector *,
00050 const cpl_bivector *) ;
00051
00054
00074
00075 cpl_polynomial * irplib_wlxcorr_best_poly(
00076 const cpl_vector * spectrum,
00077 const cpl_bivector * lines_catalog,
00078 int degree,
00079 const cpl_polynomial * guess_poly,
00080 const cpl_vector * wl_error,
00081 int nsamples,
00082 double slitw,
00083 double fwhm,
00084 int doplot,
00085 double * xc,
00086 cpl_table ** wlres)
00087 {
00088 int ntests ;
00089 int spec_sz ;
00090 cpl_polynomial ** candidates ;
00091 double xpos, wl_pos ;
00092 cpl_vector * init_pts_wl ;
00093 cpl_vector * init_pts_x ;
00094 cpl_vector * pts_wl ;
00095 cpl_vector * vxcorrs ;
00096 int best_ind ;
00097 cpl_polynomial * poly_sol ;
00098 cpl_table * spc_table ;
00099 double xc_cur ;
00100 cpl_vector * vxc ;
00101 cpl_bivector * gen_init ;
00102 double * pwl_error ;
00103 int i, j, k, l ;
00104
00105
00106 if (!spectrum || !lines_catalog || !guess_poly || !xc || !wl_error)
00107 return NULL;
00108 if (degree < 1 || degree > 3) return NULL ;
00109 if (cpl_vector_get_size(wl_error) != degree+1) return NULL ;
00110
00111
00112 ntests = (int)pow(nsamples, (degree + 1)) ;
00113 spec_sz = cpl_vector_get_size(spectrum) ;
00114 pwl_error = cpl_vector_get_data(wl_error) ;
00115
00116
00117 if (doplot) {
00118 irplib_wlxcorr_catalog_plot(lines_catalog,
00119 cpl_polynomial_eval_1d(guess_poly, 1, NULL),
00120 cpl_polynomial_eval_1d(guess_poly, spec_sz, NULL)) ;
00121 irplib_vector_plot(
00122 "set grid;set xlabel 'Position (Pixel)';set ylabel 'Intensity (ADU/sec)';",
00123 "t 'Extracted spectrum' w lines", "", spectrum);
00124 }
00125
00126
00127 init_pts_x = cpl_vector_new(degree + 1) ;
00128 init_pts_wl = cpl_vector_new(degree + 1) ;
00129 for (i=0 ; i<degree + 1 ; i++) {
00130 xpos = i*spec_sz/degree ;
00131 cpl_vector_set(init_pts_x, i, xpos) ;
00132 cpl_vector_set(init_pts_wl, i,
00133 cpl_polynomial_eval_1d(guess_poly, xpos, NULL)) ;
00134 }
00135
00136
00137 candidates = cpl_malloc(ntests * sizeof(cpl_polynomial*)) ;
00138 pts_wl = cpl_vector_new(degree + 1) ;
00139 for (i=0 ; i<nsamples ; i++) {
00140 wl_pos = cpl_vector_get(init_pts_wl, 0)
00141 - pwl_error[0]/2 + i*pwl_error[0]/nsamples ;
00142 cpl_vector_set(pts_wl, 0, wl_pos) ;
00143 for (j=0 ; j<nsamples ; j++) {
00144 wl_pos = cpl_vector_get(init_pts_wl, 1)
00145 - pwl_error[1]/2 + j*pwl_error[1]/nsamples ;
00146 cpl_vector_set(pts_wl, 1, wl_pos) ;
00147 if (degree == 1) {
00148 candidates[j+i*nsamples] =
00149 cpl_polynomial_fit_1d_create(init_pts_x, pts_wl, degree,
00150 NULL);
00151 } else {
00152 for (k=0 ; k<nsamples ; k++) {
00153 wl_pos = cpl_vector_get(init_pts_wl, 2)
00154 - pwl_error[2]/2 + k*pwl_error[2]/nsamples ;
00155 cpl_vector_set(pts_wl, 2, wl_pos) ;
00156 if (degree == 2) {
00157 candidates[k+(j+i*nsamples)*nsamples] =
00158 cpl_polynomial_fit_1d_create(init_pts_x,
00159 pts_wl, degree, NULL);
00160 } else {
00161 for (l=0 ; l<nsamples ; l++) {
00162 wl_pos = cpl_vector_get(init_pts_wl, 3)
00163 - pwl_error[3]/2 + l*pwl_error[3]/nsamples ;
00164 cpl_vector_set(pts_wl, 3, wl_pos) ;
00165 if (degree == 3) {
00166 candidates[l+(k+(j+i*nsamples)*nsamples)*nsamples]=
00167 cpl_polynomial_fit_1d_create(init_pts_x,
00168 pts_wl, degree, NULL);
00169 } else {
00170
00171 }
00172 }
00173 }
00174 }
00175 }
00176 }
00177 }
00178 cpl_vector_delete(pts_wl) ;
00179 cpl_vector_delete(init_pts_x) ;
00180 cpl_vector_delete(init_pts_wl) ;
00181
00182
00183 best_ind = 0 ;
00184 *xc = -1 ;
00185 vxc = cpl_vector_new(1) ;
00186 vxcorrs = cpl_vector_new(ntests) ;
00187 for (i=0 ; i<ntests ; i++) {
00188
00189 if ((gen_init=irplib_wlxcorr_gen_signal(lines_catalog, slitw, fwhm,
00190 candidates[i], spec_sz, 0, NULL)) != NULL) {
00191
00192 cpl_vector_correlate(vxc, cpl_bivector_get_y(gen_init),spectrum);
00193 cpl_bivector_delete(gen_init);
00194 xc_cur = cpl_vector_get(vxc, 0);
00195 if (xc_cur > *xc) {
00196 *xc = xc_cur ;
00197 best_ind = i ;
00198 }
00199 cpl_vector_set(vxcorrs, i, xc_cur) ;
00200 } else {
00201 cpl_msg_error(cpl_func, "Cannot generate the signal - abort") ;
00202 cpl_vector_delete(vxcorrs) ;
00203 cpl_vector_delete(vxc) ;
00204 for (i=0 ; i<ntests ; i++) {
00205 cpl_polynomial_delete(candidates[i]) ;
00206 }
00207 cpl_free(candidates) ;
00208 return NULL ;
00209 }
00210 }
00211 cpl_vector_delete(vxc) ;
00212
00213
00214 if (doplot) {
00215 irplib_vector_plot("set grid;", "t 'Correlation values' w lines", "",
00216 vxcorrs) ;
00217 }
00218
00219
00220
00221 cpl_vector_sort(vxcorrs, -1) ;
00222
00223
00224
00225
00226 if (doplot) {
00227 irplib_vector_plot("set grid;", "t 'Sorted correlation values' w lines",
00228 "", vxcorrs) ;
00229 }
00230 cpl_vector_delete(vxcorrs) ;
00231
00232
00233 poly_sol = cpl_polynomial_duplicate(candidates[best_ind]) ;
00234
00235 for (i=0 ; i<ntests ; i++) {
00236 cpl_polynomial_delete(candidates[i]) ;
00237 }
00238 cpl_free(candidates) ;
00239
00240
00241 if ((spc_table = irplib_wlxcorr_gen_spc_table(spectrum, lines_catalog,
00242 slitw, fwhm, guess_poly, poly_sol)) == NULL) {
00243 cpl_msg_error(cpl_func, "Cannot generate infos table") ;
00244 cpl_polynomial_delete(poly_sol) ;
00245 return NULL ;
00246 }
00247
00248
00249 if (doplot > 0) {
00250 irplib_wlxcorr_plot_spc_table(spc_table, "XC") ;
00251 }
00252
00253 if (wlres != NULL) *wlres = spc_table ;
00254 else cpl_table_delete(spc_table) ;
00255
00256 return poly_sol ;
00257 }
00258
00259
00272
00273 cpl_table * irplib_wlxcorr_gen_spc_table(
00274 const cpl_vector * spectrum,
00275 const cpl_bivector * lines_catalog,
00276 double slitw,
00277 double fwhm,
00278 const cpl_polynomial * guess_poly,
00279 const cpl_polynomial * corr_poly)
00280 {
00281 cpl_bivector * gen_init ;
00282 cpl_bivector * gen_corr ;
00283 int nsamples ;
00284 cpl_table * spc_table ;
00285 double * pgen ;
00286
00287
00288 if (spectrum == NULL) return NULL ;
00289 if (lines_catalog == NULL) return NULL ;
00290 if (guess_poly == NULL) return NULL ;
00291 if (corr_poly == NULL) return NULL ;
00292
00293
00294 nsamples = cpl_vector_get_size(spectrum) ;
00295
00296
00297 if ((gen_init=irplib_wlxcorr_gen_signal(lines_catalog, slitw, fwhm,
00298 guess_poly, nsamples, 0, NULL)) == NULL) {
00299 cpl_msg_error(cpl_func, "Cannot get the emission spectrum") ;
00300 return NULL ;
00301 }
00302
00303
00304 if ((gen_corr=irplib_wlxcorr_gen_signal(lines_catalog, slitw, fwhm,
00305 corr_poly, nsamples, 0, NULL)) == NULL) {
00306 cpl_msg_error(cpl_func, "Cannot get the emission spectrum") ;
00307 cpl_bivector_delete(gen_init) ;
00308 return NULL ;
00309 }
00310
00311
00312 spc_table = cpl_table_new(nsamples);
00313 cpl_table_new_column(spc_table, IRPLIB_COL_XC_WAVELENGTH, CPL_TYPE_DOUBLE);
00314 cpl_table_new_column(spc_table, IRPLIB_COL_XC_CAT_INIT, CPL_TYPE_DOUBLE);
00315 cpl_table_new_column(spc_table, IRPLIB_COL_XC_CAT_FINAL, CPL_TYPE_DOUBLE);
00316 cpl_table_new_column(spc_table, IRPLIB_COL_XC_OBS, CPL_TYPE_DOUBLE);
00317
00318
00319 pgen = cpl_bivector_get_x_data(gen_corr) ;
00320 cpl_table_copy_data_double(spc_table, IRPLIB_COL_XC_WAVELENGTH, pgen) ;
00321 pgen = cpl_bivector_get_y_data(gen_corr) ;
00322 cpl_table_copy_data_double(spc_table, IRPLIB_COL_XC_CAT_FINAL, pgen) ;
00323 pgen = cpl_vector_get_data(spectrum) ;
00324 cpl_table_copy_data_double(spc_table, IRPLIB_COL_XC_OBS, pgen) ;
00325 pgen = cpl_bivector_get_y_data(gen_init) ;
00326 cpl_table_copy_data_double(spc_table, IRPLIB_COL_XC_CAT_INIT, pgen);
00327 cpl_bivector_delete(gen_init);
00328 cpl_bivector_delete(gen_corr);
00329
00330 return spc_table ;
00331 }
00332
00333
00341
00342 cpl_bivector * irplib_wlxcorr_cat_extract(
00343 const cpl_bivector * lines_catalog,
00344 double wave_min,
00345 double wave_max)
00346 {
00347 int wave_min_id, wave_max_id ;
00348 double wave_cur ;
00349 cpl_vector * sub_cat_wl ;
00350 cpl_vector * sub_cat_int ;
00351 cpl_bivector * sub_cat ;
00352 int i ;
00353
00354 wave_min_id = wave_max_id = -1 ;
00355 for (i=0 ; i<cpl_bivector_get_size(lines_catalog) ; i++) {
00356 wave_cur = cpl_vector_get(cpl_bivector_get_x(lines_catalog), i) ;
00357 if ((wave_min_id<0) && (wave_cur>wave_min)) wave_min_id = i ;
00358 if (wave_cur<wave_max) wave_max_id = i ;
00359 }
00360 if (wave_min_id<0 || wave_max_id<0) {
00361 cpl_msg_error(cpl_func, "Cannot extract from the catalog spectrum") ;
00362 return NULL ;
00363 }
00364 if (wave_min_id >= wave_max_id) return NULL ;
00365 sub_cat_wl = cpl_vector_extract(cpl_bivector_get_x(lines_catalog),
00366 wave_min_id, wave_max_id, 1) ;
00367 sub_cat_int = cpl_vector_extract(cpl_bivector_get_y(lines_catalog),
00368 wave_min_id, wave_max_id, 1) ;
00369 sub_cat = cpl_bivector_wrap_vectors(sub_cat_wl, sub_cat_int) ;
00370
00371 return sub_cat ;
00372 }
00373
00374
00390
00391 cpl_bivector * irplib_wlxcorr_gen_signal(
00392 const cpl_bivector * lines_catalog,
00393 double slitw,
00394 double fwhm,
00395 const cpl_polynomial * poly,
00396 int nsamples,
00397 int search_hs,
00398 int * nb_lines)
00399 {
00400 int size, nlines ;
00401 cpl_vector * conv_kernel ;
00402 cpl_bivector * gen_spectrum ;
00403 cpl_bivector * sub_cat ;
00404 double wave_min, wave_max ;
00405 cpl_vector * wl_limits ;
00406 int ind ;
00407 double new_val, wl ;
00408 int i ;
00409
00410
00411 size = nsamples + 2 * search_hs ;
00412
00413
00414 if (!lines_catalog || !poly) return NULL ;
00415 if (size <= 1) return NULL ;
00416
00417
00418 gen_spectrum = cpl_bivector_new(size) ;
00419 cpl_vector_fill_polynomial(cpl_bivector_get_x(gen_spectrum), poly,
00420 -search_hs+1, 1) ;
00421
00422
00423 wave_min = cpl_vector_get(cpl_bivector_get_x(gen_spectrum), 0) ;
00424 wave_max = cpl_vector_get(cpl_bivector_get_x(gen_spectrum), size-1) ;
00425 sub_cat = irplib_wlxcorr_cat_extract(lines_catalog, wave_min, wave_max) ;
00426 if (sub_cat == NULL) {
00427 cpl_bivector_delete(gen_spectrum) ;
00428 return NULL ;
00429 }
00430 nlines = cpl_bivector_get_size(sub_cat) ;
00431 if (nb_lines) *nb_lines = nlines ;
00432
00433
00434 wl_limits = cpl_vector_new(size + 1);
00435 cpl_vector_fill_polynomial(wl_limits, poly, -search_hs+0.5, 1) ;
00436 if (nlines < size) {
00437
00438 cpl_vector_fill(cpl_bivector_get_y(gen_spectrum), 0.0) ;
00439 for (i=0 ; i<nlines ; i++) {
00440 wl = cpl_vector_get(cpl_bivector_get_x(sub_cat), i) ;
00441 ind = 0 ;
00442 while (wl>cpl_vector_get(wl_limits, ind) && ind<size-1)
00443 ind++ ;
00444 new_val = cpl_vector_get(cpl_bivector_get_y(gen_spectrum), ind) ;
00445 new_val += cpl_vector_get(cpl_bivector_get_y(sub_cat), i) ;
00446 cpl_vector_set(cpl_bivector_get_y(gen_spectrum), ind, new_val) ;
00447 }
00448 } else {
00449
00450 if (irplib_wlxcorr_signal_resample(cpl_bivector_get_y(gen_spectrum),
00451 wl_limits, lines_catalog)) {
00452 cpl_msg_error(cpl_func, "Cannot resample the signal") ;
00453 cpl_bivector_delete(gen_spectrum) ;
00454 cpl_vector_delete(wl_limits) ;
00455 cpl_bivector_delete(sub_cat) ;
00456 return NULL ;
00457 }
00458 }
00459 cpl_bivector_delete(sub_cat) ;
00460 cpl_vector_delete(wl_limits) ;
00461
00462
00463 if ((conv_kernel = irplib_wlxcorr_convolve_create_kernel(slitw,
00464 fwhm)) == NULL) {
00465 cpl_msg_error(cpl_func, "Cannot create convolution kernel") ;
00466 cpl_bivector_delete(gen_spectrum) ;
00467 return NULL ;
00468 }
00469
00470
00471 if (irplib_wlxcorr_convolve(cpl_bivector_get_y(gen_spectrum),conv_kernel)) {
00472 cpl_msg_error(cpl_func, "Cannot smoothe the signal");
00473 cpl_bivector_delete(gen_spectrum) ;
00474 return NULL ;
00475 }
00476 cpl_vector_delete(conv_kernel) ;
00477
00478 return gen_spectrum ;
00479 }
00480
00481
00491
00492 int irplib_wlxcorr_plot_solution(
00493 const cpl_polynomial * init,
00494 const cpl_polynomial * comp,
00495 const cpl_polynomial * sol,
00496 int pix_start,
00497 int pix_stop)
00498 {
00499 int nsamples, nplots ;
00500 cpl_vector ** vectors ;
00501 cpl_bivector * bivector ;
00502 double diff ;
00503 int i ;
00504
00505
00506 if (init == NULL || comp == NULL) return -1 ;
00507
00508
00509 nsamples = pix_stop - pix_start + 1 ;
00510 if (sol != NULL) nplots = 3 ;
00511 else nplots = 2 ;
00512
00513
00514 vectors = cpl_malloc((nplots+1)*sizeof(cpl_vector*)) ;
00515 for (i=0 ; i<nplots+1 ; i++) vectors[i] = cpl_vector_new(nsamples) ;
00516
00517
00518
00519 for (i=0 ; i<nsamples ; i++) {
00520 cpl_vector_set(vectors[0], i, pix_start+i) ;
00521 cpl_vector_set(vectors[1], i,
00522 cpl_polynomial_eval_1d(init, (double)(pix_start+i), NULL)) ;
00523 cpl_vector_set(vectors[2], i,
00524 cpl_polynomial_eval_1d(comp, (double)(pix_start+i), NULL)) ;
00525 if (sol != NULL)
00526 cpl_vector_set(vectors[3], i,
00527 cpl_polynomial_eval_1d(sol, (double)(pix_start+i), NULL)) ;
00528 }
00529
00530
00531 irplib_vectors_plot("set grid;set xlabel 'Position (pixels)';",
00532 "t '1-Initial / 2-Computed / 3-Solution' w lines",
00533 "", (const cpl_vector **)vectors, nplots+1);
00534
00535
00536 for (i=0 ; i<nplots+1 ; i++) cpl_vector_delete(vectors[i]) ;
00537 cpl_free(vectors) ;
00538
00539
00540 nplots -- ;
00541 vectors = cpl_malloc((nplots+1)*sizeof(cpl_vector*)) ;
00542 for (i=0 ; i<nplots+1 ; i++) vectors[i] = cpl_vector_new(nsamples) ;
00543
00544
00545
00546 for (i=0 ; i<nsamples ; i++) {
00547 cpl_vector_set(vectors[0], i, pix_start+i) ;
00548 diff = cpl_polynomial_eval_1d(comp, (double)(pix_start+i), NULL) -
00549 cpl_polynomial_eval_1d(init, (double)(pix_start+i), NULL) ;
00550 cpl_vector_set(vectors[1], i, diff) ;
00551 if (sol != NULL) {
00552 diff = cpl_polynomial_eval_1d(sol, (double)(pix_start+i), NULL) -
00553 cpl_polynomial_eval_1d(init, (double)(pix_start+i), NULL) ;
00554 cpl_vector_set(vectors[2], i, diff) ;
00555 }
00556 }
00557
00558
00559 if (sol == NULL) {
00560 bivector = cpl_bivector_wrap_vectors(vectors[0], vectors[1]) ;
00561 irplib_bivector_plot(
00562 "set grid;set xlabel 'Position (pixels)';set ylabel 'Wavelength difference';",
00563 "t 'Computed-Initial wavelenth' w lines", "", bivector);
00564 cpl_bivector_unwrap_vectors(bivector) ;
00565 } else {
00566 irplib_vectors_plot("set grid;set xlabel 'Position (pixels)';",
00567 "t '1-Computed - Initial / 2--Solution - Initial' w lines",
00568 "", (const cpl_vector **)vectors, nplots+1);
00569 }
00570
00571
00572 for (i=0 ; i<nplots+1 ; i++) cpl_vector_delete(vectors[i]) ;
00573 cpl_free(vectors) ;
00574
00575
00576 return 0 ;
00577 }
00578
00579
00586
00587 int irplib_wlxcorr_plot_spc_table(
00588 const cpl_table * spc_table,
00589 const char * title)
00590 {
00591 char title_loc[1024] ;
00592 cpl_vector ** vectors ;
00593 cpl_vector ** sub_vectors ;
00594 cpl_vector * tmp_vec ;
00595 int nsamples ;
00596 double hsize_nm, max, mean1, mean3 ;
00597 int start_ind, stop_ind, nblines, hsize_pix ;
00598 int i, j ;
00599
00600
00601 if (spc_table == NULL) return -1 ;
00602
00603
00604 nsamples = cpl_table_get_nrow(spc_table) ;
00605 hsize_nm = 0.2 ;
00606 hsize_pix = 10 ;
00607 nblines = 0 ;
00608 sprintf(title_loc,
00609 "t '%s - 1-Initial catalog/2-Corrected catalog/3-Observed' w lines",
00610 title) ;
00611 title_loc[1023] = (char)0 ;
00612
00613 vectors = cpl_malloc(4*sizeof(cpl_vector*)) ;
00614 vectors[0] = cpl_vector_wrap(nsamples,
00615 cpl_table_get_data_double(spc_table, IRPLIB_COL_XC_WAVELENGTH)) ;
00616 vectors[1] = cpl_vector_wrap(nsamples,
00617 cpl_table_get_data_double(spc_table, IRPLIB_COL_XC_CAT_INIT)) ;
00618 vectors[2] = cpl_vector_wrap(nsamples,
00619 cpl_table_get_data_double(spc_table, IRPLIB_COL_XC_CAT_FINAL)) ;
00620 vectors[3] = cpl_vector_wrap(nsamples,
00621 cpl_table_get_data_double(spc_table, IRPLIB_COL_XC_OBS)) ;
00622
00623
00624 mean1 = cpl_vector_get_mean(vectors[1]) ;
00625 mean3 = cpl_vector_get_mean(vectors[3]) ;
00626 if (fabs(mean3) > 1)
00627 cpl_vector_multiply_scalar(vectors[3], fabs(mean1/mean3)) ;
00628
00629 irplib_vectors_plot("set grid;set xlabel 'Wavelength (nm)';", title_loc,
00630 "", (const cpl_vector **)vectors, 4);
00631
00632
00633 if (fabs(mean3) > 1)
00634 cpl_vector_multiply_scalar(vectors[3], mean3/mean1) ;
00635
00636
00637 sprintf(title_loc,
00638 "t '%s - 1-Initial catalog/2-Corrected catalog/3-Observed (ZOOMED)' w lines",
00639 title) ;
00640 title_loc[1023] = (char)0 ;
00641 tmp_vec = cpl_vector_duplicate(vectors[2]) ;
00642 for (i=0 ; i<nblines ; i++) {
00643
00644 if ((max = cpl_vector_get_max(tmp_vec)) <= 0.0) break ;
00645 for (j=0 ; i<nsamples ; j++) {
00646 if (cpl_vector_get(tmp_vec, j) == max) break ;
00647 }
00648 if (j-hsize_pix < 0) start_ind = 0 ;
00649 else start_ind = j-hsize_pix ;
00650 if (j+hsize_pix > nsamples-1) stop_ind = nsamples-1 ;
00651 else stop_ind = j+hsize_pix ;
00652 for (j=start_ind ; j<=stop_ind ; j++) cpl_vector_set(tmp_vec, j, 0.0) ;
00653
00654 sub_vectors = cpl_malloc(4*sizeof(cpl_vector*)) ;
00655 sub_vectors[0]=cpl_vector_extract(vectors[0],start_ind,stop_ind,1);
00656 sub_vectors[1]=cpl_vector_extract(vectors[1],start_ind,stop_ind,1);
00657 sub_vectors[2]=cpl_vector_extract(vectors[2],start_ind,stop_ind,1);
00658 sub_vectors[3]=cpl_vector_extract(vectors[3],start_ind,stop_ind,1);
00659
00660 irplib_vectors_plot("set grid;set xlabel 'Wavelength (nm)';", title_loc,
00661 "", (const cpl_vector **)sub_vectors, 4);
00662
00663 cpl_vector_delete(sub_vectors[0]) ;
00664 cpl_vector_delete(sub_vectors[1]) ;
00665 cpl_vector_delete(sub_vectors[2]) ;
00666 cpl_vector_delete(sub_vectors[3]) ;
00667 cpl_free(sub_vectors) ;
00668 }
00669 cpl_vector_delete(tmp_vec) ;
00670
00671 cpl_vector_unwrap(vectors[0]) ;
00672 cpl_vector_unwrap(vectors[1]) ;
00673 cpl_vector_unwrap(vectors[2]) ;
00674 cpl_vector_unwrap(vectors[3]) ;
00675 cpl_free(vectors) ;
00676
00677 return 0 ;
00678 }
00679
00680
00688
00689 int irplib_wlxcorr_catalog_plot(
00690 const cpl_bivector * cat,
00691 double wmin,
00692 double wmax)
00693 {
00694 int start, stop ;
00695 cpl_bivector * subcat ;
00696 cpl_vector * subcat_x ;
00697 cpl_vector * subcat_y ;
00698 double * pwave ;
00699 int nvals, nvals_tot ;
00700 int i ;
00701
00702
00703 if (cat == NULL) return -1 ;
00704 if (wmax <= wmin) return -1 ;
00705
00706
00707 nvals_tot = cpl_bivector_get_size(cat) ;
00708
00709
00710 pwave = cpl_bivector_get_x_data(cat) ;
00711 if (pwave[0] >= wmin) start = 0 ;
00712 else start = -1 ;
00713 if (pwave[nvals_tot-1] <= wmax) stop = nvals_tot-1 ;
00714 else stop = -1 ;
00715 i=0 ;
00716 while ((pwave[i] < wmin) && (i<nvals_tot-1)) i++ ;
00717 start = i ;
00718 i= nvals_tot-1 ;
00719 while ((pwave[i] > wmax) && (i>0)) i-- ;
00720 stop = i ;
00721
00722 if (start>=stop) {
00723 cpl_msg_error(cpl_func, "Cannot plot the catalog") ;
00724 return -1 ;
00725 }
00726 nvals = start - stop + 1 ;
00727
00728
00729 subcat_x = cpl_vector_extract(cpl_bivector_get_x(cat), start, stop, 1) ;
00730 subcat_y = cpl_vector_extract(cpl_bivector_get_y(cat), start, stop, 1) ;
00731 subcat = cpl_bivector_wrap_vectors(subcat_x, subcat_y) ;
00732
00733
00734 irplib_bivector_plot(
00735 "set grid;set xlabel 'Wavelength (nm)';set ylabel 'Emission';",
00736 "t 'Catalog Spectrum' w lines", "", subcat);
00737 cpl_bivector_unwrap_vectors(subcat) ;
00738 cpl_vector_delete(subcat_x) ;
00739 cpl_vector_delete(subcat_y) ;
00740
00741 return 0 ;
00742 }
00743
00744
00758
00759 cpl_vector * irplib_wlxcorr_convolve_create_kernel(
00760 double slitw,
00761 double fwhm)
00762 {
00763 double sigma ;
00764 int ihtophat, gausshlen, convolen ;
00765 cpl_vector * self ;
00766 cpl_vector * tophat ;
00767 cpl_image * iself ;
00768 int i ;
00769
00770
00771 sigma = fwhm / (2.0 * sqrt(2.0*log(2.0)));
00772 ihtophat = (int)slitw/2 ;
00773 gausshlen = 1 + 5 * sigma + ihtophat ;
00774
00775
00776 convolen = 1 + 10 * sigma + 8 * ihtophat ;
00777
00778
00779 if ((slitw <= 0.0) || (fwhm <= 0.0) ||(convolen < 2 * gausshlen))
00780 return NULL ;
00781
00782
00783 self = cpl_vector_new(gausshlen) ;
00784 tophat = cpl_vector_new(convolen) ;
00785
00786
00787 iself = cpl_image_wrap_double(gausshlen, 1, cpl_vector_get_data(self));
00788
00789
00790 cpl_image_fill_gaussian(iself, 1.0, 1.0, sqrt(2.0*atan(1.0)*4.0),sigma,1.0);
00791 cpl_image_unwrap(iself) ;
00792
00793
00794 cpl_vector_fill(tophat, 0.0);
00795
00796 for (i = convolen/2-ihtophat; i < 1+convolen/2+ihtophat; i++)
00797 cpl_vector_set(tophat, i, 1.0/(1.0+2.0*ihtophat));
00798
00799
00800 if (irplib_wlxcorr_convolve(tophat, self)) {
00801 cpl_msg_error(cpl_func, "Cannot convolve") ;
00802 cpl_vector_delete(self) ;
00803 cpl_vector_delete(tophat) ;
00804 return NULL ;
00805 }
00806
00807
00808
00809 for (i = 0 ; i < gausshlen; i++)
00810 cpl_vector_set(self, i, cpl_vector_get(tophat, i + convolen/2));
00811 cpl_vector_delete(tophat) ;
00812
00813 return self;
00814 }
00815
00816
00825
00826 int irplib_wlxcorr_convolve(
00827 cpl_vector * smoothed,
00828 const cpl_vector * conv_kernel)
00829 {
00830 int nsamples ;
00831 int ihwidth ;
00832 cpl_vector * raw ;
00833 double * psmoothe ;
00834 double * praw ;
00835 double * psymm ;
00836 int i, j ;
00837
00838
00839 if ((!smoothed) || (!conv_kernel)) return -1 ;
00840
00841
00842 nsamples = cpl_vector_get_size(smoothed) ;
00843 ihwidth = cpl_vector_get_size(conv_kernel) - 1 ;
00844 if (ihwidth >= nsamples) return -1 ;
00845 psymm = cpl_vector_get_data(conv_kernel) ;
00846 psmoothe = cpl_vector_get_data(smoothed) ;
00847
00848
00849 raw = cpl_vector_duplicate(smoothed) ;
00850 praw = cpl_vector_get_data(raw) ;
00851
00852
00853 for (i=0 ; i<ihwidth ; i++) {
00854 psmoothe[i] = praw[i] * psymm[0];
00855 for (j=1 ; j <= ihwidth ; j++) {
00856 const int k = i-j < 0 ? 0 : i-j;
00857 psmoothe[i] += (praw[k]+praw[i+j]) * psymm[j];
00858 }
00859 }
00860
00861 for (i=ihwidth ; i<nsamples-ihwidth ; i++) {
00862 psmoothe[i] = praw[i] * psymm[0];
00863 for (j=1 ; j<=ihwidth ; j++)
00864 psmoothe[i] += (praw[i-j]+praw[i+j]) * psymm[j];
00865 }
00866 for (i=nsamples-ihwidth ; i<nsamples ; i++) {
00867 psmoothe[i] = praw[i] * psymm[0];
00868 for (j=1 ; j<=ihwidth ; j++) {
00869 const int k = i+j > nsamples-1 ? nsamples - 1 : i+j;
00870 psmoothe[i] += (praw[k]+praw[i-j]) * psymm[j];
00871 }
00872 }
00873 cpl_vector_delete(raw) ;
00874 return 0 ;
00875 }
00876
00879
00889
00890 static int irplib_wlxcorr_signal_resample(
00891 cpl_vector * resampled,
00892 const cpl_vector * xbounds,
00893 const cpl_bivector * hires)
00894 {
00895 cpl_vector * xhires ;
00896 cpl_vector * yhires ;
00897 double * pxhires ;
00898 double * pyhires ;
00899 double * pxbounds ;
00900 cpl_vector * ybounds ;
00901 cpl_bivector * boundary ;
00902 double * pybounds ;
00903 double * presampled ;
00904 int nsamples ;
00905 int i, itt ;
00906
00907
00908 if ((!resampled) || (!xbounds) || (!hires)) return -1 ;
00909
00910
00911 nsamples = cpl_vector_get_size(resampled) ;
00912
00913
00914 presampled = cpl_vector_get_data(resampled) ;
00915 pxbounds = cpl_vector_get_data(xbounds) ;
00916 xhires = cpl_bivector_get_x(hires) ;
00917 yhires = cpl_bivector_get_y(hires) ;
00918 pxhires = cpl_vector_get_data(xhires) ;
00919 pyhires = cpl_vector_get_data(yhires) ;
00920
00921
00922 ybounds = cpl_vector_new(cpl_vector_get_size(xbounds)) ;
00923 boundary = cpl_bivector_wrap_vectors((cpl_vector*)xbounds,ybounds) ;
00924 pybounds = cpl_vector_get_data(ybounds) ;
00925
00926
00927 if (cpl_bivector_get_size(boundary) != nsamples + 1) {
00928 cpl_bivector_unwrap_vectors(boundary) ;
00929 cpl_vector_delete(ybounds) ;
00930 return -1 ;
00931 }
00932
00933
00934 itt = cpl_vector_find(xhires, pxbounds[0]);
00935
00936
00937 if (cpl_bivector_interpolate_linear(boundary, hires)) {
00938 cpl_msg_error(cpl_func, "Cannot interpolate the signal") ;
00939 cpl_bivector_unwrap_vectors(boundary) ;
00940 cpl_vector_delete(ybounds) ;
00941 return -1 ;
00942 }
00943
00944
00945
00946 while (pxhires[itt] < pxbounds[0]) itt++;
00947
00948 for (i=0; i < nsamples; i++) {
00949
00950
00951
00952
00953 double xlow = pxbounds[i];
00954 double x = pxhires[itt];
00955
00956 if (x > pxbounds[i+1]) x = pxbounds[i+1];
00957
00958
00959 presampled[i] = pybounds[i] * (x - xlow);
00960
00961
00962 while ((pxhires[itt] < pxbounds[i+1]) &&
00963 (itt < cpl_bivector_get_size(hires))) {
00964 const double xprev = x;
00965 x = pxhires[itt+1];
00966 if (x > pxbounds[i+1]) x = pxbounds[i+1];
00967 presampled[i] += pyhires[itt] * (x - xlow);
00968 xlow = xprev;
00969 itt++;
00970 }
00971
00972
00973
00974 presampled[i] += pybounds[i+1] * (pxbounds[i+1] - xlow);
00975
00976
00977
00978 presampled[i] /= 2 * (pxbounds[i+1] - pxbounds[i]);
00979 }
00980 cpl_bivector_unwrap_vectors(boundary) ;
00981 cpl_vector_delete(ybounds) ;
00982 return 0 ;
00983 }
00984