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 #include <sinfo_skycor.h>
00036 #include <math.h>
00037 #include <sinfo_new_cube_ops.h>
00038 #include <cpl.h>
00039 #include <irplib_utils.h>
00040 #include <irplib_error.h>
00041 #include "sinfo_pfits.h"
00042 #include "sinfo_functions.h"
00043
00044 #include "sinfo_msg.h"
00045 #include "sinfo_error.h"
00046 #include "sinfo_globals.h"
00047 #include "sinfo_utils_wrappers.h"
00048 #include "sinfo_utl_cube2spectrum.h"
00049 #include "sinfo_pro_types.h"
00050
00051
00052
00053
00054 #define BAND_H_WAVE_MIN 1.445 //not used
00055 #define BAND_H_WAVE_MAX 1.820 //not used
00056
00057 #define BAND_K_WAVE_MIN 1.945 //not used
00058 #define BAND_K_WAVE_MAX 2.460 //not used
00059
00060 #define BAND_J_WAVE_MIN 1.445 //not used
00061 #define BAND_J_WAVE_MAX 1.82 //not used
00062
00063 #define SINFO_FIT_BKG_TEMP 280.
00064 #define SKY_THRES 0.95
00065 #define SKY_LINE_MAX_CUT 4
00066 #define SKY_LINE_MIN_CUT 4
00067 #define HISTO_NBINS 100
00068 #define HISTO_MIN_SIZE 10
00069 #define HISTO_Y_CUT 10
00070
00071 #define HISTO_X_LEFT_CUT 1.0
00072 #define HISTO_X_RIGHT_CUT 0.5
00073 #define HISTO_DIST_TEMPC_MIN_FCT 5.
00074 #define HISTO_DIST_TEMPC_MAX_FCT 0.25
00075
00076 #define XCOR_YSHIFT_KS_CLIP 3
00077
00078
00079 #define HPLANK 6.62606876e-34; // J s
00080 #define CLIGHT 2.99792458e+08; // m / s
00081 #define KBOLTZ 1.3806503e-23; // J / K
00082 #define AMOEBA_FTOL 1.e-5
00083 #define NBOUND 14
00084 #define NROT 25
00085 #define N_ITER_FIT_LM 15
00086 #define N_ITER_FIT_AMOEBA 10
00087
00088 double sinfo_scale_fct=1;
00089 static cpl_vector* sa_vx=NULL;
00090 static cpl_vector* sa_vy=NULL;
00091
00092 static cpl_vector* sa_ox=NULL;
00093 static cpl_vector* sa_oy=NULL;
00094 static cpl_vector* sa_sy=NULL;
00095
00096
00097
00098
00099 static cpl_imagelist*
00100 sinfo_cube_zshift_simple(cpl_imagelist* inp,
00101 const float shift);
00102
00103 static int
00104 sinfo_scales_obj_sky_cubes(cpl_frame* obj_frm,
00105 cpl_frame* sky_frm,
00106 cpl_table* bkg,
00107 cpl_table* rscale,
00108 cpl_imagelist** obj_cor);
00109
00110 static void
00111 sinfo_optimise_sky_sub(const double wtol,
00112 const double line_hw,
00113 const int method,
00114 const int do_rot,
00115 cpl_table* lrange,
00116 cpl_table* lambda,
00117 cpl_table* lr41,
00118 cpl_table* lr52,
00119 cpl_table* lr63,
00120 cpl_table* lr74,
00121 cpl_table* lr02,
00122 cpl_table* lr85,
00123 cpl_table* lr20,
00124 cpl_table* lr31,
00125 cpl_table* lr42,
00126 cpl_table* lr53,
00127 cpl_table* lr64,
00128 cpl_table* lr75,
00129 cpl_table* lr86,
00130 cpl_table* lr97,
00131 cpl_table* lr00,
00132 cpl_table** int_obj,
00133 cpl_table** int_sky,
00134 cpl_table** rscale);
00135
00136 static void
00137 sinfo_shift_sky(cpl_frame** sky_frm,
00138 cpl_table** int_sky,
00139 const double zshift);
00140
00141 static double
00142 sinfo_xcorr(cpl_table* int_obj,
00143 cpl_table* int_sky,
00144 cpl_table* lambda,
00145 const double dispersion,
00146 const double line_hw);
00147
00148 static cpl_table*
00149 sinfo_interpolate_sky(const cpl_table* inp,const cpl_table* lambdas);
00150
00151 static int
00152 sinfo_check_screw_values(cpl_table** int_obj,
00153 cpl_table** int_sky,
00154 cpl_table* grange,
00155 const double wtol);
00156
00157 static int
00158 sinfo_choose_good_sky_pixels(cpl_frame* obj_frm,
00159 cpl_image* r_img,
00160 cpl_image* g_img,
00161 const double min_frac,
00162 cpl_image** mask);
00163
00164
00165 static int
00166 sinfo_sum_spectra(const cpl_frame* obj,
00167 const cpl_frame* sky,
00168 cpl_image* mask,
00169 cpl_table* wrange,
00170 cpl_table** int_obj,
00171 cpl_table** int_sky);
00172
00173 int
00174 sinfo_thermal_background2(cpl_table* int_sky,
00175 cpl_table* lambda,
00176 cpl_table* lrange,
00177 cpl_table** bkg);
00178
00179 static int
00180 sinfo_thermal_background(cpl_table* int_sky,
00181 cpl_table* lambda,
00182 cpl_table* lrange,
00183 const double temp,
00184 const int niter,
00185 cpl_table** bkg,
00186 int* success_fit);
00187
00188 static int
00189 sinfo_object_flag_sky_pixels(cpl_frame* obj_frm,
00190 cpl_table* lambda,
00191 cpl_table* mrange,
00192 cpl_imagelist* flag_data,
00193 const double tol,
00194 cpl_image** g_img,
00195 cpl_image** r_img,
00196 cpl_image** image);
00197
00198 static int
00199 sinfo_object_flag_low_values(cpl_frame* obj_frm,
00200 const double cnt,
00201 const double sig,
00202 cpl_imagelist** flag_data);
00203
00204 static int
00205 sinfo_object_estimate_noise(cpl_frame* obj_frm, const int obj_noise_fit,
00206 double* centre, double* noise);
00207
00208
00209
00210 static int
00211 sinfo_set_ranges(cpl_frame* obj_frm,
00212 cpl_frame* sky_frm,
00213 cpl_parameterlist* cfg,
00214 cpl_table** lambda,
00215 cpl_table** lr41,
00216 cpl_table** lr52,
00217 cpl_table** lr63,
00218 cpl_table** lr74,
00219 cpl_table** lr02,
00220 cpl_table** lr85,
00221 cpl_table** lr20,
00222 cpl_table** lr31,
00223 cpl_table** lr42,
00224 cpl_table** lr53,
00225 cpl_table** lr64,
00226 cpl_table** lr75,
00227 cpl_table** lr86,
00228 cpl_table** lr97,
00229 cpl_table** lr00,
00230 cpl_table** lrange,
00231 cpl_table** grange,
00232 cpl_table** lambdas,
00233 cpl_table** mrange,
00234 int* sky_interp_sw,
00235 double* dispersion);
00236
00237
00238 static cpl_table*
00239 sinfo_table_extract_rest(cpl_table* inp,
00240 cpl_table* low,
00241 cpl_table* med,
00242 const double wtol);
00243
00244 static int
00245 sinfo_get_sub_regions(cpl_table* sky,
00246 cpl_table* x1,
00247 cpl_table* pos,
00248 const double wtol,
00249 const int npixw,
00250 cpl_table** sub_regions);
00251
00252
00253 static int
00254 sinfo_get_obj_sky_wav_sub(cpl_table* obj,
00255 cpl_table* sky,
00256 cpl_table* wav,
00257 cpl_table* sel,
00258 const double wtol,
00259 cpl_table** sub_obj,
00260 cpl_table** sub_sky,
00261 cpl_table** sub_wav);
00262
00263
00264 static cpl_table*
00265 sinfo_find_rot_waves(const double w_rot[],
00266 const int npix_w,
00267 const double w_step,
00268 cpl_table* range);
00269 static int
00270 sinfo_compute_line_ratio(cpl_table* obj,
00271 cpl_table* sky,
00272 const double wtol,
00273 const int meth,
00274 const cpl_table* sel_regions,
00275 cpl_table* cont_regions,
00276 double* r);
00277
00278 static cpl_table*
00279 sinfo_table_subtract_continuum(cpl_table* lin,cpl_table* cnt);
00280
00281 static double
00282 sinfo_fit_bkg(double p[]);
00283
00284 static double
00285 sinfo_fit_sky(double p[]);
00286
00287 static int
00288 sinfo_get_line_ratio_amoeba(cpl_table* obj,
00289 cpl_table* sky,
00290 double* r);
00291
00292 static cpl_table*
00293 sinfo_table_interpol(cpl_table* obj_lin,
00294 cpl_table* obj_cnt,
00295 cpl_table* sky_lin,
00296 cpl_table* sky_cnt,
00297 const double r);
00298
00299
00300 static int
00301 sinfo_get_line_ratio(cpl_table* obj_lin,
00302 cpl_table* obj_cnt,
00303 cpl_table* sky_lin,
00304 cpl_table* sky_cnt,
00305 const int method,
00306 double* r);
00307
00308 static cpl_table*
00309 sinfo_table_shift_simple(cpl_table* inp,
00310 const char* col,
00311 const double shift);
00312
00313 static int
00314 sinfo_table_set_column_invalid(cpl_table** int_sky,const char* col);
00315
00316 static int
00317 sinfo_table_set(cpl_table** out,
00318 const cpl_table* ref,
00319 const double val,
00320 const double tol);
00321
00322 static int
00323 sinfo_table_threshold(cpl_table** t,
00324 const char* column,
00325 const double low_cut,
00326 const double hig_cut,
00327 const double low_ass,
00328 const double hig_ass);
00329
00330
00331
00332
00333 static double sinfo_fac(const double x, const double t);
00334
00335 static int sinfo_fitbkg(const double x[],
00336 const double a[],
00337 double *result);
00338 static int sinfo_fitbkg_derivative(const double x[],
00339 const double a[],
00340 double result[]);
00341
00342
00343 static int
00344 sinfo_convolve_kernel(cpl_table** t, const int rad);
00345 int
00346 sinfo_convolve_kernel2(cpl_table** t, const int rad);
00347
00348 int
00349 sinfo_convolve_gauss(cpl_table** t, const int rad, const double fwhm);
00350 int
00351 sinfo_convolve_exp(cpl_table** t, const int rad, const double fwhm);
00352
00353 static int
00354 sinfo_table_sky_obj_flag_nan(cpl_table** s,cpl_table** o, cpl_table** w);
00355
00356
00357
00358
00359 static double
00360 sinfo_gaussian(double area,double sigma,double x,double x0,double off);
00361
00362 int
00363 sinfo_table_smooth_column(cpl_table** t, const char* c, const int r);
00364
00365 static int
00366 sinfo_image_flag_nan(cpl_image** im);
00367
00368 static int
00369 sinfo_table_flag_nan(cpl_table** t,const char* label);
00370
00371
00372 static int
00373 sinfo_cnt_mask_thresh_and_obj_finite(const cpl_image* mask,
00374 const double t,
00375 const cpl_image* obj);
00376
00377
00378
00379
00380 static cpl_table*
00381 sinfo_interpolate(const cpl_table* inp,
00382 const cpl_table* lambdas,
00383 const char* name,
00384 const char* method);
00385
00386 static cpl_table*
00387 sinfo_image2table(const cpl_image* im);
00388
00389 static int
00390 sinfo_table_extract_finite(const cpl_table* in1,
00391 const cpl_table* in2,
00392 cpl_table** ou1,
00393 cpl_table** ou2);
00394
00395 static cpl_table*
00396 sinfo_slice_z(const cpl_imagelist* cin,const int i,const int j);
00397
00398
00399
00400 static cpl_imagelist*
00401 sinfo_imagelist_select_range(const cpl_imagelist* inp,
00402 const cpl_table* full,
00403 const cpl_table* good,
00404 const double tol);
00405
00406
00407
00408 static cpl_table*
00409 sinfo_table_select_range(cpl_table* inp,
00410 cpl_table* ref,
00411 const double tol);
00412
00413 static int
00414 sinfo_table_fill_column_over_range(cpl_table** inp,
00415 const cpl_table* ref,
00416 const char* col,
00417 const double val,
00418 const double tol);
00419
00420 static double
00421 sinfo_table_column_interpolate(const cpl_table* t,
00422 const char* name,
00423 const double x);
00424
00425 static int
00426 sinfo_table_get_index_of_val(cpl_table* t,
00427 const char* name,
00428 double val,
00429 cpl_type type);
00430
00431 static int
00432 sinfo_table_get_index_of_max(cpl_table* t,
00433 const char* name,
00434 cpl_type type);
00435
00436
00437
00438
00439 static cpl_table*
00440 sinfo_where_tab_min_max(cpl_table* t,
00441 const char* col,
00442 cpl_table_select_operator op1,
00443 const double v1,
00444 cpl_table_select_operator op2,
00445 const double v2);
00446
00447
00448 static int
00449 sinfo_table_column_dindgen(cpl_table** t, const char* label);
00450
00451
00452 static int
00453 sinfo_histogram(const cpl_table* data,
00454 const int nbins,
00455 const double min,
00456 const double max,
00457 cpl_table** histo);
00458
00459
00460
00461
00462
00463
00464
00465
00466
00467
00468
00469
00474
00483 sinfo_skycor_qc*
00484 sinfo_skycor_qc_new(void)
00485 {
00486 sinfo_skycor_qc * sqc;
00487 sqc= cpl_malloc(sizeof(sinfo_skycor_qc));
00488
00489 sqc->th_fit=0;
00490
00491 return sqc;
00492
00493 }
00500 void
00501 sinfo_skycor_qc_delete(sinfo_skycor_qc** sqc)
00502 {
00503 if((*sqc) != NULL) {
00504 cpl_free(*sqc);
00505 *sqc=NULL;
00506 }
00507 }
00508
00509
00510
00521
00522 int
00523 sinfo_skycor(cpl_parameterlist * config,
00524 cpl_frame* obj_frm,
00525 cpl_frame* sky_frm,
00526 sinfo_skycor_qc* sqc,
00527 cpl_imagelist** obj_cor,
00528 cpl_table** int_obj)
00529 {
00530
00531 cpl_table* bkg=NULL;
00532
00533 cpl_table* lambda=NULL;
00534 cpl_table* lr41=NULL;
00535 cpl_table* lr52=NULL;
00536 cpl_table* lr63=NULL;
00537 cpl_table* lr74=NULL;
00538 cpl_table* lr02=NULL;
00539 cpl_table* lr85=NULL;
00540 cpl_table* lr20=NULL;
00541 cpl_table* lr31=NULL;
00542 cpl_table* lr42=NULL;
00543 cpl_table* lr53=NULL;
00544 cpl_table* lr64=NULL;
00545 cpl_table* lr75=NULL;
00546 cpl_table* lr86=NULL;
00547 cpl_table* lr97=NULL;
00548 cpl_table* lr00=NULL;
00549 cpl_table* lrange=NULL;
00550 cpl_table* mrange=NULL;
00551 cpl_table* grange=NULL;
00552 cpl_table* lambdas=NULL;
00553
00554 cpl_table* int_sky=NULL;
00555
00556 cpl_image* mask=NULL;
00557 cpl_image* gpix=NULL;
00558 cpl_image* ratio=NULL;
00559 cpl_image* ima_sky=NULL;
00560 cpl_imagelist* fdata=NULL;
00561 cpl_table* rscale=NULL;
00562 cpl_parameter* p=NULL;
00563
00564 double min_frac=SINFO_MIN_FRAC;
00565 int th_fit=0;
00566 double dispersion=0;
00567 double noise=0;
00568 double fit_temp=SINFO_FIT_BKG_TEMP;
00569
00570 double centre=0;
00571 int sky_interp_sw=0;
00572 double wshift=0;
00573 double line_hw=SINFO_LINE_HALF_WIDTH;
00574 int method=0;
00575 int do_rot=0;
00576 int obj_noise_fit=0;
00577 int niter=3;
00578
00579
00580
00581 check_nomsg(p=cpl_parameterlist_find(config,
00582 "sinfoni.sinfo_utl_skycor.min_frac"));
00583 check_nomsg(min_frac=cpl_parameter_get_double(p));
00584
00585 check_nomsg(p=cpl_parameterlist_find(config,
00586 "sinfoni.sinfo_utl_skycor.line_half_width"));
00587 check_nomsg(line_hw=cpl_parameter_get_double(p));
00588
00589 check_nomsg(p=cpl_parameterlist_find(config,
00590 "sinfoni.sinfo_utl_skycor.scale_method"));
00591 check_nomsg(method=cpl_parameter_get_int(p));
00592
00593
00594
00595 check_nomsg(p=cpl_parameterlist_find(config,
00596 "sinfoni.sinfo_utl_skycor.rot_cor"));
00597 check_nomsg(do_rot=cpl_parameter_get_bool(p));
00598
00599
00600 check_nomsg(p=cpl_parameterlist_find(config,
00601 "sinfoni.sinfo_utl_skycor.fit_obj_noise"));
00602 check_nomsg(obj_noise_fit=cpl_parameter_get_bool(p));
00603
00604
00605 check_nomsg(p=cpl_parameterlist_find(config,
00606 "sinfoni.sinfo_utl_skycor.niter"));
00607 check_nomsg(niter=cpl_parameter_get_int(p));
00608
00609
00610 sinfo_msg("Set wavelength ranges");
00611 ck0(sinfo_set_ranges(obj_frm,sky_frm,config,&lambda,
00612 &lr41,&lr52,&lr63,&lr74,&lr02,&lr85,&lr20,&lr31,&lr42,
00613 &lr53,&lr64,&lr75,&lr86,&lr97,&lr00,
00614 &lrange,&grange,&lambdas,&mrange,&sky_interp_sw,
00615 &dispersion),"Setting wavelength ranges");
00616
00617
00618
00619
00620
00621
00622
00623
00624
00625
00626
00627
00628
00629
00630
00631
00632
00633
00634
00635
00636
00637
00638
00639
00640 sinfo_msg("Estimate noise");
00641 ck0(sinfo_object_estimate_noise(obj_frm,obj_noise_fit,¢re,&noise),
00642 "Estimating noise");
00643
00644 sinfo_msg("Flag object low_levels");
00645 ck0(sinfo_object_flag_low_values(obj_frm,centre,noise,&fdata),
00646 "Flagging low pix");
00647
00648
00649
00650
00651 sinfo_msg("Flag sky pixels");
00652 ck0(sinfo_object_flag_sky_pixels(obj_frm,lambda,mrange,fdata,dispersion,
00653 &gpix,&ratio,&ima_sky),
00654 "Flagging sky pixels");
00655
00656
00657
00658
00659
00660
00661
00662
00663
00664
00665 sinfo_msg("Choose good sky (with > 95%% good spectral) pixels");
00666 ck0(sinfo_choose_good_sky_pixels(obj_frm,ratio,gpix,min_frac,&mask),
00667 "Choosing good sky pixels");
00668
00669
00670
00671
00672
00673
00674
00675 sinfo_msg("Sum obj and sky spectra");
00676 ck0(sinfo_sum_spectra(obj_frm,sky_frm,mask,lambda,int_obj,&int_sky),
00677 "summing obj & sky spectra");
00678
00679
00680
00681
00682
00683 sinfo_msg("Computes thermal background (Implement amoeba fit)");
00684 ck0(sinfo_thermal_background(int_sky,lambda,lrange,fit_temp,niter,
00685 &bkg,&th_fit),
00686 "getting termal bkg");
00687
00688 sqc->th_fit=th_fit;
00689
00690
00691
00692
00693
00694
00695
00696
00697 sinfo_msg("Subtracts thermal background from integrated OH spectrum");
00698
00699
00700 check_nomsg(cpl_table_duplicate_column(int_sky,"BKG",bkg,"INT2"));
00701 check_nomsg(cpl_table_subtract_columns(int_sky,"INT","BKG"));
00702 check_nomsg(cpl_table_erase_column(int_sky,"BKG"));
00703
00704
00705
00706
00707
00708 sinfo_msg("Checks for screw values at ends of spectrum");
00709
00710
00711 sinfo_check_screw_values(int_obj,&int_sky,grange,dispersion);
00712
00713
00714 if(sky_interp_sw == 1) {
00715 sinfo_msg("Interpolate sky if necessary");
00716 sinfo_interpolate_sky(int_sky,lambdas);
00717 }
00718
00719
00720 sinfo_msg("Crosscorrelate obj & sky to check for lambda offset");
00721
00722
00723 check_nomsg(wshift=sinfo_xcorr(*int_obj,int_sky,lambda,dispersion,line_hw));
00724
00725
00726
00727 sinfo_msg("Dispersion=%f",dispersion);
00728 sinfo_msg("Shift sky of %f pixels toward object",wshift/dispersion);
00729
00730 check_nomsg(sinfo_shift_sky(&sky_frm,&int_sky,wshift/dispersion));
00731
00732
00733
00734
00735
00736
00737 sinfo_msg("Optimise sky subtraction");
00738 check_nomsg(sinfo_optimise_sky_sub(dispersion,line_hw,method,do_rot,
00739 lrange,lambda,
00740 lr41,lr52,lr63,lr74,lr02,lr85,
00741 lr20,lr31,lr42,lr53,lr64,lr75,
00742 lr86,lr97,lr00,int_obj,&int_sky,
00743 &rscale));
00744
00745
00746
00747
00748
00749
00750 sinfo_msg("Apply same scaling to whole cubes");
00751 ck0_nomsg(sinfo_scales_obj_sky_cubes(obj_frm,sky_frm,bkg,rscale,obj_cor));
00752
00753
00754 cleanup:
00755 sinfo_free_table(&rscale);
00756 sinfo_free_imagelist(&fdata);
00757
00758 sinfo_free_table(&bkg);
00759 sinfo_free_table(&lambda);
00760 sinfo_free_table(&lrange);
00761 sinfo_free_table(&mrange);
00762 sinfo_free_table(&grange);
00763 sinfo_free_table(&lambdas);
00764 sinfo_free_image(&mask);
00765
00766 sinfo_free_table(&lr41);
00767 sinfo_free_table(&lr52);
00768 sinfo_free_table(&lr63);
00769 sinfo_free_table(&lr74);
00770 sinfo_free_table(&lr02);
00771 sinfo_free_table(&lr85);
00772 sinfo_free_table(&lr20);
00773 sinfo_free_table(&lr31);
00774 sinfo_free_table(&lr42);
00775 sinfo_free_table(&lr53);
00776 sinfo_free_table(&lr64);
00777 sinfo_free_table(&lr75);
00778 sinfo_free_table(&lr86);
00779 sinfo_free_table(&lr97);
00780 sinfo_free_table(&lr00);
00781
00782 sinfo_free_image(&gpix);
00783 sinfo_free_image(&ratio);
00784 sinfo_free_image(&ima_sky);
00785
00786 sinfo_free_table(&int_sky);
00787
00788 if (cpl_error_get_code() != CPL_ERROR_NONE) {
00789
00790 irplib_error_dump(CPL_MSG_ERROR,CPL_MSG_ERROR);
00791 return -1;
00792 } else {
00793 return 0;
00794 }
00795
00796 }
00797
00798
00799
00800
00801
00802
00833
00834
00835
00836 int
00837 sinfo_set_ranges(cpl_frame* obj_frm,
00838 cpl_frame* sky_frm,
00839 cpl_parameterlist* cfg,
00840 cpl_table** lambda,
00841 cpl_table** lr41,
00842 cpl_table** lr52,
00843 cpl_table** lr63,
00844 cpl_table** lr74,
00845 cpl_table** lr02,
00846 cpl_table** lr85,
00847 cpl_table** lr20,
00848 cpl_table** lr31,
00849 cpl_table** lr42,
00850 cpl_table** lr53,
00851 cpl_table** lr64,
00852 cpl_table** lr75,
00853 cpl_table** lr86,
00854 cpl_table** lr97,
00855 cpl_table** lr00,
00856 cpl_table** lrange,
00857 cpl_table** grange,
00858 cpl_table** lambdas,
00859 cpl_table** mrange,
00860 int* sky_interp_sw,
00861 double* dispersion)
00862
00863 {
00864
00865 cpl_propertylist* plist=NULL;
00866 double crval_obj=0;
00867 double cdelt_obj=0;
00868 double crpix_obj=0;
00869 int xsize_obj=0;
00870 int ysize_obj=0;
00871 int zsize_obj=0;
00872
00873
00874 double crval_sky=0;
00875 double cdelt_sky=0;
00876 double crpix_sky=0;
00877 int xsize_sky=0;
00878 int ysize_sky=0;
00879 int zsize_sky=0;
00880
00881 int nrow=0;
00882
00883 const double w_j_min=1.100;
00884 const double w_j_max=1.400;
00885 const double w_h_min=1.445;
00886 const double w_h_max=1.820;
00887 const double w_k_min=1.945;
00888 const double w_k_max=2.460;
00889
00890 double ws=0;
00891 double we=0;
00892 double mean=0;
00893
00894 cpl_parameter* p=NULL;
00895
00896
00897
00898 double w_bound[NBOUND]={1.067,1.125,1.196,1.252,1.289,
00899 1.400,1.472,1.5543,1.6356,1.7253,
00900 1.840,1.9570,2.095,2.300};
00901
00902 cpl_table* tmp_tbl=NULL;
00903 cpl_table* add1=NULL;
00904
00905
00906
00907
00908 cknull_nomsg(plist=cpl_propertylist_load(cpl_frame_get_filename(obj_frm),0));
00909 check_nomsg(crval_obj=sinfo_pfits_get_crval3(plist));
00910 check_nomsg(cdelt_obj=sinfo_pfits_get_cdelt3(plist));
00911 check_nomsg(crpix_obj=sinfo_pfits_get_crpix3(plist));
00912 check_nomsg(xsize_obj=sinfo_pfits_get_naxis1(plist));
00913 check_nomsg(ysize_obj=sinfo_pfits_get_naxis2(plist));
00914 check_nomsg(zsize_obj=sinfo_pfits_get_naxis3(plist));
00915
00916 sinfo_free_propertylist(&plist);
00917 *dispersion=cdelt_obj;
00918
00919
00920 check_nomsg(*lambda=cpl_table_new(zsize_obj));
00921 cpl_table_new_column(*lambda,"WAVE",CPL_TYPE_DOUBLE);
00922 cpl_table_new_column(*lambda,"INDEX",CPL_TYPE_DOUBLE);
00923 check_nomsg(sinfo_table_column_dindgen(lambda,"INDEX"));
00924 check_nomsg(sinfo_table_column_dindgen(lambda,"WAVE"));
00925
00926 check_nomsg(cpl_table_add_scalar(*lambda,"WAVE",1.));
00927 check_nomsg(cpl_table_subtract_scalar(*lambda,"WAVE",crpix_obj));
00928 check_nomsg(cpl_table_multiply_scalar(*lambda,"WAVE",cdelt_obj));
00929 check_nomsg(cpl_table_add_scalar(*lambda,"WAVE",crval_obj));
00930
00931
00932
00933
00934 cknull_nomsg(*lr41=sinfo_where_tab_min_max(*lambda,
00935 "WAVE",
00936 CPL_NOT_LESS_THAN,
00937 w_j_min,
00938 CPL_LESS_THAN,
00939 w_bound[0]));
00940
00941 cknull_nomsg(*lr52=sinfo_where_tab_min_max(*lambda,
00942 "WAVE",
00943 CPL_NOT_LESS_THAN,
00944 w_bound[0],
00945 CPL_LESS_THAN,
00946 w_bound[1]));
00947
00948 cknull_nomsg(*lr63=sinfo_where_tab_min_max(*lambda,
00949 "WAVE",
00950 CPL_NOT_LESS_THAN,
00951 w_bound[1],
00952 CPL_LESS_THAN,
00953 w_bound[2]));
00954
00955
00956 cknull_nomsg(*lr74=sinfo_where_tab_min_max(*lambda,
00957 "WAVE",
00958 CPL_NOT_LESS_THAN,
00959 w_bound[2],
00960 CPL_LESS_THAN,
00961 w_bound[3]));
00962
00963 cknull_nomsg(*lr20=sinfo_where_tab_min_max(*lambda,
00964 "WAVE",
00965 CPL_NOT_LESS_THAN,
00966 w_bound[3],
00967 CPL_LESS_THAN,
00968 w_bound[4]));
00969
00970 cknull_nomsg(*lr02=sinfo_where_tab_min_max(*lambda,
00971 "WAVE",
00972 CPL_NOT_LESS_THAN,
00973 w_bound[4],
00974 CPL_LESS_THAN,
00975 w_bound[5]));
00976
00977
00978 cknull_nomsg(*lr85=sinfo_where_tab_min_max(*lambda,
00979 "WAVE",
00980 CPL_NOT_LESS_THAN,
00981 w_bound[5],
00982 CPL_LESS_THAN,
00983 w_bound[6]));
00984
00985 cknull_nomsg(*lr31=sinfo_where_tab_min_max(*lambda,
00986 "WAVE",
00987 CPL_NOT_LESS_THAN,
00988 w_bound[6],
00989 CPL_LESS_THAN,
00990 w_bound[7]));
00991
00992 cknull_nomsg(*lr42=sinfo_where_tab_min_max(*lambda,
00993 "WAVE",
00994 CPL_NOT_LESS_THAN,
00995 w_bound[7],
00996 CPL_LESS_THAN,
00997 w_bound[8]));
00998
00999 cknull_nomsg(*lr53=sinfo_where_tab_min_max(*lambda,
01000 "WAVE",
01001 CPL_NOT_LESS_THAN,
01002 w_bound[8],
01003 CPL_LESS_THAN,
01004 w_bound[9]));
01005
01006 cknull_nomsg(*lr64=sinfo_where_tab_min_max(*lambda,
01007 "WAVE",
01008 CPL_NOT_LESS_THAN,
01009 w_bound[0],
01010 CPL_LESS_THAN,
01011 w_bound[10]));
01012
01013 cknull_nomsg(*lr75=sinfo_where_tab_min_max(*lambda,
01014 "WAVE",
01015 CPL_NOT_LESS_THAN,
01016 w_bound[10],
01017 CPL_LESS_THAN,
01018 w_bound[11]));
01019
01020 cknull_nomsg(*lr86=sinfo_where_tab_min_max(*lambda,
01021 "WAVE",
01022 CPL_NOT_LESS_THAN,
01023 w_bound[11],
01024 CPL_LESS_THAN,
01025 w_bound[12]));
01026
01027 cknull_nomsg(*lr97=sinfo_where_tab_min_max(*lambda,
01028 "WAVE",
01029 CPL_NOT_LESS_THAN,
01030 w_bound[12],
01031 CPL_LESS_THAN,
01032 w_bound[13]));
01033
01034 cknull_nomsg(*lr00=sinfo_where_tab_min_max(*lambda,
01035 "WAVE",
01036 CPL_NOT_LESS_THAN,
01037 w_bound[13],
01038 CPL_LESS_THAN,
01039 w_k_max));
01040
01041
01042 cknull_nomsg(*lrange=sinfo_where_tab_min_max(*lambda,
01043 "WAVE",
01044 CPL_NOT_LESS_THAN,
01045 w_j_min,
01046 CPL_NOT_GREATER_THAN,
01047 w_j_max));
01048
01049
01050
01051 cknull_nomsg(add1=sinfo_where_tab_min_max(*lambda,
01052 "WAVE",
01053 CPL_NOT_LESS_THAN,
01054 w_h_min,
01055 CPL_NOT_GREATER_THAN,
01056 w_h_max));
01057
01058 check_nomsg(nrow=cpl_table_get_nrow(*lrange));
01059 check_nomsg(cpl_table_insert(*lrange,add1,nrow));
01060 sinfo_free_table(&add1);
01061
01062 cknull_nomsg(add1=sinfo_where_tab_min_max(*lambda,
01063 "WAVE",
01064 CPL_NOT_LESS_THAN,
01065 w_k_min,
01066 CPL_NOT_GREATER_THAN,
01067 w_k_max));
01068
01069
01070 check_nomsg(nrow=cpl_table_get_nrow(*lrange));
01071 check_nomsg(cpl_table_insert(*lrange,add1,nrow));
01072 sinfo_free_table(&add1);
01073
01074
01075
01076 cknull_nomsg(*grange=sinfo_where_tab_min_max(*lambda,
01077 "WAVE",
01078 CPL_NOT_LESS_THAN,
01079 1.10,
01080 CPL_NOT_GREATER_THAN,
01081 1.35));
01082
01083
01084
01085
01086 cknull_nomsg(add1=sinfo_where_tab_min_max(*lambda,
01087 "WAVE",
01088 CPL_NOT_LESS_THAN,
01089 1.5,
01090 CPL_NOT_GREATER_THAN,
01091 1.7));
01092
01093 check_nomsg(nrow=cpl_table_get_nrow(*grange));
01094 check_nomsg(cpl_table_insert(*grange,add1,nrow));
01095 sinfo_free_table(&add1);
01096
01097
01098
01099 cknull_nomsg(add1=sinfo_where_tab_min_max(*lambda,
01100 "WAVE",
01101 CPL_NOT_LESS_THAN,
01102 2.0,
01103 CPL_NOT_GREATER_THAN,
01104 2.3));
01105
01106 check_nomsg(nrow=cpl_table_get_nrow(*grange));
01107 check_nomsg(cpl_table_insert(*grange,add1,nrow));
01108 sinfo_free_table(&add1);
01109
01110
01111
01112 cknull_nomsg(plist=cpl_propertylist_load(cpl_frame_get_filename(sky_frm),0));
01113 check_nomsg(crval_sky=sinfo_pfits_get_crval3(plist));
01114 check_nomsg(cdelt_sky=sinfo_pfits_get_cdelt3(plist));
01115 check_nomsg(crpix_sky=sinfo_pfits_get_crpix3(plist));
01116 check_nomsg(xsize_sky=sinfo_pfits_get_naxis1(plist));
01117 check_nomsg(ysize_sky=sinfo_pfits_get_naxis2(plist));
01118 check_nomsg(zsize_sky=sinfo_pfits_get_naxis3(plist));
01119 sinfo_free_propertylist(&plist);
01120
01121
01122 check_nomsg(*lambdas=cpl_table_new(zsize_sky));
01123 cpl_table_new_column(*lambdas,"WAVE",CPL_TYPE_DOUBLE);
01124 check_nomsg(sinfo_table_column_dindgen(lambdas,"WAVE"));
01125
01126 check_nomsg(cpl_table_add_scalar(*lambdas,"WAVE",1.));
01127 check_nomsg(cpl_table_subtract_scalar(*lambdas,"WAVE",crpix_sky));
01128 check_nomsg(cpl_table_multiply_scalar(*lambdas,"WAVE",cdelt_sky));
01129 check_nomsg(cpl_table_add_scalar(*lambdas,"WAVE",crval_sky));
01130
01131 check_nomsg(p=cpl_parameterlist_find(cfg,"sinfoni.sinfo_utl_skycor.mask_ws"));
01132 check_nomsg(ws=cpl_parameter_get_double(p));
01133 check_nomsg(p=cpl_parameterlist_find(cfg,"sinfoni.sinfo_utl_skycor.mask_we"));
01134 check_nomsg(we=cpl_parameter_get_double(p));
01135 if((ws != SINFO_MASK_WAVE_MIN) || (we != SINFO_MASK_WAVE_MAX)) {
01136 cknull_nomsg(*mrange=sinfo_where_tab_min_max(*lambda,"WAVE",
01137 CPL_NOT_LESS_THAN,ws,
01138 CPL_NOT_GREATER_THAN,we));
01139 } else {
01140 check_nomsg(*mrange=cpl_table_duplicate(*lrange));
01141 }
01142
01143
01144 check_nomsg(cpl_table_duplicate_column(*lambda,"WDIFF",*lambdas,"WAVE"));
01145 check_nomsg(cpl_table_subtract_columns(*lambda,"WDIFF","WAVE"));
01146 check_nomsg(mean=cpl_table_get_column_mean(*lambda,"WDIFF"));
01147 check_nomsg(nrow=cpl_table_get_nrow(*lambda));
01148 sinfo_msg_warning("diff %f",nrow*mean);
01149 if((fabs(nrow*mean) > 0) || (zsize_obj != zsize_sky)) {
01150 sinfo_msg("We have to interpolate sky frame - this is not good");
01151 *sky_interp_sw=1;
01152 }
01153
01154
01155 return 0;
01156
01157 cleanup:
01158 sinfo_free_table(&add1);
01159 sinfo_free_table(&tmp_tbl);
01160 sinfo_free_propertylist(&plist);
01161 return -1;
01162
01163 }
01164
01165
01176
01177
01178 static int
01179 sinfo_table_column_dindgen(cpl_table** t, const char* label)
01180 {
01181
01182 int sz=0;
01183 int i=0;
01184
01185 cknull(*t,"Null input vector");
01186 check(sz=cpl_table_get_nrow(*t),"Getting size of a vector");
01187 for(i=0;i<sz;i++) {
01188 cpl_table_set(*t,label,i,(double)i);
01189 }
01190
01191 return 0;
01192 cleanup:
01193 return -1;
01194
01195 }
01196
01197
01208
01209
01210
01211 int
01212 sinfo_sum_spectra(const cpl_frame* obj_frm,
01213 const cpl_frame* sky_frm,
01214 cpl_image* mask,
01215 cpl_table* wrange,
01216 cpl_table** int_obj,
01217 cpl_table** int_sky)
01218 {
01219
01220
01221
01222 cpl_image* obj_slice=NULL;
01223 cpl_image* sky_slice=NULL;
01224 cpl_image* gslice=NULL;
01225 cpl_image* pos_tmp=NULL;
01226 cpl_image* msk_tmp=NULL;
01227 cpl_imagelist* obj=NULL;
01228 cpl_imagelist* sky=NULL;
01229
01230
01231 cpl_table* loop=NULL;
01232 cpl_table* opos_tbl=NULL;
01233 cpl_table* spos_tbl=NULL;
01234 cpl_table* tmp_tbl=NULL;
01235 cpl_table* loop_tbl=NULL;
01236
01237 double med=0;
01238 double sdv=0;
01239 double avg=0;
01240
01241 int zsize=0;
01242 int i=0;
01243 int pos_i=0;
01244
01245
01246
01247 cknull_nomsg(obj=cpl_imagelist_load(cpl_frame_get_filename(obj_frm),
01248 CPL_TYPE_DOUBLE,0));
01249 cknull_nomsg(sky=cpl_imagelist_load(cpl_frame_get_filename(sky_frm),
01250 CPL_TYPE_DOUBLE,0));
01251
01252 check_nomsg(zsize=cpl_imagelist_get_size(obj));
01253 check_nomsg(*int_obj = cpl_table_new(zsize));
01254 check_nomsg(*int_sky = cpl_table_new(zsize));
01255 check_nomsg(cpl_table_duplicate_column(*int_obj,"WAVE",wrange,"WAVE"));
01256 check_nomsg(cpl_table_duplicate_column(*int_sky,"WAVE",wrange,"WAVE"));
01257 check_nomsg(cpl_table_new_column(*int_obj,"INT",CPL_TYPE_DOUBLE));
01258 check_nomsg(cpl_table_new_column(*int_sky,"INT",CPL_TYPE_DOUBLE));
01259 check_nomsg(cpl_table_fill_column_window_double(*int_obj,"INT",0,zsize,0));
01260 check_nomsg(cpl_table_fill_column_window_double(*int_sky,"INT",0,zsize,0));
01261
01262
01263
01264 cknull_nomsg(loop_tbl=sinfo_image2table(mask));
01265 check_nomsg(cpl_table_and_selected_double(loop_tbl,"VALUE",
01266 CPL_GREATER_THAN,0.5));
01267 check_nomsg(loop=cpl_table_extract_selected(loop_tbl));
01268 sinfo_free_table(&loop_tbl);
01269 sinfo_free_table(&loop);
01270
01271
01272 for (i=0;i<zsize;i++) {
01273 check_nomsg(obj_slice = cpl_imagelist_get(obj,i));
01274
01275
01276 pos_i=sinfo_cnt_mask_thresh_and_obj_finite(mask,0.5,obj_slice);
01277 if (pos_i >= 1) {
01278 if ((pos_i) < 3 ) {
01279
01280
01281
01282 check_nomsg(cpl_table_set_double(*int_obj,"INT",i,
01283 cpl_image_get_mean(obj_slice)));
01284 } else {
01285
01286
01287
01288
01289 check_nomsg(gslice = cpl_image_duplicate(obj_slice));
01290 check_nomsg(sinfo_image_flag_nan(&gslice));
01291
01292
01293
01294
01295
01296
01297
01298 check_nomsg(med = cpl_image_get_median(gslice));
01299 check_nomsg(sdv = cpl_image_get_stdev(gslice));
01300
01301
01302 check_nomsg(cpl_image_threshold(gslice,med-3*sdv,med+3*sdv,0,0));
01303 check_nomsg(avg= cpl_image_get_mean(gslice));
01304 sinfo_free_image(&gslice);
01305 check_nomsg(cpl_table_set_double(*int_obj,"INT",i,avg));
01306
01307 }
01308 }
01309
01310
01311 check_nomsg(sky_slice = cpl_imagelist_get(sky,i));
01312
01313 pos_i=sinfo_cnt_mask_thresh_and_obj_finite(mask,0.5,sky_slice);
01314 if (pos_i >= 1) {
01315 if ((pos_i) < 3) {
01316
01317
01318
01319 check_nomsg(cpl_table_set_double(*int_sky,"INT",i,
01320 cpl_image_get_mean(sky_slice)));
01321 } else {
01322
01323
01324
01325 check_nomsg(gslice = cpl_image_duplicate(sky_slice));
01326 check_nomsg(sinfo_image_flag_nan(&gslice));
01327
01328
01329 check_nomsg(med = cpl_image_get_median(gslice));
01330 check_nomsg(sdv = cpl_image_get_stdev(gslice));
01331
01332 check_nomsg(cpl_image_threshold(gslice,med-3*sdv,med+3*sdv,0,0));
01333 check_nomsg(avg= cpl_image_get_mean(gslice));
01334 sinfo_free_image(&gslice);
01335 check_nomsg(cpl_table_set_double(*int_sky,"INT",i,avg));
01336
01337
01338
01339
01340
01341
01342
01343 }
01344 }
01345 }
01346
01347 sinfo_free_imagelist(&obj);
01348 sinfo_free_imagelist(&sky);
01349
01350
01351 return 0;
01352
01353 cleanup:
01354 sinfo_free_image(&gslice);
01355 sinfo_free_image(&pos_tmp);
01356 sinfo_free_image(&msk_tmp);
01357 sinfo_free_table(&tmp_tbl);
01358 sinfo_free_table(&opos_tbl);
01359 sinfo_free_table(&spos_tbl);
01360 sinfo_free_table(&loop_tbl);
01361 sinfo_free_table(&loop);
01362 sinfo_free_imagelist(&obj);
01363 sinfo_free_imagelist(&sky);
01364
01365 return -1;
01366 }
01367
01368
01369
01370
01371
01372
01383
01384
01385
01386 static int
01387 sinfo_cnt_mask_thresh_and_obj_finite(const cpl_image* mask,
01388 const double t,
01389 const cpl_image* obj)
01390 {
01391
01392 int cnt=0;
01393 int sxm=0;
01394 int sym=0;
01395 int sxo=0;
01396 int syo=0;
01397 int i=0;
01398 double* pm=NULL;
01399 double* po=NULL;
01400
01401 check_nomsg(sxm=cpl_image_get_size_x(mask));
01402 check_nomsg(sym=cpl_image_get_size_y(mask));
01403 check_nomsg(sxo=cpl_image_get_size_x(obj));
01404 check_nomsg(syo=cpl_image_get_size_y(obj));
01405 if( sxm != sxo || sym != syo) {
01406 goto cleanup;
01407 }
01408 check_nomsg(pm=cpl_image_get_data_double(mask));
01409 check_nomsg(po=cpl_image_get_data_double(obj));
01410
01411 for(i=0;i<sxm*sym;i++) {
01412
01413 if( (pm[i] > t) && (!isnan(po[i]))) { cnt++; }
01414
01415 }
01416
01417 return cnt;
01418 cleanup:
01419 return -1;
01420
01421 }
01422
01423
01424
01425
01426
01427
01438
01439
01440
01441 static int
01442 sinfo_image_flag_nan(cpl_image** im)
01443 {
01444
01445 int cnt=0;
01446 int sx=0;
01447 int sy=0;
01448 int i=0;
01449 int j=0;
01450
01451 double* pi=NULL;
01452
01453 check_nomsg(sx=cpl_image_get_size_x(*im));
01454 check_nomsg(sy=cpl_image_get_size_y(*im));
01455 check_nomsg(pi=cpl_image_get_data_double(*im));
01456
01457 for(j=0;j<sy;j++) {
01458 for(i=0;i<sx;i++) {
01459 if(isnan(pi[j*sx+i])) {
01460 check_nomsg(cpl_image_reject(*im,i+1,j+1));
01461 cnt++;
01462 }
01463 }
01464 }
01465
01466 return cnt;
01467 cleanup:
01468 return -1;
01469
01470 }
01471
01472
01473
01474
01486
01487
01488 int
01489 sinfo_object_estimate_noise(cpl_frame* obj_frm,
01490 const int obj_noise_fit,
01491 double* centre,
01492 double* noise)
01493 {
01494
01495 const int nbins=HISTO_NBINS;
01496
01497 int xsz=0;
01498 int ysz=0;
01499 int zsz=0;
01500 int n=0;
01501 int i=0;
01502 int k=0;
01503 int r=0;
01504
01505 int max_h=0;
01506 int min_x=0;
01507 int max_x=0;
01508 int status=0;
01509 int max_pos=0;
01510 int min_pos=0;
01511 int min_xi_sz=0;
01512 int ndist=0;
01513
01514 double avg_d=0;
01515 double std_d=0;
01516 double hmin=0;
01517 double hmax=0;
01518 double kappa=3;
01519
01520 double* pdata=NULL;
01521 double* disth=NULL;
01522 double* distx=NULL;
01523
01524 double peak=0;
01525 double tempc=0;
01526 double value=0;
01527 double thres=0;
01528 double val=0;
01529 double x0=0;
01530 double sigma=0;
01531 double area=0;
01532 double offset=0;
01533
01534
01535
01536 cpl_propertylist* plist=NULL;
01537 cpl_imagelist* obj_cub=NULL;
01538 cpl_table* data_tbl=NULL;
01539 cpl_table* histo=NULL;
01540 cpl_image* img=NULL;
01541 cpl_table* dist=NULL;
01542 cpl_table* min_xi=NULL;
01543 cpl_table* tmp_tbl1=NULL;
01544 cpl_table* tmp_tbl2=NULL;
01545 cpl_vector* vx=NULL;
01546 cpl_vector* vy=NULL;
01547 cpl_vector* sx=NULL;
01548 cpl_vector* sy=NULL;
01549 int counter=0;
01550
01551
01552 cknull_nomsg(plist=cpl_propertylist_load(cpl_frame_get_filename(obj_frm),0));
01553 check_nomsg(xsz=sinfo_pfits_get_naxis1(plist));
01554 check_nomsg(ysz=sinfo_pfits_get_naxis2(plist));
01555 check_nomsg(zsz=sinfo_pfits_get_naxis3(plist));
01556 sinfo_free_propertylist(&plist);
01557
01558 cknull_nomsg(obj_cub=cpl_imagelist_load(cpl_frame_get_filename(obj_frm),
01559 CPL_TYPE_DOUBLE,0));
01560
01561 n=xsz*ysz*zsz;
01562 check_nomsg(data_tbl=cpl_table_new(n));
01563 check_nomsg(cpl_table_new_column(data_tbl,"DATA",CPL_TYPE_DOUBLE));
01564
01565
01566 for(k=0;k<zsz;k++) {
01567 check_nomsg(img=cpl_imagelist_get(obj_cub,k));
01568 check_nomsg(pdata=cpl_image_get_data(img));
01569 for(i=0;i<xsz*ysz;i++) {
01570 if(!isnan(pdata[i])) {
01571 cpl_table_set_double(data_tbl,"DATA",r,pdata[i]);
01572 r++;
01573 }
01574 }
01575 }
01576 sinfo_free_imagelist(&obj_cub);
01577
01578 check_nomsg(cpl_table_erase_invalid(data_tbl));
01579 check_nomsg(avg_d=cpl_table_get_column_mean(data_tbl,"DATA"));
01580 check_nomsg(std_d=cpl_table_get_column_stdev(data_tbl,"DATA"));
01581
01582
01583 hmin=avg_d-kappa*std_d;
01584 hmax=avg_d+kappa*std_d;
01585 sinfo_msg("mean=%f stdv=%f",avg_d,std_d);
01586 sinfo_msg("hmin=%f hmax=%f",hmin,hmax);
01587 sinfo_msg("Computes histogram");
01588 ck0(sinfo_histogram(data_tbl,nbins,hmin,hmax,&histo),"building histogram");
01589
01590 value=(double)(hmax-hmin)/nbins/2.;
01591
01592
01593
01594 while(min_xi_sz < HISTO_MIN_SIZE && counter < 10) {
01595 counter++;
01596 check_nomsg(max_h=cpl_table_get_column_max(histo,"HY"));
01597
01598 check_nomsg(max_pos=sinfo_table_get_index_of_max(histo,"HY",CPL_TYPE_INT));
01599
01600
01601
01602
01603
01604
01605
01606
01607 min_x=max_pos-1;
01608 max_x=max_pos+2;
01609
01610
01611 sinfo_free_table(&tmp_tbl1);
01612
01613
01614 check_nomsg(tmp_tbl1=sinfo_extract_table_rows(histo,"HL",
01615 CPL_LESS_THAN,
01616 cpl_table_get(histo,"HL",max_pos,&status)));
01617 thres=cpl_table_get_column_max(tmp_tbl1,"HY")/HISTO_Y_CUT;
01618
01619
01620
01621 sinfo_free_table(&min_xi);
01622 check_nomsg(min_xi=sinfo_extract_table_rows(tmp_tbl1,"HY",
01623 CPL_GREATER_THAN,thres));
01624
01625
01626
01627
01628
01629 min_xi_sz=cpl_table_get_nrow(min_xi);
01630 val=cpl_table_get(min_xi,"HL",0,&status);
01631
01632 check_nomsg(min_pos=sinfo_table_get_index_of_val(histo,"HL",val,
01633 CPL_TYPE_DOUBLE));
01634
01635
01636
01637
01638
01639 if (min_xi_sz > 0) {
01640 min_x = min_pos-HISTO_X_LEFT_CUT*(max_pos-min_pos);
01641 max_x = max_pos+HISTO_X_RIGHT_CUT*(max_pos-min_pos);
01642 }
01643
01644
01645 check_nomsg(hmin=sinfo_table_column_interpolate(histo,"HL",min_x));
01646 check_nomsg(hmax=sinfo_table_column_interpolate(histo,"HL",max_x));
01647
01648
01649 sinfo_free_table(&histo);
01650 ck0(sinfo_histogram(data_tbl,nbins,hmin,hmax,&histo),"building histogram");
01651
01652 check_nomsg(cpl_table_add_scalar(histo,"HL",(hmax-hmin)/nbins/2));
01653
01654
01655
01656
01657 }
01658 sinfo_free_table(&data_tbl);
01659 sinfo_free_table(&min_xi);
01660
01661
01662
01663 check_nomsg(peak=cpl_table_get_column_max(histo,"HY"));
01664
01665 sinfo_free_table(&tmp_tbl1);
01666
01667 check_nomsg(tmp_tbl1=sinfo_extract_table_rows(histo,"HY",CPL_EQUAL_TO,peak));
01668
01669
01670
01671
01672 check_nomsg(*centre=cpl_table_get_column_mean(tmp_tbl1,"HL"));
01673
01674
01675 sinfo_free_table(&tmp_tbl1);
01676 check_nomsg(tmp_tbl1=sinfo_extract_table_rows(histo,"HY",
01677 CPL_GREATER_THAN,
01678 peak/HISTO_Y_CUT));
01679 sinfo_free_table(&tmp_tbl2);
01680 check_nomsg(tmp_tbl2=sinfo_extract_table_rows(tmp_tbl1,"HY",
01681 CPL_LESS_THAN,peak));
01682 sinfo_free_table(&tmp_tbl1);
01683
01684 check_nomsg(tempc=*centre-cpl_table_get_column_min(tmp_tbl2,"HL"));
01685
01686 sinfo_free_table(&tmp_tbl2);
01687
01688 check_nomsg(dist=sinfo_where_tab_min_max(histo,"HL",
01689 CPL_GREATER_THAN,*centre-HISTO_DIST_TEMPC_MIN_FCT*tempc,
01690 CPL_NOT_GREATER_THAN,*centre+HISTO_DIST_TEMPC_MAX_FCT*tempc));
01691
01692 offset=cpl_table_get_column_min(histo,"HY");
01693 sinfo_free_table(&histo);
01694
01695
01696 check_nomsg(ndist=cpl_table_get_nrow(dist));
01697 check_nomsg(cpl_table_cast_column(dist,"HY","HYdouble",CPL_TYPE_DOUBLE));
01698 check_nomsg(disth=cpl_table_get_data_double(dist,"HYdouble"));
01699 check_nomsg(distx=cpl_table_get_data_double(dist,"HL"));
01700
01701
01702
01703
01704
01705 *noise=tempc/2;
01706
01707
01708
01709 if(obj_noise_fit == 1) {
01710 check_nomsg(vy=cpl_vector_wrap(ndist,disth));
01711 check_nomsg(vx=cpl_vector_wrap(ndist,distx));
01712 check_nomsg(sx=cpl_vector_new(ndist));
01713 check_nomsg(cpl_vector_fill(sx,1.));
01714 check_nomsg(sy=cpl_vector_duplicate(sx));
01715 x0=*centre;
01716 sigma=tempc/2;
01717
01718 check_nomsg(cpl_vector_fit_gaussian(vx,NULL,
01719 vy,NULL,
01720 CPL_FIT_ALL,
01721 &x0,&sigma,&area,&offset,
01722 NULL,NULL,NULL));
01723
01724
01725
01726
01727
01728 *noise=sigma;
01729 sinfo_unwrap_vector(&vx);
01730 sinfo_unwrap_vector(&vy);
01731 sinfoni_free_vector(&sx);
01732 sinfoni_free_vector(&sy);
01733 }
01734 sinfo_free_table(&dist);
01735
01736
01737 return 0;
01738
01739 cleanup:
01740 sinfo_free_imagelist(&obj_cub);
01741 sinfo_free_propertylist(&plist);
01742 sinfo_free_table(&min_xi);
01743 sinfo_free_table(&tmp_tbl1);
01744 sinfo_free_table(&tmp_tbl2);
01745 sinfo_free_table(&histo);
01746 sinfo_free_table(&dist);
01747 sinfo_free_table(&data_tbl);
01748 sinfoni_free_vector(&sx);
01749 sinfoni_free_vector(&sy);
01750 sinfo_unwrap_vector(&vx);
01751 sinfo_unwrap_vector(&vy);
01752
01753 return -1;
01754
01755 }
01756
01757
01770 static cpl_table*
01771 sinfo_where_tab_min_max(cpl_table* t,
01772 const char* col,
01773 cpl_table_select_operator op1,
01774 const double v1,
01775 cpl_table_select_operator op2,
01776 const double v2)
01777 {
01778
01779 cpl_table* res=NULL;
01780 cpl_table* tmp=NULL;
01781
01782 check_nomsg(cpl_table_and_selected_double(t,col,op1,v1));
01783 check_nomsg(tmp=cpl_table_extract_selected(t));
01784 check_nomsg(cpl_table_and_selected_double(tmp,col,op2,v2));
01785 check_nomsg(res=cpl_table_extract_selected(tmp));
01786 check_nomsg(cpl_table_select_all(t));
01787 sinfo_free_table(&tmp);
01788
01789 return res;
01790
01791 cleanup:
01792 sinfo_free_table(&tmp);
01793 sinfo_free_table(&res);
01794
01795 return NULL;
01796
01797 }
01798
01822
01823
01824 static int
01825 sinfo_histogram(const cpl_table* data,
01826 const int nbins,
01827 const double min,
01828 const double max,
01829 cpl_table** histo)
01830 {
01831 cpl_table* tmp_tbl1=NULL;
01832 cpl_table* tmp_tbl2=NULL;
01833 cpl_table* dat=NULL;
01834 int ntot=0;
01835 int i=0;
01836 int* phy=NULL;
01837 double* pdt=NULL;
01838
01839
01840 double vtmp=0;
01841 double vstp=0;
01842 double vmax=0;
01843 double vmin=0;
01844
01845 int h=0;
01846 check_nomsg(dat=cpl_table_duplicate(data));
01847 check_nomsg(cpl_table_cast_column(dat,"DATA","DDATA",CPL_TYPE_DOUBLE));
01848
01849
01850
01851
01852
01853 check_nomsg(cpl_table_and_selected_double(dat,"DDATA",
01854 CPL_NOT_GREATER_THAN,max));
01855
01856
01857
01858 check_nomsg(tmp_tbl1=cpl_table_extract_selected(dat));
01859
01860
01861
01862
01863
01864
01865
01866
01867
01868 check_nomsg(cpl_table_and_selected_double(tmp_tbl1,"DDATA",
01869 CPL_GREATER_THAN,min));
01870 check_nomsg(tmp_tbl2=cpl_table_extract_selected(tmp_tbl1));
01871
01872
01873
01874
01875
01876 sinfo_free_table(&tmp_tbl1);
01877
01878
01879
01880
01881
01882
01883
01884 check_nomsg(ntot=cpl_table_get_nrow(tmp_tbl2));
01885
01886
01887 check_nomsg(vmin=cpl_table_get_column_min(tmp_tbl2,"DDATA"));
01888 check_nomsg(vmax=cpl_table_get_column_max(tmp_tbl2,"DDATA"));
01889 vstp=(vmax-vmin)/(nbins-1);
01890
01891 check_nomsg(*histo=cpl_table_new(nbins));
01892 check_nomsg(cpl_table_new_column(*histo,"HX",CPL_TYPE_DOUBLE));
01893 check_nomsg(cpl_table_new_column(*histo,"HL",CPL_TYPE_DOUBLE));
01894 check_nomsg(cpl_table_new_column(*histo,"HY",CPL_TYPE_INT));
01895
01896
01897 check_nomsg(cpl_table_fill_column_window(*histo,"HL",0,nbins,0));
01898 check_nomsg(cpl_table_fill_column_window(*histo,"HY",0,nbins,0));
01899
01900 check_nomsg(phy=cpl_table_get_data_int(*histo,"HY"));
01901
01902 check_nomsg(pdt=cpl_table_get_data_double(dat,"DATA"));
01903
01904 for(i=0;i<nbins;i++) {
01905 cpl_table_set_double(*histo,"HX",i,(double)i);
01906 vtmp=vmin+i*vstp;
01907 cpl_table_set_double(*histo,"HL",i,vtmp);
01908 }
01909 h=0;
01910
01911 for(i=0;i<ntot;i++) {
01912 h=floor((pdt[i]-vmin)/vstp);
01913 if((h<nbins) && (h>-1)) {
01914 phy[h]++;
01915 }
01916 }
01917
01918
01919 sinfo_free_table(&tmp_tbl2);
01920 sinfo_free_table(&dat);
01921
01922
01923 return 0;
01924 cleanup:
01925 sinfo_free_table(&tmp_tbl1);
01926 sinfo_free_table(&tmp_tbl2);
01927 sinfo_free_table(&dat);
01928
01929 return -1;
01930
01931 }
01932
01933
01934
01944
01945
01946 int
01947 sinfo_object_flag_low_values(cpl_frame* obj_frm,
01948 const double cnt,
01949 const double sig,
01950 cpl_imagelist** flag_data)
01951 {
01952
01953 int xsz=0;
01954 int ysz=0;
01955 int zsz=0;
01956 int n=0;
01957 int i=0;
01958 int k=0;
01959 int r=0;
01960
01961 cpl_propertylist* plist=NULL;
01962 cpl_table* data_tbl=NULL;
01963 cpl_table* flag_tbl=NULL;
01964 cpl_image* img=NULL;
01965 cpl_imagelist* obj_cub=NULL;
01966
01967 double* pdata=NULL;
01968 double* pflag=NULL;
01969
01970
01971 cknull_nomsg(plist=cpl_propertylist_load(cpl_frame_get_filename(obj_frm),0));
01972 check_nomsg(xsz=sinfo_pfits_get_naxis1(plist));
01973 check_nomsg(ysz=sinfo_pfits_get_naxis2(plist));
01974 check_nomsg(zsz=sinfo_pfits_get_naxis3(plist));
01975 sinfo_free_propertylist(&plist);
01976
01977 cknull_nomsg(obj_cub=cpl_imagelist_load(cpl_frame_get_filename(obj_frm),
01978 CPL_TYPE_DOUBLE,0));
01979
01980 n=xsz*ysz*zsz;
01981 check_nomsg(data_tbl=cpl_table_new(n));
01982 check_nomsg(cpl_table_new_column(data_tbl,"DATA",CPL_TYPE_DOUBLE));
01983
01984 for(k=0;k<zsz;k++) {
01985 check_nomsg(img=cpl_imagelist_get(obj_cub,k));
01986 check_nomsg(pdata=cpl_image_get_data_double(img));
01987 for(i=0;i<xsz*ysz;i++) {
01988 if(!isnan(pdata[i])) {
01989 check_nomsg(cpl_table_set_double(data_tbl,"DATA",r,pdata[i]));
01990 r++;
01991 }
01992 }
01993 }
01994
01995 check_nomsg(cpl_table_erase_invalid(data_tbl));
01996 sinfo_msg("Background level: %f Noise: %f",cnt,sig);
01997 check_nomsg(cpl_table_and_selected_double(data_tbl,"DATA",
01998 CPL_LESS_THAN,cnt+2*sig));
01999 check_nomsg(flag_tbl=cpl_table_extract_selected(data_tbl));
02000 sinfo_free_table(&data_tbl);
02001
02002 sinfo_free_table(&flag_tbl);
02003
02004 check_nomsg(*flag_data=cpl_imagelist_new());
02005 for(i=0;i<zsz;i++) {
02006 check_nomsg(img=cpl_image_new(xsz,ysz,CPL_TYPE_DOUBLE));
02007 check_nomsg(cpl_image_add_scalar(img,0));
02008 check_nomsg(cpl_imagelist_set(*flag_data,cpl_image_duplicate(img),i));
02009 sinfo_free_image(&img);
02010 }
02011 for(k=0;k<zsz;k++) {
02012 check_nomsg(img=cpl_imagelist_get(*flag_data,k));
02013 pflag=cpl_image_get_data_double(cpl_imagelist_get(*flag_data,k));
02014 pdata=cpl_image_get_data_double(cpl_imagelist_get(obj_cub,k));
02015 for(i=0;i<xsz*ysz;i++) {
02016 if((!isnan(pdata[i])) && pdata[i] < (cnt+2*sig)) {
02017 pflag[i]=1;
02018 }
02019 }
02020 }
02021
02022 sinfo_free_imagelist(&obj_cub);
02023
02024
02025
02026
02027 return 0;
02028
02029 cleanup:
02030 sinfo_free_propertylist(&plist);
02031 sinfo_free_imagelist(&obj_cub);
02032 sinfo_free_table(&data_tbl);
02033 sinfo_free_table(&flag_tbl);
02034
02035 return -1;
02036 }
02037
02038
02052
02053
02054
02055 int
02056 sinfo_object_flag_sky_pixels(cpl_frame* obj_frm,
02057 cpl_table* lambda,
02058 cpl_table* mrange,
02059 cpl_imagelist* flag_data,
02060 const double tol,
02061 cpl_image** g_img,
02062 cpl_image** r_img,
02063 cpl_image** image)
02064 {
02065
02066 int xsz=0;
02067 int ysz=0;
02068 int zsz=0;
02069 int i=0;
02070 int j=0;
02071 int gpix_i=0;
02072 double tot=0;
02073 double all_pix=0;
02074 double flag_pix=0;
02075 double ratio=0;
02076
02077 double* pr_img=NULL;
02078 double* pg_img=NULL;
02079 double* pimage=NULL;
02080 cpl_propertylist* plist=NULL;
02081 cpl_imagelist* osel=NULL;
02082 cpl_imagelist* fsel=NULL;
02083 cpl_table* gpix=NULL;
02084 cpl_table* gspec=NULL;
02085 cpl_table* fspec=NULL;
02086 cpl_table* ospec=NULL;
02087
02088 cpl_imagelist* obj_cub=NULL;
02089
02090
02091 cknull_nomsg(plist=cpl_propertylist_load(cpl_frame_get_filename(obj_frm),0));
02092
02093 check_nomsg(xsz=sinfo_pfits_get_naxis1(plist));
02094 check_nomsg(ysz=sinfo_pfits_get_naxis2(plist));
02095 check_nomsg(zsz=sinfo_pfits_get_naxis3(plist));
02096 sinfo_free_propertylist(&plist);
02097 cknull_nomsg(obj_cub=cpl_imagelist_load(cpl_frame_get_filename(obj_frm),
02098 CPL_TYPE_DOUBLE,0));
02099
02100
02101
02102 check_nomsg(*r_img=cpl_image_new(xsz,ysz,CPL_TYPE_DOUBLE));
02103 check_nomsg(*g_img=cpl_image_new(xsz,ysz,CPL_TYPE_DOUBLE));
02104 check_nomsg(*image=cpl_image_new(xsz,ysz,CPL_TYPE_DOUBLE));
02105
02106 cknull_nomsg(pr_img=cpl_image_get_data_double(*r_img));
02107 cknull_nomsg(pg_img=cpl_image_get_data_double(*g_img));
02108 cknull_nomsg(pimage=cpl_image_get_data_double(*image));
02109
02110
02111
02112
02113
02114
02115
02116
02117
02118
02119 cknull_nomsg(osel=sinfo_imagelist_select_range(obj_cub,lambda,
02120 mrange,tol));
02121
02122 sinfo_free_imagelist(&obj_cub);
02123
02124 cknull_nomsg(fsel=sinfo_imagelist_select_range(flag_data,lambda,
02125 mrange,tol));
02126
02127
02128
02129
02130
02131
02132 for(j=0;j<ysz;j++) {
02133 for(i=0;i<xsz;i++) {
02134
02135 cknull_nomsg(ospec=sinfo_slice_z(osel,i,j));
02136 cknull_nomsg(fspec=sinfo_slice_z(fsel,i,j));
02137
02138 check_nomsg(gpix_i=sinfo_table_extract_finite(ospec,fspec,&gpix,&gspec));
02139
02140 if(gpix_i > 0) {
02141
02142 all_pix=(double)gpix_i;
02143
02144
02145
02146
02147
02148
02149
02150
02151
02152
02153
02154
02155
02156
02157
02158 check_nomsg(flag_pix=cpl_table_and_selected_double(gspec,"VALUE",
02159 CPL_GREATER_THAN,0.5));
02160
02161
02162
02163 ratio=flag_pix/all_pix;
02164
02165 if(all_pix > cpl_table_get_nrow(mrange)/2) {
02166 pg_img[i+j*xsz]=1;
02167 pr_img[i+j*xsz]=ratio;
02168 }
02169
02170 check_nomsg(pimage[i+j*xsz]=cpl_table_get_column_mean(gpix,"VALUE"));
02171
02172 }
02173 sinfo_free_table(&ospec);
02174 sinfo_free_table(&fspec);
02175 sinfo_free_table(&gpix);
02176 sinfo_free_table(&gspec);
02177
02178 }
02179 }
02180 sinfo_free_imagelist(&osel);
02181 sinfo_free_imagelist(&fsel);
02182
02183
02184
02185
02186
02187
02188
02189 check_nomsg(tot=cpl_image_get_flux(*g_img));
02190 if(tot < 1) {
02191 sinfo_msg_error("no good spaxel");
02192 goto cleanup;
02193 }
02194
02195 return 0;
02196
02197
02198 cleanup:
02199 sinfo_free_propertylist(&plist);
02200 sinfo_free_imagelist(&obj_cub);
02201 sinfo_free_imagelist(&osel);
02202 sinfo_free_imagelist(&fsel);
02203 sinfo_free_table(&ospec);
02204 sinfo_free_table(&fspec);
02205 sinfo_free_table(&gpix);
02206 sinfo_free_table(&gspec);
02207
02208 return -1;
02209
02210
02211 }
02212
02221 int
02222 sinfo_choose_good_sky_pixels(cpl_frame* obj_frm,
02223 cpl_image* r_img,
02224 cpl_image* g_img,
02225 const double min_frac,
02226 cpl_image** mask)
02227 {
02228
02229 int xsz=0;
02230 int ysz=0;
02231 int zsz=0;
02232 int r2i=0;
02233 int status=0;
02234 double tot=0;
02235 double thres=SKY_THRES;
02236 double cum_x_max=0;
02237
02238 cpl_image* r2img=NULL;
02239 cpl_propertylist* plist=NULL;
02240 cpl_table* cum=NULL;
02241 cpl_table* hcum=NULL;
02242 cpl_table* thres_tbl=NULL;
02243
02244 cknull_nomsg(plist=cpl_propertylist_load(cpl_frame_get_filename(obj_frm),0));
02245 check_nomsg(xsz=sinfo_pfits_get_naxis1(plist));
02246 check_nomsg(ysz=sinfo_pfits_get_naxis2(plist));
02247 check_nomsg(zsz=sinfo_pfits_get_naxis3(plist));
02248 sinfo_free_propertylist(&plist);
02249
02250
02251 check_nomsg(*mask=cpl_image_new(xsz,ysz,CPL_TYPE_DOUBLE));
02252
02253
02254 check_nomsg(r2img=cpl_image_duplicate(r_img));
02255 check_nomsg(cpl_image_threshold(r2img,thres,thres,0,1));
02256 check_nomsg(r2i=cpl_image_get_flux(r2img));
02257
02258 if(r2i>0) {
02259 sinfo_free_image(&(*mask));
02260 check_nomsg(*mask=cpl_image_duplicate(r2img));
02261 }
02262 sinfo_free_image(&r2img);
02263 check_nomsg(r2img=cpl_image_duplicate(r_img));
02264 check_nomsg(cpl_image_threshold(r2img,thres,CPL_PIXEL_MAXVAL,
02265 0,CPL_PIXEL_MAXVAL));
02266 sinfo_free_image(&r2img);
02267
02268 check_nomsg(tot=cpl_image_get_flux(g_img));
02269
02270
02271 sinfo_msg("%2.2d spaxels (%4.1f %% of good pixels) are designated as sky",
02272 r2i,100.*r2i/tot);
02273
02274
02275 if (1.*r2i/tot < min_frac) {
02276 sinfo_msg("this is too small - will increase it to %4.1f %%",
02277 100.*min_frac);
02278 check_nomsg(cum=cpl_table_new(xsz*ysz));
02279 check_nomsg(cpl_table_new_column(cum,"X",CPL_TYPE_DOUBLE));
02280 sinfo_table_column_dindgen(&cum,"X");
02281 check_nomsg(cpl_table_add_scalar(cum,"X",1.));
02282
02283
02284 hcum = sinfo_image2table(r_img);
02285 check_nomsg(sinfo_sort_table_1(hcum,"VALUE",FALSE));
02286
02287
02288 check_nomsg(cpl_table_duplicate_column(cum,"H",hcum,"VALUE"));
02289 check_nomsg(cum_x_max=cpl_table_get_column_max(cum,"X"));
02290 check_nomsg(cpl_table_duplicate_column(cum,"R",cum,"X"));
02291 check_nomsg(cpl_table_divide_scalar(cum,"R",cum_x_max));
02292 check_nomsg(cpl_table_and_selected_double(cum,"R",
02293 CPL_NOT_LESS_THAN,
02294 (1.-min_frac)));
02295 check_nomsg(thres_tbl=cpl_table_extract_selected(cum));
02296
02297 check_nomsg(thres = cpl_table_get(thres_tbl,"R",0,&status));
02298
02299 sinfo_free_image(&(*mask));
02300
02301
02302 check_nomsg(*mask=cpl_image_duplicate(r_img));
02303 check_nomsg(cpl_image_threshold(*mask,thres,thres,0,1));
02304 }
02305 sinfo_free_table(&cum);
02306 sinfo_free_table(&hcum);
02307 sinfo_free_table(&thres_tbl);
02308
02309 return 0;
02310 cleanup:
02311
02312 sinfo_free_propertylist(&plist);
02313 sinfo_free_image(&r2img);
02314 sinfo_free_table(&cum);
02315 sinfo_free_table(&hcum);
02316 sinfo_free_table(&thres_tbl);
02317
02318 return -1;
02319
02320 }
02321
02322
02331
02332
02333 static double
02334 sinfo_fit_bkg(double p[])
02335
02336 {
02337 double* px=NULL;
02338 double* py=NULL;
02339 double* pv=NULL;
02340 cpl_vector* vtmp=NULL;
02341 double max=0;
02342 int i=0;
02343 int np=0;
02344
02345 double chi2=0;
02346
02347 check_nomsg(px= cpl_vector_get_data(sa_vx));
02348 check_nomsg(py= cpl_vector_get_data(sa_vy));
02349 check_nomsg(np= cpl_vector_get_size(sa_vx));
02350 check_nomsg(vtmp=cpl_vector_duplicate(sa_vy));
02351 check_nomsg(pv=cpl_vector_get_data(vtmp));
02352
02353 for(i=0;i<np;i++) {
02354 pv[i]=sinfo_fac(px[i],p[2]);
02355
02356 }
02357 check_nomsg(max=cpl_vector_get_max(vtmp));
02358 if(max> 0) {
02359 check_nomsg(cpl_vector_divide_scalar(vtmp,max));
02360 check_nomsg(cpl_vector_multiply_scalar(vtmp,p[1]));
02361 check_nomsg(cpl_vector_add_scalar(vtmp,p[0]));
02362 }
02363
02364
02365 for(i=0;i<np;i++) {
02366 chi2+=(py[i]-pv[i])*(py[i]-pv[i]);
02367 }
02368 sinfoni_free_vector(&vtmp);
02369 return chi2;
02370 cleanup:
02371 sinfoni_free_vector(&vtmp);
02372 return -1;
02373
02374 }
02375
02376
02377
02389
02390
02391 int
02392 sinfo_thermal_background2(cpl_table* int_sky,
02393 cpl_table* lambda,
02394 cpl_table* lrange,
02395 cpl_table** bkg)
02396 {
02397
02398 int n2=0;
02399 int i=0;
02400 int j=0;
02401 int nrow=0;
02402
02403 cpl_table* tmp1=NULL;
02404 cpl_table* tmp2=NULL;
02405
02406 double max=0;
02407 double wmin=0;
02408 double wmax=0;
02409 double p0[3];
02410 const int MP=4;
02411 const int NP=3;
02412 double y[MP];
02413 double** ap=NULL;
02414 int nfunc=0;
02415 int status=0;
02416 int row=0;
02417 double bkg_min=0;
02418 double bkg_max=0;
02419 double p0_min=0;
02420 double p0_max=0;
02421 double p1_min=0;
02422 double p1_max=0;
02423 double p2_min=0;
02424 double p2_max=0;
02425 double* pw=NULL;
02426 double* pf=NULL;
02427
02428 ap=(double**) cpl_calloc(MP,sizeof(double*));
02429
02430 for(i=0;i<MP;i++) {
02431 ap[i]=cpl_calloc(NP,sizeof(double));
02432 }
02433
02434 cknull(int_sky,"Null input table sky");
02435 cknull(lambda,"Null input table lambda");
02436 cknull(lrange,"Null input table lrange");
02437
02438
02439
02440 check_nomsg(wmin=cpl_table_get_column_min(lrange,"WAVE"));
02441 check_nomsg(wmax=cpl_table_get_column_max(lrange,"WAVE"));
02442 check_nomsg(cpl_table_and_selected_double(int_sky,"WAVE",
02443 CPL_NOT_LESS_THAN,wmin));
02444 check_nomsg(cpl_table_and_selected_double(int_sky,"WAVE",
02445 CPL_NOT_GREATER_THAN,wmax));
02446 check_nomsg(tmp1=cpl_table_extract_selected(int_sky));
02447
02448 check_nomsg(row=sinfo_table_get_index_of_val(tmp1,"WAVE",
02449 wmax,CPL_TYPE_DOUBLE));
02450 check_nomsg(max=cpl_table_get_double(tmp1,"INT",row,&status));
02451 check_nomsg(sinfo_table_flag_nan(&tmp1,"INT"));
02452 check_nomsg(cpl_table_erase_invalid(tmp1));
02453 check_nomsg(cpl_table_and_selected_double(tmp1,"INT",CPL_NOT_EQUAL_TO,0));
02454 check_nomsg(tmp2=cpl_table_extract_selected(tmp1));
02455
02456 sinfo_free_table(&tmp1);
02457 check_nomsg(n2=cpl_table_get_nrow(tmp2));
02458 check_nomsg(sa_vx=cpl_vector_wrap(n2,
02459 cpl_table_get_data_double(tmp2,"WAVE")));
02460 check_nomsg(sa_vy=cpl_vector_wrap(n2,
02461 cpl_table_get_data_double(tmp2,"INT")));
02462
02463
02464 for(i=0;i<MP;i++) {
02465 for(j=0;j<NP;j++) {
02466 ap[i][j]=0;
02467 }
02468 }
02469
02470 check_nomsg(bkg_min=cpl_table_get_column_min(tmp2,"INT"));
02471 check_nomsg(bkg_max=cpl_table_get_double(tmp2,"INT",row,&status));
02472
02473
02474
02475 p0_min=bkg_min*0.9;
02476 p0_max=bkg_min*1.1;
02477 p1_min=bkg_max*0.9;
02478 p1_max=bkg_max*1.1;
02479 p1_min=200;
02480 p2_max=300;
02481
02482 ap[0][0]=p0_min; ap[0][1]=p1_min; ap[0][2]=p2_min;
02483 ap[1][0]=p0_max; ap[1][1]=p1_min; ap[1][2]=p2_min;
02484 ap[2][0]=p0_min; ap[2][1]=p1_max; ap[2][2]=p2_min;
02485 ap[3][0]=p0_min; ap[3][1]=p1_min; ap[3][2]=p2_max;
02486
02487 sinfo_msg("Before amoeba fit");
02488 for(i=0;i<MP;i++) {
02489 for(j=0;j<NP;j++) {
02490 sinfo_msg("ap[%d][%d]=%g",i,j,ap[i][j]);
02491 }
02492 }
02493
02494
02495
02496
02497 for(i=0;i<MP;i++) {
02498 p0[0]=ap[i][0];
02499 p0[1]=ap[i][1];
02500 p0[2]=ap[i][2];
02501 y[i]=sinfo_fit_bkg(p0);
02502 }
02503
02504 sinfo_msg("p0=%g %g %g",p0[0],p0[1],p0[2]);
02505 for(i=0;i<N_ITER_FIT_AMOEBA;i++) {
02506 check_nomsg(sinfo_fit_amoeba(ap,y,NP,AMOEBA_FTOL,sinfo_fit_bkg,&nfunc));
02507 sinfo_msg("After amoeba fit");
02508 sinfo_msg("iter=%d ap=%g %g %g",i,ap[0][0],ap[0][1],ap[0][2]);
02509 }
02510 sinfo_unwrap_vector(&sa_vx);
02511 sinfo_unwrap_vector(&sa_vy);
02512 sinfo_free_table(&tmp2);
02513
02514
02515 sinfo_msg("After amoeba fit");
02516 for(i=0;i<MP;i++) {
02517 for(j=0;j<NP;j++) {
02518 sinfo_msg("ap[%d][%d]=%g",i,j,ap[i][j]);
02519 }
02520 sinfo_msg("y[%d]=%g",i,y[i]);
02521 }
02522
02523
02524
02525 check_nomsg(nrow=cpl_table_get_nrow(lambda));
02526 check_nomsg(*bkg=cpl_table_new(nrow));
02527 check_nomsg(cpl_table_duplicate_column(*bkg,"WAVE",lambda,"WAVE"));
02528 check_nomsg(cpl_table_new_column(*bkg,"INT2",CPL_TYPE_DOUBLE));
02529 cpl_table_fill_column_window(*bkg,"INT2",0,nrow,0.);
02530 check_nomsg(pw=cpl_table_get_data_double(*bkg,"WAVE"));
02531 check_nomsg(pf=cpl_table_get_data_double(*bkg,"INT2"));
02532
02533 for(i=0;i<nrow;i++) {
02534 pf[i]=sinfo_fac(pw[i],ap[0][2]);
02535 }
02536 check_nomsg(max=cpl_table_get_column_max(*bkg,"INT2"));
02537
02538 if(max != 0) {
02539 check_nomsg(cpl_table_divide_scalar(*bkg,"INT2",max));
02540 }
02541 check_nomsg(cpl_table_multiply_scalar(*bkg,"INT2",ap[0][1]));
02542 check_nomsg(cpl_table_add_scalar(*bkg,"INT2",ap[0][0]));
02543
02544 sinfo_new_destroy_2Ddoublearray(&ap,MP);
02545
02546
02547 return 0;
02548
02549 cleanup:
02550 sinfo_new_destroy_2Ddoublearray(&ap,MP);
02551 sinfo_free_table(&tmp1);
02552 sinfo_free_table(&tmp2);
02553 sinfo_unwrap_vector(&sa_vx);
02554 sinfo_unwrap_vector(&sa_vy);
02555 return -1;
02556
02557 }
02558
02559
02560
02561
02573
02574
02575 int
02576 sinfo_thermal_background(cpl_table* int_sky,
02577 cpl_table* lambda,
02578 cpl_table* lrange,
02579 const double temp,
02580 const int niter,
02581 cpl_table** bkg,
02582 int* success_fit)
02583 {
02584
02585 int npix=0;
02586 int i=0;
02587 int row=0;
02588 const int ndim=3;
02589 int ia[ndim];
02590
02591 int NPOINTS=0;
02592
02593
02594 double temp1=0;
02595 double temp2=0;
02596
02597
02598
02599
02600 double max_tmp2=0;
02601 double* ptmp1=NULL;
02602 double thermal=0;
02603 double* pw=NULL;
02604 double p0[3];
02605 double wmin=0;
02606 double wmax=0;
02607 double ga0=0;
02608 double ga1=0;
02609
02610 double ga2=0;
02611 double mse=0;
02612 double chired=0;
02613
02614
02615 cpl_vector *a = cpl_vector_new(ndim);
02616 cpl_table* xlr=NULL;
02617 cpl_table* ylr=NULL;
02618 cpl_table* wlr=NULL;
02619 cpl_table* tmp=NULL;
02620 cpl_table* temp2_tbl=NULL;
02621
02622 cpl_vector* y=NULL;
02623 cpl_vector* sy=NULL;
02624
02625 cpl_matrix* x_matrix=NULL;
02626 double bkg_min=0;
02627 double bkg_max=0;
02628 int status=0;
02629 double avg=0;
02630 double sdv=0;
02631 double med=0;
02632
02633
02634
02635 check_nomsg(wmin=cpl_table_get_column_min(lrange,"WAVE"));
02636 check_nomsg(wmax=cpl_table_get_column_max(lrange,"WAVE"));
02637
02638 bkg_min=sinfo_fac(wmin,temp);
02639 bkg_max=sinfo_fac(wmax,temp);
02640
02641
02642
02643
02644 check_nomsg(cpl_table_and_selected_double(lambda,"WAVE",
02645 CPL_NOT_LESS_THAN,wmin));
02646 check_nomsg(tmp=cpl_table_extract_selected(lambda));
02647
02648 check_nomsg(cpl_table_and_selected_double(tmp,"WAVE",
02649 CPL_NOT_GREATER_THAN,wmax));
02650 check_nomsg(xlr=cpl_table_extract_selected(tmp));
02651 sinfo_free_table(&tmp);
02652
02653
02654 check_nomsg(cpl_table_and_selected_double(int_sky,"WAVE",
02655 CPL_NOT_LESS_THAN,wmin));
02656 check_nomsg(tmp=cpl_table_extract_selected(int_sky));
02657 check_nomsg(cpl_table_and_selected_double(tmp,"WAVE",
02658 CPL_NOT_GREATER_THAN,wmax));
02659
02660
02661
02662 check_nomsg(cpl_table_and_selected_double(tmp,"INT",CPL_GREATER_THAN,-2));
02663 check_nomsg(ylr=cpl_table_extract_selected(tmp));
02664
02665 sinfo_free_table(&tmp);
02666 check_nomsg(tmp=cpl_table_duplicate(ylr));
02667 sinfo_free_table(&ylr);
02668
02669 check_nomsg(avg=cpl_table_get_column_mean(tmp,"INT"));
02670 check_nomsg(sdv=cpl_table_get_column_stdev(tmp,"INT"));
02671 check_nomsg(cpl_table_and_selected_double(tmp,"INT",
02672 CPL_LESS_THAN,avg+10*sdv));
02673
02674 check_nomsg(ylr=cpl_table_extract_selected(tmp));
02675 sinfo_free_table(&tmp);
02676
02677
02678
02679
02680
02681
02682 check_nomsg(cpl_table_and_selected_double(ylr,"INT",CPL_NOT_EQUAL_TO,0));
02683
02684 check_nomsg(wlr=cpl_table_extract_selected(ylr));
02685
02686
02687 check_nomsg(p0[0]=cpl_table_get_column_min(wlr,"INT"));
02688 check_nomsg(row=sinfo_table_get_index_of_val(ylr,"WAVE",
02689 wmax,CPL_TYPE_DOUBLE));
02690 check_nomsg(p0[1]=cpl_table_get_double(ylr,"INT",row,&status));
02691 p0[2]=temp;
02692
02693
02694 ga0=p0[0];
02695 ga1=p0[1]/bkg_max;
02696
02697 ga2=p0[2];
02698
02699
02700 check_nomsg(sinfo_table_flag_nan(&wlr,"INT"));
02701 check_nomsg(cpl_table_erase_invalid(wlr));
02702
02703
02704
02705
02706
02707 check_nomsg(NPOINTS=cpl_table_get_nrow(ylr));
02708
02709 check_nomsg(x_matrix = cpl_matrix_wrap(NPOINTS,1,
02710 cpl_table_get_data_double(ylr,"WAVE")));
02711 check_nomsg(y=cpl_vector_wrap(NPOINTS,cpl_table_get_data_double(ylr,"INT")));
02712
02713
02714
02715 check_nomsg(cpl_vector_set(a, 0, ga0));
02716 check_nomsg(cpl_vector_set(a, 1, ga1));
02717 check_nomsg(cpl_vector_set(a, 2, ga2));
02718
02719 check_nomsg(sy=cpl_vector_duplicate(y));
02720 check_nomsg(cpl_vector_power(sy,2));
02721 check_nomsg(cpl_vector_power(sy,0.5));
02722
02723
02724 ia[0] = 1;
02725 ia[1] = 1;
02726 ia[2] = 1;
02727
02728
02729 for(i=0;i<niter;i++) {
02730
02731
02732
02733
02734
02735
02736
02737 if(CPL_ERROR_NONE != sinfo_fit_lm(x_matrix,NULL,y,sy,a,ia,sinfo_fitbkg,
02738 sinfo_fitbkg_derivative,
02739 &mse,&chired,NULL)) {
02740 sinfo_msg_warning("Thermal background fit failed");
02741 irplib_error_reset();
02742 cpl_error_reset();
02743 *success_fit=1;
02744
02745 goto recover;
02746 }
02747
02748 bkg_max=sinfo_fac(wmax,cpl_vector_get(a,2));
02749
02750
02751 sinfo_msg("after fit: a=%g %g %g chired=%g",
02752 cpl_vector_get(a,0),
02753 cpl_vector_get(a,1),
02754 cpl_vector_get(a,2),
02755 chired);
02756
02757
02758
02759 }
02760
02761 sinfo_msg("Last fit: a=%g %g %g chired=%g",
02762 cpl_vector_get(a,0),
02763 cpl_vector_get(a,1),
02764 cpl_vector_get(a,2),
02765 chired);
02766
02767 sinfo_unwrap_vector(&y);
02768 sinfoni_free_vector(&sy);
02769 sinfo_unwrap_matrix(&x_matrix);
02770 sinfo_free_table(&xlr);
02771 sinfo_free_table(&ylr);
02772 sinfo_free_table(&wlr);
02773
02774 ga0=cpl_vector_get(a,0);
02775 ga1=cpl_vector_get(a,1);
02776 ga2=cpl_vector_get(a,2);
02777
02778 check_nomsg(npix=cpl_table_get_nrow(lrange));
02779 check_nomsg(pw=cpl_table_get_data_double(lrange,"WAVE"));
02780 check_nomsg(temp2_tbl=cpl_table_new(npix));
02781 check_nomsg(cpl_table_new_column(temp2_tbl,"TEMP2",CPL_TYPE_DOUBLE));
02782
02783 for(i=0;i<npix;i++) {
02784 temp2=sinfo_fac(pw[i],ga2);
02785 check_nomsg(cpl_table_set_double(temp2_tbl,"TEMP2",i,temp2));
02786 }
02787 check_nomsg(max_tmp2=cpl_table_get_column_max(temp2_tbl,"TEMP2"));
02788 sinfo_free_table(&temp2_tbl);
02789
02790
02791
02792 check_nomsg(npix=cpl_table_get_nrow(lambda));
02793 check_nomsg(pw=cpl_table_get_data_double(lambda,"WAVE"));
02794 check_nomsg(*bkg=cpl_table_new(npix));
02795 check_nomsg(cpl_table_new_column(*bkg,"WAVE",CPL_TYPE_DOUBLE));
02796 check_nomsg(cpl_table_new_column(*bkg,"TEMP1",CPL_TYPE_DOUBLE));
02797 check_nomsg(cpl_table_new_column(*bkg,"INT",CPL_TYPE_DOUBLE));
02798 check_nomsg(cpl_table_new_column(*bkg,"INT2",CPL_TYPE_DOUBLE));
02799
02800 for(i=0;i<npix;i++) {
02801 check_nomsg(cpl_table_set_double(*bkg,"WAVE",i,pw[i]));
02802 temp1=sinfo_fac(pw[i],ga2);
02803 check_nomsg(cpl_table_set_double(*bkg,"TEMP1",i,temp1));
02804 }
02805
02806 check_nomsg(ptmp1=cpl_table_get_data_double(*bkg,"TEMP1"));
02807 bkg_max=sinfo_fac(wmax,ga2);
02808
02809 for(i=0;i<npix;i++) {
02810 thermal=ga0+ptmp1[i]/max_tmp2*ga1;
02811 check_nomsg(cpl_table_set_double(*bkg,"INT",i,thermal));
02812 thermal=ga0+ga1*sinfo_fac(pw[i],ga2);
02813 check_nomsg(cpl_table_set_double(*bkg,"INT2",i,thermal));
02814 }
02815 sinfoni_free_vector(&a);
02816
02817 return 0;
02818
02819 recover:
02820 sinfo_msg_warning("Recover fit of thermal background");
02821 check_nomsg(npix=cpl_table_get_nrow(lambda));
02822 check_nomsg(*bkg=cpl_table_new(npix));
02823 check_nomsg(cpl_table_new_column(*bkg,"TEMP1",CPL_TYPE_DOUBLE));
02824 check_nomsg(cpl_table_new_column(*bkg,"INT",CPL_TYPE_DOUBLE));
02825 check_nomsg(cpl_table_new_column(*bkg,"INT2",CPL_TYPE_DOUBLE));
02826 med=cpl_table_get_column_median(ylr,"INT");
02827 for(i=0;i<npix;i++) {
02828 check_nomsg(cpl_table_set_double(*bkg,"INT",i,med));
02829 check_nomsg(cpl_table_set_double(*bkg,"INT2",i,med));
02830 }
02831
02832 sinfoni_free_vector(&a);
02833 sinfo_unwrap_vector(&y);
02834 sinfoni_free_vector(&sy);
02835 sinfo_unwrap_matrix(&x_matrix);
02836 sinfo_free_table(&xlr);
02837 sinfo_free_table(&ylr);
02838 sinfo_free_table(&wlr);
02839 sinfo_free_table(&tmp);
02840
02841 return 0;
02842
02843
02844 cleanup:
02845 sinfoni_free_vector(&a);
02846 sinfo_unwrap_vector(&y);
02847 sinfoni_free_vector(&sy);
02848 sinfo_unwrap_matrix(&x_matrix);
02849
02850 sinfo_free_table(&xlr);
02851 sinfo_free_table(&ylr);
02852 sinfo_free_table(&wlr);
02853 sinfo_free_table(&tmp);
02854 sinfo_free_table(&temp2_tbl);
02855
02856 return -1;
02857
02858 }
02868 static cpl_table*
02869 sinfo_slice_z(const cpl_imagelist* cin,const int i,const int j)
02870 {
02871 int sx=0;
02872 int sy=0;
02873 int sz=0;
02874 int k=0;
02875 cpl_table* tout=NULL;
02876 cpl_image* img=NULL;
02877 double* pim=NULL;
02878
02879 cknull(cin,"null input imagelist");
02880 check_nomsg(sz=cpl_imagelist_get_size(cin));
02881 check_nomsg(img=cpl_imagelist_get(cin,0));
02882 check_nomsg(sx=cpl_image_get_size_x(img));
02883 check_nomsg(sy=cpl_image_get_size_y(img));
02884 check_nomsg(tout=cpl_table_new(sz));
02885 check_nomsg(cpl_table_new_column(tout,"VALUE",CPL_TYPE_DOUBLE));
02886 for(k=0;k<sz;k++) {
02887 check_nomsg(img=cpl_imagelist_get(cin,k));
02888 check_nomsg(pim=cpl_image_get_data_double(img));
02889 check_nomsg(cpl_table_set(tout,"VALUE",k,pim[j*sx+i]));
02890 }
02891
02892 return tout;
02893 cleanup:
02894 sinfo_free_table(&tout);
02895
02896 return NULL;
02897
02898 }
02899
02910 double
02911 sinfo_xcorr(cpl_table* int_obj,
02912 cpl_table* int_sky,
02913 cpl_table* lambda,
02914 const double dispersion,
02915 const double line_hw)
02916 {
02917
02918 cpl_table* z=NULL;
02919 cpl_table* tmp_sky=NULL;
02920
02921 cpl_table* z_diff=NULL;
02922 cpl_table* z_pos=NULL;
02923
02924 int z_ext=0;
02925 double z_mean=0;
02926 double z_sdv=0;
02927 int nrow=0;
02928 int i=0;
02929
02930 double g_lam=0;
02931 double g_err=0;
02932
02933 double sky_max=0;
02934
02935 double* pint=NULL;
02936 int* ppos=NULL;
02937 int status=0;
02938 int zsize=0;
02939 int iq=0;
02940
02941 double g_diff=0;
02942 int jz=0;
02943 int z1=0;
02944 int nfit=0;
02945 int npos=0;
02946 cpl_table* z_good=NULL;
02947 cpl_table* w_tbl=NULL;
02948 cpl_table* o_tbl=NULL;
02949 cpl_table* s_tbl=NULL;
02950 cpl_vector* vw=NULL;
02951 cpl_vector* vs=NULL;
02952 cpl_vector* vo=NULL;
02953 cpl_vector* sx=NULL;
02954 cpl_vector* sy=NULL;
02955
02956
02957
02958 double zfit=0;
02959
02960
02961 double mse=0;
02962
02963 double ws=0;
02964 double wc_s=0;
02965 double sig_s=0;
02966 double bkg_s=0;
02967 double amp_s=0;
02968 double area_s=0;
02969
02970 double wo=0;
02971 double wc_o=0;
02972 double sig_o=0;
02973 double bkg_o=0;
02974 double amp_o=0;
02975 double area_o=0;
02976
02977 cpl_polynomial* cfit=NULL;
02978 int pows[2];
02979 cpl_vector* vx=NULL;
02980 cpl_vector* vy=NULL;
02981
02982
02983
02984
02985
02986
02987 zsize=cpl_table_get_nrow(int_obj);
02988 check_nomsg(z = cpl_table_duplicate(int_sky));
02989 ck0_nomsg(sinfo_table_flag_nan(&z,"INT"));
02990
02991 check_nomsg(z_mean=cpl_table_get_column_mean(z,"INT"));
02992 if(z_mean < 0) {
02993 check_nomsg(cpl_table_multiply_scalar(z,"INT",-1));
02994 }
02995
02996
02997
02998
02999 check_nomsg(tmp_sky=cpl_table_duplicate(int_sky));
03000 ck0_nomsg(sinfo_table_flag_nan(&tmp_sky,"INT"));
03001 check_nomsg(sky_max=cpl_table_get_column_max(tmp_sky,"INT"));
03002 sinfo_free_table(&tmp_sky);
03003
03004
03005 check_nomsg(nrow=cpl_table_get_nrow(z));
03006 check_nomsg(pint=cpl_table_get_data_double(z,"INT"));
03007 check_nomsg(sky_max=cpl_table_get_column_max(z,"INT"));
03008 for(i=0;i<nrow;i++) {
03009 if(pint[i]<sky_max/SKY_LINE_MIN_CUT) {
03010 pint[i]=0;
03011 }
03012 }
03013
03014
03015
03016
03017
03018 check_nomsg(z_diff=cpl_table_duplicate(z));
03019 check_nomsg(cpl_table_duplicate_column(z_diff,"INT1",z,"INT"));
03020 check_nomsg(cpl_table_duplicate_column(z_diff,"INT2",z,"INT"));
03021 check_nomsg(cpl_table_shift_column(z_diff,"INT1",-1));
03022 check_nomsg(cpl_table_duplicate_column(z_diff,"DIFF",z_diff,"INT2"));
03023 check_nomsg(cpl_table_subtract_columns(z_diff,"DIFF","INT1"));
03024
03025 check_nomsg(cpl_table_erase_window(z_diff,nrow-2,2));
03026
03027
03028
03029 check_nomsg(cpl_table_new_column(z_diff,"POS",CPL_TYPE_INT));
03030 check_nomsg(cpl_table_fill_column_window_int(z_diff,"POS",0,nrow,0));
03031
03032 check_nomsg(pint=cpl_table_get_data_double(z_diff,"DIFF"));
03033 check_nomsg(ppos=cpl_table_get_data_int(z_diff,"POS"));
03034 check_nomsg(nrow=cpl_table_get_nrow(z_diff));
03035 for(i=1;i<nrow;i++) {
03036 if(!isnan(pint[i]) && (pint[i]>0 && pint[i-1]<0)) {
03037 ppos[i]=i;
03038 }
03039 }
03040
03041
03042 check_nomsg(cpl_table_select_all(z_diff));
03043 check_nomsg(cpl_table_and_selected_int(z_diff,"POS",CPL_GREATER_THAN,0));
03044 check_nomsg(z_pos=cpl_table_extract_selected(z_diff));
03045 sinfo_free_table(&z_diff);
03046
03047
03048
03049
03050
03051
03052
03053
03054
03055
03056 g_lam = 0.;
03057 g_diff = 0.;
03058 g_err = 0.;
03059 check_nomsg(npos=cpl_table_get_nrow(z_pos));
03060 z_ext = line_hw ;
03061 check_nomsg(cpl_table_new_column(z_pos,"STATUS",CPL_TYPE_INT));
03062 check_nomsg(cpl_table_fill_column_window_int(z_pos,"STATUS",0,npos,0));
03063 check_nomsg(cpl_table_new_column(z_pos,"SIGS",CPL_TYPE_DOUBLE));
03064 check_nomsg(cpl_table_new_column(z_pos,"WAVES",CPL_TYPE_DOUBLE));
03065 check_nomsg(cpl_table_new_column(z_pos,"BKGS",CPL_TYPE_DOUBLE));
03066 check_nomsg(cpl_table_new_column(z_pos,"AREAS",CPL_TYPE_DOUBLE));
03067 check_nomsg(cpl_table_new_column(z_pos,"AMPS",CPL_TYPE_DOUBLE));
03068
03069
03070 check_nomsg(cpl_table_new_column(z_pos,"SIGO",CPL_TYPE_DOUBLE));
03071 check_nomsg(cpl_table_new_column(z_pos,"WAVEO",CPL_TYPE_DOUBLE));
03072 check_nomsg(cpl_table_new_column(z_pos,"BKGO",CPL_TYPE_DOUBLE));
03073 check_nomsg(cpl_table_new_column(z_pos,"AREAO",CPL_TYPE_DOUBLE));
03074 check_nomsg(cpl_table_new_column(z_pos,"AMPO",CPL_TYPE_DOUBLE));
03075
03076 check_nomsg(cpl_table_new_column(z_pos,"WAVEC",CPL_TYPE_DOUBLE));
03077 check_nomsg(cpl_table_new_column(z_pos,"WDIF",CPL_TYPE_DOUBLE));
03078 check_nomsg(cpl_table_new_column(z_pos,"ERR",CPL_TYPE_DOUBLE));
03079
03080
03081 nfit=2*z_ext+1;
03082
03083
03084
03085
03086 for (jz=0;jz<npos;jz++) {
03087 check_nomsg(z1 = cpl_table_get_int(z_pos,"POS",jz,&status));
03088
03089
03090 if((z1-z_ext) > 0 && (z1+z_ext) < zsize) {
03091 check_nomsg(cpl_table_select_all(int_sky));
03092 check_nomsg(cpl_table_select_all(int_obj));
03093 check_nomsg(cpl_table_select_all(lambda));
03094 check_nomsg(cpl_table_and_selected_window(int_sky,z1-z_ext,nfit));
03095 check_nomsg(s_tbl=cpl_table_extract_selected(int_sky));
03096 check_nomsg(cpl_table_and_selected_window(lambda,z1-z_ext,nfit));
03097 check_nomsg(w_tbl=cpl_table_extract_selected(lambda));
03098 check_nomsg(cpl_table_and_selected_window(int_obj,z1-z_ext,nfit));
03099 check_nomsg(o_tbl=cpl_table_extract_selected(int_obj));
03100
03101
03102 check_nomsg(vw=cpl_vector_wrap(nfit,
03103 cpl_table_get_data_double(w_tbl,"WAVE")));
03104 check_nomsg(vs=cpl_vector_wrap(nfit,
03105 cpl_table_get_data_double(s_tbl,"INT")));
03106 check_nomsg(vo=cpl_vector_wrap(nfit,
03107 cpl_table_get_data_double(o_tbl,"INT")));
03108
03109 check_nomsg(sx=cpl_vector_new(nfit));
03110 check_nomsg(cpl_vector_fill(sx,10.));
03111 check_nomsg(sy=cpl_vector_duplicate(sx));
03112
03113 check_nomsg(ws=cpl_table_get_double(lambda,"WAVE",z1,&status));
03114 check_nomsg(amp_s=cpl_table_get_double(z_pos,"INT",jz,&status));
03115 wc_s=ws;
03116 sig_s=z_ext;
03117 area_s=2*z_ext;
03118 bkg_s=0;
03119 if(wc_s < 2.35) {
03120
03121 if(CPL_ERROR_NONE == cpl_vector_fit_gaussian(vw,NULL,
03122 vs,NULL,
03123 CPL_FIT_ALL,
03124 &wc_s,&sig_s,
03125 &area_s,&bkg_s,
03126 NULL,NULL,NULL)) {
03127 amp_s=sinfo_gaussian(area_s,sig_s,ws,wc_s,bkg_s);
03128
03129
03130
03131
03132
03133
03134
03135 check_nomsg(cpl_table_set_double(z_pos,"WAVES",jz,wc_s));
03136 check_nomsg(cpl_table_set_double(z_pos,"SIGS",jz,sig_s));
03137 check_nomsg(cpl_table_set_double(z_pos,"AREAS",jz,area_s));
03138 check_nomsg(cpl_table_set_double(z_pos,"BKGS",jz,bkg_s));
03139 check_nomsg(cpl_table_set_double(z_pos,"AMPS",jz,amp_s));
03140 } else {
03141 cpl_error_reset();
03142 check_nomsg(cpl_table_set_double(z_pos,"WAVES",jz,wc_s));
03143 check_nomsg(cpl_table_set_double(z_pos,"SIGS",jz,sig_s));
03144 check_nomsg(cpl_table_set_double(z_pos,"AREAS",jz,area_s));
03145 check_nomsg(cpl_table_set_double(z_pos,"BKGS",jz,bkg_s));
03146 check_nomsg(cpl_table_set_double(z_pos,"AMPS",jz,amp_s));
03147 check_nomsg(cpl_table_set_int(z_pos,"STATUS",jz,-1));
03148 }
03149 check_nomsg(wo=cpl_table_get_double(lambda,"WAVE",z1,&status));
03150 check_nomsg(amp_o=cpl_table_get_double(z_pos,"INT",jz,&status));
03151 wc_o=wo;
03152 sig_o=z_ext;
03153 area_o=2*z_ext;
03154 bkg_o=0;
03155 if(CPL_ERROR_NONE == cpl_vector_fit_gaussian(vw,NULL,
03156 vo,sy,
03157 CPL_FIT_ALL,
03158 &wc_o,&sig_o,
03159 &area_o,&bkg_o,
03160 NULL,NULL,NULL)) {
03161
03162 amp_o=sinfo_gaussian(area_o,sig_o,wo,wc_o,bkg_o);
03163
03164
03165
03166
03167
03168
03169
03170 check_nomsg(cpl_table_set_double(z_pos,"WAVEO",jz,wc_o));
03171 check_nomsg(cpl_table_set_double(z_pos,"SIGO",jz,sig_o));
03172 check_nomsg(cpl_table_set_double(z_pos,"AREAO",jz,area_o));
03173 check_nomsg(cpl_table_set_double(z_pos,"BKGO",jz,bkg_o));
03174 check_nomsg(cpl_table_set_double(z_pos,"AMPO",jz,amp_o));
03175 } else {
03176 cpl_error_reset();
03177 check_nomsg(cpl_table_set_double(z_pos,"WAVEO",jz,wc_o));
03178 check_nomsg(cpl_table_set_double(z_pos,"SIGO",jz,sig_o));
03179 check_nomsg(cpl_table_set_double(z_pos,"AREAO",jz,area_o));
03180 check_nomsg(cpl_table_set_double(z_pos,"BKGO",jz,bkg_o));
03181 check_nomsg(cpl_table_set_double(z_pos,"AMPO",jz,amp_o));
03182 check_nomsg(cpl_table_set_int(z_pos,"STATUS",jz,-2));
03183
03184
03185
03186
03187
03188
03189
03190
03191
03192
03193 }
03194 check_nomsg(cpl_table_set_double(z_pos,"ERR",
03195 jz,sqrt(sig_s*sig_s+sig_o*sig_o)));
03196 check_nomsg(cpl_table_set_double(z_pos,"WDIF",jz,wc_s-wc_o));
03197 check_nomsg(cpl_table_set_double(z_pos,"WAVEC",jz,(wc_o+wc_s)/2));
03198
03199 } else {
03200 check_nomsg(cpl_table_set_int(z_pos,"STATUS",jz,-3));
03201 }
03202 sinfo_unwrap_vector(&vw);
03203 sinfo_unwrap_vector(&vs);
03204 sinfo_unwrap_vector(&vo);
03205 sinfoni_free_vector(&sx);
03206 sinfoni_free_vector(&sy);
03207 sinfo_free_table(&w_tbl);
03208 sinfo_free_table(&s_tbl);
03209 sinfo_free_table(&o_tbl);
03210 }
03211 }
03212
03213
03214 check_nomsg(cpl_table_duplicate_column(z_pos,"YDIF",z_pos,"WDIF"));
03215 check_nomsg(cpl_table_divide_scalar(z_pos,"YDIF",dispersion));
03216
03217
03218
03219
03220 check_nomsg(cpl_table_and_selected_int(z_pos,"STATUS",CPL_GREATER_THAN,-1));
03221
03222
03223
03224 check_nomsg(z_good=cpl_table_extract_selected(z_pos));
03225 check_nomsg(npos=cpl_table_get_nrow(z_good));
03226 sinfo_free_table(&z_pos);
03227 if(npos == 0) {
03228 return 0;
03229 }
03230 check_nomsg(z_pos=cpl_table_duplicate(z_good));
03231 check_nomsg(z_mean = cpl_table_get_column_median(z_pos,"WDIF"));
03232 check_nomsg(z_sdv = cpl_table_get_column_stdev(z_pos,"WDIF"));
03233
03234
03235
03236 check_nomsg(cpl_table_duplicate_column(z_pos,"CHECK",
03237 z_pos,"WDIF"));
03238
03239 cpl_table_erase_column(z_pos,"AMPO");
03240 cpl_table_erase_column(z_pos,"SIGO");
03241 cpl_table_erase_column(z_pos,"AREAO");
03242 cpl_table_erase_column(z_pos,"BKGO");
03243 cpl_table_erase_column(z_pos,"WAVEO");
03244 cpl_table_erase_column(z_pos,"AMPS");
03245 cpl_table_erase_column(z_pos,"SIGS");
03246 cpl_table_erase_column(z_pos,"AREAS");
03247 cpl_table_erase_column(z_pos,"BKGS");
03248 cpl_table_erase_column(z_pos,"WAVES");
03249 cpl_table_erase_column(z_pos,"STATUS");
03250 cpl_table_erase_column(z_pos,"INT");
03251 cpl_table_erase_column(z_pos,"INT1");
03252 cpl_table_erase_column(z_pos,"INT2");
03253 cpl_table_erase_column(z_pos,"ERR");
03254 cpl_table_erase_column(z_pos,"POS");
03255 cpl_table_erase_column(z_pos,"DIFF");
03256
03257
03258
03259
03260
03261 sinfo_msg("ks-clip1");
03262 sinfo_msg("iter mean (um) sdv (um) mean (pix) sdv (pix)");
03263
03264
03265
03266 for (iq = 0;iq<XCOR_YSHIFT_KS_CLIP;iq++) {
03267
03268 sinfo_msg(" %d %3.2g %3.2g %5.4g %5.4g",
03269 iq,z_mean,z_sdv,z_mean/dispersion,z_sdv/dispersion);
03270
03271
03272 check_nomsg(cpl_table_subtract_scalar(z_pos,"CHECK",z_mean));
03273 check_nomsg(cpl_table_duplicate_column(z_pos,"CHECKW",z_pos,"CHECK"));
03274 check_nomsg(cpl_table_multiply_columns(z_pos,"CHECKW","CHECK"));
03275 check_nomsg(cpl_table_power_column(z_pos,"CHECKW",0.5));
03276 check_nomsg(cpl_table_add_scalar(z_pos,"CHECK",z_mean));
03277 check_nomsg(cpl_table_and_selected_double(z_pos,"CHECKW",
03278 CPL_NOT_GREATER_THAN,2*z_sdv));
03279 sinfo_free_table(&z_good);
03280 check_nomsg(z_good=cpl_table_extract_selected(z_pos));
03281
03282 check_nomsg(cpl_table_select_all(z_pos));
03283
03284
03285 check_nomsg(z_mean = cpl_table_get_column_median(z_good,"WDIF"));
03286 check_nomsg(z_sdv = cpl_table_get_column_stdev(z_good,"WDIF"));
03287 sinfo_free_table(&z_good);
03288 check_nomsg(cpl_table_erase_column(z_pos,"CHECKW"));
03289
03290 }
03291
03292
03293
03294
03295
03296
03297
03298
03299
03300
03301
03302
03303 cpl_table_select_all(z_pos);
03304 check_nomsg(cpl_table_new_column(z_pos,"ZFIT",CPL_TYPE_DOUBLE));
03305 check_nomsg(nfit=cpl_table_get_nrow(z_pos));
03306 check_nomsg(cpl_table_fill_column_window(z_pos,"ZFIT",0,nfit,0));
03307
03308 check_nomsg(z_good=cpl_table_duplicate(z_pos));
03309
03310
03311 sinfo_msg("ks-clip2");
03312 sinfo_msg("iter mean (um) sdv (um) mean (pix) sdv (pix)");
03313 check_nomsg(cpl_table_select_all(z_good));
03314
03315 for(iq=0;iq<XCOR_YSHIFT_KS_CLIP;iq++) {
03316
03317 check_nomsg(nfit=cpl_table_get_nrow(z_good));
03318
03319 if(nfit>0) {
03320 check_nomsg(vx=cpl_vector_wrap(nfit,
03321 cpl_table_get_data_double(z_good,"WAVE")));
03322 check_nomsg(vy=cpl_vector_wrap(nfit,
03323 cpl_table_get_data_double(z_good,"WDIF")));
03324 check_nomsg(cfit=cpl_polynomial_fit_1d_create(vx,vy,0,&mse));
03325 pows[0]=0;
03326 pows[1]=0;
03327 check_nomsg(zfit=cpl_polynomial_get_coeff(cfit,pows));
03328 sinfo_free_polynomial(&cfit);
03329
03330
03331
03332
03333 check_nomsg(cpl_table_fill_column_window(z_good,"ZFIT",0,nfit,zfit));
03334 check_nomsg(cpl_table_duplicate_column(z_good,"WRES",z_good,"WDIF"));
03335 check_nomsg(cpl_table_subtract_columns(z_good,"WRES","ZFIT"));
03336
03337 check_nomsg(z_sdv=cpl_table_get_column_stdev(z_good,"WRES"));
03338 check_nomsg(z_mean=cpl_table_get_column_mean(z_good,"WDIF"));
03339
03340 sinfo_msg(" %d %3.2g %3.2g %5.4g %5.4g",
03341 iq,z_mean,z_sdv,z_mean/dispersion,z_sdv/dispersion);
03342
03343 check_nomsg(nfit=cpl_table_get_nrow(z_pos));
03344 check_nomsg(cpl_table_fill_column_window(z_pos,"ZFIT",0,nfit,zfit));
03345 check_nomsg(cpl_table_duplicate_column(z_pos,"WRES",z_pos,"WDIF"));
03346 check_nomsg(cpl_table_subtract_columns(z_pos,"WRES","ZFIT"));
03347
03348 check_nomsg(cpl_table_multiply_columns(z_pos,"WRES","WRES"));
03349 check_nomsg(cpl_table_power_column(z_pos,"WRES",0.5));
03350
03351
03352
03353
03354
03355
03356
03357 check_nomsg(cpl_table_and_selected_double(z_pos,"WRES",
03358 CPL_NOT_GREATER_THAN,3*z_sdv));
03359
03360 check_nomsg(sinfo_free_table(&z_good));
03361 check_nomsg(z_good=cpl_table_extract_selected(z_pos));
03362
03363
03364 check_nomsg(cpl_table_select_all(z_pos));
03365 check_nomsg(cpl_table_select_all(z_good));
03366 check_nomsg(cpl_table_erase_column(z_good,"WRES"));
03367 check_nomsg(cpl_table_erase_column(z_pos,"WRES"));
03368
03369 sinfo_unwrap_vector(&vx);
03370 sinfo_unwrap_vector(&vy);
03371
03372 }
03373
03374 }
03375
03376
03377
03378 sinfo_unwrap_vector(&vx);
03379 sinfo_unwrap_vector(&vy);
03380 sinfo_free_polynomial(&cfit);
03381 sinfo_free_table(&z);
03382 sinfo_free_table(&z_pos);
03383 sinfo_free_table(&z_good);
03384
03385
03386 return zfit;
03387 cleanup:
03388
03389 sinfo_free_table(&z_good);
03390 sinfo_free_table(&z);
03391 sinfo_free_table(&z_diff);
03392 sinfo_free_table(&tmp_sky);
03393 sinfo_free_table(&z_pos);
03394 sinfo_unwrap_vector(&vw);
03395 sinfo_unwrap_vector(&vs);
03396 sinfo_unwrap_vector(&vo);
03397 sinfoni_free_vector(&sx);
03398 sinfoni_free_vector(&sy);
03399 sinfo_unwrap_vector(&vx);
03400 sinfo_unwrap_vector(&vy);
03401 sinfo_free_table(&w_tbl);
03402 sinfo_free_table(&s_tbl);
03403 sinfo_free_table(&o_tbl);
03404 sinfo_free_polynomial(&cfit);
03405
03406 return 0;
03407
03408
03409 }
03419 static int
03420 sinfo_table_flag_nan(cpl_table** t,const char* label)
03421 {
03422
03423 int sz=0;
03424 int i=0;
03425 double* pt=NULL;
03426
03427 check_nomsg(sz=cpl_table_get_nrow(*t));
03428 check_nomsg(pt=cpl_table_get_data_double(*t,label));
03429 for(i=0;i<sz;i++) {
03430 if(isnan(pt[i])) {
03431 check_nomsg(cpl_table_set_invalid(*t ,label,i));
03432 }
03433 }
03434
03435 return 0;
03436
03437 cleanup:
03438
03439 return -1;
03440 }
03441
03451 static int
03452 sinfo_table_sky_obj_flag_nan(cpl_table** s,cpl_table** o, cpl_table** w)
03453 {
03454
03455 int no=0;
03456 int ns=0;
03457 int nw=0;
03458 int ni=0;
03459
03460 int i=0;
03461 double* po=NULL;
03462 double* ps=NULL;
03463 double* pw=NULL;
03464
03465 check_nomsg(no=cpl_table_get_nrow(*o));
03466 check_nomsg(ns=cpl_table_get_nrow(*s));
03467 check_nomsg(nw=cpl_table_get_nrow(*w));
03468
03469 if(no != ns || ns != nw || no != nw) {
03470 sinfo_msg_error("different input tables sizes");
03471 goto cleanup;
03472 }
03473 check_nomsg(po=cpl_table_get_data_double(*o,"INT"));
03474 check_nomsg(ps=cpl_table_get_data_double(*s,"INT"));
03475 check_nomsg(pw=cpl_table_get_data_double(*w,"WAVE"));
03476 for(i=0;i<no;i++) {
03477 if(isnan(po[i]) || isnan(ps[i]) || isnan(pw[i])) {
03478 check_nomsg(cpl_table_set_invalid(*o ,"INT",i));
03479 check_nomsg(cpl_table_set_invalid(*s ,"INT",i));
03480 check_nomsg(cpl_table_set_invalid(*w ,"WAVE",i));
03481 ni++;
03482 }
03483 }
03484
03485 return no-ni;
03486
03487 cleanup:
03488
03489 return -1;
03490 }
03491
03492
03493
03494
03495
03496
03497
03498
03499
03500
03501
03502
03503
03504
03505
03506
03507
03508
03509
03510
03511
03512
03513
03514
03515
03516
03517
03518
03519
03520
03521
03522
03523
03524
03525
03526
03527
03528
03529
03530
03531
03558 void
03559 sinfo_optimise_sky_sub(const double wtol,
03560 const double line_hw,
03561 const int method,
03562 const int do_rot,
03563 cpl_table* lrange,
03564 cpl_table* lambda,
03565 cpl_table* lr41,
03566 cpl_table* lr52,
03567 cpl_table* lr63,
03568 cpl_table* lr74,
03569 cpl_table* lr02,
03570 cpl_table* lr85,
03571 cpl_table* lr20,
03572 cpl_table* lr31,
03573 cpl_table* lr42,
03574 cpl_table* lr53,
03575 cpl_table* lr64,
03576 cpl_table* lr75,
03577 cpl_table* lr86,
03578 cpl_table* lr97,
03579 cpl_table* lr00,
03580 cpl_table** int_obj,
03581 cpl_table** int_sky,
03582 cpl_table** rscale)
03583
03584 {
03585
03586 int npixw=2*line_hw;
03587 cpl_array* do_hk=NULL;
03588 cpl_array* rfit=NULL;
03589 int i=0;
03590 cpl_table* sky_lr=NULL;
03591 cpl_table* obj_lr=NULL;
03592 cpl_table* wav_lr=NULL;
03593 double sky_med=0;
03594 double sky_sdv=0;
03595 int lr41_i=0;
03596 int lr52_i=0;
03597 int lr63_i=0;
03598 int lr74_i=0;
03599 int lr02_i=0;
03600 int lr85_i=0;
03601 int lr20_i=0;
03602 int lr31_i=0;
03603 int lr42_i=0;
03604 int lr53_i=0;
03605 int lr64_i=0;
03606 int lr75_i=0;
03607 int lr86_i=0;
03608 int lr97_i=0;
03609 int lr00_i=0;
03610
03611 int xxx1_i=0;
03612 int status=0;
03613 int finite_pix_i=0;
03614 double sky_thresh=0.;
03615
03616 cpl_table* rat_sky=NULL;
03617
03618 cpl_table* xxx1=NULL;
03619 cpl_table* xxx2=NULL;
03620 cpl_table* xxx1_sub=NULL;
03621 cpl_table* line_regions=NULL;
03622 cpl_table* cont_regions=NULL;
03623 int line_i=0;
03624 int cont_i=0;
03625 double fmed=0;
03626 double fsdv=0;
03627 cpl_table* fline_res=NULL;
03628 int fclip_i=0;
03629 int fline_i=0;
03630 cpl_table* rscale0=NULL;
03631 double r=0;
03632 cpl_table* obj_cont=NULL;
03633 cpl_table* sky_cont=NULL;
03634 cpl_table* obj_line=NULL;
03635 cpl_table* sky_line=NULL;
03636
03637
03638
03639 int low_pos_i=0;
03640 int med_pos_i=0;
03641 int hi_pos_i=0;
03642
03643 cpl_table* finite_pix=NULL;
03644 cpl_table* tmp_tbl=NULL;
03645
03646 cpl_table* low_scale=NULL;
03647 cpl_table* med_scale=NULL;
03648 cpl_table* hi_regions=NULL;
03649
03650 cpl_table* low_regions=NULL;
03651 cpl_table* med_regions=NULL;
03652
03653
03654 cpl_table* low_pos=NULL;
03655 cpl_table* med_pos=NULL;
03656 cpl_table* hi_pos=NULL;
03657 cpl_table* llr_xxx1=NULL;
03658
03659 double rhi=0;
03660 double rmed=0;
03661 double rlow=0;
03662
03663 double min_lrange=0;
03664 double max_lrange=0;
03665
03666 int nrow=0;
03667
03668
03669 double w_rot_low[NROT]={1.00852,1.03757,1.09264,1.15388,1.22293,
03670 1.30216,1.45190,1.52410,1.60308,1.69037,
03671 1.78803,2.02758,2.18023,1.02895,1.08343,
03672 1.14399,1.21226,1.29057,1.43444,1.50555,
03673 1.58333,1.66924,1.76532,2.00082,2.15073};
03674
03675
03676 double w_rot_med[NROT]={1.00282,1.02139,1.04212,1.07539,1.09753,
03677 1.13542,1.15917,1.20309,1.22870,1.28070,
03678 1.30853,1.41861,1.46048,1.48877,1.53324,
03679 1.56550,1.61286,1.65024,1.70088,1.74500,
03680 1.79940,1.97719,2.04127,2.12496,2.19956};
03681
03682
03683
03684 check_nomsg(do_hk = cpl_array_new(NBOUND+1,CPL_TYPE_INT));
03685 check_nomsg(rfit = cpl_array_new(NBOUND+1,CPL_TYPE_DOUBLE));
03686
03687 lr41_i=cpl_table_get_nrow(lr41);
03688 lr52_i=cpl_table_get_nrow(lr52);
03689 lr63_i=cpl_table_get_nrow(lr63);
03690 lr74_i=cpl_table_get_nrow(lr74);
03691 lr02_i=cpl_table_get_nrow(lr02);
03692 lr85_i=cpl_table_get_nrow(lr85);
03693 lr20_i=cpl_table_get_nrow(lr20);
03694 lr31_i=cpl_table_get_nrow(lr31);
03695 lr42_i=cpl_table_get_nrow(lr42);
03696 lr53_i=cpl_table_get_nrow(lr53);
03697 lr64_i=cpl_table_get_nrow(lr64);
03698 lr75_i=cpl_table_get_nrow(lr75);
03699 lr86_i=cpl_table_get_nrow(lr86);
03700 lr97_i=cpl_table_get_nrow(lr97);
03701 check_nomsg(lr00_i=cpl_table_get_nrow(lr00));
03702
03703 cpl_array_set_int(do_hk,0,lr41_i);
03704 cpl_array_set_int(do_hk,1,lr52_i);
03705 cpl_array_set_int(do_hk,2,lr63_i);
03706 cpl_array_set_int(do_hk,3,lr74_i);
03707 cpl_array_set_int(do_hk,4,lr02_i);
03708 cpl_array_set_int(do_hk,5,lr85_i);
03709 cpl_array_set_int(do_hk,6,lr20_i);
03710 cpl_array_set_int(do_hk,7,lr31_i);
03711 cpl_array_set_int(do_hk,8,lr42_i);
03712 cpl_array_set_int(do_hk,9,lr53_i);
03713 cpl_array_set_int(do_hk,10,lr64_i);
03714 cpl_array_set_int(do_hk,11,lr75_i);
03715 cpl_array_set_int(do_hk,12,lr86_i);
03716 cpl_array_set_int(do_hk,13,lr97_i);
03717 check_nomsg(cpl_array_set_int(do_hk,14,lr00_i));
03718
03719 check_nomsg(rscale0=cpl_table_duplicate(*int_sky));
03720 check_nomsg(cpl_table_new_column(rscale0,"RATIO",CPL_TYPE_DOUBLE));
03721
03722
03723 for (i=0;i<NBOUND+1;i++) {
03724 if (cpl_array_get_int(do_hk,i,&status) > 0) {
03725
03726
03727 switch(i) {
03728
03729
03730 case 0:
03731 ck0_nomsg(sinfo_get_obj_sky_wav_sub(*int_obj,*int_sky,lambda,lr41,wtol,
03732 &obj_lr,&sky_lr,&wav_lr));
03733 break;
03734
03735 case 1:
03736 ck0_nomsg(sinfo_get_obj_sky_wav_sub(*int_obj,*int_sky,lambda,lr52,wtol,
03737 &obj_lr,&sky_lr,&wav_lr));
03738 break;
03739
03740 case 2:
03741 ck0_nomsg(sinfo_get_obj_sky_wav_sub(*int_obj,*int_sky,lambda,lr63,wtol,
03742 &obj_lr,&sky_lr,&wav_lr));
03743 break;
03744
03745 case 3:
03746 ck0_nomsg(sinfo_get_obj_sky_wav_sub(*int_obj,*int_sky,lambda,lr74,wtol,
03747 &obj_lr,&sky_lr,&wav_lr));
03748 break;
03749
03750 case 4:
03751 ck0_nomsg(sinfo_get_obj_sky_wav_sub(*int_obj,*int_sky,lambda,lr02,wtol,
03752 &obj_lr,&sky_lr,&wav_lr));
03753 break;
03754
03755 case 5:
03756 ck0_nomsg(sinfo_get_obj_sky_wav_sub(*int_obj,*int_sky,lambda,lr85,wtol,
03757 &obj_lr,&sky_lr,&wav_lr));
03758 break;
03759 case 6:
03760 ck0_nomsg(sinfo_get_obj_sky_wav_sub(*int_obj,*int_sky,lambda,lr20,wtol,
03761 &obj_lr,&sky_lr,&wav_lr));
03762 break;
03763 case 7:
03764 ck0_nomsg(sinfo_get_obj_sky_wav_sub(*int_obj,*int_sky,lambda,lr31,wtol,
03765 &obj_lr,&sky_lr,&wav_lr));
03766 break;
03767 case 8:
03768 ck0_nomsg(sinfo_get_obj_sky_wav_sub(*int_obj,*int_sky,lambda,lr42,wtol,
03769 &obj_lr,&sky_lr,&wav_lr));
03770 break;
03771 case 9:
03772 ck0_nomsg(sinfo_get_obj_sky_wav_sub(*int_obj,*int_sky,lambda,lr53,wtol,
03773 &obj_lr,&sky_lr,&wav_lr));
03774 break;
03775 case 10:
03776 ck0_nomsg(sinfo_get_obj_sky_wav_sub(*int_obj,*int_sky,lambda,lr64,wtol,
03777 &obj_lr,&sky_lr,&wav_lr));
03778 break;
03779 case 11:
03780 ck0_nomsg(sinfo_get_obj_sky_wav_sub(*int_obj,*int_sky,lambda,lr75,wtol,
03781 &obj_lr,&sky_lr,&wav_lr));
03782 break;
03783 case 12:
03784 ck0_nomsg(sinfo_get_obj_sky_wav_sub(*int_obj,*int_sky,lambda,lr86,wtol,
03785 &obj_lr,&sky_lr,&wav_lr));
03786 break;
03787 case 13:
03788 ck0_nomsg(sinfo_get_obj_sky_wav_sub(*int_obj,*int_sky,lambda,lr97,wtol,
03789 &obj_lr,&sky_lr,&wav_lr));
03790 break;
03791 case 14:
03792 ck0_nomsg(sinfo_get_obj_sky_wav_sub(*int_obj,*int_sky,lambda,lr00,
03793 wtol,&obj_lr,&sky_lr,&wav_lr));
03794 break;
03795 default:
03796 sinfo_msg_error("case not supported");
03797 goto cleanup;
03798 }
03799 if(sky_lr == NULL || obj_lr == NULL || wav_lr == NULL) {
03800 finite_pix_i=0;
03801 sinfo_msg("no good pix left");
03802 } else {
03803
03804 check_nomsg(cpl_table_save(sky_lr,NULL,NULL,"out_skylr0.fits",0));
03805 check_nomsg(cpl_table_save(obj_lr,NULL,NULL,"out_objlr0.fits",0));
03806
03807
03808 check_nomsg(finite_pix_i=sinfo_table_sky_obj_flag_nan(&sky_lr,
03809 &obj_lr,
03810 &wav_lr));
03811
03812
03813 }
03814
03815
03816 if (finite_pix_i > npixw) {
03817
03818
03819 check_nomsg(cpl_table_erase_invalid(obj_lr));
03820 check_nomsg(cpl_table_erase_invalid(sky_lr));
03821 check_nomsg(cpl_table_erase_invalid(wav_lr));
03822
03823
03824 check_nomsg(sky_med=cpl_table_get_column_median(sky_lr,"INT"));
03825 check_nomsg(sky_sdv=cpl_table_get_column_stdev(sky_lr,"INT"));
03826 check_nomsg(cpl_table_select_all(sky_lr));
03827 sky_thresh=sky_med+sky_sdv;
03828
03829 check_nomsg(xxx1_i=cpl_table_and_selected_double(sky_lr,"INT",
03830 CPL_GREATER_THAN,sky_thresh));
03831 check_nomsg(xxx1 = cpl_table_extract_selected(sky_lr));
03832 check_nomsg(cpl_table_select_all(sky_lr));
03833
03834 if (xxx1_i > 0) {
03835
03836
03837
03838 check_nomsg(xxx2 = cpl_table_duplicate(sky_lr));
03839
03840
03841
03842 ck0_nomsg(sinfo_table_threshold(&xxx2,"INT",sky_thresh,
03843 sky_thresh,0.,10.));
03844
03845
03846
03847
03848
03849
03850
03851 check_nomsg(sinfo_convolve_kernel(&xxx2,npixw/2));
03852
03853
03854
03855 check_nomsg(line_i=cpl_table_and_selected_double(xxx2,"CNV",
03856 CPL_GREATER_THAN,0));
03857
03858 check_nomsg(line_regions=cpl_table_extract_selected(xxx2));
03859 check_nomsg(cpl_table_erase_column(line_regions,"INT"));
03860 check_nomsg(cpl_table_erase_column(line_regions,"CNV"));
03861
03862 check_nomsg(cpl_table_select_all(xxx2));
03863
03864
03865 check_nomsg(cont_i=cpl_table_and_selected_double(xxx2,"CNV",
03866 CPL_EQUAL_TO,0));
03867 check_nomsg(cont_regions=cpl_table_extract_selected(xxx2));
03868 check_nomsg(cpl_table_erase_column(cont_regions,"INT"));
03869 check_nomsg(cpl_table_erase_column(cont_regions,"CNV"));
03870 check_nomsg(cpl_table_select_all(xxx2));
03871 sinfo_free_table(&xxx2);
03872
03873
03874 if (line_i >= 3 && cont_i >= 3) {
03875
03876
03877 if (i == 0) sinfo_msg("optimising 4-1 transitions");
03878 if (i == 1) sinfo_msg("optimising 5-2 transitions");
03879 if (i == 2) sinfo_msg("optimising 6-3 transitions");
03880 if (i == 3) sinfo_msg("optimising 7-4 transitions");
03881 if (i == 4) sinfo_msg("optimising 0-2 transitions");
03882 if (i == 5) sinfo_msg("optimising 8-5 transitions");
03883 if (i == 6) sinfo_msg("optimising 2-0 transitions");
03884 if (i == 7) sinfo_msg("optimising 3-1 transitions");
03885 if (i == 8) sinfo_msg("optimising 4-2 transitions");
03886 if (i == 9) sinfo_msg("optimising 5-3 transitions");
03887 if (i == 10) sinfo_msg("optimising 6-4 transitions");
03888 if (i == 11) sinfo_msg("optimising 7-5 transitions");
03889 if (i == 12) sinfo_msg("optimising 8-6 transitions");
03890 if (i == 13) sinfo_msg("optimising 9-7 transitions");
03891 if (i == 14) sinfo_msg("optimising final bit");
03892
03893
03894
03895
03896
03897 sinfo_free_table(&obj_cont);
03898 sinfo_free_table(&sky_cont);
03899 sinfo_free_table(&sky_line);
03900 sinfo_free_table(&obj_line);
03901
03902 cknull_nomsg(obj_line=sinfo_table_select_range(obj_lr,line_regions,wtol));
03903 cknull_nomsg(sky_line=sinfo_table_select_range(sky_lr,line_regions,wtol));
03904 cknull_nomsg(obj_cont=sinfo_table_select_range(obj_lr,cont_regions,wtol));
03905 cknull_nomsg(sky_cont=sinfo_table_select_range(sky_lr,cont_regions,wtol));
03906
03907
03908
03909
03910
03911
03912
03913
03914
03915
03916
03917
03918
03919
03920
03921
03922
03923
03924
03925
03926
03927 sinfo_free_table(&fline_res);
03928
03929
03930 ck0_nomsg(sinfo_get_line_ratio(obj_line,obj_cont,
03931 sky_line,sky_cont,method,&r));
03932 sinfo_msg("1st Line ratio %g",r);
03933
03934
03935 if(cpl_table_get_nrow(obj_cont) > 0) {
03936 check_nomsg(fline_res=sinfo_table_interpol(obj_line,obj_cont,
03937 sky_line,sky_cont,
03938 r));
03939 } else {
03940 check_nomsg(fline_res=cpl_table_duplicate(obj_line));
03941 }
03942
03943
03944 cpl_table_select_all(fline_res);
03945 check_nomsg(fmed = cpl_table_get_column_median(fline_res,"INT"));
03946 check_nomsg(fsdv = cpl_table_get_column_stdev(fline_res,"INT"));
03947
03948 check_nomsg(cpl_table_duplicate_column(fline_res,"AINT",
03949 fline_res,"INT"));
03950 check_nomsg(cpl_table_multiply_columns(fline_res,"AINT","INT"));
03951 check_nomsg(cpl_table_power_column(fline_res,"AINT",0.5));
03952 check_nomsg(fclip_i=cpl_table_and_selected_double(fline_res,"AINT",
03953 CPL_GREATER_THAN,
03954 fmed+3*fsdv));
03955
03956 check_nomsg(cpl_table_select_all(fline_res));
03957
03958
03959 if (fclip_i > 0) {
03960
03961
03962
03963 check_nomsg(line_i=cpl_table_and_selected_double(fline_res,
03964 "AINT",
03965 CPL_NOT_GREATER_THAN,
03966 fmed+3*fsdv));
03967
03968 sinfo_free_table(&line_regions);
03969
03970 check_nomsg(line_regions=cpl_table_extract_selected(fline_res));
03971 check_nomsg(cpl_table_erase_column(line_regions,"INT"));
03972 check_nomsg(cpl_table_erase_column(line_regions,"AINT"));
03973
03974 sinfo_free_table(&obj_line);
03975 sinfo_free_table(&sky_line);
03976
03977
03978
03979
03980
03981
03982
03983
03984
03985 obj_line=sinfo_table_select_range(obj_lr,line_regions,wtol);
03986 sky_line=sinfo_table_select_range(sky_lr,line_regions,wtol);
03987 fline_i=cpl_table_get_nrow(line_regions);
03988
03989
03990 if(fline_i>=3) {
03991
03992
03993 ck0_nomsg(sinfo_get_line_ratio(obj_line,obj_cont,
03994 sky_line,sky_cont,method,&r));
03995
03996 sinfo_msg("2nd Line ratio %g",r);
03997
03998 } else {
03999 irplib_error_reset();
04000 cpl_error_reset();
04001 }
04002
04003 sinfo_free_table(&sky_line);
04004 sinfo_free_table(&obj_line);
04005 }
04006
04007 sinfo_msg("use %d pixels for line and %d for continuum estimation",
04008 cpl_table_get_nrow(line_regions),cpl_table_get_nrow(cont_regions));
04009
04010 sinfo_msg("OH spectrum scaling = %f ",r);
04011 check_nomsg(cpl_array_set_double(rfit,i,r));
04012 ck0_nomsg(sinfo_table_set(&rscale0,wav_lr,r,wtol));
04013
04014 }
04015 }
04016 }
04017
04018 }
04019
04020 sinfo_free_table(&xxx1);
04021 sinfo_free_table(&xxx2);
04022 sinfo_free_table(&sky_lr);
04023 sinfo_free_table(&obj_lr);
04024 sinfo_free_table(&wav_lr);
04025
04026 sinfo_free_table(&line_regions);
04027 sinfo_free_table(&cont_regions);
04028
04029 }
04030
04031 sinfo_free_array(&do_hk);
04032 sinfo_free_array(&rfit);
04033
04034
04035
04036
04037 check_nomsg(cpl_table_select_all(rscale0));
04038
04039
04040
04041
04042 check_nomsg(*rscale = cpl_table_extract_selected(rscale0));
04043 sinfo_free_table(&rscale0);
04044
04045
04046 check_nomsg(rat_sky=cpl_table_duplicate(*int_sky));
04047 check_nomsg(cpl_table_duplicate_column(rat_sky,"RATIO",*rscale,"RATIO"));
04048 check_nomsg(cpl_table_multiply_columns(rat_sky,"INT","RATIO"));
04049
04050
04051
04052
04053
04054
04055
04056
04057
04058 if (do_rot == 1) {
04059
04060
04061 check_nomsg(min_lrange=cpl_table_get_column_min(lrange,"WAVE"));
04062 check_nomsg(max_lrange=cpl_table_get_column_max(lrange,"WAVE"));
04063
04064 check_nomsg(finite_pix_i=sinfo_table_sky_obj_flag_nan(&rat_sky,
04065 int_obj,
04066 &lambda));
04067
04068
04069 check_nomsg(finite_pix=cpl_table_duplicate(lambda));
04070
04071
04072 check_nomsg(cpl_table_erase_invalid(finite_pix));
04073
04074
04075 if (finite_pix_i > npixw) {
04076
04077
04078
04079
04080 check_nomsg(cpl_table_and_selected_double(finite_pix,"WAVE",
04081 CPL_GREATER_THAN,
04082 min_lrange));
04083
04084 check_nomsg(cpl_table_and_selected_double(finite_pix,"WAVE",
04085 CPL_LESS_THAN,
04086 max_lrange));
04087
04088
04089
04090 check_nomsg(tmp_tbl=cpl_table_extract_selected(finite_pix));
04091 sinfo_free_table(&finite_pix);
04092 check_nomsg(finite_pix=cpl_table_duplicate(tmp_tbl));
04093 sinfo_free_table(&tmp_tbl);
04094 sinfo_free_table(&sky_lr);
04095 sinfo_free_table(&obj_lr);
04096 sinfo_free_table(&wav_lr);
04097
04098
04099 cknull(sky_lr=sinfo_table_select_range(rat_sky,finite_pix,wtol),
04100 "extracting sky sub range");
04101 cknull(obj_lr=sinfo_table_select_range(*int_obj,finite_pix,wtol),
04102 "extracting obj sub range");
04103 cknull(wav_lr=sinfo_table_select_range(lambda,finite_pix,wtol),
04104 "extracting sky sub range");
04105
04106
04107
04108
04109
04110
04111
04112
04113
04114 if(1 == cpl_table_has_valid(sky_lr,"INT")) {
04115 check_nomsg(sky_med=cpl_table_get_column_median(sky_lr,"INT"));
04116 check_nomsg(sky_sdv=cpl_table_get_column_stdev(sky_lr,"INT"));
04117 sky_thresh=sky_med+sky_sdv;
04118
04119 check_nomsg(xxx1_i=cpl_table_and_selected_double(sky_lr,"INT",
04120 CPL_GREATER_THAN,sky_thresh));
04121 check_nomsg(xxx1=cpl_table_extract_selected(sky_lr));
04122 check_nomsg(cpl_table_select_all(sky_lr));
04123 } else {
04124 xxx1_i=0;
04125 }
04126 if (xxx1_i > 0) {
04127 sinfo_msg("xxx1_i=%d",xxx1_i);
04128
04129 sinfo_msg("wav_lr wmin=%g wmax=%g",
04130 cpl_table_get_column_min(wav_lr,"WAVE"),
04131 cpl_table_get_column_max(wav_lr,"WAVE"));
04132
04133 cknull_nomsg(llr_xxx1=sinfo_table_select_range(wav_lr,xxx1,wtol));
04134
04135
04136 cknull(low_pos=sinfo_find_rot_waves(w_rot_low,npixw,wtol,llr_xxx1),
04137 "Determining low positions");
04138
04139
04140 check_nomsg(low_pos_i=cpl_table_get_nrow(low_pos));
04141
04142 cknull(med_pos=sinfo_find_rot_waves(w_rot_med,npixw,wtol,llr_xxx1),
04143 "Determining med positions");
04144 check_nomsg(med_pos_i=cpl_table_get_nrow(med_pos));
04145
04146
04147
04148
04149
04150
04151
04152
04153
04154
04155
04156
04157
04158
04159
04160 cknull(hi_pos=sinfo_table_extract_rest(xxx1,low_pos,med_pos,wtol),
04161 "determining hi position");
04162 check_nomsg(hi_pos_i=cpl_table_get_nrow(hi_pos));
04163
04164
04165
04166
04167 check_nomsg(xxx2 = cpl_table_duplicate(sky_lr));
04168 check_nomsg(nrow=cpl_table_get_nrow(sky_lr));
04169
04170
04171
04172
04173
04174
04175
04176
04177 ck0_nomsg(sinfo_table_threshold(&xxx2,"INT",sky_thresh,
04178 sky_thresh,0.,10.));
04179 sinfo_msg("sky_thresh=%g %g %f",sky_thresh,sky_med,sky_sdv);
04180
04181 check_nomsg(sinfo_convolve_kernel(&xxx2,npixw/2));
04182
04183
04184 check_nomsg(cont_i=cpl_table_and_selected_double(xxx2,"CNV",
04185 CPL_EQUAL_TO,0));
04186 check_nomsg(cont_regions=cpl_table_extract_selected(xxx2));
04187
04188 sinfo_free_table(&xxx2);
04189 check_nomsg(cpl_table_erase_column(cont_regions,"INT"));
04190 check_nomsg(cpl_table_erase_column(cont_regions,"CNV"));
04191
04192 check(low_pos_i=sinfo_get_sub_regions(sky_lr,xxx1,low_pos,wtol,
04193 npixw,&low_regions),"failed determining low regions");
04194
04195 check(med_pos_i=sinfo_get_sub_regions(sky_lr,xxx1,med_pos,wtol,
04196 npixw,&med_regions),"failed determining med regions");
04197
04198
04199 check(hi_pos_i=sinfo_get_sub_regions(sky_lr,xxx1,hi_pos,wtol,
04200 npixw,&hi_regions),"failed determining hi regions");
04201
04202
04203
04204
04205
04206
04207
04208
04209
04210 sinfo_msg("hi_pos_i : %d med_pos_i : %d low_pos_i : %d cont_i: %d",
04211 hi_pos_i, med_pos_i, low_pos_i, cont_i);
04212
04213
04214 if (hi_pos_i >= 3 && med_pos_i >= 3 && low_pos_i >= 3 && cont_i >= 3) {
04215
04216
04217 ck0_nomsg(sinfo_compute_line_ratio(obj_lr, sky_lr,wtol, method,
04218 hi_regions,cont_regions,&rhi));
04219 sinfo_msg("high rotational OH scaling %g",rhi);
04220
04221
04222 ck0_nomsg(sinfo_compute_line_ratio(obj_lr, sky_lr,wtol, method,
04223 med_regions,cont_regions,&rmed));
04224
04225 sinfo_msg("P1(3.5) & R1(1.5) rotational OH scaling %g ",rmed);
04226
04227
04228 ck0_nomsg(sinfo_compute_line_ratio(obj_lr, sky_lr,wtol, method,
04229 low_regions,cont_regions,&rlow));
04230 sinfo_msg("P1(2.5) & Q1(1.5) rotational OH scaling %g",rlow);
04231
04232 cknull(low_scale=sinfo_find_rot_waves(w_rot_low,npixw,wtol,lambda),
04233 "Determining low scale");
04234
04235
04236
04237 cknull(med_scale=sinfo_find_rot_waves(w_rot_low,npixw,wtol,lambda),
04238 "Determining low scale");
04239 check_nomsg(cpl_table_multiply_scalar(*rscale,"RATIO",rhi));
04240 ck0_nomsg(sinfo_table_fill_column_over_range(rscale,med_scale,
04241 "RATIO",rmed/rhi,wtol));
04242 ck0_nomsg(sinfo_table_fill_column_over_range(rscale,low_scale,
04243 "RATIO",rlow/rhi,wtol));
04244
04245 }
04246 }
04247 }
04248 }
04249
04250
04251
04252 check_nomsg(cpl_table_duplicate_column(*int_sky,"INTC",*int_sky,"INT"));
04253
04254
04255
04256 check_nomsg(cpl_table_duplicate_column(*int_sky,"RATIO",*rscale,"RATIO"));
04257 check_nomsg(cpl_table_multiply_columns(*int_sky,"INTC","RATIO"));
04258
04259 check_nomsg(cpl_table_duplicate_column(*int_obj,"INTC",*int_obj,"INT"));
04260 check_nomsg(cpl_table_duplicate_column(*int_obj,"SKYC",*int_sky,"INTC"));
04261 check_nomsg(cpl_table_subtract_columns(*int_obj,"INTC","SKYC"));
04262
04263 check_nomsg(cpl_table_erase_column(*int_sky,"INT"));
04264 check_nomsg(cpl_table_name_column(*int_sky,"INTC","INT"));
04265
04266
04267 cleanup:
04268 sinfo_free_table(&llr_xxx1);
04269 sinfo_free_table(&hi_pos);
04270 sinfo_free_table(&low_pos);
04271 sinfo_free_table(&med_pos);
04272 sinfo_free_table(&low_regions);
04273 sinfo_free_table(&med_regions);
04274 sinfo_free_table(&hi_regions);
04275 sinfo_free_table(&low_scale);
04276 sinfo_free_table(&med_scale);
04277
04278
04279 sinfo_free_table(&finite_pix);
04280 sinfo_free_table(&xxx1_sub);
04281 sinfo_free_table(&tmp_tbl);
04282 sinfo_free_table(&rat_sky);
04283 sinfo_free_table(&fline_res);
04284 sinfo_free_table(&sky_cont);
04285 sinfo_free_table(&obj_cont);
04286 sinfo_free_table(&obj_line);
04287 sinfo_free_table(&sky_line);
04288 sinfo_free_table(&rscale0);
04289 sinfo_free_table(&xxx1);
04290 sinfo_free_table(&xxx2);
04291 sinfo_free_table(&line_regions);
04292 sinfo_free_table(&cont_regions);
04293 sinfo_free_table(&sky_lr);
04294 sinfo_free_table(&obj_lr);
04295 sinfo_free_table(&wav_lr);
04296 sinfo_free_array(&rfit);
04297 sinfo_free_array(&do_hk);
04298 return;
04299
04300 }
04309 static int
04310 sinfo_table_get_index_of_max(cpl_table* t,const char* name,cpl_type type)
04311 {
04312
04313 int i=0;
04314 int result=0;
04315 int nrow=0;
04316 int* pi=NULL;
04317 float* pf=NULL;
04318 double* pd=NULL;
04319 double max=0;
04320
04321
04322 if(t == NULL) {
04323 cpl_error_set(cpl_func, CPL_ERROR_NULL_INPUT);
04324 return result;
04325 }
04326 max=cpl_table_get_column_max(t,name);
04327 nrow=cpl_table_get_nrow(t);
04328 switch(type) {
04329
04330 case CPL_TYPE_INT:
04331 pi=cpl_table_get_data_int(t,name);
04332 for(i=0;i<nrow;i++) {
04333 if(pi[i]==(int)max) result=i;
04334 }
04335 break;
04336 case CPL_TYPE_FLOAT:
04337 pf=cpl_table_get_data_float(t,name);
04338 for(i=0;i<nrow;i++) {
04339 if(pf[i]==(float)max) result=i;
04340 }
04341 break;
04342 case CPL_TYPE_DOUBLE:
04343 pd=cpl_table_get_data_double(t,name);
04344 for(i=0;i<nrow;i++) {
04345 if(pd[i]==max) result=i;
04346 }
04347 break;
04348 default:
04349 sinfo_msg_error("Wrong column type");
04350 cpl_error_set(cpl_func, CPL_ERROR_TYPE_MISMATCH);
04351 return result;
04352
04353 }
04354 return result;
04355 }
04356
04357
04358
04368 static int
04369 sinfo_table_get_index_of_val(cpl_table* t,
04370 const char* name,
04371 double val,
04372 cpl_type type)
04373 {
04374
04375 int i=0;
04376 int result=0;
04377 int nrow=0;
04378 int* pi=NULL;
04379 float* pf=NULL;
04380 double* pd=NULL;
04381
04382 if(t == NULL) {
04383 cpl_error_set(cpl_func, CPL_ERROR_NULL_INPUT);
04384 return result;
04385 }
04386
04387 nrow=cpl_table_get_nrow(t);
04388 switch(type) {
04389
04390 case CPL_TYPE_INT:
04391 pi=cpl_table_get_data_int(t,name);
04392 for(i=0;i<nrow;i++) {
04393 if(pi[i]==(int)val) result=i;
04394 }
04395 break;
04396 case CPL_TYPE_FLOAT:
04397 pf=cpl_table_get_data_float(t,name);
04398 for(i=0;i<nrow;i++) {
04399 if(pf[i]==(float)val) result=i;
04400 }
04401 break;
04402 case CPL_TYPE_DOUBLE:
04403 pd=cpl_table_get_data_double(t,name);
04404 for(i=0;i<nrow;i++) {
04405 if(pd[i]==val) result=i;
04406 }
04407 break;
04408 default:
04409 sinfo_msg_error("Wrong column type");
04410 cpl_error_set(cpl_func, CPL_ERROR_TYPE_MISMATCH);
04411 return result;
04412
04413 }
04414 return result;
04415 }
04416
04428 static double
04429 sinfo_table_column_interpolate(const cpl_table* t,
04430 const char* name,
04431 const double x)
04432 {
04433
04434 double val1=0;
04435 double val2=0;
04436 int x1=0;
04437 int x2=0;
04438 double m=0;
04439 double y=0;
04440 int status=0;
04441 int nrow=0;
04442 nrow=cpl_table_get_nrow(t);
04443 if ((1<x) && (x<nrow-1)) {
04444 x1=x-1;
04445 x2=x+1;
04446 } else if (x<2) {
04447 x1=0;
04448 x2=1;
04449 } else {
04450 x1=nrow-2;
04451 x2=nrow-1;
04452 }
04453 check_nomsg(val1=cpl_table_get(t,name,x1,&status));
04454 check_nomsg(val2=cpl_table_get(t,name,x2,&status));
04455
04456 m=(val2-val1)/(x2-x1);
04457 y=val1+m*(x-x1);
04458
04459 return y;
04460
04461 cleanup:
04462
04463 return -1;
04464
04465
04466 }
04467
04476 static cpl_imagelist*
04477 sinfo_imagelist_select_range(const cpl_imagelist* inp,
04478 const cpl_table* full,
04479 const cpl_table* good,
04480 const double tol)
04481 {
04482 cpl_imagelist* out=NULL;
04483 int osz=0;
04484 int isz=0;
04485 int ksz=0;
04486 int k=0;
04487 int i=0;
04488 int status=0;
04489
04490 double wave_chk=0;
04491 double wave_sel=0;
04492
04493 cpl_image* img=NULL;
04494
04495
04496
04497
04498
04499 check_nomsg(osz=cpl_table_get_nrow(good));
04500 check_nomsg(ksz=cpl_imagelist_get_size(inp));
04501 check_nomsg(isz=cpl_table_get_nrow(good));
04502 check_nomsg(out=cpl_imagelist_new());
04503
04504
04505 for(k=0;k<ksz;k++) {
04506 check_nomsg(img=cpl_imagelist_get(inp,k));
04507 check_nomsg(wave_chk=cpl_table_get(full,"WAVE",k,&status));
04508 if(i<isz) {
04509 check_nomsg(wave_sel=cpl_table_get(good,"WAVE",i,&status));
04510 }
04511
04512 if(fabs(wave_chk - wave_sel) < tol) {
04513 check_nomsg(cpl_imagelist_set(out,cpl_image_duplicate(img),i));
04514 i++;
04515 }
04516 }
04517 if(i==0) {
04518 sinfo_msg_error("No lines selected");
04519 goto cleanup;
04520 }
04521 return out;
04522
04523 cleanup:
04524
04525 return NULL;
04526
04527 }
04528
04538 static int
04539 sinfo_table_extract_finite(const cpl_table* in1,
04540 const cpl_table* in2,
04541 cpl_table** ou1,
04542 cpl_table** ou2)
04543 {
04544
04545 int size1=0;
04546 int size2=0;
04547 int i=0;
04548 int ninv1=0;
04549 int ninv2=0;
04550 double* pout1=NULL;
04551 double* pout2=NULL;
04552
04553 cknull(in1,"null input image");
04554 cknull(in2,"null input image");
04555 cknull_nomsg(*ou1=cpl_table_duplicate(in1));
04556 cknull_nomsg(*ou2=cpl_table_duplicate(in2));
04557
04558 check_nomsg(size1=cpl_table_get_nrow(*ou1));
04559 check_nomsg(size2=cpl_table_get_nrow(*ou2));
04560
04561 check_nomsg(pout1=cpl_table_get_data_double(*ou1,"VALUE"));
04562 check_nomsg(pout2=cpl_table_get_data_double(*ou2,"VALUE"));
04563 for(i=0;i<size1;i++) {
04564 if (isnan(pout1[i])) {
04565 check_nomsg(cpl_table_set_invalid(*ou1,"VALUE",i));
04566 check_nomsg(cpl_table_set_invalid(*ou2,"VALUE",i));
04567 }
04568 }
04569 ninv1=cpl_table_count_invalid(*ou1,"VALUE");
04570 ninv2=cpl_table_count_invalid(*ou2,"VALUE");
04571 if(ninv1==size1) {
04572 goto cleanup;
04573 }
04574 if(ninv2==size2) {
04575 goto cleanup;
04576 }
04577 check_nomsg(cpl_table_erase_invalid(*ou1));
04578 check_nomsg(cpl_table_erase_invalid(*ou2));
04579 return (size1-ninv1);
04580
04581 cleanup:
04582 return 0;
04583
04584 }
04585
04592 static cpl_table*
04593 sinfo_image2table(const cpl_image* im)
04594 {
04595 cpl_table* out=NULL;
04596 int sx=0;
04597 int sy=0;
04598 double* pim=NULL;
04599 double* pval=NULL;
04600 int i=0;
04601 int j=0;
04602 int k=0;
04603
04604 cknull(im,"input image is NULL");
04605
04606 check_nomsg(sx=cpl_image_get_size_x(im));
04607 check_nomsg(sy=cpl_image_get_size_y(im));
04608 check_nomsg(pim=cpl_image_get_data_double(im));
04609 check_nomsg(out=cpl_table_new(sx*sy));
04610 check_nomsg(cpl_table_new_column(out,"VALUE",CPL_TYPE_DOUBLE));
04611 check_nomsg(pval=cpl_table_get_data_double(out,"VALUE"));
04612
04613 for(j=0;j<sy;j++) {
04614 for(i=0;i<sx;i++) {
04615
04616
04617
04618
04619 cpl_table_set_double(out,"VALUE",k++,pim[j*sx+i]);
04620 }
04621 }
04622
04623 return out;
04624 cleanup:
04625 sinfo_free_table(&out);
04626 return NULL;
04627
04628 }
04637 int
04638 sinfo_check_screw_values(cpl_table** int_obj,
04639 cpl_table** int_sky,
04640 cpl_table* grange,
04641 const double wtol)
04642 {
04643
04644 cpl_table* xsky=NULL;
04645 cpl_table* xobj=NULL;
04646 cpl_table* gsky_pos=NULL;
04647 cpl_table* gobj_pos=NULL;
04648 cpl_table* tbl1=NULL;
04649
04650 double sky_min=0;
04651 double sky_max=0;
04652 double gsky_min=0;
04653 double gsky_max=0;
04654 double obj_min=0;
04655 double obj_max=0;
04656 double gobj_min=0;
04657 double gobj_max=0;
04658
04659 int gsky_pos_i=0;
04660 int gobj_pos_i=0;
04661 int ns1=0;
04662 int ns2=0;
04663 int no1=0;
04664 int no2=0;
04665
04666 cknull(*int_sky,"Null input sky spectrum");
04667 cknull(*int_obj,"Null input obj spectrum");
04668 cknull(grange,"Null input wavelength range");
04669 cknull_nomsg(xsky=sinfo_table_select_range(*int_sky,grange,wtol));
04670 check_nomsg(sky_min=cpl_table_get_column_min(xsky,"INT"));
04671 check_nomsg(sky_max=cpl_table_get_column_max(xsky,"INT"));
04672
04673 gsky_max = (sky_max>0) ? sky_max : 0;
04674 gsky_min = (sky_min<0) ? sky_min : 0;
04675
04676 check_nomsg(cpl_table_select_all(*int_sky));
04677 check_nomsg(ns1=cpl_table_and_selected_double(*int_sky,"INT",
04678 CPL_GREATER_THAN,gsky_max));
04679
04680 ck0_nomsg(sinfo_table_set_column_invalid(int_sky,"INT"));
04681 check_nomsg(tbl1=cpl_table_extract_selected(*int_sky));
04682 check_nomsg(cpl_table_select_all(*int_sky));
04683 check_nomsg(ns2=cpl_table_and_selected_double(*int_sky,"INT",
04684 CPL_LESS_THAN,gsky_min));
04685
04686 ck0_nomsg(sinfo_table_set_column_invalid(int_sky,"INT"));
04687
04688 check_nomsg(gsky_pos=cpl_table_extract_selected(*int_sky));
04689 check_nomsg(cpl_table_insert(gsky_pos,tbl1,cpl_table_get_nrow(gsky_pos)));
04690 sinfo_free_table(&tbl1);
04691
04692 check_nomsg(gsky_pos_i=cpl_table_get_nrow(gsky_pos));
04693
04694
04695
04696 sinfo_free_table(&gsky_pos);
04697
04698 sinfo_free_table(&xsky);
04699 sinfo_msg("gsky_min=%f gsky_max=%f",gsky_min,gsky_max);
04700
04701 cknull_nomsg(xobj=sinfo_table_select_range(*int_obj,grange,wtol));
04702 check_nomsg(obj_min=cpl_table_get_column_min(xobj,"INT"));
04703 check_nomsg(obj_max=cpl_table_get_column_max(xobj,"INT"));
04704
04705 gobj_max = (obj_max>0) ? obj_max : 0;
04706 gobj_min = (obj_min<0) ? obj_min : 0;
04707
04708 check_nomsg(cpl_table_select_all(*int_obj));
04709 check_nomsg(no1=cpl_table_and_selected_double(*int_obj,"INT",
04710 CPL_GREATER_THAN,gobj_max));
04711
04712 ck0_nomsg(sinfo_table_set_column_invalid(int_obj,"INT"));
04713 check_nomsg(tbl1=cpl_table_extract_selected(*int_obj));
04714 check_nomsg(cpl_table_select_all(*int_obj));
04715 check_nomsg(no2=cpl_table_and_selected_double(*int_obj,"INT",
04716 CPL_LESS_THAN,gobj_min));
04717
04718 ck0_nomsg(sinfo_table_set_column_invalid(int_obj,"INT"));
04719
04720 check_nomsg(gobj_pos=cpl_table_extract_selected(*int_obj));
04721
04722 check_nomsg(cpl_table_insert(gobj_pos,tbl1,cpl_table_get_nrow(gobj_pos)));
04723
04724 sinfo_free_table(&tbl1);
04725 check_nomsg(gobj_pos_i=cpl_table_get_nrow(gobj_pos));
04726
04727
04728 sinfo_free_table(&gobj_pos);
04729 sinfo_msg("gobj_min=%f gobj_max=%f",gobj_min,gobj_max);
04730
04731 sinfo_free_table(&xobj);
04732
04733 return 0;
04734 cleanup:
04735 sinfo_free_table(&xsky);
04736 sinfo_free_table(&xobj);
04737 sinfo_free_table(&gsky_pos);
04738 sinfo_free_table(&gobj_pos);
04739 sinfo_free_table(&tbl1);
04740
04741 return -1;
04742
04743
04744 }
04745
04746
04747
04757 static int
04758 sinfo_table_fill_column_over_range(cpl_table** inp,
04759 const cpl_table* ref,
04760 const char* col,
04761 const double val,
04762 const double tol)
04763 {
04764
04765 int i=0;
04766 int k=0;
04767 int nref=0;
04768 int ninp=0;
04769
04770 double* pwin=NULL;
04771 double* pcin=NULL;
04772 double* pwrf=NULL;
04773
04774 cknull(inp,"null input table");
04775 cknull(ref,"null reference table");
04776
04777 check_nomsg(ninp=cpl_table_get_nrow(*inp));
04778 check_nomsg(nref=cpl_table_get_nrow(ref));
04779 check_nomsg(pwin=cpl_table_get_data_double(*inp,"WAVE"));
04780 check_nomsg(pcin=cpl_table_get_data_double(*inp,col));
04781 check_nomsg(pwrf=cpl_table_get_data_double(ref,"WAVE"));
04782
04783 k=0;
04784 i=0;
04785
04786
04787
04788
04789
04790 if(pwin[0]<=pwrf[0]) {
04791
04792 for(i=0;i<ninp;i++) {
04793 if(k<nref) {
04794
04795
04796
04797
04798 if(fabs(pwin[i] - pwrf[k])< tol) {
04799 pcin[i]=val;
04800 k++;
04801 }
04802 }
04803 }
04804 } else {
04805
04806
04807
04808 for(k=0;k<nref;k++) {
04809 if(i<ninp) {
04810
04811
04812
04813
04814 if(fabs(pwin[i] - pwrf[k])< tol) {
04815 pcin[i]=val;
04816 i++;
04817 }
04818 }
04819 }
04820 }
04821
04822 return 0;
04823
04824 cleanup:
04825 return -1;
04826
04827 }
04828
04829
04838 static cpl_table*
04839 sinfo_table_select_range(cpl_table* inp, cpl_table* ref,const double tol)
04840 {
04841
04842 cpl_table* out=NULL;
04843 int ninp=0;
04844 int nref=0;
04845 int nout=0;
04846
04847 int i=0;
04848 int k=0;
04849 double* pou;
04850 double* prf;
04851 double wmin=0;
04852 double wmax=0;
04853 cpl_table* tmp=NULL;
04854 cknull(inp,"null input table");
04855 cknull(ref,"null reference table");
04856
04857 check_nomsg(ninp=cpl_table_get_nrow(inp));
04858 check_nomsg(nref=cpl_table_get_nrow(ref));
04859 if(ninp != nref) {
04860
04861 check_nomsg(wmin=cpl_table_get_column_min(ref,"WAVE"));
04862 check_nomsg(wmax=cpl_table_get_column_max(ref,"WAVE"));
04863 cpl_table_select_all(inp);
04864 check_nomsg(ninp=cpl_table_and_selected_double(inp,"WAVE",
04865 CPL_NOT_LESS_THAN,wmin));
04866 check_nomsg(tmp=cpl_table_extract_selected(inp));
04867 check_nomsg(ninp=cpl_table_and_selected_double(tmp,"WAVE",
04868 CPL_NOT_GREATER_THAN,wmax));
04869 check_nomsg(out=cpl_table_extract_selected(tmp));
04870 sinfo_free_table(&tmp);
04871 } else {
04872 check_nomsg(out=cpl_table_duplicate(inp));
04873 }
04874
04875 check_nomsg(nout=cpl_table_get_nrow(out));
04876 if(nout == 0) {
04877 sinfo_msg("nout=%d",nout);
04878 goto cleanup;
04879 }
04880 tmp=cpl_table_duplicate(out);
04881 sinfo_free_table(&out);
04882 check_nomsg(pou=cpl_table_get_data_double(tmp,"WAVE"));
04883 check_nomsg(prf=cpl_table_get_data_double(ref,"WAVE"));
04884
04885 check_nomsg(cpl_table_new_column(tmp,"FLAG",CPL_TYPE_INT));
04886 check_nomsg(cpl_table_fill_column_window(tmp,"FLAG",0,nout,-1));
04887
04888 k=0;
04889 i=0;
04890
04891
04892
04893
04894
04895 if(pou[0]<=prf[0]) {
04896
04897 for(i=0;i<nout;i++) {
04898 if(k<nref) {
04899 if(fabs(pou[i] - prf[k])< tol) {
04900 check_nomsg(cpl_table_set_int(tmp,"FLAG",i,1));
04901 k++;
04902 }
04903 }
04904 }
04905 } else {
04906
04907
04908
04909 for(k=0;k<nref;k++) {
04910 if(i<nout) {
04911
04912
04913
04914
04915 if(fabs(pou[i] - prf[k])< tol) {
04916 check_nomsg(cpl_table_set_int(tmp,"FLAG",i,1));
04917 i++;
04918 }
04919 }
04920 }
04921 }
04922
04923 check_nomsg(nout=cpl_table_and_selected_int(tmp,"FLAG",CPL_GREATER_THAN,0));
04924 check_nomsg(out=cpl_table_extract_selected(tmp));
04925 sinfo_free_table(&tmp);
04926 check_nomsg(cpl_table_erase_column(out,"FLAG"));
04927 if(nout==0) {
04928 goto cleanup;
04929 }
04930
04931
04932 return out;
04933
04934 cleanup:
04935 sinfo_free_table(&tmp);
04936 sinfo_free_table(&out);
04937 return NULL;
04938
04939 }
04940
04950 cpl_table*
04951 sinfo_interpolate_sky(const cpl_table* inp,const cpl_table* lambdas)
04952 {
04953
04954 cpl_table* new=NULL;
04955 new = sinfo_interpolate(inp,lambdas,"WAVE","lsquadratic");;
04956
04957 return new;
04958 }
04959
04960
04970 static cpl_table*
04971 sinfo_interpolate(const cpl_table* inp,
04972 const cpl_table* lambdas,
04973 const char* name,
04974 const char* method)
04975 {
04976
04977 cpl_table* out=NULL;
04978
04979
04980 cknull_nomsg(lambdas);
04981 cknull_nomsg(name);
04982 cknull_nomsg(method);
04983
04984 out=cpl_table_duplicate(inp);
04985 return out;
04986
04987 cleanup:
04988
04989 return NULL;
04990
04991 }
04992
04993
04994
05007 static double
05008 sinfo_gaussian(double area,double sigma,double x,double x0,double off)
05009 {
05010 return area/sqrt(PI_NUMB*sigma*sigma)*
05011 exp(-(x-x0)*(x-x0)/(2*sigma*sigma))+off;
05012 }
05013
05020 int
05021 sinfo_table_smooth_column(cpl_table** t, const char* c, const int r)
05022 {
05023 int nrow=0;
05024 int i=0;
05025 int j=0;
05026 double* pval=NULL;
05027 double sum;
05028 check_nomsg(nrow=cpl_table_get_nrow(*t));
05029 check_nomsg(pval=cpl_table_get_data_double(*t,c));
05030 for(i=r;i<nrow;i++) {
05031 sum=0;
05032 for(j=-r;j<=r;j++) {
05033 sum+=pval[i+j];
05034 }
05035 pval[i]=sum/(2*r+1);
05036 }
05037 return 0;
05038 cleanup:
05039 return -1;
05040 }
05041
05050 void
05051 sinfo_shift_sky(cpl_frame** sky_frm,
05052 cpl_table** int_sky,
05053 const double zshift)
05054 {
05055
05056 int xsz=0;
05057 int ysz=0;
05058 int zsz=0;
05059 cpl_propertylist* plist=NULL;
05060 cpl_imagelist* sky_cub=NULL;
05061 cpl_imagelist* sky_shift=NULL;
05062 cpl_table* int_sky_dup=NULL;
05063
05064 double min=0;
05065 double max=0;
05066
05067
05068 cknull_nomsg(plist=cpl_propertylist_load(
05069 cpl_frame_get_filename(*sky_frm),0));
05070
05071 check_nomsg(xsz=sinfo_pfits_get_naxis1(plist));
05072 check_nomsg(ysz=sinfo_pfits_get_naxis2(plist));
05073 check_nomsg(zsz=sinfo_pfits_get_naxis3(plist));
05074 sinfo_free_propertylist(&plist);
05075
05076 cknull_nomsg(sky_cub=cpl_imagelist_load(cpl_frame_get_filename(*sky_frm),
05077 CPL_TYPE_FLOAT,0));
05078
05079 check_nomsg(min=cpl_table_get_column_min(*int_sky,"INT"));
05080 check_nomsg(max=cpl_table_get_column_max(*int_sky,"INT"));
05081 int_sky_dup=cpl_table_duplicate(*int_sky);
05082 sinfo_free_table(&(*int_sky));
05083
05084
05085
05086
05087
05088
05089
05090
05091
05092
05093
05094
05095
05096 check_nomsg(*int_sky=sinfo_table_shift_simple(int_sky_dup,"INT",zshift));
05097
05098
05099
05100
05101
05102
05103 sinfo_free_table(&int_sky_dup);
05104
05105
05106
05107
05108
05109
05110
05111
05112
05113
05114
05115
05116
05117
05118
05119 check_nomsg(sky_shift=sinfo_cube_zshift_simple(sky_cub,(float)zshift));
05120
05121
05122
05123
05124
05125 sinfo_free_imagelist(&sky_shift);
05126
05127
05128
05129
05130
05131
05132
05133
05134
05135
05136
05137
05138 sinfo_free_imagelist(&sky_cub);
05139
05140 return;
05141
05142 cleanup:
05143 sinfo_free_table(&int_sky_dup);
05144 sinfo_free_propertylist(&plist);
05145 sinfo_free_imagelist(&sky_shift);
05146 sinfo_free_imagelist(&sky_cub);
05147 return;
05148 }
05149
05150
05157 static int
05158 sinfo_convolve_kernel(cpl_table** t, const int rad)
05159 {
05160 int i=0;
05161 int j=0;
05162 int np=0;
05163 double val=0;
05164 double* pidata=NULL;
05165 double* pcdata=NULL;
05166 double dw=0;
05167 double dr=0;
05168 double ws=0;
05169 double we=0;
05170 cknull(*t,"null input table");
05171 check_nomsg(cpl_table_new_column(*t,"CNV",CPL_TYPE_DOUBLE));
05172 check_nomsg(pidata=cpl_table_get_data_double(*t,"INT"));
05173 check_nomsg(pcdata=cpl_table_get_data_double(*t,"CNV"));
05174 check_nomsg(ws=cpl_table_get_column_min(*t,"WAVE"));
05175 check_nomsg(we=cpl_table_get_column_max(*t,"WAVE"));
05176 check_nomsg(np=cpl_table_get_nrow(*t));
05177 dw=(we-ws)/(np-1);
05178 dr=(we-ws)/(rad-1);
05179
05180 for(i=0;i<rad;i++) {
05181 pcdata[i]=0.;
05182 }
05183 for(i=np-rad;i<np;i++) {
05184 pcdata[i]=0.;
05185 }
05186 for(i=rad;i<np-rad;i++) {
05187 val=0;
05188 for(j=-rad;j<rad;j++) {
05189 val+=pidata[i+j];
05190 }
05191
05192 check_nomsg(cpl_table_set_double(*t,"CNV",i,val));
05193 }
05194 return 0;
05195
05196 cleanup:
05197 return -1;
05198
05199 }
05200
05201
05202
05210 int
05211 sinfo_convolve_kernel2(cpl_table** t, const int rad)
05212 {
05213 int i=0;
05214 int j=0;
05215 int np=0;
05216 double val=0;
05217 double* pidata=NULL;
05218 double* pcdata=NULL;
05219 double dw=0;
05220 double dr=0;
05221 double ws=0;
05222 double we=0;
05223 cknull(*t,"null input table");
05224 check_nomsg(cpl_table_new_column(*t,"CNV",CPL_TYPE_DOUBLE));
05225 check_nomsg(pidata=cpl_table_get_data_double(*t,"INT"));
05226 check_nomsg(pcdata=cpl_table_get_data_double(*t,"CNV"));
05227 check_nomsg(ws=cpl_table_get_column_min(*t,"WAVE"));
05228 check_nomsg(we=cpl_table_get_column_max(*t,"WAVE"));
05229 check_nomsg(np=cpl_table_get_nrow(*t));
05230 dw=(we-ws)/(np-1);
05231 dr=(we-ws)/(rad-1);
05232
05233 for(i=0;i<rad;i++) {
05234 pcdata[i]=0.;
05235 }
05236 for(i=np-rad;i<np;i++) {
05237 pcdata[i]=0.;
05238 }
05239 for(i=0;i<np-rad;i++) {
05240 val=0;
05241 for(j=0;j<rad;j++) {
05242 val+=pidata[i+j];
05243 }
05244
05245 check_nomsg(cpl_table_set_double(*t,"CNV",i,val));
05246 }
05247 return 0;
05248
05249 cleanup:
05250 return -1;
05251
05252 }
05253
05254
05255
05263 int
05264 sinfo_convolve_exp(cpl_table** t, const int rad, const double fwhm)
05265 {
05266 int i=0;
05267 int j=0;
05268 int np=0;
05269 double ln2=0.693147180560;
05270 double k=ln2/fwhm;
05271 double val=0;
05272 double* pidata=NULL;
05273 double* pcdata=NULL;
05274 double dw=0;
05275 double dr=0;
05276 double ws=0;
05277 double we=0;
05278 cknull(*t,"null input table");
05279 check_nomsg(cpl_table_new_column(*t,"CNV",CPL_TYPE_DOUBLE));
05280 check_nomsg(pidata=cpl_table_get_data_double(*t,"INT"));
05281 check_nomsg(pcdata=cpl_table_get_data_double(*t,"CNV"));
05282 check_nomsg(ws=cpl_table_get_column_min(*t,"WAVE"));
05283 check_nomsg(we=cpl_table_get_column_max(*t,"WAVE"));
05284 check_nomsg(np=cpl_table_get_nrow(*t));
05285 dw=(we-ws)/(np-1);
05286 dr=(we-ws)/(rad-1);
05287
05288 for(i=0;i<rad;i++) {
05289 pcdata[i]=0.;
05290 }
05291 for(i=np-rad;i<np;i++) {
05292 pcdata[i]=0.;
05293 }
05294 for(i=rad;i<np-rad;i++) {
05295 val=0;
05296 for(j=-rad;j<rad;j++) {
05297 val+=pidata[i+j]*k*pow(2.0,-2.0*fabs(i-rad)/fwhm);
05298 }
05299
05300 check_nomsg(cpl_table_set_double(*t,"CNV",i,val));
05301 }
05302 return 0;
05303
05304 cleanup:
05305 return -1;
05306
05307 }
05308
05309
05318 int
05319 sinfo_convolve_gauss(cpl_table** t, const int rad, const double fwhm)
05320 {
05321 int i=0;
05322 int j=0;
05323 int np=0;
05324 double val=0;
05325 double sigma=fwhm/2.3548;
05326 double sigma2=sigma*sigma;
05327 double* pidata=NULL;
05328 double* pcdata=NULL;
05329 double dw=0;
05330 double dr=0;
05331 double ws=0;
05332 double we=0;
05333 double tx=0;
05334
05335 cknull(*t,"null input table");
05336 check_nomsg(cpl_table_new_column(*t,"CNV",CPL_TYPE_DOUBLE));
05337 check_nomsg(pidata=cpl_table_get_data_double(*t,"INT"));
05338 check_nomsg(pcdata=cpl_table_get_data_double(*t,"CNV"));
05339 check_nomsg(ws=cpl_table_get_column_min(*t,"WAVE"));
05340 check_nomsg(we=cpl_table_get_column_max(*t,"WAVE"));
05341 check_nomsg(np=cpl_table_get_nrow(*t));
05342 dw=(we-ws)/(np-1);
05343 dr=(we-ws)/(rad-1);
05344
05345 for(i=0;i<rad;i++) {
05346 pcdata[i]=0.;
05347 }
05348 for(i=np-rad;i<np;i++) {
05349 pcdata[i]=0.;
05350 }
05351 for(i=rad;i<np-rad;i++) {
05352 val=0;
05353 for(j=-rad;j<rad;j++) {
05354 tx=i-rad;
05355 val+=pidata[i+j]*exp(-tx*tx/2.0/sigma2)/(sigma*sqrt(2.0*PI_NUMB));
05356 }
05357
05358 check_nomsg(cpl_table_set_double(*t,"CNV",i,val));
05359 }
05360 return 0;
05361
05362 cleanup:
05363 return -1;
05364
05365 }
05366
05378 int
05379 sinfo_scales_obj_sky_cubes(cpl_frame* obj_frm,
05380 cpl_frame* sky_frm,
05381 cpl_table* bkg,
05382 cpl_table* rscale,
05383 cpl_imagelist** obj_cor)
05384 {
05385
05386
05387 int i=0;
05388 int j=0;
05389 int k=0;
05390 int xsz=0;
05391 int ysz=0;
05392 int zsz=0;
05393
05394 double* podata=NULL;
05395 double* psdata=NULL;
05396 double* pbdata=NULL;
05397 double* pcdata=NULL;
05398 double* pscale=NULL;
05399
05400
05401 cpl_image* imgo=NULL;
05402 cpl_image* imgs=NULL;
05403 cpl_image* imgc=NULL;
05404 cpl_propertylist* plist=NULL;
05405 cpl_imagelist* obj_cub=NULL;
05406 cpl_imagelist* sky_cub=NULL;
05407
05408
05409
05410
05411 cknull_nomsg(plist=cpl_propertylist_load(cpl_frame_get_filename(obj_frm),0));
05412 check_nomsg(xsz=sinfo_pfits_get_naxis1(plist));
05413 check_nomsg(ysz=sinfo_pfits_get_naxis2(plist));
05414 check_nomsg(zsz=sinfo_pfits_get_naxis3(plist));
05415 sinfo_free_propertylist(&plist);
05416 cknull_nomsg(obj_cub=cpl_imagelist_load(cpl_frame_get_filename(obj_frm),
05417 CPL_TYPE_DOUBLE,0));
05418
05419
05420
05421 cknull_nomsg(plist=cpl_propertylist_load(cpl_frame_get_filename(sky_frm),0));
05422 check_nomsg(xsz=sinfo_pfits_get_naxis1(plist));
05423 check_nomsg(ysz=sinfo_pfits_get_naxis2(plist));
05424 check_nomsg(zsz=sinfo_pfits_get_naxis3(plist));
05425 sinfo_free_propertylist(&plist);
05426 cknull_nomsg(sky_cub=cpl_imagelist_load(cpl_frame_get_filename(sky_frm),
05427 CPL_TYPE_DOUBLE,0));
05428
05429 check_nomsg(*obj_cor=cpl_imagelist_duplicate(obj_cub));
05430
05431 for(k=0;k<zsz;k++) {
05432 check_nomsg(imgo=cpl_imagelist_get(obj_cub,k));
05433 check_nomsg(imgc=cpl_imagelist_get(*obj_cor,k));
05434 check_nomsg(imgs=cpl_imagelist_get(sky_cub,k));
05435
05436 check_nomsg(podata=cpl_image_get_data_double(imgo));
05437 check_nomsg(pcdata=cpl_image_get_data_double(imgc));
05438 check_nomsg(psdata=cpl_image_get_data_double(imgs));
05439 check_nomsg(pbdata=cpl_table_get_data_double(bkg,"INT2"));
05440 check_nomsg(pscale=cpl_table_get_data_double(rscale,"RATIO"));
05441
05442 for (j=0;j<ysz; j++) {
05443 for (i=0;i<xsz; i++) {
05444 if(!isnan(podata[i+j*xsz]) &&
05445 !isnan(psdata[i+j*xsz]) &&
05446 !isnan(pbdata[k]) &&
05447 !isnan(pscale[k])) {
05448 pcdata[i+j*xsz] = podata[i+j*xsz]-
05449 (psdata[i+j*xsz]-pbdata[k])*pscale[k];
05450 }
05451 }
05452 }
05453 }
05454
05455 sinfo_free_imagelist(&sky_cub);
05456 sinfo_free_imagelist(&obj_cub);
05457
05458 return 0;
05459 cleanup:
05460
05461
05462 sinfo_free_imagelist(&sky_cub);
05463 sinfo_free_imagelist(&obj_cub);
05464 return -1;
05465 }
05466
05467
05484 static int
05485 sinfo_fitbkg(const double x[],
05486 const double a[],
05487 double *result)
05488 {
05489
05490
05491 double fac = sinfo_fac(x[0],a[2]);
05492
05493
05494
05495
05496
05497
05498
05499
05500 *result = a[0]+a[1]*fac;
05501
05502 return 0;
05503 }
05504
05528 static int
05529 sinfo_fitbkg_derivative(const double x[],
05530 const double a[],
05531 double d[])
05532 {
05533 double c=14387.7;
05534
05535
05536
05537
05538
05539 double fac = sinfo_fac(x[0],a[2]);
05540 double f2=0;
05541 double f1=0;
05542 double da=0.001;
05543 f1=a[0]+a[1]*fac;
05544
05545
05546 f2=a[0]+a[1]*sinfo_fac(x[0],a[2]+da*a[2]);
05547 d[0]=1.;
05548 d[1]=fac;
05549 d[2]=a[1]*fac*fac*x[0]*x[0]*x[0]*x[0]*c/(a[2]*a[2])*exp(c/(x[0]*a[2]));
05550
05551
05552
05553
05554 return 0;
05555 }
05556
05557
05558
05572 static double
05573 sinfo_fac(const double x, const double t)
05574 {
05575
05576 double c=14387.7;
05577
05578
05579 return pow(x,-5.)/(exp(c/(x*fabs(t)))-1.);
05580 }
05581
05591 static int
05592 sinfo_table_threshold(cpl_table** t,
05593 const char* column,
05594 const double low_cut,
05595 const double hig_cut,
05596 const double low_ass,
05597 const double hig_ass)
05598 {
05599
05600 int nrow=0;
05601 int i=0;
05602 double* pdata=NULL;
05603 cknull(*t,"null input table!");
05604
05605 check_nomsg(nrow=cpl_table_get_nrow(*t));
05606 check_nomsg(pdata=cpl_table_get_data_double(*t,column));
05607
05608 for(i=0;i<nrow;i++) {
05609
05610 if(pdata[i]<low_cut) {
05611 pdata[i]=low_ass;
05612 }
05613 if (pdata[i] >= hig_cut) {
05614 pdata[i]=hig_ass;
05615 }
05616
05617 }
05618
05619 return 0;
05620
05621 cleanup:
05622
05623 return -1;
05624 }
05625
05633 static int
05634 sinfo_table_set_column_invalid(cpl_table** t,const char* col)
05635 {
05636 int nrow=0;
05637 int i=0;
05638 cknull(*t,"input table is NULL");
05639 check_nomsg(nrow=cpl_table_get_nrow(*t));
05640 for(i=0;i<nrow;i++) {
05641 if( cpl_table_is_selected(*t,i) ) {
05642 cpl_table_set_invalid(*t,col,i);
05643 }
05644 }
05645 return 0;
05646 cleanup:
05647 return -1;
05648
05649 }
05650
05651
05652
05653
05654 static int
05655 sinfo_table_set(cpl_table** inp,
05656 const cpl_table* ref,
05657 const double val,
05658 const double tol)
05659 {
05660
05661 int ninp=0;
05662 int nref=0;
05663 double* piw=NULL;
05664 double* prw=NULL;
05665 double* pir=NULL;
05666 int i=0;
05667 int k=0;
05668 cknull(*inp,"NULL input table");
05669 cknull(ref,"NULL reference table");
05670
05671 check_nomsg(ninp=cpl_table_get_nrow(*inp));
05672 check_nomsg(nref=cpl_table_get_nrow(ref));
05673
05674 check_nomsg(prw=cpl_table_get_data_double(ref,"WAVE"));
05675 check_nomsg(piw=cpl_table_get_data_double(*inp,"WAVE"));
05676 check_nomsg(pir=cpl_table_get_data_double(*inp,"RATIO"));
05677
05678
05679 for(i=0;i<ninp;i++) {
05680
05681 if(fabs(piw[i]-prw[k]) < tol) {
05682 check_nomsg(cpl_table_set_double(*inp,"RATIO",i,val));
05683 k++;
05684 }
05685 }
05686 return 0;
05687
05688 cleanup:
05689
05690 return -1;
05691
05692 }
05693
05694
05695
05696 static cpl_table*
05697 sinfo_table_shift_simple(cpl_table* inp,
05698 const char* col,
05699 const double shift)
05700 {
05701
05702 int nrow=0;
05703 cpl_table* out=NULL;
05704 int is=(int)shift;
05705 double ds=shift-is;
05706 double* pi=NULL;
05707 double* po=NULL;
05708 double m=0;
05709 int i=0;
05710 cknull(inp,"null input table");
05711
05712 check_nomsg(nrow=cpl_table_get_nrow(inp));
05713 check_nomsg(out=cpl_table_duplicate(inp));
05714 check_nomsg(cpl_table_fill_column_window(out,col,0,nrow,0));
05715 check_nomsg(pi=cpl_table_get_data_double(inp,col));
05716 check_nomsg(po=cpl_table_get_data_double(out,col));
05717
05718
05719 for(i=0;i<nrow;i++) {
05720 if((i+is)>0 && (i+is+1) < nrow) {
05721 m=pi[i+is+1]-pi[i+is];
05722 po[i]=pi[i+is]+m*ds;
05723 }
05724 }
05725 return out;
05726 cleanup:
05727 sinfo_free_table(&out);
05728 return NULL;
05729
05730 }
05731
05732
05733
05734
05735 static cpl_imagelist*
05736 sinfo_cube_zshift_simple(cpl_imagelist* inp,
05737 const float shift)
05738 {
05739
05740 int nx=0;
05741 int ny=0;
05742 int nz=0;
05743
05744 int i=0;
05745 int j=0;
05746 int k=0;
05747 int ks=(int)shift;
05748
05749 float ds=shift-ks;
05750 float* pu=NULL;
05751 float* pl=NULL;
05752 float* po=NULL;
05753
05754 float int2=0;
05755 float int1=0;
05756 float m=0;
05757
05758 cpl_imagelist* out=NULL;
05759 cpl_image* imgu=NULL;
05760 cpl_image* imgl=NULL;
05761 cpl_image* imgo=NULL;
05762
05763
05764 cknull(inp,"null input cube");
05765
05766 check_nomsg(nz=cpl_imagelist_get_size(inp));
05767 check_nomsg(out=cpl_imagelist_duplicate(inp));
05768 check_nomsg(imgo=cpl_imagelist_get(out,0));
05769 check_nomsg(nx=cpl_image_get_size_x(imgo));
05770 check_nomsg(ny=cpl_image_get_size_y(imgo));
05771
05772 for(k=0;k<nz;k++) {
05773 if((k+ks)>0 && (k+ks+1) < nz) {
05774
05775 check_nomsg(imgu=cpl_imagelist_get(inp,k+ks+1));
05776 check_nomsg(imgl=cpl_imagelist_get(inp,k+ks));
05777 check_nomsg(imgo=cpl_imagelist_get(out,k));
05778
05779 check_nomsg(pu=cpl_image_get_data_float(imgu));
05780 check_nomsg(pl=cpl_image_get_data_float(imgl));
05781 check_nomsg(po=cpl_image_get_data_float(imgo));
05782
05783
05784 for(j=0;j<ny;j++) {
05785 for(i=0;i<nx;i++) {
05786 int2=pu[nx*j+i];
05787 int1=pl[nx*j+i];
05788 m=int2-int1;
05789 po[nx*j+i]=int1+m*ds;
05790 }
05791 }
05792 }
05793
05794
05795 }
05796 return out;
05797 cleanup:
05798 sinfo_free_imagelist(&out);
05799 return NULL;
05800
05801 }
05802
05803
05814 static int
05815 sinfo_get_line_ratio(cpl_table* obj_lin,
05816 cpl_table* obj_cnt,
05817 cpl_table* sky_lin,
05818 cpl_table* sky_cnt,
05819 const int method,
05820 double* r)
05821 {
05822
05823 int nobj;
05824 int nsky;
05825 int i=0;
05826
05827 cpl_table* obj_dif=NULL;
05828 cpl_table* sky_dif=NULL;
05829
05830 double* poi=NULL;
05831 double* psi=NULL;
05832 double* pvd=NULL;
05833 double* pvn=NULL;
05834 double* pvr=NULL;
05835
05836 cpl_vector* num=NULL;
05837 cpl_vector* den=NULL;
05838 cpl_vector* rat=NULL;
05839 cpl_vector* wav=NULL;
05840 double mnum=0;
05841 double mden=0;
05842 double tnum=0;
05843 double tden=0;
05844 int pows[2];
05845 cpl_polynomial* cfit=NULL;
05846 double mse=0;
05847
05848 cknull(obj_lin,"null obj line table");
05849 cknull(sky_lin,"null sky line table");
05850
05851 cknull(obj_cnt,"null obj cont table");
05852 cknull(sky_cnt,"null sky cont table");
05853
05854
05855 cknull_nomsg(obj_dif=sinfo_table_subtract_continuum(obj_lin,obj_cnt));
05856 cknull_nomsg(sky_dif=sinfo_table_subtract_continuum(sky_lin,sky_cnt));
05857
05858 check_nomsg(nobj=cpl_table_get_nrow(obj_dif));
05859 check_nomsg(nsky=cpl_table_get_nrow(sky_dif));
05860
05861
05862
05863 if(nobj != nsky) {
05864 sinfo_msg_error("obj and sky table must have the same no of rows!");
05865 sinfo_msg_error("nobj=%d nsky=%d",nobj,nsky);
05866 goto cleanup;
05867 }
05868
05869
05870 if(method == 0) {
05871 ck0_nomsg(sinfo_get_line_ratio_amoeba(obj_dif,sky_dif,r));
05872 sinfo_free_table(&obj_dif);
05873 sinfo_free_table(&sky_dif);
05874 return 0;
05875 }
05876
05877
05878 check_nomsg(poi=cpl_table_get_data_double(obj_dif,"INT"));
05879 check_nomsg(psi=cpl_table_get_data_double(sky_dif,"INT"));
05880
05881 check_nomsg(num=cpl_vector_new(nobj));
05882 check_nomsg(den=cpl_vector_new(nobj));
05883 check_nomsg(rat=cpl_vector_new(nobj));
05884 check_nomsg(cpl_vector_fill(num,0));
05885 check_nomsg(cpl_vector_fill(den,0));
05886 check_nomsg(cpl_vector_fill(rat,0));
05887 check_nomsg(pvd=cpl_vector_get_data(den));
05888 check_nomsg(pvn=cpl_vector_get_data(num));
05889 check_nomsg(pvr=cpl_vector_get_data(rat));
05890
05891 for(i=0;i<nobj;i++) {
05892 if(!isnan(psi[i]) && !isnan(poi[i]) && !isinf(psi[i]) && !isinf(poi[i]) ) {
05893 pvn[i]=psi[i]*poi[i];
05894 pvd[i]=psi[i]*psi[i];
05895 if(psi[i] != 0) {
05896 pvr[i]=poi[i]/psi[i];
05897 }
05898 }
05899 }
05900 sinfo_free_table(&sky_dif);
05901
05902 check_nomsg(mnum=cpl_vector_get_median(num));
05903 check_nomsg(mden=cpl_vector_get_median(den));
05904 check_nomsg(tnum=cpl_vector_get_mean(num)*nobj);
05905 check_nomsg(tden=cpl_vector_get_mean(den)*nobj);
05906
05907
05908
05909 sinfoni_free_vector(&num);
05910 sinfoni_free_vector(&den);
05911 if(method == 1) {
05912 *r=tnum/tden;
05913 } else if (method == 2) {
05914 *r=mnum/mden;
05915 } else if (method == 3) {
05916 *r=cpl_vector_get_median(rat);
05917 } else if (method == 4) {
05918 *r=cpl_vector_get_mean(rat);
05919 } else if (method == 5) {
05920
05921 check_nomsg(wav=cpl_vector_wrap(nobj,
05922 cpl_table_get_data_double(obj_dif,"WAVE")));
05923 check_nomsg(cfit=cpl_polynomial_fit_1d_create(wav,rat,0,&mse));
05924 sinfo_unwrap_vector(&wav);
05925 pows[0]=0;
05926 pows[1]=0;
05927 check_nomsg(*r=cpl_polynomial_get_coeff(cfit,pows));
05928 sinfo_free_polynomial(&cfit);
05929
05930 }
05931
05932 sinfo_free_table(&obj_dif);
05933 sinfoni_free_vector(&rat);
05934 return 0;
05935 cleanup:
05936 sinfo_free_table(&obj_dif);
05937 sinfo_free_table(&sky_dif);
05938 sinfoni_free_vector(&num);
05939 sinfoni_free_vector(&den);
05940 sinfoni_free_vector(&rat);
05941 sinfo_unwrap_vector(&wav);
05942
05943 return -1;
05944 }
05945
05946
05947
05948
05958 static int
05959 sinfo_get_line_ratio_amoeba(cpl_table* obj,
05960 cpl_table* sky,
05961 double* r)
05962 {
05963
05964
05965 int i=0;
05966 const int MP=2;
05967 const int NP=1;
05968 double y[MP];
05969 double p0[NP];
05970 double** ap=NULL;
05971 int nfunc=0;
05972 int np=0;
05973 check_nomsg(np=cpl_table_get_nrow(obj));
05974 check_nomsg(sa_ox=cpl_vector_wrap(np,cpl_table_get_data_double(obj,"WAVE")));
05975 check_nomsg(sa_oy=cpl_vector_wrap(np,cpl_table_get_data_double(obj,"INT")));
05976 check_nomsg(sa_sy=cpl_vector_wrap(np,cpl_table_get_data_double(sky,"INT")));
05977
05978
05979
05980 ap=(double**) cpl_calloc(MP,sizeof(double*));
05981 for(i=0;i<MP;i++) {
05982 ap[i]=cpl_calloc(NP,sizeof(double));
05983 }
05984
05985 ap[0][0]=-1.;
05986 ap[1][0]=+1.;
05987
05988
05989
05990 for(i=0;i<MP;i++) {
05991 p0[0]=ap[i][0];
05992 y[i]=sinfo_fit_sky(p0);
05993 }
05994
05995
05996 check_nomsg(sinfo_fit_amoeba(ap,y,NP,AMOEBA_FTOL,sinfo_fit_sky,&nfunc));
05997
05998 sinfo_msg("After amoeba fit");
05999 sinfo_msg("ap[0][0]=%g ap[0][1]=%g",ap[0][0],ap[1][0]);
06000
06001 *r=ap[0][0];
06002
06003 sinfo_unwrap_vector(&sa_ox);
06004 sinfo_unwrap_vector(&sa_oy);
06005 sinfo_unwrap_vector(&sa_sy);
06006 sinfo_new_destroy_2Ddoublearray(&ap,MP);
06007
06008
06009 return 0;
06010
06011 cleanup:
06012 sinfo_unwrap_vector(&sa_ox);
06013 sinfo_unwrap_vector(&sa_oy);
06014 sinfo_unwrap_vector(&sa_sy);
06015 sinfo_new_destroy_2Ddoublearray(&ap,MP);
06016
06017 return -1;
06018
06019 }
06020
06021
06022
06031
06032
06033 static double
06034 sinfo_fit_sky(double p[])
06035
06036 {
06037 double* ps=NULL;
06038 double* po=NULL;
06039 double* pv=NULL;
06040 cpl_vector* vtmp=NULL;
06041 int i=0;
06042 int np=0;
06043 int pows[2];
06044 double mse=0;
06045 double cont=0;
06046 cpl_polynomial* pfit=NULL;
06047
06048 double rms=0;
06049
06050
06051
06052 check_nomsg(pfit=cpl_polynomial_fit_1d_create(sa_ox,sa_oy,0,&mse));
06053 pows[0]=0;
06054 pows[1]=0;
06055 check_nomsg(cont=cpl_polynomial_get_coeff(pfit,pows));
06056 check_nomsg(sinfo_free_polynomial(&pfit));
06057 check_nomsg(cpl_vector_subtract_scalar(sa_oy,cont));
06058
06059
06060
06061 check_nomsg(pfit=cpl_polynomial_fit_1d_create(sa_ox,sa_sy,0,&mse));
06062 pows[0]=0;
06063 pows[1]=0;
06064 check_nomsg(cont=cpl_polynomial_get_coeff(pfit,pows));
06065 check_nomsg(sinfo_free_polynomial(&pfit));
06066 check_nomsg(cpl_vector_subtract_scalar(sa_sy,cont));
06067
06068
06069 check_nomsg(po= cpl_vector_get_data(sa_oy));
06070 check_nomsg(ps= cpl_vector_get_data(sa_sy));
06071
06072 check_nomsg(np=cpl_vector_get_size(sa_oy));
06073 check_nomsg(vtmp=cpl_vector_new(np));
06074 check_nomsg(pv= cpl_vector_get_data(vtmp));
06075
06076
06077 for(i=0;i<np;i++) {
06078 pv[i]=po[i]-ps[i]*p[0];
06079 }
06080
06081 check_nomsg(rms=cpl_vector_get_stdev(vtmp));
06082 sinfoni_free_vector(&vtmp);
06083 return rms;
06084 cleanup:
06085 sinfoni_free_vector(&vtmp);
06086 return -1;
06087
06088 }
06089
06090
06091
06103 static cpl_table*
06104 sinfo_table_interpol(cpl_table* obj_lin,
06105 cpl_table* obj_cnt,
06106 cpl_table* sky_lin,
06107 cpl_table* sky_cnt,
06108 const double r)
06109 {
06110
06111 cpl_table* out=NULL;
06112 cpl_table* obj_dif=NULL;
06113 cpl_table* sky_dif=NULL;
06114 cknull(obj_lin,"null line table");
06115 cknull(obj_cnt,"null cont table");
06116
06117 cknull_nomsg(obj_dif=sinfo_table_subtract_continuum(obj_lin,obj_cnt));
06118 cknull_nomsg(sky_dif=sinfo_table_subtract_continuum(sky_lin,sky_cnt));
06119
06120 check_nomsg(out=cpl_table_duplicate(obj_dif));
06121 check_nomsg(cpl_table_duplicate_column(out,"CSKY",sky_dif,"INT"));
06122 check_nomsg(cpl_table_multiply_scalar(out,"CSKY",r));
06123 check_nomsg(cpl_table_subtract_columns(out,"INT","CSKY"));
06124
06125 sinfo_free_table(&obj_dif);
06126 sinfo_free_table(&sky_dif);
06127
06128 return out;
06129
06130 cleanup:
06131 sinfo_free_table(&obj_dif);
06132 sinfo_free_table(&sky_dif);
06133 sinfo_free_table(&out);
06134 return NULL;
06135
06136 }
06137
06138
06139
06140
06141
06142
06151 static cpl_table*
06152 sinfo_table_subtract_continuum(cpl_table* lin,
06153 cpl_table* cnt)
06154
06155 {
06156
06157 cpl_table* out=NULL;
06158 int nlin=0;
06159 int ncnt=0;
06160 int i=0;
06161 double* poi=NULL;
06162 cpl_vector* vx=NULL;
06163 cpl_vector* vy=NULL;
06164 cpl_polynomial* cfit=NULL;
06165 int pows[2];
06166 double mse=0;
06167 double yfit=0;
06168
06169 cknull(lin,"null line table");
06170 cknull(cnt,"null cont table");
06171 check_nomsg(out=cpl_table_duplicate(lin));
06172 check_nomsg(cpl_table_new_column(out,"CONT",CPL_TYPE_DOUBLE));
06173 check_nomsg(nlin=cpl_table_get_nrow(lin));
06174 check_nomsg(ncnt=cpl_table_get_nrow(cnt));
06175
06176 check_nomsg(cpl_table_fill_column_window(out,"CONT",0,nlin,0));
06177
06178
06179 check_nomsg(vx=cpl_vector_wrap(ncnt,cpl_table_get_data_double(cnt,"WAVE")));
06180 check_nomsg(vy=cpl_vector_wrap(ncnt,cpl_table_get_data_double(cnt,"INT")));
06181 check_nomsg(cfit=cpl_polynomial_fit_1d_create(vx,vy,0,&mse));
06182 sinfo_unwrap_vector(&vx);
06183 sinfo_unwrap_vector(&vy);
06184
06185 pows[0]=0;
06186 pows[1]=0;
06187 check_nomsg(yfit=cpl_polynomial_get_coeff(cfit,pows));
06188 sinfo_free_polynomial(&cfit);
06189
06190
06191 check_nomsg(poi=cpl_table_get_data_double(out,"CONT"));
06192 for(i=0;i<nlin;i++) {
06193 poi[i]=yfit;
06194 }
06195
06196 check_nomsg(cpl_table_subtract_columns(out,"INT","CONT"));
06197 check_nomsg(cpl_table_erase_column(out,"CONT"));
06198
06199
06200
06201 return out;
06202
06203 cleanup:
06204 sinfo_unwrap_vector(&vx);
06205 sinfo_unwrap_vector(&vy);
06206 sinfo_free_polynomial(&cfit);
06207 sinfo_free_table(&out);
06208 return NULL;
06209
06210 }
06211
06212
06213 static int
06214 sinfo_compute_line_ratio(cpl_table* obj,
06215 cpl_table* sky,
06216 const double wtol,
06217 const int meth,
06218 const cpl_table* sel_regions,
06219 cpl_table* cont_regions,
06220 double* r)
06221 {
06222 cpl_table* line_regions=NULL;
06223 cpl_table* obj_cnt=NULL;
06224 cpl_table* sky_cnt=NULL;
06225 cpl_table* obj_lin=NULL;
06226 cpl_table* sky_lin=NULL;
06227 cpl_table* lres=NULL;
06228 double fmed=0;
06229 double fsdv=0;
06230 double fthresh=0;
06231 int fclip_i=0;
06232 int line_i=0;
06233
06234
06235
06236 check_nomsg(line_regions = cpl_table_duplicate(sel_regions));
06237
06238
06239 check_nomsg(obj_lin=sinfo_table_select_range(obj,line_regions,wtol));
06240 check_nomsg(sky_lin=sinfo_table_select_range(sky,line_regions,wtol));
06241 check_nomsg(obj_cnt=sinfo_table_select_range(obj,cont_regions,wtol));
06242 check_nomsg(sky_cnt=sinfo_table_select_range(sky,cont_regions,wtol));
06243
06244 ck0_nomsg(sinfo_get_line_ratio(obj_lin,obj_cnt,sky_lin,sky_cnt,meth,r));
06245
06246
06247
06248
06249
06250
06251
06252
06253
06254 check_nomsg(lres=sinfo_table_interpol(obj_lin,obj_cnt,sky_lin,sky_cnt,*r));
06255
06256 check_nomsg(fmed = cpl_table_get_column_median(lres,"INT"));
06257 check_nomsg(fsdv = cpl_table_get_column_stdev(lres,"INT"));
06258 fthresh=fmed+3*fsdv;
06259
06260 check_nomsg(cpl_table_duplicate_column(lres,"AINT",lres,"INT"));
06261 check_nomsg(cpl_table_multiply_columns(lres,"AINT","INT"));
06262 check_nomsg(cpl_table_power_column(lres,"AINT",0.5));
06263 check_nomsg(fclip_i=cpl_table_and_selected_double(lres,"AINT",
06264 CPL_GREATER_THAN,
06265 fthresh));
06266 check_nomsg(cpl_table_select_all(lres));
06267
06268
06269 if (fclip_i > 0) {
06270
06271 check_nomsg(line_i=cpl_table_and_selected_double(lres,"AINT",
06272 CPL_LESS_THAN,
06273 fthresh));
06274 sinfo_free_table(&line_regions);
06275 check_nomsg(line_regions=cpl_table_extract_selected(lres));
06276 sinfo_free_table(&lres);
06277
06278 check_nomsg(cpl_table_erase_column(line_regions,"INT"));
06279 check_nomsg(cpl_table_erase_column(line_regions,"AINT"));
06280
06281
06282 if (line_i >= 3) {
06283
06284 sinfo_free_table(&obj_lin);
06285 sinfo_free_table(&sky_lin);
06286 check_nomsg(obj_lin=sinfo_table_select_range(obj,line_regions,wtol));
06287 check_nomsg(sky_lin=sinfo_table_select_range(sky,line_regions,wtol));
06288
06289 sinfo_free_table(&line_regions);
06290
06291
06292
06293 ck0_nomsg(sinfo_get_line_ratio(obj_lin,obj_cnt,sky_lin,sky_cnt,meth,r));
06294 }
06295 }
06296 *r=fabs(*r);
06297
06298 sinfo_free_table(&obj_cnt);
06299 sinfo_free_table(&sky_cnt);
06300 sinfo_free_table(&sky_lin);
06301 sinfo_free_table(&obj_lin);
06302 sinfo_free_table(&lres);
06303 sinfo_free_table(&line_regions);
06304
06305
06306 return 0;
06307
06308
06309 cleanup:
06310
06311
06312 sinfo_free_table(&obj_cnt);
06313 sinfo_free_table(&sky_cnt);
06314 sinfo_free_table(&sky_lin);
06315 sinfo_free_table(&obj_lin);
06316
06317 sinfo_free_table(&lres);
06318 sinfo_free_table(&line_regions);
06319
06320 return -1;
06321
06322 }
06334 static cpl_table*
06335 sinfo_find_rot_waves(
06336 const double w_rot[],
06337 const int npix_w,
06338 const double w_step,
06339 cpl_table* range
06340 )
06341 {
06342 int i=0;
06343 int x_i=0;
06344 int r_start=0;
06345 double w_min=0;
06346 double w_max=0;
06347
06348 cpl_table* w_sel=NULL;
06349 cpl_table* res=NULL;
06350
06351 check_nomsg(res = cpl_table_new(0));
06352
06353 check_nomsg(cpl_table_copy_structure(res,range));
06354
06355 for (i=0; i< NROT; i++) {
06356
06357
06358
06359
06360 w_min=w_rot[i]-npix_w*w_step;
06361 w_max=w_rot[i]+npix_w*w_step;
06362
06363 check_nomsg(cpl_table_and_selected_double(range,"WAVE",
06364 CPL_GREATER_THAN,w_min));
06365 check_nomsg(cpl_table_and_selected_double(range,"WAVE",
06366 CPL_LESS_THAN,w_max));
06367 sinfo_free_table(&w_sel);
06368 check_nomsg(w_sel=cpl_table_extract_selected(range));
06369 check_nomsg(x_i=cpl_table_get_nrow(w_sel));
06370
06371 if (x_i > 0) {
06372 check_nomsg(r_start=cpl_table_get_nrow(res));
06373
06374 check_nomsg(cpl_table_insert(res,w_sel,r_start));
06375 }
06376 check_nomsg(cpl_table_select_all(range));
06377 }
06378
06379
06380 sinfo_free_table(&w_sel);
06381
06382
06383 return res;
06384
06385 cleanup:
06386 sinfo_free_table(&w_sel);
06387 sinfo_free_table(&res);
06388 return NULL;
06389
06390 }
06391
06404 static int
06405 sinfo_get_obj_sky_wav_sub(cpl_table* obj,
06406 cpl_table* sky,
06407 cpl_table* wav,
06408 cpl_table* sel,
06409 const double wtol,
06410 cpl_table** sub_obj,
06411 cpl_table** sub_sky,
06412 cpl_table** sub_wav)
06413
06414 {
06415 cknull_nomsg(*sub_obj = sinfo_table_select_range(obj,sel,wtol));
06416 cknull_nomsg(*sub_sky = sinfo_table_select_range(sky,sel,wtol));
06417 cknull_nomsg(*sub_wav = sinfo_table_select_range(wav,sel,wtol));
06418 return 0;
06419
06420 cleanup:
06421 sinfo_free_table(&(*sub_obj));
06422 sinfo_free_table(&(*sub_sky));
06423 sinfo_free_table(&(*sub_wav));
06424
06425 return -1;
06426
06427 }
06428
06429 static int
06430 sinfo_get_sub_regions(cpl_table* sky,
06431 cpl_table* x1,
06432 cpl_table* pos,
06433 const double wtol,
06434 const int npixw,
06435 cpl_table** res)
06436 {
06437
06438 cpl_table* x1_sub=NULL;
06439 cpl_table* x2=NULL;
06440
06441 int nrow=0;
06442 int np=0;
06443
06444 cknull(sky,"Null input sky table");
06445 cknull(x1 ,"Null input x1 table");
06446 cknull(pos,"Null input pos table");
06447
06448 check_nomsg(x2=cpl_table_duplicate(sky));
06449 check_nomsg(nrow=cpl_table_get_nrow(sky));
06450 check_nomsg(cpl_table_fill_column_window(x2,"INT",0,nrow,0));
06451
06452
06453
06454
06455
06456
06457 x1_sub=sinfo_table_select_range(x1,pos,wtol);
06458
06459 if(x1_sub != NULL) {
06460 ck0_nomsg(sinfo_table_fill_column_over_range(&x2,x1_sub,"INT",10.,wtol));
06461 sinfo_free_table(&x1_sub);
06462 check_nomsg(sinfo_convolve_kernel(&x2,npixw/2));
06463 check_nomsg(np=cpl_table_and_selected_double(x2,"CNV",CPL_GREATER_THAN,0));
06464 check_nomsg(*res=cpl_table_extract_selected(x2));
06465 sinfo_free_table(&x2);
06466 check_nomsg(cpl_table_erase_column(*res,"INT"));
06467 check_nomsg(cpl_table_erase_column(*res,"CNV"));
06468
06469 } else {
06470
06471 irplib_error_reset();
06472 cpl_error_reset();
06473 sinfo_free_table(&x1_sub);
06474 sinfo_free_table(&x2);
06475
06476 return np;
06477 }
06478
06479 return np;
06480 cleanup:
06481
06482 sinfo_free_table(&x1_sub);
06483 sinfo_free_table(&x2);
06484 return -1;
06485
06486 }
06487
06488 static cpl_table*
06489 sinfo_table_extract_rest(cpl_table* inp,
06490 cpl_table* low,
06491 cpl_table* med,
06492 const double wtol)
06493 {
06494
06495 cpl_table* out=NULL;
06496 double* pinp=NULL;
06497 double* plow=NULL;
06498 double* pmed=NULL;
06499 int nlow=0;
06500 int nmed=0;
06501
06502 int nrow=0;
06503 int i=0;
06504 int k=0;
06505 cpl_table* tmp=NULL;
06506
06507 cknull(inp,"null input table");
06508
06509
06510 check_nomsg(tmp=cpl_table_duplicate(inp));
06511 check_nomsg(nrow=cpl_table_get_nrow(tmp));
06512 check_nomsg(cpl_table_new_column(tmp,"SEL",CPL_TYPE_INT));
06513 check_nomsg(cpl_table_fill_column_window_int(tmp,"SEL",0,nrow,0));
06514
06515 check_nomsg(pinp=cpl_table_get_data_double(inp,"WAVE"));
06516 check_nomsg(plow=cpl_table_get_data_double(low,"WAVE"));
06517 check_nomsg(pmed=cpl_table_get_data_double(med,"WAVE"));
06518 nlow=cpl_table_get_nrow(low);
06519
06520
06521
06522 if(nlow > 0) {
06523 k=0;
06524 for(i=0;i<nrow;i++) {
06525 if(fabs(pinp[i]-plow[k]) < wtol) {
06526 cpl_table_set_int(tmp,"SEL",k,-1);
06527 k++;
06528 }
06529 }
06530 }
06531 nmed=cpl_table_get_nrow(med);
06532
06533 k=0;
06534 if(nmed > 0) {
06535 for(i=0;i<nrow;i++) {
06536 if(fabs(pinp[i]-pmed[k]) < wtol) {
06537 cpl_table_set_int(tmp,"SEL",k,-1);
06538 k++;
06539 }
06540 }
06541 }
06542
06543 check_nomsg(cpl_table_and_selected_int(tmp,"SEL",CPL_GREATER_THAN,-1));
06544 check_nomsg(out=cpl_table_extract_selected(tmp));
06545 sinfo_free_table(&tmp);
06546 check_nomsg(cpl_table_erase_column(out,"SEL"));
06547
06548 return out;
06549
06550 cleanup:
06551 sinfo_free_table(&tmp);
06552 return NULL;
06553
06554 }
06555