/* $Id: sinfo_skycor.c,v 1.50 2012-05-04 08:11:35 amodigli Exp $ * * This file is part of the SINFONI Pipeline * Copyright (C) 2002,2003 European Southern Observatory * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by * the Free Software Foundation; either version 2 of the License, or * (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program; if not, write to the Free Software * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ /* * $Author: amodigli $ * $Date: 2012-05-04 08:11:35 $ * $Revision: 1.50 $ * $Name: not supported by cvs2svn $ */ #ifdef HAVE_CONFIG_H #include #endif /*----------------------------------------------------------------------------- Includes ----------------------------------------------------------------------------*/ #include #include #include #include #include #include #include "sinfo_pfits.h" #include "sinfo_fit.h" #include "sinfo_functions.h" #include "sinfo_msg.h" #include "sinfo_error.h" #include "sinfo_globals.h" #include "sinfo_utils_wrappers.h" #include "sinfo_utl_cube2spectrum.h" #include "sinfo_pro_types.h" /*----------------------------------------------------------------------------- Defines ----------------------------------------------------------------------------*/ #define BAND_H_WAVE_MIN 1.445 //not used #define BAND_H_WAVE_MAX 1.820 //not used #define BAND_K_WAVE_MIN 1.945 //not used #define BAND_K_WAVE_MAX 2.460 //not used #define BAND_J_WAVE_MIN 1.445 //not used #define BAND_J_WAVE_MAX 1.82 //not used #define SINFO_FIT_BKG_TEMP 280. #define SKY_THRES 0.95 #define SKY_LINE_MAX_CUT 4 /* this should be 4 */ #define SKY_LINE_MIN_CUT 4 /* this should be 4 */ #define XCOR_YSHIFT_KS_CLIP 5 /* this should be 5 */ #define HPLANK 6.62606876e-34; // J s #define CLIGHT 2.99792458e+08; // m / s #define KBOLTZ 1.3806503e-23; // J / K #define AMOEBA_FTOL 1.e-5 #define NBOUND 14 #define NROT 25 #define N_ITER_FIT_LM 15 #define N_ITER_FIT_AMOEBA 10 double sinfo_scale_fct=1; static cpl_vector* sa_vx=NULL; static cpl_vector* sa_vy=NULL; static cpl_vector* sa_ox=NULL; static cpl_vector* sa_oy=NULL; static cpl_vector* sa_sy=NULL; /*----------------------------------------------------------------------------- Functions prototypes -----------------------------------------------------------------------------*/ cpl_vector* sinfo_sky_background_estimate(cpl_vector *spectrum, int msize, int fsize); static int sinfo_scales_obj_sky_cubes(cpl_imagelist* obj_cub, cpl_imagelist* sky_cub, cpl_table* bkg, cpl_table* rscale, cpl_imagelist** obj_cor); static int sinfo_sub_thr_bkg_from_obj_cube(cpl_imagelist* obj_cub, cpl_table* int_sky, cpl_imagelist** obj_cor); static cpl_vector* sinfo_filter_min(const cpl_vector* vi, const int size); static cpl_vector* sinfo_filter_max(const cpl_vector* vi, const int size); static cpl_vector* sinfo_filter_smo(const cpl_vector* vi, const int size); static cpl_imagelist* sinfo_cube_zshift_simple(cpl_imagelist* inp, const float shift); static void sinfo_optimise_sky_sub(const double wtol, const double line_hw, const int method, const int do_rot, cpl_table* lrange, cpl_table* lambda, cpl_table* lr41, cpl_table* lr52, cpl_table* lr63, cpl_table* lr74, cpl_table* lr02, cpl_table* lr85, cpl_table* lr20, cpl_table* lr31, cpl_table* lr42, cpl_table* lr53, cpl_table* lr64, cpl_table* lr75, cpl_table* lr86, cpl_table* lr97, cpl_table* lr00, cpl_table** int_obj, cpl_table** int_sky, cpl_table** rscale); static void sinfo_shift_sky(cpl_frame** sky_frm, cpl_table** int_sky, const double zshift); static double sinfo_xcorr(cpl_table* int_obj, cpl_table* int_sky, cpl_table* lambda, const double dispersion, const double line_hw); static cpl_table* sinfo_interpolate_sky(const cpl_table* inp,const cpl_table* lambdas); static int sinfo_check_screw_values(cpl_table** int_obj, cpl_table** int_sky, cpl_table* grange, const double wtol); static int sinfo_choose_good_sky_pixels(cpl_frame* obj_frm, cpl_image* r_img, cpl_image* g_img, const double min_frac, cpl_image** mask); static int sinfo_sum_spectra(const cpl_frame* obj, const cpl_frame* sky, cpl_image* mask, cpl_table* wrange, const int llx, const int lly, const int urx, const int ury, cpl_table** int_obj, cpl_table** int_sky); int sinfo_thermal_background2(cpl_table* int_sky, cpl_table* lambda, cpl_table* lrange, cpl_table** bkg); static int sinfo_thermal_background(cpl_table* int_sky, cpl_table* lambda, cpl_table* lrange, const double temp, const int niter, const int filter_width, const double wtol, cpl_table** bkg, int* success_fit); static int sinfo_object_flag_sky_pixels(cpl_frame* obj_frm, cpl_table* lambda, cpl_table* mrange, cpl_imagelist* flag_data, const double tol, cpl_image** g_img, cpl_image** r_img, cpl_image** image); static int sinfo_object_flag_low_values(cpl_frame* obj_frm, const double cnt, const double sig, cpl_imagelist** flag_data); static int sinfo_object_estimate_noise(cpl_frame* obj_frm, const int obj_noise_fit, double* centre, double* noise); static int sinfo_set_ranges(cpl_frame* obj_frm, cpl_frame* sky_frm, cpl_parameterlist* cfg, cpl_table** lambda, cpl_table** lr41, cpl_table** lr52, cpl_table** lr63, cpl_table** lr74, cpl_table** lr02, cpl_table** lr85, cpl_table** lr20, cpl_table** lr31, cpl_table** lr42, cpl_table** lr53, cpl_table** lr64, cpl_table** lr75, cpl_table** lr86, cpl_table** lr97, cpl_table** lr00, cpl_table** lrange, cpl_table** grange, cpl_table** lambdas, cpl_table** mrange, int* sky_interp_sw, double* dispersion); static cpl_table* sinfo_table_extract_rest(cpl_table* inp, cpl_table* low, cpl_table* med, const double wtol); static int sinfo_get_sub_regions(cpl_table* sky, cpl_table* x1, cpl_table* pos, const double wtol, const int npixw, cpl_table** sub_regions); static int sinfo_get_obj_sky_wav_sub(cpl_table* obj, cpl_table* sky, cpl_table* wav, cpl_table* sel, const double wtol, cpl_table** sub_obj, cpl_table** sub_sky, cpl_table** sub_wav); static cpl_table* sinfo_find_rot_waves(const double w_rot[], const int npix_w, const double w_step, cpl_table* range); static int sinfo_compute_line_ratio(cpl_table* obj, cpl_table* sky, const double wtol, const int meth, const cpl_table* sel_regions, cpl_table* cont_regions, double* r); static cpl_table* sinfo_table_subtract_continuum(cpl_table* lin,cpl_table* cnt); static double sinfo_fit_bkg(double p[]); static double sinfo_fit_sky(double p[]); static int sinfo_get_line_ratio_amoeba(cpl_table* obj, cpl_table* sky, double* r); static cpl_table* sinfo_table_interpol(cpl_table* obj_lin, cpl_table* obj_cnt, cpl_table* sky_lin, cpl_table* sky_cnt, const double r); static int sinfo_get_line_ratio(cpl_table* obj_lin, cpl_table* obj_cnt, cpl_table* sky_lin, cpl_table* sky_cnt, const int method, double* r); static cpl_table* sinfo_table_shift_simple(cpl_table* inp, const char* col, const double shift); /* static int sinfo_table_set_column_invalid(cpl_table** int_sky,const char* col); */ static int sinfo_table_set(cpl_table** out, const cpl_table* ref, const double val, const double tol); static int sinfo_table_threshold(cpl_table** t, const char* column, const double low_cut, const double hig_cut, const double low_ass, const double hig_ass); static double sinfo_fac(const double x, const double t); static int sinfo_fitbkg(const double x[], const double a[], double *result); static int sinfo_fitbkg_derivative(const double x[], const double a[], double result[]); static int sinfo_convolve_kernel(cpl_table** t, const int rad); int sinfo_convolve_kernel2(cpl_table** t, const int rad); int sinfo_convolve_gauss(cpl_table** t, const int rad, const double fwhm); int sinfo_convolve_exp(cpl_table** t, const int rad, const double fwhm); static int sinfo_table_sky_obj_flag_nan(cpl_table** s,cpl_table** o, cpl_table** w); static int sinfo_table_set_nan_out_min_max(cpl_table** s, const char* c, const double min, const double max); static double sinfo_gaussian_amp(double area,double sigma,double x,double x0,double off); static double sinfo_gaussian_area(double amp,double sigma,double x,double x0,double off); int sinfo_table_smooth_column(cpl_table** t, const char* c, const int r); static int sinfo_image_flag_nan(cpl_image** im); static int sinfo_table_flag_nan(cpl_table** t,const char* label); static int sinfo_cnt_mask_thresh_and_obj_finite(const cpl_image* mask, const double t, const cpl_image* obj); static cpl_table* sinfo_interpolate(const cpl_table* inp, const cpl_table* lambdas, const char* name, const char* method); static cpl_table* sinfo_image2table(const cpl_image* im); static int sinfo_table_extract_finite(const cpl_table* in1, const cpl_table* in2, cpl_table** ou1, cpl_table** ou2); static cpl_table* sinfo_slice_z(const cpl_imagelist* cin,const int i,const int j); static cpl_imagelist* sinfo_imagelist_select_range(const cpl_imagelist* inp, const cpl_table* full, const cpl_table* good, const double tol); static cpl_table* sinfo_table_select_range(cpl_table* inp, cpl_table* ref, const double tol); static int sinfo_table_fill_column_over_range(cpl_table** inp, const cpl_table* ref, const char* col, const double val, const double tol); static int sinfo_table_column_dindgen(cpl_table** t, const char* label); /**@{*/ /*---------------------------------------------------------------------------*/ /** * @addtogroup sinfo_utl_skycor Functions to correct sky residuals on \ science cubes */ /*---------------------------------------------------------------------------*/ /** @name sinfo_skycor_qc_new @brief structure to init skycor_qc @return pointer to structure */ sinfo_skycor_qc* sinfo_skycor_qc_new(void) { sinfo_skycor_qc * sqc; sqc= cpl_malloc(sizeof(sinfo_skycor_qc)); sqc->th_fit=0; return sqc; } /** @name sinfo_skycor_qc_delete @brief function to free a skycor_qc structure @return void */ void sinfo_skycor_qc_delete(sinfo_skycor_qc** sqc) { if((*sqc) != NULL) { cpl_free(*sqc); *sqc=NULL; } } /** @brief Execute the plugin instance given by the interface @param config parameter configuration @param obj_frm input object frm @param sky_frm input sky frm @param sqc QC parameters @param obj_cor corrected object cube @param int_obj corrected object spectrum @return 0 if everything is ok, -1 else */ /*---------------------------------------------------------------------------*/ int sinfo_skycor(cpl_parameterlist * config, cpl_frame* obj_frm, cpl_frame* sky_frm, sinfo_skycor_qc* sqc, cpl_imagelist** obj_cor, cpl_table** int_obj) { cpl_table* bkg=NULL; cpl_table* lambda=NULL; cpl_table* lr41=NULL; cpl_table* lr52=NULL; cpl_table* lr63=NULL; cpl_table* lr74=NULL; cpl_table* lr02=NULL; cpl_table* lr85=NULL; cpl_table* lr20=NULL; cpl_table* lr31=NULL; cpl_table* lr42=NULL; cpl_table* lr53=NULL; cpl_table* lr64=NULL; cpl_table* lr75=NULL; cpl_table* lr86=NULL; cpl_table* lr97=NULL; cpl_table* lr00=NULL; cpl_table* lrange=NULL; cpl_table* mrange=NULL; cpl_table* grange=NULL; cpl_table* lambdas=NULL; cpl_table* int_sky=NULL; cpl_image* mask=NULL; cpl_image* gpix=NULL; cpl_image* ratio=NULL; cpl_image* ima_sky=NULL; cpl_imagelist* fdata=NULL; cpl_table* rscale=NULL; cpl_parameter* p=NULL; int th_fit=0; double dispersion=0; double noise=0; //double temp=252.69284; double centre=0; int sky_interp_sw=0; double wshift=0; double pshift=0; int method=0; int do_rot=0; int obj_noise_fit=0; int niter=3; double min_frac=0.8; double line_hw=7; double fit_temp=280; int filter_width=SINFO_SKY_BKG_FILTER_WIDTH; int llx=0; int lly=0; int urx=64; int ury=64; cpl_imagelist* obj_cub=NULL; cpl_imagelist* sky_cub=NULL; int sub_thr_bkg=0; check_nomsg(p=cpl_parameterlist_find(config, "sinfoni.sinfo_utl_skycor.min_frac")); check_nomsg(min_frac=cpl_parameter_get_double(p)); check_nomsg(p=cpl_parameterlist_find(config, "sinfoni.sinfo_utl_skycor.line_half_width")); check_nomsg(line_hw=cpl_parameter_get_double(p)); check_nomsg(p=cpl_parameterlist_find(config, "sinfoni.sinfo_utl_skycor.sky_bkg_filter_width")); check_nomsg(filter_width=cpl_parameter_get_int(p)); check_nomsg(p=cpl_parameterlist_find(config, "sinfoni.sinfo_utl_skycor.scale_method")); check_nomsg(method=cpl_parameter_get_int(p)); check_nomsg(p=cpl_parameterlist_find(config, "sinfoni.sinfo_utl_skycor.rot_cor")); check_nomsg(do_rot=cpl_parameter_get_bool(p)); check_nomsg(p=cpl_parameterlist_find(config, "sinfoni.sinfo_utl_skycor.sub_thr_bkg_from_obj")); check_nomsg(sub_thr_bkg=cpl_parameter_get_bool(p)); check_nomsg(p=cpl_parameterlist_find(config, "sinfoni.sinfo_utl_skycor.fit_obj_noise")); check_nomsg(obj_noise_fit=cpl_parameter_get_bool(p)); check_nomsg(p=cpl_parameterlist_find(config, "sinfoni.sinfo_utl_skycor.niter")); check_nomsg(niter=cpl_parameter_get_int(p)); check_nomsg(p=cpl_parameterlist_find(config, "sinfoni.sinfo_utl_skycor.pshift")); check_nomsg(pshift=cpl_parameter_get_double(p)); check_nomsg(p=cpl_parameterlist_find(config, "sinfoni.sinfo_utl_skycor.llx")); check_nomsg(llx=cpl_parameter_get_int(p)); check_nomsg(p=cpl_parameterlist_find(config, "sinfoni.sinfo_utl_skycor.lly")); check_nomsg(lly=cpl_parameter_get_int(p)); check_nomsg(p=cpl_parameterlist_find(config, "sinfoni.sinfo_utl_skycor.urx")); check_nomsg(urx=cpl_parameter_get_int(p)); check_nomsg(p=cpl_parameterlist_find(config, "sinfoni.sinfo_utl_skycor.ury")); check_nomsg(ury=cpl_parameter_get_int(p)); // set wavelength ranges sinfo_msg("Set wavelength ranges"); ck0(sinfo_set_ranges(obj_frm,sky_frm,config,&lambda, &lr41,&lr52,&lr63,&lr74,&lr02,&lr85,&lr20,&lr31,&lr42, &lr53,&lr64,&lr75,&lr86,&lr97,&lr00, &lrange,&grange,&lambdas,&mrange,&sky_interp_sw, &dispersion),"Setting wavelength ranges"); //check_nomsg(cpl_table_save(grange, NULL, NULL, "out_grange.fits", //CPL_IO_DEFAULT)); /* sinfo_msg("lr20=%d",cpl_table_get_nrow(lr20)); sinfo_msg("lr31=%d",cpl_table_get_nrow(lr31)); sinfo_msg("lr42=%d",cpl_table_get_nrow(lr42)); sinfo_msg("lr53=%d",cpl_table_get_nrow(lr53)); sinfo_msg("lr64=%d",cpl_table_get_nrow(lr64)); sinfo_msg("lr75=%d",cpl_table_get_nrow(lr75)); sinfo_msg("lr86=%d",cpl_table_get_nrow(lr86)); sinfo_msg("lr97=%d",cpl_table_get_nrow(lr97)); sinfo_msg("lr00=%d",cpl_table_get_nrow(lr00)); sinfo_msg("min_lrange=%f",cpl_table_get_column_min(lrange,"INDEX")); sinfo_msg("min_grange=%f",cpl_table_get_column_min(grange,"INDEX")); sinfo_msg("min_srange=%f",cpl_table_get_column_min(lambdas,"WAVE")); sinfo_msg("min_mrange=%f",cpl_table_get_column_min(mrange,"INDEX")); sinfo_msg("max_lrange=%f",cpl_table_get_column_max(lrange,"INDEX")); sinfo_msg("max_grange=%f",cpl_table_get_column_max(grange,"INDEX")); sinfo_msg("max_srange=%f",cpl_table_get_column_max(lambdas,"WAVE")); sinfo_msg("max_mrange=%f",cpl_table_get_column_max(mrange,"INDEX")); */ sinfo_msg("Estimate noise"); ck0(sinfo_object_estimate_noise(obj_frm,obj_noise_fit,¢re,&noise), "Estimating noise"); sinfo_msg("Background=%f Noise=%f",centre,noise); sinfo_msg("Flag object low_levels"); ck0(sinfo_object_flag_low_values(obj_frm,centre,noise,&fdata), "Flagging low pix"); //cpl_imagelist_save(fdata,"out_fdata.fits", // CPL_BPP_IEEE_FLOAT,NULL,CPL_IO_DEFAULT); sinfo_msg("Flag sky pixels"); ck0(sinfo_object_flag_sky_pixels(obj_frm,lambda,mrange,fdata,dispersion, &gpix,&ratio,&ima_sky), "Flagging sky pixels"); //cpl_image_save(gpix,"out_gpix.fits",CPL_BPP_IEEE_FLOAT, // NULL,CPL_IO_DEFAULT); //cpl_image_save(ratio,"out_ratio.fits",CPL_BPP_IEEE_FLOAT, // NULL,CPL_IO_DEFAULT); //cpl_image_save(ima_sky,"out_ima_sky.fits", // CPL_BPP_IEEE_FLOAT,NULL,CPL_IO_DEFAULT); // choose pixels which seems to be good sky pixels // (95% spectral pixels are flagged) sinfo_msg("Choose good sky (with > 95%% good spectral) pixels"); ck0(sinfo_choose_good_sky_pixels(obj_frm,ratio,gpix,min_frac,&mask), "Choosing good sky pixels"); //cpl_image_save(mask,"out_mask.fits",CPL_BPP_IEEE_FLOAT,NULL,CPL_IO_DEFAULT); // threshold ratio for fraction 'minfract' of spatial pixels to be 'sky' //sinfo_msg("To do: flag_threshold_sky_pixels"); // sum spectra of flagged pixels in object and sky frames sinfo_msg("Sum obj and sky spectra"); ck0(sinfo_sum_spectra(obj_frm,sky_frm,mask,lambda,llx,lly,urx,ury,int_obj, &int_sky),"summing obj & sky spectra"); //check_nomsg(cpl_table_save(*int_obj, NULL, NULL, "out_int_obj.fits", // CPL_IO_DEFAULT)); //check_nomsg(cpl_table_save(int_sky, NULL, NULL, "out_int_sky.fits", //CPL_IO_DEFAULT)); // Computes thermal background sinfo_msg("Computes thermal background"); /* ck0(sinfo_thermal_background2(int_sky,lambda,lrange,&bkg), "getting termal bkg"); */ ck0(sinfo_thermal_background(int_sky,lambda,lrange,fit_temp,niter, filter_width,dispersion,&bkg,&th_fit), "getting termal bkg"); check_nomsg(cpl_table_duplicate_column(*int_obj,"INT_SKY_ORG",int_sky,"INT")); check_nomsg(cpl_table_duplicate_column(*int_obj,"INT_BKG_FIT",bkg,"INT2")); check_nomsg(cpl_table_duplicate_column(*int_obj,"INT_BKG_SMO",int_sky, "INT_BKG_SMO")); sqc->th_fit=th_fit; //check_nomsg(cpl_table_save(bkg,NULL,NULL,"out_thermal_background.fits", //CPL_IO_DEFAULT)); /* ck0(sinfo_pro_save_tbl(bkg,set,set,"out_thermal_background.fits", "THERMAL_BACKGROUND",NULL,cpl_func,config), "Error saving %s","THERMAL_BACKGROUND"); */ sinfo_msg("Subtracts thermal background from integrated OH spectrum"); //sinfo_msg("nrow=%d %d",cpl_table_get_nrow(int_sky), // cpl_table_get_nrow(bkg)); check_nomsg(cpl_table_duplicate_column(int_sky,"BKG",bkg,"INT2")); check_nomsg(cpl_table_duplicate_column(int_sky,"INT0",int_sky,"INT")); check_nomsg(cpl_table_subtract_columns(int_sky,"INT","BKG")); //check_nomsg(cpl_table_duplicate_column(int_obj,"INT", // int_obj,"INT_OBJ_COR")); if(sub_thr_bkg == 1) { check_nomsg(cpl_table_duplicate_column(*int_obj,"THR_BKG",bkg,"INT2")); check_nomsg(cpl_table_subtract_columns(*int_obj,"INT","THR_BKG")); } //check_nomsg(cpl_table_save(*int_obj, NULL, NULL, "out_int_obj.fits", //CPL_IO_DEFAULT)); //check_nomsg(cpl_table_save(int_sky, NULL, NULL, "out_int_sky.fits", //CPL_IO_DEFAULT)); //check_nomsg(cpl_table_erase_column(int_sky,"BKG")); //check_nomsg(cpl_table_save(grange, NULL, NULL, "out_grange_6.fits", //CPL_IO_DEFAULT)); //check_nomsg(cpl_table_save(int_sky, NULL, NULL, "out_int_sky_sub.fits", //CPL_IO_DEFAULT)); // check for screw values at ends of spectrum sinfo_msg("Checks for screw values at ends of spectrum"); sinfo_check_screw_values(int_obj,&int_sky,grange,dispersion); //check_nomsg(cpl_table_save(grange, NULL, NULL, "out_grange_7.fits", //CPL_IO_DEFAULT)); //check_nomsg(cpl_table_save(*int_obj, NULL, NULL, "out_int_obj_chk.fits", //CPL_IO_DEFAULT)); //check_nomsg(cpl_table_save(int_sky, NULL, NULL, "out_int_sky_chk.fits", //CPL_IO_DEFAULT)); if(sky_interp_sw == 1) { sinfo_msg("Interpolate sky if necessary"); sinfo_interpolate_sky(int_sky,lambdas); } sinfo_msg("Crosscorrelate obj & sky to check for lambda offset"); //check_nomsg(cpl_table_save(*int_obj, NULL, NULL, "out_int_obj_chk.fits", //CPL_IO_DEFAULT)); //check_nomsg(cpl_table_save(int_sky, NULL, NULL, "out_int_sky_chk.fits", //CPL_IO_DEFAULT)); check_nomsg(wshift=sinfo_xcorr(*int_obj,int_sky,lambda,dispersion,line_hw)); //wshift=-1.7164495*dispersion; sinfo_msg("Dispersion=%f",dispersion); if(pshift == 0) { pshift=wshift/dispersion; } sinfo_msg("Shift sky of %f pixels toward object",pshift); check_nomsg(sinfo_shift_sky(&sky_frm,&int_sky,pshift)); //check_nomsg(cpl_table_save(*int_obj, NULL, NULL, "out_pip3_int_obj.fits", //CPL_IO_DEFAULT)); //check_nomsg(cpl_table_save(int_sky, NULL, NULL, "out_pip3_int_sky.fits", //CPL_IO_DEFAULT)); //DEBUG sinfo_msg("Optimise sky subtraction"); check_nomsg(sinfo_optimise_sky_sub(dispersion,line_hw,method,do_rot, lrange,lambda, lr41,lr52,lr63,lr74,lr02,lr85, lr20,lr31,lr42,lr53,lr64,lr75, lr86,lr97,lr00,int_obj,&int_sky, &rscale)); //check_nomsg(cpl_table_save(*int_obj, NULL, NULL, "out_int_obj.fits", //CPL_IO_DEFAULT)); //check_nomsg(cpl_table_save(int_sky, NULL, NULL, "out_int_sky.fits", //CPL_IO_DEFAULT)); sinfo_msg("Apply same scaling to whole cubes"); cknull_nomsg(obj_cub=cpl_imagelist_load(cpl_frame_get_filename(obj_frm), CPL_TYPE_DOUBLE,0)); cknull_nomsg(sky_cub=cpl_imagelist_load(cpl_frame_get_filename(sky_frm), CPL_TYPE_DOUBLE,0)); if(sub_thr_bkg == 1) { ck0_nomsg(sinfo_sub_thr_bkg_from_obj_cube(obj_cub,int_sky,obj_cor)); } else { check_nomsg(*obj_cor=cpl_imagelist_duplicate(obj_cub)); } ck0_nomsg(sinfo_scales_obj_sky_cubes(*obj_cor,sky_cub,bkg,rscale,obj_cor)); check_nomsg(cpl_table_name_column(*int_obj,"INT","INT_OBJ_ORG")); check_nomsg(cpl_table_name_column(*int_obj,"INTC","INT_OBJ_COR")); check_nomsg(cpl_table_name_column(*int_obj,"SKYC","INT_SKY_COR")); cleanup: sinfo_free_table(&rscale); sinfo_free_imagelist(&fdata); sinfo_free_table(&bkg); sinfo_free_table(&lambda); sinfo_free_table(&lrange); sinfo_free_table(&mrange); sinfo_free_table(&grange); sinfo_free_table(&lambdas); sinfo_free_image(&mask); sinfo_free_table(&lr41); sinfo_free_table(&lr52); sinfo_free_table(&lr63); sinfo_free_table(&lr74); sinfo_free_table(&lr02); sinfo_free_table(&lr85); sinfo_free_table(&lr20); sinfo_free_table(&lr31); sinfo_free_table(&lr42); sinfo_free_table(&lr53); sinfo_free_table(&lr64); sinfo_free_table(&lr75); sinfo_free_table(&lr86); sinfo_free_table(&lr97); sinfo_free_table(&lr00); sinfo_free_image(&gpix); sinfo_free_image(&ratio); sinfo_free_image(&ima_sky); //sinfo_free_table(&int_obj); sinfo_free_table(&int_sky); sinfo_free_imagelist(&obj_cub); sinfo_free_imagelist(&sky_cub); if (cpl_error_get_code() != CPL_ERROR_NONE) { return -1; } else { return 0; } } /*-------------------------------------------------------------------------*/ /** @name sinfo_set_ranges @memo defines several wavelength ranges @param obj_frm input object frame @param int_sky input sky spectra table @param obj_cor output object cube corrected for thermal background @return int 0, -1 @doc Returns in case of succes -1 else. */ /*--------------------------------------------------------------------------*/ static int sinfo_sub_thr_bkg_from_obj_cube(cpl_imagelist* obj_cub, cpl_table* int_sky, cpl_imagelist** obj_cor) { double* pthr_bkg=NULL; int zsz=0; int k=0; cpl_image* imgo=NULL; check_nomsg(pthr_bkg=cpl_table_get_data_double(int_sky,"BKG")); check_nomsg(zsz=cpl_imagelist_get_size(obj_cub)); check_nomsg(*obj_cor=cpl_imagelist_duplicate(obj_cub)); for(k=0;k 0) || (zsize_obj != zsize_sky)) { sinfo_msg("We have to interpolate sky frame - this is not good"); *sky_interp_sw=1; } return 0; cleanup: sinfo_free_table(&add1); sinfo_free_table(&tmp_tbl); sinfo_free_propertylist(&plist); return -1; } /*-------------------------------------------------------------------------*/ /** @name sinfo_table_dindgen @memo Fill a table column with values from 0 to table_size-1 @param t table to be filled @param label name of table column @return int 0, -1 @doc Returns in case of succes -1 else. */ /*--------------------------------------------------------------------------*/ static int sinfo_table_column_dindgen(cpl_table** t, const char* label) { int sz=0; int i=0; cknull(*t,"Null input vector"); check(sz=cpl_table_get_nrow(*t),"Getting size of a vector"); for(i=0;i 0.5); //TO BE REMOVED: loop_tbl is not used cknull_nomsg(loop_tbl=sinfo_image2table(mask)); check_nomsg(cpl_table_and_selected_double(loop_tbl,"VALUE", CPL_GREATER_THAN,0.5)); check_nomsg(loop=cpl_table_extract_selected(loop_tbl)); sinfo_free_table(&loop_tbl); sinfo_free_table(&loop); //Determines object spectrum for (i=0;i 0.5 && finite(obj_slice),pos_i); pos_i=sinfo_cnt_mask_thresh_and_obj_finite(mask,0.5,obj_slice); if (pos_i >= 1) { if ((pos_i) < 3 ) { //int_obj[i] = mean(obj_slice[pos]); //TODO here obj_slice should be considered only on pos: // one should do a selection/thresholding check_nomsg(cpl_table_set_double(*int_obj,"INT",i, cpl_image_get_mean_window(obj_slice, llx,lly, urx,ury))); } else { // select only poisitions where mask>0.5 and obj is finite // gslice = obj_slice[pos]; //sinfo_msg("obj pos_i=%d",pos_i); check_nomsg(gslice = cpl_image_duplicate(obj_slice)); check_nomsg(sinfo_image_flag_nan(&gslice)); /* sinfo_msg("obj: min=%f max=%f", cpl_image_get_min(obj_slice), cpl_image_get_max(obj_slice)); */ //check_nomsg(cpl_image_threshold(gslice,SINFO_DBL_MIN,3.0e6,1,0)); //check_nomsg(cpl_image_multiply(gslice,mask)); if(cpl_image_count_rejected(gslice) < 2048) { //2048=64*64/2 check_nomsg(med = cpl_image_get_median_window(gslice,llx,lly,urx,ury)); check_nomsg(sdv = cpl_image_get_stdev_window(gslice,llx,lly,urx,ury)); //sinfo_msg("med=%f sdv=%f",med,sdv); //avg = mean(gslice[where(gslice < med+3*sdv && gslice > med-3*sdv)]); check_nomsg(cpl_image_threshold(gslice,med-3*sdv,med+3*sdv,0,0)); check_nomsg(avg= cpl_image_get_mean_window(gslice,llx,lly,urx,ury)); check_nomsg(cpl_table_set_double(*int_obj,"INT",i,avg)); } else { check_nomsg(cpl_table_set_invalid(*int_obj,"INT",i)); } sinfo_free_image(&gslice); //sinfo_msg("sky int=%f",avg); } } //Determines sky spectrum check_nomsg(sky_slice = cpl_imagelist_get(sky,i)); //pos = where(mask > 0.5 and finite(sky_slice),pos_i); pos_i=sinfo_cnt_mask_thresh_and_obj_finite(mask,0.5,sky_slice); if (pos_i >= 1) { if ((pos_i) < 3) { //int_obj[i] = mean(obj_slice[pos]); //TODO here obj_slice should be considered only on pos: // one should do a selection/thresholding check_nomsg(cpl_table_set_double(*int_sky,"INT",i, cpl_image_get_mean_window(sky_slice, llx,lly, urx,ury))); } else { //sinfo_msg("pos_i=%d",pos_i); // select only poisitions where mask>0.5 and obj is finite // gslice = obj_slice[pos]; check_nomsg(gslice = cpl_image_duplicate(sky_slice)); check_nomsg(sinfo_image_flag_nan(&gslice)); //check_nomsg(cpl_image_threshold(gslice,SINFO_DBL_MIN,3.0e6,1,0)); //check_nomsg(cpl_image_multiply(gslice,mask)); if(cpl_image_count_rejected(gslice) < 2048) { //2048=64*64/2 check_nomsg(med = cpl_image_get_median_window(gslice,llx,lly,urx,ury)); check_nomsg(sdv = cpl_image_get_stdev_window(gslice,llx,lly,urx,ury)); //avg = mean(gslice[where(gslice < med+3*sdv && gslice > med-3*sdv)]); check_nomsg(cpl_image_threshold(gslice,med-3*sdv,med+3*sdv,0,0)); check_nomsg(avg= cpl_image_get_mean_window(gslice,llx,lly,urx,ury)); check_nomsg(cpl_table_set_double(*int_sky,"INT",i,avg)); } else { check_nomsg(cpl_table_set_invalid(*int_sky,"INT",i)); } sinfo_free_image(&gslice); /* if(i<100) { sinfo_msg("sky: wave=%f int=%f", cpl_table_get_double(*int_sky,"WAVE",i,&status),avg); } */ } } } sinfo_free_imagelist(&obj); sinfo_free_imagelist(&sky); return 0; cleanup: sinfo_free_image(&gslice); sinfo_free_image(&pos_tmp); sinfo_free_image(&msk_tmp); sinfo_free_table(&tmp_tbl); sinfo_free_table(&opos_tbl); sinfo_free_table(&spos_tbl); sinfo_free_table(&loop_tbl); sinfo_free_table(&loop); sinfo_free_imagelist(&obj); sinfo_free_imagelist(&sky); return -1; } /*-------------------------------------------------------------------------*/ /** @name sinfo_cnt_mask_thresh_and_obj_finite @memo count pixel values satisfying conditions @param mask input mask frame @param t input mask min threshold @param obj input mask frame @return No of points satisfying condition or -1 @doc */ /*--------------------------------------------------------------------------*/ static int sinfo_cnt_mask_thresh_and_obj_finite(const cpl_image* mask, const double t, const cpl_image* obj) { int cnt=0; int sxm=0; int sym=0; int sxo=0; int syo=0; int i=0; const double* pm=NULL; const double* po=NULL; check_nomsg(sxm=cpl_image_get_size_x(mask)); check_nomsg(sym=cpl_image_get_size_y(mask)); check_nomsg(sxo=cpl_image_get_size_x(obj)); check_nomsg(syo=cpl_image_get_size_y(obj)); if( sxm != sxo || sym != syo) { goto cleanup; } check_nomsg(pm=cpl_image_get_data_double_const(mask)); check_nomsg(po=cpl_image_get_data_double_const(obj)); for(i=0;i t) && (!irplib_isnan(po[i]))) { cnt++; } } return cnt; cleanup: return -1; } /*-------------------------------------------------------------------------*/ /** @name sinfo_flag_nan @memo reject as bad pix image NANs @param im input image frame @return if success The number of bad bixels else -1 */ /*--------------------------------------------------------------------------*/ static int sinfo_image_flag_nan(cpl_image** im) { int cnt=0; int sx=0; int sy=0; int i=0; int j=0; double* pi=NULL; check_nomsg(sx=cpl_image_get_size_x(*im)); check_nomsg(sy=cpl_image_get_size_y(*im)); check_nomsg(pi=cpl_image_get_data_double(*im)); for(j=0;j 0) { min_x = min_pos-HISTO_X_LEFT_CUT*(max_pos-min_pos); max_x = max_pos+HISTO_X_RIGHT_CUT*(max_pos-min_pos); } //sinfo_msg("min_x=%d max_x=%d",min_x,max_x); check_nomsg(hmin=sinfo_table_column_interpolate(histo,"HL",min_x)); check_nomsg(hmax=sinfo_table_column_interpolate(histo,"HL",max_x)); //sinfo_msg("hmin=%f hmax=%f min_xi_sz=%d",hmin,hmax,min_xi_sz); //cpl_table_save(histo, NULL, NULL, "out_histo.fits", CPL_IO_DEFAULT); sinfo_free_table(&histo); ck0(sinfo_histogram(data_tbl,nbins,hmin,hmax,&histo),"building histogram"); //cpl_table_save(histo, NULL, NULL, "out_histo1.fits", CPL_IO_DEFAULT); check_nomsg(cpl_table_add_scalar(histo,"HL",(hmax-hmin)/nbins/2)); //cpl_table_save(histo, NULL, NULL, "out_histo2.fits", CPL_IO_DEFAULT); } sinfo_free_table(&data_tbl); sinfo_free_table(&min_xi); //cpl_table_save(histo, NULL, NULL, "out_histo.fits", CPL_IO_DEFAULT); check_nomsg(peak=cpl_table_get_column_max(histo,"HY")); //sinfo_msg("peak=%f",peak); sinfo_free_table(&tmp_tbl1); check_nomsg(tmp_tbl1=sinfo_extract_table_rows(histo,"HY",CPL_EQUAL_TO,peak)); //cpl_table_save(tmp_tbl1, NULL, NULL, "out_tmp_tbl1.fits", CPL_IO_DEFAULT); check_nomsg(*centre=cpl_table_get_column_mean(tmp_tbl1,"HL")); //sinfo_msg("Background level=%f",*centre); sinfo_free_table(&tmp_tbl1); check_nomsg(tmp_tbl1=sinfo_extract_table_rows(histo,"HY", CPL_GREATER_THAN, peak/HISTO_Y_CUT)); sinfo_free_table(&tmp_tbl2); check_nomsg(tmp_tbl2=sinfo_extract_table_rows(tmp_tbl1,"HY", CPL_LESS_THAN,peak)); sinfo_free_table(&tmp_tbl1); check_nomsg(tempc=*centre-cpl_table_get_column_min(tmp_tbl2,"HL")); //sinfo_msg("min HX %f",cpl_table_get_column_min(tmp_tbl2,"HL")); sinfo_free_table(&tmp_tbl2); //sinfo_msg("Tempc=%f",tempc); check_nomsg(dist=sinfo_where_tab_min_max(histo,"HL", CPL_GREATER_THAN,*centre-HISTO_DIST_TEMPC_MIN_FCT*tempc, CPL_NOT_GREATER_THAN,*centre+HISTO_DIST_TEMPC_MAX_FCT*tempc)); offset=cpl_table_get_column_min(histo,"HY"); sinfo_free_table(&histo); check_nomsg(ndist=cpl_table_get_nrow(dist)); check_nomsg(cpl_table_cast_column(dist,"HY","HYdouble",CPL_TYPE_DOUBLE)); check_nomsg(disth=cpl_table_get_data_double(dist,"HYdouble")); check_nomsg(distx=cpl_table_get_data_double(dist,"HL")); //cpl_table_save(dist, NULL, NULL, "out_dist.fits", CPL_IO_DEFAULT); //TODO //gaussfit(distx,disty,dista,nterms=3); //*noise=dista[2]; *noise=tempc/2; /* THIS DOES NOT WORK */ //sinfo_msg("FWHM/2=%f",*noise); if(obj_noise_fit == 1) { check_nomsg(vy=cpl_vector_wrap(ndist,disth)); check_nomsg(vx=cpl_vector_wrap(ndist,distx)); check_nomsg(sx=cpl_vector_new(ndist)); check_nomsg(cpl_vector_fill(sx,1.)); check_nomsg(sy=cpl_vector_duplicate(sx)); x0=*centre; sigma=tempc/2; check_nomsg(cpl_vector_fit_gaussian(vx,NULL, vy,NULL, CPL_FIT_ALL, &x0,&sigma,&area,&offset, NULL,NULL,NULL)); //sinfo_msg("Gauss fit parameters:" // "x0=%f sigma=%f area=%f offset=%f mse=%f chired=%f", // x0,sigma,area,offset,mse,chired); //sinfo_msg("Background level=%f",*centre); //sinfo_msg("Noise=%f",sigma); *noise=sigma; sinfo_unwrap_vector(&vx); sinfo_unwrap_vector(&vy); sinfo_free_my_vector(&sx); sinfo_free_my_vector(&sy); } sinfo_free_table(&dist); //*noise=18.7448; //*noise=20.585946; return 0; cleanup: sinfo_free_imagelist(&obj_cub); sinfo_free_propertylist(&plist); sinfo_free_table(&min_xi); sinfo_free_table(&tmp_tbl1); sinfo_free_table(&tmp_tbl2); sinfo_free_table(&histo); sinfo_free_table(&dist); sinfo_free_table(&data_tbl); sinfo_free_my_vector(&sx); sinfo_free_my_vector(&sy); sinfo_unwrap_vector(&vx); sinfo_unwrap_vector(&vy); return -1; } /** @name sinfo_where_tab_min_max @brief finds table rows comprised between a min and a max @param t input table @param col column involved in selection @param op1 first operation of selection @param v1 first threshold value @param op2 second operation of selection @param v2 second threshold value @return new allocated sub table satisfying conditions */ cpl_table* sinfo_where_tab_min_max(cpl_table* t, const char* col, cpl_table_select_operator op1, const double v1, cpl_table_select_operator op2, const double v2) { cpl_table* res=NULL; cpl_table* tmp=NULL; check_nomsg(cpl_table_and_selected_double(t,col,op1,v1)); check_nomsg(tmp=cpl_table_extract_selected(t)); check_nomsg(cpl_table_and_selected_double(tmp,col,op2,v2)); check_nomsg(res=cpl_table_extract_selected(tmp)); check_nomsg(cpl_table_select_all(t)); sinfo_free_table(&tmp); return res; cleanup: sinfo_free_table(&tmp); sinfo_free_table(&res); return NULL; } /*-------------------------------------------------------------------------*/ /** @name sinfo_histogram @memo computes histogram @param data input data table @param nbins input size of histogram bin @param min input min data to be considered in histogram @param max input max data to be considered in histogram @param hist output histogram: table with results: # H: histogram # X: array of data values corresponding to the center of each bin # Xmean: 1D array corresponding to the mean of the data values entering each histogram bin Returns in case of succes -1 else. @doc Compute the histogram with the IDL intrinsic function HISTOGRAM, using the input options specified by the parameters Dmin, Dmax, Bin. All the computations are performed in floating-point arithmetics. Then compute arrays of values corresponding to each histogram bin, useful for plots, fitting, etc. */ /*--------------------------------------------------------------------------*/ int sinfo_histogram(const cpl_table* data, const int nbins, const double min, const double max, cpl_table** histo) { cpl_table* tmp_tbl1=NULL; cpl_table* tmp_tbl2=NULL; cpl_table* dat=NULL; int ntot=0; int i=0; int* phy=NULL; double* pdt=NULL; /* double* phx=NULL; */ double vtmp=0; double vstp=0; double vmax=0; double vmin=0; int h=0; check_nomsg(dat=cpl_table_duplicate(data)); check_nomsg(cpl_table_cast_column(dat,"DATA","DDATA",CPL_TYPE_DOUBLE)); /* sinfo_msg("min=%f max=%f", cpl_table_get_column_min(dat,"DDATA"), cpl_table_get_column_max(dat,"DDATA")); */ check_nomsg(cpl_table_and_selected_double(dat,"DDATA", CPL_NOT_GREATER_THAN,max)); /* check_nomsg(cpl_table_and_selected_double(dat,"DDATA",CPL_LESS_THAN,max)); */ check_nomsg(tmp_tbl1=cpl_table_extract_selected(dat)); /* sinfo_msg("min=%f max=%f", cpl_table_get_column_min(tmp_tbl1,"DDATA"), cpl_table_get_column_max(tmp_tbl1,"DDATA")); */ /* check_nomsg(cpl_table_and_selected_double(tmp_tbl1,"DDATA", CPL_NOT_LESS_THAN,min)); */ check_nomsg(cpl_table_and_selected_double(tmp_tbl1,"DDATA", CPL_GREATER_THAN,min)); check_nomsg(tmp_tbl2=cpl_table_extract_selected(tmp_tbl1)); /* sinfo_msg("min=%f max=%f", cpl_table_get_column_min(tmp_tbl2,"DDATA"), cpl_table_get_column_max(tmp_tbl2,"DDATA")); */ sinfo_free_table(&tmp_tbl1); /* check_nomsg(tmp_tbl1=sinfo_extract_table_rows(dat,"DDATA", CPL_NOT_GREATER_THAN,max)); check_nomsg(tmp_tbl2=sinfo_extract_table_rows(tmp_tbl1,"DDATA", CPL_NOT_LESS_THAN,min)); */ check_nomsg(ntot=cpl_table_get_nrow(tmp_tbl2)); /* not necessry to sort: check_nomsg(sinfo_sort_table_1(tmp_tbl2,"DDATA",FALSE));*/ check_nomsg(vmin=cpl_table_get_column_min(tmp_tbl2,"DDATA")); check_nomsg(vmax=cpl_table_get_column_max(tmp_tbl2,"DDATA")); vstp=(vmax-vmin)/(nbins-1); /* sinfo_msg("vmin=%f vmax=%f step=%f",vmin,vmax,vstp); */ check_nomsg(*histo=cpl_table_new(nbins)); check_nomsg(cpl_table_new_column(*histo,"HX",CPL_TYPE_DOUBLE)); check_nomsg(cpl_table_new_column(*histo,"HL",CPL_TYPE_DOUBLE)); check_nomsg(cpl_table_new_column(*histo,"HY",CPL_TYPE_INT)); /*check_nomsg(cpl_table_fill_column_window(*histo,"HX",0,nbins,0)); */ check_nomsg(cpl_table_fill_column_window(*histo,"HL",0,nbins,0)); check_nomsg(cpl_table_fill_column_window(*histo,"HY",0,nbins,0)); check_nomsg(phy=cpl_table_get_data_int(*histo,"HY")); /*check_nomsg(phx=cpl_table_get_data_double(*histo,"HX")); */ check_nomsg(pdt=cpl_table_get_data_double(dat,"DATA")); for(i=0;i-1)) { phy[h]++; } } //cpl_table_save(*histo, NULL, NULL, "out_histo_p5.fits", CPL_IO_DEFAULT); sinfo_free_table(&tmp_tbl2); sinfo_free_table(&dat); return 0; cleanup: sinfo_free_table(&tmp_tbl1); sinfo_free_table(&tmp_tbl2); sinfo_free_table(&dat); return -1; } /*-------------------------------------------------------------------------*/ /** @name sinfo_object_flag_low_values @memo flags low values in object frame: val>=cnt+2*sig @param obj_frm input object frame @param cnt input value involved in thresold @param sig input value involved in threshold @param flag_data output imagelist @return 0 in case of succes -1 else. */ /*--------------------------------------------------------------------------*/ int sinfo_object_flag_low_values(cpl_frame* obj_frm, const double cnt, const double sig, cpl_imagelist** flag_data) { int xsz=0; int ysz=0; int zsz=0; int n=0; int i=0; int k=0; int r=0; cpl_propertylist* plist=NULL; cpl_table* data_tbl=NULL; cpl_table* flag_tbl=NULL; cpl_image* img=NULL; cpl_imagelist* obj_cub=NULL; double* pdata=NULL; double* pflag=NULL; /* Get Object relevant information */ cknull_nomsg(plist=cpl_propertylist_load(cpl_frame_get_filename(obj_frm),0)); check_nomsg(xsz=sinfo_pfits_get_naxis1(plist)); check_nomsg(ysz=sinfo_pfits_get_naxis2(plist)); check_nomsg(zsz=sinfo_pfits_get_naxis3(plist)); sinfo_free_propertylist(&plist); cknull_nomsg(obj_cub=cpl_imagelist_load(cpl_frame_get_filename(obj_frm), CPL_TYPE_DOUBLE,0)); n=xsz*ysz*zsz; check_nomsg(data_tbl=cpl_table_new(n)); check_nomsg(cpl_table_new_column(data_tbl,"DATA",CPL_TYPE_DOUBLE)); for(k=0;k 0) { // build two arrays of proper size all_pix=(double)gpix_i; /* sinfo_msg("flagspec: min=%f max=%f", cpl_table_get_column_min(fspec,"VALUE"), cpl_table_get_column_max(fspec,"VALUE")); sinfo_msg("good flagspec: min=%f max=%f", cpl_table_get_column_min(gspec,"VALUE"), cpl_table_get_column_max(gspec,"VALUE")); sinfo_msg("nfspec=%d",cpl_table_get_nrow(fspec)); check_nomsg(cpl_table_save(fspec, NULL, NULL, "out_fspec.fits",CPL_IO_DEFAULT)); check_nomsg(cpl_table_save(gspec, NULL, NULL, "out_gspec.fits", CPL_IO_DEFAULT)); */ //check_nomsg(flag_pix=cpl_table_and_selected_double(fspec,"VALUE", // CPL_GREATER_THAN,0.5)); //sinfo_msg("all_pix=%f flag_pix=%f",all_pix,flag_pix); check_nomsg(flag_pix=cpl_table_and_selected_double(gspec,"VALUE", CPL_GREATER_THAN,0.5)); //sinfo_msg("all_pix=%f flag_pix=%f",all_pix,flag_pix); // flag_pix = float(n_elements(where(fspec[gpix] > 0.5))); // compute the ratio between the two arrays ratio=flag_pix/all_pix; // considers only pixels with have at least half good pixels if(all_pix > cpl_table_get_nrow(mrange)/2) { pg_img[i+j*xsz]=1; pr_img[i+j*xsz]=ratio; } //mean(ospec(gpix)) check_nomsg(pimage[i+j*xsz]=cpl_table_get_column_mean(gpix,"VALUE")); //sinfo_msg("ix=%d iy=%d r=%f",i,j,ratio); } sinfo_free_table(&ospec); sinfo_free_table(&fspec); sinfo_free_table(&gpix); sinfo_free_table(&gspec); } /* end for over i */ } /* end for over j */ sinfo_free_imagelist(&osel); sinfo_free_imagelist(&fsel); /* cpl_image_save(*r_img,"out_r_img.fits",CPL_BPP_IEEE_FLOAT, NULL,CPL_IO_DEFAULT); cpl_image_save(*g_img,"out_g_img.fits",CPL_BPP_IEEE_FLOAT, NULL,CPL_IO_DEFAULT); cpl_image_save(*image,"out_image.fits",CPL_BPP_IEEE_FLOAT, NULL,CPL_IO_DEFAULT); */ // get total(g_arr) check_nomsg(tot=cpl_image_get_flux(*g_img)); if(tot < 1) { sinfo_msg_error("no good spaxel"); goto cleanup; } return 0; cleanup: sinfo_free_propertylist(&plist); sinfo_free_imagelist(&obj_cub); sinfo_free_imagelist(&osel); sinfo_free_imagelist(&fsel); sinfo_free_table(&ospec); sinfo_free_table(&fspec); sinfo_free_table(&gpix); sinfo_free_table(&gspec); return -1; } /** @name sinfo_choose_good_sky_pixels @param obj_frm input object frame @param r_img input image ratio good/all pixels @param g_img input image with good pixels @param min_frac input threshold setting the minimum fraction of good pixels @param mask output mask indicating pixels which satisfy the condition */ int sinfo_choose_good_sky_pixels(cpl_frame* obj_frm, cpl_image* r_img, cpl_image* g_img, const double min_frac, cpl_image** mask) { int xsz=0; int ysz=0; //int zsz=0; int r2i=0; int status=0; double tot=0; double thres=SKY_THRES; double cum_x_max=0; cpl_image* r2img=NULL; cpl_propertylist* plist=NULL; cpl_table* cum=NULL; cpl_table* hcum=NULL; cpl_table* thres_tbl=NULL; cknull_nomsg(plist=cpl_propertylist_load(cpl_frame_get_filename(obj_frm),0)); check_nomsg(xsz=sinfo_pfits_get_naxis1(plist)); check_nomsg(ysz=sinfo_pfits_get_naxis2(plist)); //check_nomsg(zsz=sinfo_pfits_get_naxis3(plist)); sinfo_free_propertylist(&plist); // choose pixels which seem to be sky (ie 95% of spectral pixels are flagged) check_nomsg(*mask=cpl_image_new(xsz,ysz,CPL_TYPE_DOUBLE)); //r2 = where(r_img >= thres,r2i); // count good pixels: set to 0 what < thres and to 1 what > thres check_nomsg(r2img=cpl_image_duplicate(r_img)); check_nomsg(cpl_image_threshold(r2img,thres,thres,0,1)); check_nomsg(r2i=cpl_image_get_flux(r2img)); if(r2i>0) { sinfo_free_image(&(*mask)); check_nomsg(*mask=cpl_image_duplicate(r2img)); } sinfo_free_image(&r2img); check_nomsg(r2img=cpl_image_duplicate(r_img)); check_nomsg(cpl_image_threshold(r2img,thres,SINFO_DBL_MAX,0,SINFO_DBL_MAX)); sinfo_free_image(&r2img); check_nomsg(tot=cpl_image_get_flux(g_img)); sinfo_msg("%2.2d spaxels (%4.1f %% of good pixels) are designated as sky", r2i,100.*r2i/tot); //threshold ratio for fraction 'minfrac' of spatial pixels to be 'sky' if (1.*r2i/tot < min_frac) { sinfo_msg("this is too small - will increase it to %4.1f %%", 100.*min_frac); check_nomsg(cum=cpl_table_new(xsz*ysz)); check_nomsg(cpl_table_new_column(cum,"X",CPL_TYPE_DOUBLE)); sinfo_table_column_dindgen(&cum,"X"); check_nomsg(cpl_table_add_scalar(cum,"X",1.)); //hcum = r_img(sort(r_img)); hcum = sinfo_image2table(r_img); check_nomsg(sinfo_sort_table_1(hcum,"VALUE",FALSE)); //thresh = hcum[where(xcum/max(xcum) >= 1.-min_frac)]; check_nomsg(cpl_table_duplicate_column(cum,"H",hcum,"VALUE")); check_nomsg(cum_x_max=cpl_table_get_column_max(cum,"X")); check_nomsg(cpl_table_duplicate_column(cum,"R",cum,"X")); check_nomsg(cpl_table_divide_scalar(cum,"R",cum_x_max)); check_nomsg(cpl_table_and_selected_double(cum,"R", CPL_NOT_LESS_THAN, (1.-min_frac))); check_nomsg(thres_tbl=cpl_table_extract_selected(cum)); check_nomsg(thres = cpl_table_get(thres_tbl,"R",0,&status)); //*mask[where(r_img >= thresh)] = 1; sinfo_free_image(&(*mask)); check_nomsg(*mask=cpl_image_duplicate(r_img)); check_nomsg(cpl_image_threshold(*mask,thres,thres,0,1)); } sinfo_free_table(&cum); sinfo_free_table(&hcum); sinfo_free_table(&thres_tbl); return 0; cleanup: sinfo_free_propertylist(&plist); sinfo_free_image(&r2img); sinfo_free_table(&cum); sinfo_free_table(&hcum); sinfo_free_table(&thres_tbl); return -1; } /*-------------------------------------------------------------------------*/ /** @name sinfo_fit_bkg @memo Computes chi2 of difference INT(sky)-thermal_background @param sa pointer to sinfo_amoeba structure @param p input fit parameters Returns chi2= sum_i[int(sky)_i-therma_i] */ /*--------------------------------------------------------------------------*/ static double sinfo_fit_bkg(double p[]) { double* px=NULL; double* py=NULL; double* pv=NULL; cpl_vector* vtmp=NULL; double max=0; int i=0; int np=0; double chi2=0; check_nomsg(px= cpl_vector_get_data(sa_vx)); check_nomsg(py= cpl_vector_get_data(sa_vy)); check_nomsg(np= cpl_vector_get_size(sa_vx)); check_nomsg(vtmp=cpl_vector_duplicate(sa_vy)); check_nomsg(pv=cpl_vector_get_data(vtmp)); for(i=0;i 0) { check_nomsg(cpl_vector_divide_scalar(vtmp,max)); check_nomsg(cpl_vector_multiply_scalar(vtmp,p[1])); check_nomsg(cpl_vector_add_scalar(vtmp,p[0])); } for(i=0;i0) { for(i=0;i pi[j]) { min=pi[j]; } } po[i]=min; } // To prevent border effects: for (i = 0; i < start; i++) { po[i] = po[start]; } for (i = end; i < length; i++) { po[i] = po[end-1]; } return vo; cleanup: return NULL; } /** @name sinfo_filter_max @memo apply a filter max of a given size to a vector @param vi input vector @param size size of filter @return a new allocated vector resulting from the filter proces @doc The following static function passes a max filter of given box size on the data buffer. The box size must be a positive odd integer. */ static cpl_vector* sinfo_filter_max(const cpl_vector* vi, const int size) { cpl_vector* vo=NULL; double max=0; int start=size/2; int end=0; int length=0; int i=0; int j=0; const double* pi=NULL; double* po=NULL; cknull(vi,"null input vector"); pi=cpl_vector_get_data_const(vi); length=cpl_vector_get_size(vi); end=length-size/2; vo=cpl_vector_new(length); po=cpl_vector_get_data(vo); for(i=start; i < end; i++) { max=pi[i-start]; for(j=i-start+1;j * * CPL_ERROR_NULL_INPUT * * The input spectrum is a NULL pointer. * * * * CPL_ERROR_ILLEGAL_INPUT * * Either msize is less than 3, or fsize is less * than msize, or fsize is greater than length/2. * * * * @enderror * * This function fills the array @em back with the estimated values * of the background along the input spectrum. The algorithm is * based on the assumption that there is at least one background * value at any position of the min-filter box running along the * spectrum. A min-filter is passed on the spectrum, and the result * is smoothed by averaging on a running box of size @em fsize. * The min-filter is run between the positions @em msize / 2 * and @em length - @em msize / 2, and the min value found at * such positions is then repeated up to the spectrum ends. Similarly, * the running average is limited to the interval from @em fsize / 2 * and @em length - @em fsize / 2, leaving the most external values * untouched. After this, a max filter and a smoothing using boxes * with double the specified sizes are run, as a way to eliminate * the contamination from occasional cold pixels. Finally, the * min filter and the smoothing are applied again to obviate the * slight background over-estimation introduced by the max filter. * * It is required that the @em back array is at least long as the * array @em spectrum. Moreover @em msize must be greater than 1, * and @em fsize greater than, or equal to, @em msize. Likewise, * @em length must be greater than twice @em fsize. If such conditions * are not met, or if the input arrays are @c NULL pointers, this * function will set an error code, and leave the @em back array * untouched. If either @em msize or @em fsize are even numbers, * they are made odd by adding 1. Suggested values for @em msize * and @em fsize are 15 pixels for typical arc lamp spectra. */ cpl_vector* sinfo_sky_background_estimate(cpl_vector *spectrum, int msize, int fsize) { cpl_vector * minf=NULL; cpl_vector * maxf=NULL; cpl_vector * smof=NULL; cpl_vector * back=NULL; double* pb=NULL; double* ps=NULL; int i=0; int length=0; cknull(spectrum,"null input data"); if (msize % 2 == 0) msize++; if (fsize % 2 == 0) fsize++; check_nomsg(length=cpl_vector_get_size(spectrum)); if (msize < 3 || fsize < msize || length < 2*fsize) return NULL; cknull_nomsg(minf = sinfo_filter_min(spectrum, msize)); cknull_nomsg(smof = sinfo_filter_smo(minf, fsize)); cpl_vector_delete(minf); cknull_nomsg(maxf = sinfo_filter_max(smof,2*msize+1)); cpl_vector_delete(smof); cknull_nomsg(smof = sinfo_filter_smo(maxf, 2*fsize+1)); cpl_vector_delete(maxf); cknull_nomsg(minf = sinfo_filter_min(smof, 2*msize+1)); cpl_vector_delete(smof); cknull_nomsg(smof = sinfo_filter_smo(minf, 2*fsize+1)); cpl_vector_delete(minf); cknull_nomsg(back=cpl_vector_new(length)); cknull_nomsg(pb=cpl_vector_get_data(back)); cknull_nomsg(ps=cpl_vector_get_data(smof)); for (i = 0; i < length; i++) { pb[i] = ps[i]; } cpl_vector_delete(smof); return back; cleanup: return NULL; } /** @name sinfo_slice_z @memo extract a cube slice (spectrum) along z at i,j @param cin input imagelist @param i inpout position (x) @param j inpout position (y) @return if success a table containg the spectrum along z at i,j else NULL */ static cpl_table* sinfo_slice_z(const cpl_imagelist* cin,const int i,const int j) { int sx=0; //int sy=0; int sz=0; int k=0; cpl_table* tout=NULL; const cpl_image* img=NULL; const double* pim=NULL; cknull(cin,"null input imagelist"); check_nomsg(sz=cpl_imagelist_get_size(cin)); check_nomsg(img=cpl_imagelist_get_const(cin,0)); check_nomsg(sx=cpl_image_get_size_x(img)); //check_nomsg(sy=cpl_image_get_size_y(img)); check_nomsg(tout=cpl_table_new(sz)); check_nomsg(cpl_table_new_column(tout,"VALUE",CPL_TYPE_DOUBLE)); for(k=0;k0 && pint[i-1]<0)) { ppos[i]=i; } } //check_nomsg(cpl_table_save(z_diff,NULL,NULL,"out_z_diff.fits", // CPL_IO_DEFAULT)); check_nomsg(cpl_table_select_all(z_diff)); check_nomsg(cpl_table_and_selected_int(z_diff,"POS",CPL_GREATER_THAN,0)); check_nomsg(z_pos=cpl_table_extract_selected(z_diff)); sinfo_free_table(&z_diff); //check_nomsg(cpl_table_save(z_pos,NULL,NULL, // "out_z_pos.fits",CPL_IO_DEFAULT)); //Do a gaussian fit in a range of size 2*zext centered at //each line maximum position (fit the line) to get in corresponding arrays: // 1) line lambda position of object and sky // 2) line object -sky intensity // 3) line object-sky intensity error /* g_lam = 0.; g_diff = 0.; g_err = 0.; */ check_nomsg(npos=cpl_table_get_nrow(z_pos)); z_ext = line_hw ; check_nomsg(cpl_table_new_column(z_pos,"STATUS_S",CPL_TYPE_INT)); check_nomsg(cpl_table_new_column(z_pos,"STATUS_O",CPL_TYPE_INT)); check_nomsg(cpl_table_fill_column_window_int(z_pos,"STATUS_S",0,npos,0)); check_nomsg(cpl_table_fill_column_window_int(z_pos,"STATUS_O",0,npos,0)); check_nomsg(cpl_table_new_column(z_pos,"SIGS",CPL_TYPE_DOUBLE)); check_nomsg(cpl_table_new_column(z_pos,"WAVES",CPL_TYPE_DOUBLE)); check_nomsg(cpl_table_new_column(z_pos,"BKGS",CPL_TYPE_DOUBLE)); check_nomsg(cpl_table_new_column(z_pos,"AREAS",CPL_TYPE_DOUBLE)); check_nomsg(cpl_table_new_column(z_pos,"AMPS",CPL_TYPE_DOUBLE)); check_nomsg(cpl_table_new_column(z_pos,"SIGO",CPL_TYPE_DOUBLE)); check_nomsg(cpl_table_new_column(z_pos,"WAVEO",CPL_TYPE_DOUBLE)); check_nomsg(cpl_table_new_column(z_pos,"BKGO",CPL_TYPE_DOUBLE)); check_nomsg(cpl_table_new_column(z_pos,"AREAO",CPL_TYPE_DOUBLE)); check_nomsg(cpl_table_new_column(z_pos,"AMPO",CPL_TYPE_DOUBLE)); check_nomsg(cpl_table_new_column(z_pos,"WAVEC",CPL_TYPE_DOUBLE)); check_nomsg(cpl_table_new_column(z_pos,"WDIF",CPL_TYPE_DOUBLE)); check_nomsg(cpl_table_new_column(z_pos,"ERR",CPL_TYPE_DOUBLE)); nfit=2*z_ext+1; //sinfo_msg("npos=%d z_ext=%d",npos,z_ext); //sinfo_table_column_dump(z_pos,"POS",CPL_TYPE_INT); //check_nomsg(cpl_table_save(z_pos,NULL,NULL, // "out_z_pos_0.fits",CPL_IO_DEFAULT)); //check_nomsg(cpl_table_save(int_obj,NULL,NULL, // "out_int_obj_0.fits",CPL_IO_DEFAULT)); //check_nomsg(cpl_table_save(int_sky,NULL,NULL, // "out_int_sky_0.fits",CPL_IO_DEFAULT)); for (jz=0;jz 0 && (z1+z_ext) < zsize) { check_nomsg(cpl_table_select_all(int_sky)); check_nomsg(cpl_table_select_all(int_obj)); check_nomsg(cpl_table_select_all(lambda)); check_nomsg(cpl_table_and_selected_window(int_sky,z1-z_ext,nfit)); check_nomsg(s_tbl=cpl_table_extract_selected(int_sky)); check_nomsg(cpl_table_and_selected_window(lambda,z1-z_ext,nfit)); check_nomsg(w_tbl=cpl_table_extract_selected(lambda)); check_nomsg(cpl_table_and_selected_window(int_obj,z1-z_ext,nfit)); check_nomsg(o_tbl=cpl_table_extract_selected(int_obj)); check_nomsg(vw=cpl_vector_wrap(nfit, cpl_table_get_data_double(w_tbl,"WAVE"))); check_nomsg(vs=cpl_vector_wrap(nfit, cpl_table_get_data_double(s_tbl,"INT"))); check_nomsg(vo=cpl_vector_wrap(nfit, cpl_table_get_data_double(o_tbl,"INT"))); check_nomsg(sx=cpl_vector_new(nfit)); check_nomsg(cpl_vector_fill(sx,10.)); check_nomsg(sy=cpl_vector_duplicate(sx)); // Check if the object line is in emission or absorbtion o1=cpl_vector_get(vo,0); o2=cpl_vector_get(vo,nfit-1); oc=(o1+o2)*0.5; om=cpl_vector_get_median_const(vo); if(om1) { check_nomsg(z_sdv = cpl_table_get_column_stdev(z_good,"WDIF")); } else { z_sdv=0; } sinfo_free_table(&z_good); check_nomsg(cpl_table_erase_column(z_pos,"CHECKW")); } /* do a poly fit of wdif versus wave*/ /* for (iq = 0; iq<3; iq++) { // sinfo_msg("%d %f %f",iq,mean(zfit),zsdv); par1 = poly_fit(g_lam[z_good],g_diff[z_good],poly_n); z_fit = g_diff*0.; for (ii=0;ii0) { check_nomsg(vx=cpl_vector_wrap(nfit, cpl_table_get_data_double(z_good,"WAVE"))); check_nomsg(vy=cpl_vector_wrap(nfit, cpl_table_get_data_double(z_good,"WDIF"))); check_nomsg(cfit=sinfo_polynomial_fit_1d_create(vx,vy,0,&mse)); pows[0]=0; pows[1]=0; check_nomsg(zfit=cpl_polynomial_get_coeff(cfit,pows)); sinfo_free_polynomial(&cfit); //sinfo_msg("coeff 0=%g um %g pix",zfit,zfit/dispersion); //computes residuals=difference-fit and their standard deviation //and then do a kappa-sigma clip of outliers (out of 3 sigma) check_nomsg(cpl_table_fill_column_window(z_good,"ZFIT",0,nfit,zfit)); check_nomsg(cpl_table_duplicate_column(z_good,"WRES",z_good,"WDIF")); check_nomsg(cpl_table_subtract_columns(z_good,"WRES","ZFIT")); if(nfit>1) { //sinfo_msg("nfit=%d",nfit); //cpl_table_dump(z_good,0,nfit,stdout); check_nomsg(z_sdv=cpl_table_get_column_stdev(z_good,"WRES")); //sinfo_msg("z_sdv=%f",z_sdv); } else { z_sdv=0; } check_nomsg(z_mean=cpl_table_get_column_mean(z_good,"WDIF")); sinfo_msg(" %d %3.2g %3.2g %5.4g %5.4g", iq,z_mean,z_sdv,z_mean/dispersion,z_sdv/dispersion); check_nomsg(nfit=cpl_table_get_nrow(z_pos)); check_nomsg(cpl_table_fill_column_window(z_pos,"ZFIT",0,nfit,zfit)); check_nomsg(cpl_table_duplicate_column(z_pos,"WRES",z_pos,"WDIF")); check_nomsg(cpl_table_subtract_columns(z_pos,"WRES","ZFIT")); check_nomsg(cpl_table_multiply_columns(z_pos,"WRES","WRES")); check_nomsg(cpl_table_power_column(z_pos,"WRES",0.5)); //cpl_table_dump(z_pos,0,cpl_table_get_nrow(z_pos),stdout); /* sinfo_msg("min=%g max=%g ndat=%d", cpl_table_get_column_min(z_pos,"WRES"), cpl_table_get_column_max(z_pos,"WRES"), cpl_table_get_nrow(z_pos)); */ check_nomsg(cpl_table_and_selected_double(z_pos,"WRES", CPL_NOT_GREATER_THAN,3*z_sdv)); check_nomsg(sinfo_free_table(&z_good)); check_nomsg(z_good=cpl_table_extract_selected(z_pos)); check_nomsg(cpl_table_select_all(z_pos)); check_nomsg(cpl_table_select_all(z_good)); check_nomsg(cpl_table_erase_column(z_good,"WRES")); check_nomsg(cpl_table_erase_column(z_pos,"WRES")); sinfo_unwrap_vector(&vx); sinfo_unwrap_vector(&vy); } } //sinfo_msg(">>mean=%g",cpl_table_get_column_mean(z_good,"WDIF")); //check_nomsg(cpl_table_save(z_good,NULL,NULL, // "out_z_pos3.fits",CPL_IO_DEFAULT)); sinfo_unwrap_vector(&vx); sinfo_unwrap_vector(&vy); sinfo_free_polynomial(&cfit); sinfo_free_table(&z); sinfo_free_table(&z_pos); sinfo_free_table(&z_good); return zfit; cleanup: sinfo_free_table(&z_good); sinfo_free_table(&z); sinfo_free_table(&z_diff); sinfo_free_table(&tmp_sky); sinfo_free_table(&z_pos); sinfo_unwrap_vector(&vw); sinfo_unwrap_vector(&vs); sinfo_unwrap_vector(&vo); sinfo_free_my_vector(&sx); sinfo_free_my_vector(&sy); sinfo_unwrap_vector(&vx); sinfo_unwrap_vector(&vy); sinfo_free_table(&w_tbl); sinfo_free_table(&s_tbl); sinfo_free_table(&o_tbl); sinfo_free_polynomial(&cfit); return 0; } /** @name sinfo_table_flag_nan @memo flag NAN values in a table column @param t input table @param label column name @param min minimum value @param min maximum value @return 0 if success, -1 else */ static int sinfo_table_set_nan_out_min_max(cpl_table** t, const char* c, const double min, const double max) { int sz=0; int i=0; double* pt=NULL; check_nomsg(sz=cpl_table_get_nrow(*t)); check_nomsg(pt=cpl_table_get_data_double(*t,c)); for(i=0;i max) { check_nomsg(cpl_table_set_invalid(*t ,c,i)); } } return 0; cleanup: return -1; } /** @name sinfo_table_flag_nan @memo flag NAN values in a table column @param t input table @param label column name @return 0 if success, -1 else */ static int sinfo_table_flag_nan(cpl_table** t,const char* label) { int sz=0; int i=0; double* pt=NULL; check_nomsg(sz=cpl_table_get_nrow(*t)); check_nomsg(pt=cpl_table_get_data_double(*t,label)); for(i=0;i 0) { new_sky= spline(lambda[z_good]-z_fit[z_good],old_sky[z_good],lambda); new_fin= where(finite(new_sky,/infinity) || finite(old_sky,/nan),newfin_i); if (new_fin_i > 0) new_sky[new_fin] = !values.f_nan; sky_out[ix,iy,*] = new_sky; } } } sky = sky_out; } cleanup: return; } */ /** @name sinfo_optimise_sky_sub @memo optimse sky subtraction @param wtol tolerance on wavelength for selection @param line_hw line half width @param method optimization method to determine ratio line(obj)/line(sky) @param lambda wavelength range @param lr41 4-1 transition wavelength range table @param lr52 5-2 transition wavelength range table @param lr63 6-3 transition wavelength range table @param lr74 7-4 transition wavelength range table @param lr02 0-2 transition wavelength range table @param lr85 8-5 transition wavelength range table @param lr20 2-0 transition wavelength range table @param lr31 3-1 transition wavelength range table @param lr42 4-2 transition wavelength range table @param lr53 5-3 transition wavelength range table @param lr64 6-4 transition wavelength range table @param lr75 7-5 transition wavelength range table @param lr86 8-6 transition wavelength range table @param lr97 9-7 transition wavelength range table @param lr00 0-0 transition wavelength range table @param int_obj input/output(corrected) object spectrum @param int_sky input/output(corrected) sky spectrum @param rscale output computed ratio table */ void sinfo_optimise_sky_sub(const double wtol, const double line_hw, const int method, const int do_rot, cpl_table* lrange, cpl_table* lambda, cpl_table* lr41, cpl_table* lr52, cpl_table* lr63, cpl_table* lr74, cpl_table* lr02, cpl_table* lr85, cpl_table* lr20, cpl_table* lr31, cpl_table* lr42, cpl_table* lr53, cpl_table* lr64, cpl_table* lr75, cpl_table* lr86, cpl_table* lr97, cpl_table* lr00, cpl_table** int_obj, cpl_table** int_sky, cpl_table** rscale) { int npixw=2*line_hw; //full width in pixels of unresolved emission line cpl_array* do_hk=NULL; cpl_array* rfit=NULL; int i=0; cpl_table* sky_lr=NULL; cpl_table* obj_lr=NULL; cpl_table* wav_lr=NULL; double sky_med=0; double sky_sdv=0; int lr41_i=0; int lr52_i=0; int lr63_i=0; int lr74_i=0; int lr02_i=0; int lr85_i=0; int lr20_i=0; int lr31_i=0; int lr42_i=0; int lr53_i=0; int lr64_i=0; int lr75_i=0; int lr86_i=0; int lr97_i=0; int lr00_i=0; int xxx1_i=0; int status=0; int finite_pix_i=0; double sky_thresh=0.; cpl_table* rat_sky=NULL; cpl_table* xxx1=NULL; cpl_table* xxx2=NULL; cpl_table* xxx1_sub=NULL; cpl_table* line_regions=NULL; cpl_table* cont_regions=NULL; int line_i=0; int cont_i=0; double fmed=0; double fsdv=0; cpl_table* fline_res=NULL; int fclip_i=0; int fline_i=0; cpl_table* rscale0=NULL; double r=0; cpl_table* obj_cont=NULL; cpl_table* sky_cont=NULL; cpl_table* obj_line=NULL; cpl_table* sky_line=NULL; //Rotational parameters int low_pos_i=0; int med_pos_i=0; int hi_pos_i=0; cpl_table* finite_pix=NULL; cpl_table* tmp_tbl=NULL; cpl_table* low_scale=NULL; cpl_table* med_scale=NULL; cpl_table* hi_regions=NULL; cpl_table* low_regions=NULL; cpl_table* med_regions=NULL; cpl_table* low_pos=NULL; cpl_table* med_pos=NULL; cpl_table* hi_pos=NULL; cpl_table* llr_xxx1=NULL; double rhi=0; double rmed=0; double rlow=0; double min_lrange=0; double max_lrange=0; int nrow=0; double w_rot_low[NROT]={1.00852,1.03757,1.09264,1.15388,1.22293, 1.30216,1.45190,1.52410,1.60308,1.69037, 1.78803,2.02758,2.18023,1.02895,1.08343, 1.14399,1.21226,1.29057,1.43444,1.50555, 1.58333,1.66924,1.76532,2.00082,2.15073}; double w_rot_med[NROT]={1.00282,1.02139,1.04212,1.07539,1.09753, 1.13542,1.15917,1.20309,1.22870,1.28070, 1.30853,1.41861,1.46048,1.48877,1.53324, 1.56550,1.61286,1.65024,1.70088,1.74500, 1.79940,1.97719,2.04127,2.12496,2.19956}; check_nomsg(do_hk = cpl_array_new(NBOUND+1,CPL_TYPE_INT)); check_nomsg(rfit = cpl_array_new(NBOUND+1,CPL_TYPE_DOUBLE)); lr41_i=cpl_table_get_nrow(lr41); lr52_i=cpl_table_get_nrow(lr52); lr63_i=cpl_table_get_nrow(lr63); lr74_i=cpl_table_get_nrow(lr74); lr02_i=cpl_table_get_nrow(lr02); lr85_i=cpl_table_get_nrow(lr85); lr20_i=cpl_table_get_nrow(lr20); lr31_i=cpl_table_get_nrow(lr31); lr42_i=cpl_table_get_nrow(lr42); lr53_i=cpl_table_get_nrow(lr53); lr64_i=cpl_table_get_nrow(lr64); lr75_i=cpl_table_get_nrow(lr75); lr86_i=cpl_table_get_nrow(lr86); lr97_i=cpl_table_get_nrow(lr97); check_nomsg(lr00_i=cpl_table_get_nrow(lr00)); cpl_array_set_int(do_hk,0,lr41_i); cpl_array_set_int(do_hk,1,lr52_i); cpl_array_set_int(do_hk,2,lr63_i); cpl_array_set_int(do_hk,3,lr74_i); cpl_array_set_int(do_hk,4,lr02_i); cpl_array_set_int(do_hk,5,lr85_i); cpl_array_set_int(do_hk,6,lr20_i); cpl_array_set_int(do_hk,7,lr31_i); cpl_array_set_int(do_hk,8,lr42_i); cpl_array_set_int(do_hk,9,lr53_i); cpl_array_set_int(do_hk,10,lr64_i); cpl_array_set_int(do_hk,11,lr75_i); cpl_array_set_int(do_hk,12,lr86_i); cpl_array_set_int(do_hk,13,lr97_i); check_nomsg(cpl_array_set_int(do_hk,14,lr00_i)); check_nomsg(rscale0=cpl_table_duplicate(*int_sky)); check_nomsg(cpl_table_new_column(rscale0,"RATIO",CPL_TYPE_DOUBLE)); check_nomsg(nrow=cpl_table_get_nrow(rscale0)); check_nomsg(cpl_table_fill_column_window(rscale0,"RATIO",0,nrow,0)); // For each range extract proper: obj, sky, wave spectra for (i=0;i 0) { switch(i) { case 0: ck0_nomsg(sinfo_get_obj_sky_wav_sub(*int_obj,*int_sky,lambda,lr41,wtol, &obj_lr,&sky_lr,&wav_lr)); break; case 1: ck0_nomsg(sinfo_get_obj_sky_wav_sub(*int_obj,*int_sky,lambda,lr52,wtol, &obj_lr,&sky_lr,&wav_lr)); break; case 2: ck0_nomsg(sinfo_get_obj_sky_wav_sub(*int_obj,*int_sky,lambda,lr63,wtol, &obj_lr,&sky_lr,&wav_lr)); break; case 3: ck0_nomsg(sinfo_get_obj_sky_wav_sub(*int_obj,*int_sky,lambda,lr74,wtol, &obj_lr,&sky_lr,&wav_lr)); break; case 4: ck0_nomsg(sinfo_get_obj_sky_wav_sub(*int_obj,*int_sky,lambda,lr02,wtol, &obj_lr,&sky_lr,&wav_lr)); break; case 5: ck0_nomsg(sinfo_get_obj_sky_wav_sub(*int_obj,*int_sky,lambda,lr85,wtol, &obj_lr,&sky_lr,&wav_lr)); break; case 6: ck0_nomsg(sinfo_get_obj_sky_wav_sub(*int_obj,*int_sky,lambda,lr20,wtol, &obj_lr,&sky_lr,&wav_lr)); break; case 7: ck0_nomsg(sinfo_get_obj_sky_wav_sub(*int_obj,*int_sky,lambda,lr31,wtol, &obj_lr,&sky_lr,&wav_lr)); break; case 8: ck0_nomsg(sinfo_get_obj_sky_wav_sub(*int_obj,*int_sky,lambda,lr42,wtol, &obj_lr,&sky_lr,&wav_lr)); break; case 9: ck0_nomsg(sinfo_get_obj_sky_wav_sub(*int_obj,*int_sky,lambda,lr53,wtol, &obj_lr,&sky_lr,&wav_lr)); break; case 10: ck0_nomsg(sinfo_get_obj_sky_wav_sub(*int_obj,*int_sky,lambda,lr64,wtol, &obj_lr,&sky_lr,&wav_lr)); break; case 11: ck0_nomsg(sinfo_get_obj_sky_wav_sub(*int_obj,*int_sky,lambda,lr75,wtol, &obj_lr,&sky_lr,&wav_lr)); break; case 12: ck0_nomsg(sinfo_get_obj_sky_wav_sub(*int_obj,*int_sky,lambda,lr86,wtol, &obj_lr,&sky_lr,&wav_lr)); break; case 13: ck0_nomsg(sinfo_get_obj_sky_wav_sub(*int_obj,*int_sky,lambda,lr97,wtol, &obj_lr,&sky_lr,&wav_lr)); break; case 14: ck0_nomsg(sinfo_get_obj_sky_wav_sub(*int_obj,*int_sky,lambda,lr00, wtol,&obj_lr,&sky_lr,&wav_lr)); break; default: sinfo_msg_error("case not supported"); goto cleanup; } if(sky_lr == NULL || obj_lr == NULL || wav_lr == NULL) { finite_pix_i=0; sinfo_msg("no good pix left"); } else { //AMO: the following 2 seems to be critical for robustness //check_nomsg(cpl_table_save(sky_lr,NULL,NULL,"out_skylr0.fits", // CPL_IO_DEFAULT)); //check_nomsg(cpl_table_save(obj_lr,NULL,NULL,"out_objlr0.fits", // CPL_IO_DEFAULT)); //check_nomsg(cpl_table_save(wav_lr,NULL,NULL,"out_wavlr0.fits", // CPL_IO_DEFAULT)); check_nomsg(finite_pix_i=sinfo_table_sky_obj_flag_nan(&sky_lr, &obj_lr, &wav_lr)); } if (finite_pix_i > npixw) { // identify sky lines //sinfo_msg("finite_pix_i=%d",finite_pix_i); check_nomsg(cpl_table_erase_invalid(obj_lr)); check_nomsg(cpl_table_erase_invalid(sky_lr)); check_nomsg(cpl_table_erase_invalid(wav_lr)); check_nomsg(sky_med=cpl_table_get_column_median(sky_lr,"INT")); check_nomsg(sky_sdv=cpl_table_get_column_stdev(sky_lr,"INT")); check_nomsg(cpl_table_select_all(sky_lr)); sky_thresh=sky_med+sky_sdv; //sinfo_msg("sky_thresh=%f",sky_thresh); check_nomsg(xxx1_i=cpl_table_and_selected_double(sky_lr,"INT", CPL_GREATER_THAN,sky_thresh)); check_nomsg(xxx1 = cpl_table_extract_selected(sky_lr)); check_nomsg(cpl_table_select_all(sky_lr)); if (xxx1_i > 0) { //separate line and continuum regions //by convolving with a hat region of large as a line //sinfo_msg("xxx1_i=%d",xxx1_i); check_nomsg(xxx2 = cpl_table_duplicate(sky_lr)); //check_nomsg(cpl_table_save(sky_lr,NULL,NULL, // "out_skylr.fits",CPL_IO_DEFAULT)); //check_nomsg(cpl_table_save(obj_lr,NULL,NULL, // "out_objlr.fits",CPL_IO_DEFAULT)); //check_nomsg(cpl_table_save(xxx2,NULL,NULL, // "out_xxx2_0.fits",CPL_IO_DEFAULT)); ck0_nomsg(sinfo_table_threshold(&xxx2,"INT",sky_thresh, sky_thresh,0.,10.)); //check_nomsg(cpl_table_save(xxx2,NULL,NULL, // "out_xxx2_1.fits",CPL_IO_DEFAULT)); /* TODO xxx2[xxx1] = 10.; */ //sinfo_msg("npixw/2=%d",npixw); check_nomsg(sinfo_convolve_kernel(&xxx2,npixw/2)); //check_nomsg(cpl_table_save(xxx2,NULL,NULL, // "out_xxx2_2.fits",CPL_IO_DEFAULT)); // get line_regions check_nomsg(line_i=cpl_table_and_selected_double(xxx2,"CNV", CPL_GREATER_THAN,0)); check_nomsg(line_regions=cpl_table_extract_selected(xxx2)); //check_nomsg(cpl_table_save(line_regions,NULL,NULL, //"out_line_regions.fits",CPL_IO_DEFAULT)); check_nomsg(cpl_table_erase_column(line_regions,"INT")); check_nomsg(cpl_table_erase_column(line_regions,"CNV")); check_nomsg(cpl_table_select_all(xxx2)); // get cont_regions check_nomsg(cont_i=cpl_table_and_selected_double(xxx2,"CNV", CPL_EQUAL_TO,0)); check_nomsg(cont_regions=cpl_table_extract_selected(xxx2)); //check_nomsg(cpl_table_save(line_regions,NULL,NULL, //"out_cont_regions.fits",CPL_IO_DEFAULT)); check_nomsg(cpl_table_erase_column(cont_regions,"INT")); check_nomsg(cpl_table_erase_column(cont_regions,"CNV")); check_nomsg(cpl_table_select_all(xxx2)); sinfo_free_table(&xxx2); if (line_i >= 3 && cont_i >= 3) { //If we have enough line points and continuum points //sinfo_msg("line_i=%d cont_i=%d",line_i,cont_i); if (i == 0) sinfo_msg("optimising 4-1 transitions"); if (i == 1) sinfo_msg("optimising 5-2 transitions"); if (i == 2) sinfo_msg("optimising 6-3 transitions"); if (i == 3) sinfo_msg("optimising 7-4 transitions"); if (i == 4) sinfo_msg("optimising 0-2 transitions"); if (i == 5) sinfo_msg("optimising 8-5 transitions"); if (i == 6) sinfo_msg("optimising 2-0 transitions"); if (i == 7) sinfo_msg("optimising 3-1 transitions"); if (i == 8) sinfo_msg("optimising 4-2 transitions"); if (i == 9) sinfo_msg("optimising 5-3 transitions"); if (i == 10) sinfo_msg("optimising 6-4 transitions"); if (i == 11) sinfo_msg("optimising 7-5 transitions"); if (i == 12) sinfo_msg("optimising 8-6 transitions"); if (i == 13) sinfo_msg("optimising 9-7 transitions"); if (i == 14) sinfo_msg("optimising final bit"); // Fit the object profile='fline_res' of the sky line residuals // left after proper scaled sky spectrum lines (and continua) // subtraction. Then determines median and stdev to flag outliers //Free memory for each loop sinfo_free_table(&obj_cont); sinfo_free_table(&sky_cont); sinfo_free_table(&sky_line); sinfo_free_table(&obj_line); //Identify obj lines and continuum, same for sky cknull_nomsg(obj_line=sinfo_table_select_range(obj_lr,line_regions, wtol)); cknull_nomsg(sky_line=sinfo_table_select_range(sky_lr,line_regions, wtol)); cknull_nomsg(obj_cont=sinfo_table_select_range(obj_lr,cont_regions, wtol)); cknull_nomsg(sky_cont=sinfo_table_select_range(sky_lr,cont_regions, wtol)); //Here was commented //check_nomsg(cpl_table_save(line_regions,NULL,NULL, // "out_line.fits",CPL_IO_DEFAULT)); //check_nomsg(cpl_table_save(cont_regions,NULL,NULL, // "out_cont.fits",CPL_IO_DEFAULT)); //check_nomsg(cpl_table_save(obj_cont,NULL,NULL, // "out_obj_cont.fits",CPL_IO_DEFAULT)); //check_nomsg(cpl_table_save(obj_line,NULL,NULL, // "out_obj_line.fits",CPL_IO_DEFAULT)); //check_nomsg(cpl_table_save(sky_line,NULL,NULL, // "out_sky_line.fits",CPL_IO_DEFAULT)); //check_nomsg(cpl_table_save(sky_cont,NULL,NULL, // "out_sky_cont.fits",CPL_IO_DEFAULT)); sinfo_free_table(&fline_res); //FIXME: in some cases obj_cont is empty //sinfo_msg("first line ratio determination"); ck0_nomsg(sinfo_get_line_ratio(obj_line,obj_cont, sky_line,sky_cont,method,&r)); sinfo_msg("1st Line ratio %g",r); if(cpl_table_get_nrow(obj_cont) > 0) { check_nomsg(fline_res=sinfo_table_interpol(obj_line,obj_cont, sky_line,sky_cont, r)); } else { check_nomsg(fline_res=cpl_table_duplicate(obj_line)); } // check if there are outliers cpl_table_select_all(fline_res); check_nomsg(fmed = cpl_table_get_column_median(fline_res,"INT")); check_nomsg(fsdv = cpl_table_get_column_stdev(fline_res,"INT")); check_nomsg(cpl_table_duplicate_column(fline_res,"AINT", fline_res,"INT")); check_nomsg(cpl_table_multiply_columns(fline_res,"AINT","INT")); check_nomsg(cpl_table_power_column(fline_res,"AINT",0.5)); check_nomsg(fclip_i=cpl_table_and_selected_double(fline_res,"AINT", CPL_GREATER_THAN, fmed+3*fsdv)); check_nomsg(cpl_table_select_all(fline_res)); if (fclip_i > 0) { // do a k-sigma clip to select a better line region //sinfo_msg("fclip_i=%d",fclip_i); //Find again line_regions check_nomsg(line_i=cpl_table_and_selected_double(fline_res, "AINT", CPL_NOT_GREATER_THAN, fmed+3*fsdv)); // get new (better) line_regions sinfo_free_table(&line_regions); //sinfo_msg("line_i=%d",line_i); check_nomsg(line_regions=cpl_table_extract_selected(fline_res)); check_nomsg(cpl_table_erase_column(line_regions,"INT")); check_nomsg(cpl_table_erase_column(line_regions,"AINT")); sinfo_free_table(&obj_line); sinfo_free_table(&sky_line); //check_nomsg(cpl_table_save(obj_lr,NULL,NULL, // "out_obj_lr.fits",CPL_IO_DEFAULT)); //check_nomsg(cpl_table_save(line_regions,NULL,NULL, // "out_line_regions.fits", // CPL_IO_DEFAULT)); // The following 2 may return an error so we do not check and // later we reset the error obj_line=sinfo_table_select_range(obj_lr,line_regions,wtol); sky_line=sinfo_table_select_range(sky_lr,line_regions,wtol); fline_i=cpl_table_get_nrow(line_regions); //sinfo_msg("fline_i=%d",fline_i); if(fline_i>=3) { // repeat the determination of the line ratio //sinfo_msg("second line ratio determination"); ck0_nomsg(sinfo_get_line_ratio(obj_line,obj_cont, sky_line,sky_cont,method,&r)); sinfo_msg("2nd Line ratio %g",r); } else { cpl_error_reset(); } sinfo_free_table(&sky_line); sinfo_free_table(&obj_line); } cpl_msg_info(cpl_func,"use %" CPL_SIZE_FORMAT " pixels for line and %" CPL_SIZE_FORMAT " for continuum estimation", cpl_table_get_nrow(line_regions),cpl_table_get_nrow(cont_regions)); sinfo_msg("OH spectrum scaling = %f ",r); check_nomsg(cpl_array_set_double(rfit,i,r)); ck0_nomsg(sinfo_table_set(&rscale0,wav_lr,r,wtol)); } /* end if line_i */ } /* end if xxx1_i */ } /* end finite_pix_i */ } sinfo_free_table(&xxx1); sinfo_free_table(&xxx2); sinfo_free_table(&sky_lr); sinfo_free_table(&obj_lr); sinfo_free_table(&wav_lr); sinfo_free_table(&line_regions); sinfo_free_table(&cont_regions); } /* end for loop on i */ sinfo_free_array(&do_hk); sinfo_free_array(&rfit); //sinfo_msg("n scale=%d",cpl_table_get_nrow(rscale0)); //check_nomsg(cpl_table_save(rscale0,NULL,NULL, // "out_rscale0.fits",CPL_IO_DEFAULT)); check_nomsg(cpl_table_select_all(rscale0)); /* TODO: here one has to implementa an interpol function check_nomsg(range0_i=cpl_table_and_selected_double(rscale0,"RATIO", CPL_NOT_EQUAL_TO,0)); */ check_nomsg(*rscale = cpl_table_extract_selected(rscale0)); sinfo_free_table(&rscale0); check_nomsg(rat_sky=cpl_table_duplicate(*int_sky)); check_nomsg(cpl_table_duplicate_column(rat_sky,"RATIO",*rscale,"RATIO")); check_nomsg(cpl_table_duplicate_column(*int_obj,"COR_FCT_VIB", *rscale,"RATIO")); //check_nomsg(cpl_table_save(rat_sky,NULL,NULL, // "rat_sky0.fits",CPL_IO_DEFAULT)); check_nomsg(cpl_table_multiply_columns(rat_sky,"INT","RATIO")); //check_nomsg(cpl_table_save(*int_obj,NULL,NULL, // "out_obj0.fits",CPL_IO_DEFAULT)); //check_nomsg(cpl_table_save(*int_sky,NULL,NULL, // "out_sky0.fits",CPL_IO_DEFAULT)); /* check_nomsg(cpl_table_duplicate_column(*int_obj,"INTC",*int_obj,"INT")); check_nomsg(cpl_table_duplicate_column(*int_obj,"SKYC",*int_sky,"INTC")); check_nomsg(cpl_table_subtract_columns(*int_obj,"INTC","SKYC")); */ // do simple rotational correction if (do_rot == 1) { //finite_pix = where(finite(int_sky) && finite(int_obj),finite_pix_i); check_nomsg(min_lrange=cpl_table_get_column_min(lrange,"WAVE")); check_nomsg(max_lrange=cpl_table_get_column_max(lrange,"WAVE")); //sinfo_msg("min_lrange=%g max_lrange=%g",min_lrange,max_lrange); //check_nomsg(finite_pix_i=sinfo_table_sky_obj_flag_nan(&rat_sky, // int_obj, // &lambda)); check_nomsg(finite_pix_i=sinfo_table_sky_obj_flag_nan(int_sky, int_obj, &lambda)); check_nomsg(finite_pix=cpl_table_duplicate(lambda)); //TODO: lambda invalid values need to be reset to valid (?) check_nomsg(cpl_table_erase_invalid(finite_pix)); if (finite_pix_i > npixw) { //finite_pix = finite_pix[where(finite_pix > min(lrange) && // finite_pix < max(lrange))]; check_nomsg(cpl_table_and_selected_double(finite_pix,"WAVE", CPL_GREATER_THAN, min_lrange)); check_nomsg(cpl_table_and_selected_double(finite_pix,"WAVE", CPL_LESS_THAN, max_lrange)); check_nomsg(tmp_tbl=cpl_table_extract_selected(finite_pix)); sinfo_free_table(&finite_pix); check_nomsg(finite_pix=cpl_table_duplicate(tmp_tbl)); sinfo_free_table(&tmp_tbl); sinfo_free_table(&sky_lr); sinfo_free_table(&obj_lr); sinfo_free_table(&wav_lr); cknull(sky_lr=sinfo_table_select_range(rat_sky,finite_pix,wtol), "extracting sky sub range"); cknull(obj_lr=sinfo_table_select_range(*int_obj,finite_pix,wtol), "extracting obj sub range"); cknull(wav_lr=sinfo_table_select_range(lambda,finite_pix,wtol), "extracting sky sub range"); //check_nomsg(cpl_table_save(rat_sky,NULL,NULL, // "out_rat_sky.fits",CPL_IO_DEFAULT)); //check_nomsg(cpl_table_save(finite_pix,NULL,NULL, // "out_finite_pix.fits",CPL_IO_DEFAULT)); //check_nomsg(cpl_table_save(sky_lr,NULL,NULL, // "out_sky_lr.fits",CPL_IO_DEFAULT)); //check_nomsg(cpl_table_save(obj_lr,NULL,NULL, // "out_obj_lr.fits",CPL_IO_DEFAULT)); //check_nomsg(cpl_table_save(wav_lr,NULL,NULL, // "out_wav_lr.fits",CPL_IO_DEFAULT)); //The following may fail (sky_lr may be empty) so we do not check if(1 == cpl_table_has_valid(sky_lr,"INT")) { check_nomsg(sky_med=cpl_table_get_column_median(sky_lr,"INT")); check_nomsg(sky_sdv=cpl_table_get_column_stdev(sky_lr,"INT")); sky_thresh=sky_med+sky_sdv; //xxx1 = where(sky_lr > median(sky_lr)+stddev(sky_lr),xxx1_i); check_nomsg(xxx1_i=cpl_table_and_selected_double(sky_lr,"INT", CPL_GREATER_THAN,sky_thresh)); check_nomsg(xxx1=cpl_table_extract_selected(sky_lr)); check_nomsg(cpl_table_select_all(sky_lr)); } else { xxx1_i=0; } if (xxx1_i > 0) { sinfo_msg("xxx1_i=%d",xxx1_i); sinfo_msg("wav_lr wmin=%g wmax=%g", cpl_table_get_column_min(wav_lr,"WAVE"), cpl_table_get_column_max(wav_lr,"WAVE")); cknull_nomsg(llr_xxx1=sinfo_table_select_range(wav_lr,xxx1,wtol)); //check_nomsg(cpl_table_save(wav_lr,NULL,NULL, // "out_llr_xxx1.fits",CPL_IO_DEFAULT)); cknull(low_pos=sinfo_find_rot_waves(w_rot_low,npixw,wtol,llr_xxx1), "Determining low positions"); check_nomsg(low_pos_i=cpl_table_get_nrow(low_pos)); //check_nomsg(cpl_table_dump(low_pos,0,low_pos_i,stdout)); cknull(med_pos=sinfo_find_rot_waves(w_rot_med,npixw,wtol,llr_xxx1), "Determining med positions"); check_nomsg(med_pos_i=cpl_table_get_nrow(med_pos)); //check_nomsg(cpl_table_dump(med_pos,0,med_pos_i,stdout)); //TODO: //hipos = [0] //for i=0,n_elements(xxx1)-1 do begin // x1 = where(lowpos eq i,x1_i) // x2 = where(medpos eq i,x2_i) // if (x1_i eq 0 and x2_i eq 0) then hipos = [hipos,i] //endfor //hipos = hipos[1:n_elements(hipos)-1] //TODO: hi_pos=sinfo_find_rot_waves(w_rot_hi,npixw,wtol,wav_lr); cknull(hi_pos=sinfo_table_extract_rest(xxx1,low_pos,med_pos,wtol), "determining hi position"); check_nomsg(hi_pos_i=cpl_table_get_nrow(hi_pos)); //check_nomsg(cpl_table_dump(hi_pos,0,hi_pos_i,stdout)); //xxx2[xxx1] = 10.; check_nomsg(xxx2 = cpl_table_duplicate(sky_lr)); check_nomsg(nrow=cpl_table_get_nrow(sky_lr)); //check_nomsg(cpl_table_save(xxx2,NULL,NULL, // "out_xxx1.fits",CPL_IO_DEFAULT)); //check_nomsg(cpl_table_save(xxx2,NULL,NULL, // "out_xxx2_0.fits",CPL_IO_DEFAULT)); // AMO: Why the following? //check_nomsg(cpl_table_fill_column_window(xxx2,"INT",0,nrow,0)); //xxx2 = convol(xxx2,replicate(1,npixw),/edge_truncate,/center); //cont_regions = where(xxx2 == 0,cont_i); ck0_nomsg(sinfo_table_threshold(&xxx2,"INT",sky_thresh, sky_thresh,0.,10.)); sinfo_msg("sky_thresh=%g %g %f",sky_thresh,sky_med,sky_sdv); //check_nomsg(cpl_table_save(xxx2,NULL,NULL, // "out_xxx2_1.fits",CPL_IO_DEFAULT)); check_nomsg(sinfo_convolve_kernel(&xxx2,npixw/2)); //check_nomsg(cpl_table_save(xxx2,NULL,NULL, // "out_xxx2_2.fits",CPL_IO_DEFAULT)); check_nomsg(cont_i=cpl_table_and_selected_double(xxx2,"CNV", CPL_EQUAL_TO,0)); check_nomsg(cont_regions=cpl_table_extract_selected(xxx2)); sinfo_free_table(&xxx2); check_nomsg(cpl_table_erase_column(cont_regions,"INT")); check_nomsg(cpl_table_erase_column(cont_regions,"CNV")); check(low_pos_i=sinfo_get_sub_regions(sky_lr,xxx1,low_pos,wtol, npixw,&low_regions),"failed determining low regions"); check(med_pos_i=sinfo_get_sub_regions(sky_lr,xxx1,med_pos,wtol, npixw,&med_regions),"failed determining med regions"); check(hi_pos_i=sinfo_get_sub_regions(sky_lr,xxx1,hi_pos,wtol, npixw,&hi_regions),"failed determining hi regions"); /* sinfo_msg("xxx1: wmin=%g wmax=%g", cpl_table_get_column_min(xxx1,"WAVE"), cpl_table_get_column_max(xxx1,"WAVE")); sinfo_msg("low_pos: wmin=%g wmax=%g", cpl_table_get_column_min(low_pos,"WAVE"), cpl_table_get_column_max(low_pos,"WAVE")); */ sinfo_msg("hi_pos_i : %d med_pos_i : %d low_pos_i : %d cont_i: %d", hi_pos_i, med_pos_i, low_pos_i, cont_i); if (hi_pos_i >= 3 && med_pos_i >= 3 && low_pos_i >= 3 && cont_i >= 3) { //compute line ratio for hi_regions ck0_nomsg(sinfo_compute_line_ratio(obj_lr, sky_lr,wtol, method, hi_regions,cont_regions,&rhi)); sinfo_msg("high rotational OH scaling %g",rhi); //compute line ratio for med_regions ck0_nomsg(sinfo_compute_line_ratio(obj_lr, sky_lr,wtol, method, med_regions,cont_regions,&rmed)); sinfo_msg("P1(3.5) & R1(1.5) rotational OH scaling %g ",rmed); //compute line ratio for med_regions ck0_nomsg(sinfo_compute_line_ratio(obj_lr, sky_lr,wtol, method, low_regions,cont_regions,&rlow)); sinfo_msg("P1(2.5) & Q1(1.5) rotational OH scaling %g",rlow); cknull(low_scale=sinfo_find_rot_waves(w_rot_low,npixw,wtol,lambda), "Determining low scale"); cknull(med_scale=sinfo_find_rot_waves(w_rot_low,npixw,wtol,lambda), "Determining low scale"); check_nomsg(cpl_table_multiply_scalar(*rscale,"RATIO",rhi)); ck0_nomsg(sinfo_table_fill_column_over_range(rscale,med_scale, "RATIO",rmed/rhi,wtol)); ck0_nomsg(sinfo_table_fill_column_over_range(rscale,low_scale, "RATIO",rlow/rhi,wtol)); } } //xxx1_i > 0 }//finitepix > npixw }//do_rot==1 //end of new rotational bit //check_nomsg(cpl_table_save(*int_obj,NULL,NULL, // "out_obj.fits",CPL_IO_DEFAULT)); //check_nomsg(cpl_table_save(*int_sky,NULL,NULL, // "out_sky.fits",CPL_IO_DEFAULT)); check_nomsg(cpl_table_duplicate_column(*int_sky,"INTC",*int_sky,"INT")); //sinfo_msg("n sky=%d",cpl_table_get_nrow(*int_sky)); //sinfo_msg("n scale=%d",cpl_table_get_nrow(*rscale)); check_nomsg(cpl_table_duplicate_column(*int_obj,"COR_FCT_ALL", *rscale,"RATIO")); check_nomsg(cpl_table_duplicate_column(*int_sky,"RATIO",*rscale,"RATIO")); check_nomsg(cpl_table_multiply_columns(*int_sky,"INTC","RATIO")); check_nomsg(cpl_table_duplicate_column(*int_obj,"INTC",*int_obj,"INT")); //check_nomsg(cpl_table_save(*int_obj,NULL,NULL, // "out_obj1.fits",CPL_IO_DEFAULT)); //check_nomsg(cpl_table_save(*int_sky,NULL,NULL, // "out_sky1.fits",CPL_IO_DEFAULT)); check_nomsg(cpl_table_duplicate_column(*int_obj,"SKYC",*int_sky,"INTC")); check_nomsg(cpl_table_subtract_columns(*int_obj,"INTC","SKYC")); check_nomsg(cpl_table_erase_column(*int_sky,"INT")); check_nomsg(cpl_table_name_column(*int_sky,"INTC","INT")); cleanup: sinfo_free_table(&llr_xxx1); sinfo_free_table(&hi_pos); sinfo_free_table(&low_pos); sinfo_free_table(&med_pos); sinfo_free_table(&low_regions); sinfo_free_table(&med_regions); sinfo_free_table(&hi_regions); sinfo_free_table(&low_scale); sinfo_free_table(&med_scale); sinfo_free_table(&finite_pix); sinfo_free_table(&xxx1_sub); sinfo_free_table(&tmp_tbl); sinfo_free_table(&rat_sky); sinfo_free_table(&fline_res); sinfo_free_table(&sky_cont); sinfo_free_table(&obj_cont); sinfo_free_table(&obj_line); sinfo_free_table(&sky_line); sinfo_free_table(&rscale0); sinfo_free_table(&xxx1); sinfo_free_table(&xxx2); sinfo_free_table(&line_regions); sinfo_free_table(&cont_regions); sinfo_free_table(&sky_lr); sinfo_free_table(&obj_lr); sinfo_free_table(&wav_lr); sinfo_free_array(&rfit); sinfo_free_array(&do_hk); return; } /** @name sinfo_table_get_index_of_max @memo find raw index of table column's max @param t input table @param name input table column @param type CPL_TYPE of input table column @return table raw correspoinding to table max */ int sinfo_table_get_index_of_max(cpl_table* t,const char* name,cpl_type type) { int i=0; int result=0; int nrow=0; int* pi=NULL; float* pf=NULL; double* pd=NULL; double max=0; if(t == NULL) { cpl_error_set(cpl_func, CPL_ERROR_NULL_INPUT); return result; } max=cpl_table_get_column_max(t,name); nrow=cpl_table_get_nrow(t); switch(type) { case CPL_TYPE_INT: pi=cpl_table_get_data_int(t,name); for(i=0;i0) ? sky_max : 0; gsky_min = (sky_min<0) ? sky_min : 0; //gsky_pos = where(int_sky > 1.*gsky_max || int_sky < 1.*gsky_min,gskypos_i); check_nomsg(cpl_table_select_all(*int_sky)); ck0_nomsg(sinfo_table_set_nan_out_min_max(int_sky,"INT",gsky_min,gsky_max)); //check_nomsg(cpl_table_save(*int_sky,NULL,NULL, // "out_gsky_pos.fits",CPL_IO_DEFAULT)); sinfo_free_table(&xsky); //sinfo_msg("gsky_min=%f gsky_max=%f",gsky_min,gsky_max); cknull_nomsg(xobj=sinfo_table_select_range(*int_obj,grange,wtol)); check_nomsg(obj_min=cpl_table_get_column_min(xobj,"INT")); check_nomsg(obj_max=cpl_table_get_column_max(xobj,"INT")); //check_nomsg(cpl_table_save(xobj,NULL,NULL,"out_xobj.fits",CPL_IO_DEFAULT)); gobj_max = (obj_max>0) ? obj_max : 0; gobj_min = (obj_min<0) ? obj_min : 0; //gobj_pos=where(int_obj > 1.*gobj_max || int_obj < 1.*gobj_min,gobj_pos_i); check_nomsg(cpl_table_select_all(*int_obj)); ck0_nomsg(sinfo_table_set_nan_out_min_max(int_obj,"INT",gobj_min,gobj_max)); //check_nomsg(cpl_table_save(*int_obj,NULL,NULL, // "out_gobj_pos.fits",CPL_IO_DEFAULT)); //sinfo_msg("gobj_min=%f gobj_max=%f",gobj_min,gobj_max); sinfo_free_table(&xobj); return 0; cleanup: sinfo_free_table(&xsky); sinfo_free_table(&xobj); return -1; } /** @name sinfo_table_fill_column_over_range @param inp input table @param ref input table specifying the wavelength range @param col input table column @param tol input tolerance @return pointer to a new allocated table corresponding to the selection */ static int sinfo_table_fill_column_over_range(cpl_table** inp, const cpl_table* ref, const char* col, const double val, const double tol) { int i=0; int k=0; int nref=0; int ninp=0; double* pwin=NULL; double* pcin=NULL; const double* pwrf=NULL; cknull(inp,"null input table"); cknull(ref,"null reference table"); check_nomsg(ninp=cpl_table_get_nrow(*inp)); check_nomsg(nref=cpl_table_get_nrow(ref)); check_nomsg(pwin=cpl_table_get_data_double(*inp,"WAVE")); check_nomsg(pcin=cpl_table_get_data_double(*inp,col)); check_nomsg(pwrf=cpl_table_get_data_double_const(ref,"WAVE")); k=0; i=0; /* sinfo_msg("ninp=%d nref=%d",ninp,nref); sinfo_msg("winp(0)=%f wref(0)=%f tol=%f diff=%f", pwin[0],pwrf[0],tol,fabs(pwin[0]-pwrf[0])); */ if(pwin[0]<=pwrf[0]) { //sinfo_msg("case 1"); for(i=0;ipwrf[0] //sinfo_msg("case 2"); for(k=0;kprf[0] //sinfo_msg_debug("case 2"); for(k=0;k low cut threshold */ static int sinfo_table_threshold(cpl_table** t, const char* column, const double low_cut, const double hig_cut, const double low_ass, const double hig_ass) { int nrow=0; int i=0; double* pdata=NULL; cknull(*t,"null input table!"); check_nomsg(nrow=cpl_table_get_nrow(*t)); check_nomsg(pdata=cpl_table_get_data_double(*t,column)); for(i=0;i= hig_cut) { pdata[i]=hig_ass; } } return 0; cleanup: return -1; } /** @name sinfo_table_set_column_invalid @param t input/output table @param col column to check static int sinfo_table_set_column_invalid(cpl_table** t,const char* col) { int nrow=0; int i=0; cknull(*t,"input table is NULL"); check_nomsg(nrow=cpl_table_get_nrow(*t)); for(i=0;i0 && (i+is+1) < nrow) { m=pi[i+is+1]-pi[i+is]; po[i]=pi[i+is]+m*ds; } } return out; cleanup: sinfo_free_table(&out); return NULL; } static cpl_imagelist* sinfo_cube_zshift_simple(cpl_imagelist* inp, const float shift) { int nx=0; int ny=0; int nz=0; int i=0; int j=0; int k=0; int ks=(int)shift; float ds=shift-ks; float* pu=NULL; float* pl=NULL; float* po=NULL; float int2=0; float int1=0; float m=0; cpl_imagelist* out=NULL; cpl_image* imgu=NULL; cpl_image* imgl=NULL; cpl_image* imgo=NULL; cknull(inp,"null input cube"); check_nomsg(nz=cpl_imagelist_get_size(inp)); check_nomsg(out=cpl_imagelist_duplicate(inp)); check_nomsg(imgo=cpl_imagelist_get(out,0)); check_nomsg(nx=cpl_image_get_size_x(imgo)); check_nomsg(ny=cpl_image_get_size_y(imgo)); for(k=0;k0 && (k+ks+1) < nz) { check_nomsg(imgu=cpl_imagelist_get(inp,k+ks+1)); check_nomsg(imgl=cpl_imagelist_get(inp,k+ks)); check_nomsg(imgo=cpl_imagelist_get(out,k)); check_nomsg(pu=cpl_image_get_data_float(imgu)); check_nomsg(pl=cpl_image_get_data_float(imgl)); check_nomsg(po=cpl_image_get_data_float(imgo)); for(j=0;j fmed+3*fsdv,fclip_i); check_nomsg(cpl_table_duplicate_column(lres,"AINT",lres,"INT")); check_nomsg(cpl_table_multiply_columns(lres,"AINT","INT")); check_nomsg(cpl_table_power_column(lres,"AINT",0.5)); check_nomsg(fclip_i=cpl_table_and_selected_double(lres,"AINT", CPL_GREATER_THAN, fthresh)); check_nomsg(cpl_table_select_all(lres)); if (fclip_i > 0) { //line_regions = line_regions[where(abs(fline_res) < fmed+3*fsdv)]; check_nomsg(line_i=cpl_table_and_selected_double(lres,"AINT", CPL_LESS_THAN, fthresh)); sinfo_free_table(&line_regions); check_nomsg(line_regions=cpl_table_extract_selected(lres)); sinfo_free_table(&lres); check_nomsg(cpl_table_erase_column(line_regions,"INT")); check_nomsg(cpl_table_erase_column(line_regions,"AINT")); if (line_i >= 3) { sinfo_free_table(&obj_lin); sinfo_free_table(&sky_lin); check_nomsg(obj_lin=sinfo_table_select_range(obj,line_regions,wtol)); check_nomsg(sky_lin=sinfo_table_select_range(sky,line_regions,wtol)); sinfo_free_table(&line_regions); //r = amoeba(1.e-5,function_name='fitsky',p0=[1.],scale=[0.1]); ck0_nomsg(sinfo_get_line_ratio(obj_lin,obj_cnt,sky_lin,sky_cnt,meth,r)); } } *r=fabs(*r); //Free memory sinfo_free_table(&obj_cnt); sinfo_free_table(&sky_cnt); sinfo_free_table(&sky_lin); sinfo_free_table(&obj_lin); sinfo_free_table(&lres); sinfo_free_table(&line_regions); return 0; cleanup: sinfo_free_table(&obj_cnt); sinfo_free_table(&sky_cnt); sinfo_free_table(&sky_lin); sinfo_free_table(&obj_lin); sinfo_free_table(&lres); sinfo_free_table(&line_regions); return -1; } /** @name sinfo_find_rot_waves @memo find wavelength range corresponding to a given rotational level @param w_rot array with wavelengths corresponding to rotational levels @param npix_w number of pixels corresponding to line width @param w_step wavelength sampling step @param range full wavelength range @return wavelength range corresponding to a given rotational level */ static cpl_table* sinfo_find_rot_waves( const double w_rot[], const int npix_w, const double w_step, cpl_table* range ) { int i=0; int x_i=0; int r_start=0; double w_min=0; double w_max=0; cpl_table* w_sel=NULL; cpl_table* res=NULL; check_nomsg(res = cpl_table_new(0)); check_nomsg(cpl_table_copy_structure(res,range)); for (i=0; i< NROT; i++) { //x = where(lambda > l_rot_low[i]-npixw*cdelto && // lambda < l_rot_low[i]+npixw*cdelto,x_i); w_min=w_rot[i]-npix_w*w_step; w_max=w_rot[i]+npix_w*w_step; check_nomsg(cpl_table_and_selected_double(range,"WAVE", CPL_GREATER_THAN,w_min)); check_nomsg(cpl_table_and_selected_double(range,"WAVE", CPL_LESS_THAN,w_max)); sinfo_free_table(&w_sel); check_nomsg(w_sel=cpl_table_extract_selected(range)); check_nomsg(x_i=cpl_table_get_nrow(w_sel)); if (x_i > 0) { check_nomsg(r_start=cpl_table_get_nrow(res)); //sinfo_msg("i=%d x_i=%d w_min=%g w_max=%g",i,x_i,w_min,w_max); check_nomsg(cpl_table_insert(res,w_sel,r_start)); } check_nomsg(cpl_table_select_all(range)); } //res = range[1:cpl_table_get_nrow(res)-1]; sinfo_free_table(&w_sel); return res; cleanup: sinfo_free_table(&w_sel); sinfo_free_table(&res); return NULL; } /** @name sinfo_get_obj_sky_wav_sub @param obj input object spectrum @param sky input sky spectrum @param wav input full wavelength range @param sel input selection wavelength range @param wtol intup wavelength tolerance @param sub_obj output sub set of object spectrum @param sub_sky output sub set of sky spectrum @param sub_wav output sub set of wavelength range */ static int sinfo_get_obj_sky_wav_sub(cpl_table* obj, cpl_table* sky, cpl_table* wav, cpl_table* sel, const double wtol, cpl_table** sub_obj, cpl_table** sub_sky, cpl_table** sub_wav) { cknull_nomsg(*sub_obj = sinfo_table_select_range(obj,sel,wtol)); cknull_nomsg(*sub_sky = sinfo_table_select_range(sky,sel,wtol)); cknull_nomsg(*sub_wav = sinfo_table_select_range(wav,sel,wtol)); return 0; cleanup: sinfo_free_table(&(*sub_obj)); sinfo_free_table(&(*sub_sky)); sinfo_free_table(&(*sub_wav)); return -1; } static int sinfo_get_sub_regions(cpl_table* sky, cpl_table* x1, cpl_table* pos, const double wtol, const int npixw, cpl_table** res) { cpl_table* x1_sub=NULL; cpl_table* x2=NULL; int nrow=0; int np=0; cknull(sky,"Null input sky table"); cknull(x1 ,"Null input x1 table"); cknull(pos,"Null input pos table"); check_nomsg(x2=cpl_table_duplicate(sky)); check_nomsg(nrow=cpl_table_get_nrow(sky)); check_nomsg(cpl_table_fill_column_window(x2,"INT",0,nrow,0)); //x2[x1[pos]] = 10.; //x2 = convol(x2,replicate(1,npixw),/edge_truncate,/center); //res = where(x2 > 0,hi_i); //cpl_table_save(x1, NULL, NULL, "out_x1.fits", CPL_IO_DEFAULT); x1_sub=sinfo_table_select_range(x1,pos,wtol); if(x1_sub != NULL) { ck0_nomsg(sinfo_table_fill_column_over_range(&x2,x1_sub,"INT",10.,wtol)); sinfo_free_table(&x1_sub); check_nomsg(sinfo_convolve_kernel(&x2,npixw/2)); check_nomsg(np=cpl_table_and_selected_double(x2,"CNV",CPL_GREATER_THAN,0)); check_nomsg(*res=cpl_table_extract_selected(x2)); sinfo_free_table(&x2); check_nomsg(cpl_table_erase_column(*res,"INT")); check_nomsg(cpl_table_erase_column(*res,"CNV")); } else { cpl_error_reset(); sinfo_free_table(&x1_sub); sinfo_free_table(&x2); return np; } return np; cleanup: sinfo_free_table(&x1_sub); sinfo_free_table(&x2); return -1; } static cpl_table* sinfo_table_extract_rest(cpl_table* inp, cpl_table* low, cpl_table* med, const double wtol) { cpl_table* out=NULL; double* pinp=NULL; double* plow=NULL; double* pmed=NULL; int nlow=0; int nmed=0; int nrow=0; int i=0; int k=0; cpl_table* tmp=NULL; cknull(inp,"null input table"); check_nomsg(tmp=cpl_table_duplicate(inp)); check_nomsg(nrow=cpl_table_get_nrow(tmp)); check_nomsg(cpl_table_new_column(tmp,"SEL",CPL_TYPE_INT)); check_nomsg(cpl_table_fill_column_window_int(tmp,"SEL",0,nrow,0)); check_nomsg(pinp=cpl_table_get_data_double(inp,"WAVE")); check_nomsg(plow=cpl_table_get_data_double(low,"WAVE")); check_nomsg(pmed=cpl_table_get_data_double(med,"WAVE")); nlow=cpl_table_get_nrow(low); //check_nomsg(cpl_table_save(low,NULL,NULL,"out_low.fits",CPL_IO_DEFAULT)); if(nlow > 0) { k=0; for(i=0;i 0) { for(i=0;i