sinfo_skycor.c

00001 /* $Id: sinfo_skycor.c,v 1.20 2007/01/10 08:28:17 amodigli Exp $
00002  *
00003  * This file is part of the SINFONI Pipeline
00004  * Copyright (C) 2002,2003 European Southern Observatory
00005  *
00006  * This program is free software; you can redistribute it and/or modify
00007  * it under the terms of the GNU General Public License as published by
00008  * the Free Software Foundation; either version 2 of the License, or
00009  * (at your option) any later version.
00010  *
00011  * This program is distributed in the hope that it will be useful,
00012  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00013  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014  * GNU General Public License for more details.
00015  *
00016  * You should have received a copy of the GNU General Public License
00017  * along with this program; if not, write to the Free Software
00018  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00019  */
00020 
00021 /*
00022  * $Author: amodigli $
00023  * $Date: 2007/01/10 08:28:17 $
00024  * $Revision: 1.20 $
00025  * $Name:  $
00026  */
00027 
00028 #ifdef HAVE_CONFIG_H
00029 #include <config.h>
00030 #endif
00031 
00032 /*-----------------------------------------------------------------------------
00033                                 Includes
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                                 Defines
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      /* this should be 4 */
00066 #define SKY_LINE_MIN_CUT 4      /* this should be 4 */
00067 #define HISTO_NBINS 100
00068 #define HISTO_MIN_SIZE 10
00069 #define HISTO_Y_CUT 10          /* 50 is a better value this 
00070                                    affects histo right marging*/
00071 #define HISTO_X_LEFT_CUT  1.0   /* 0.2 is better */
00072 #define HISTO_X_RIGHT_CUT  0.5  /* 1.0 is better */
00073 #define HISTO_DIST_TEMPC_MIN_FCT 5.   /* 5.0 */
00074 #define HISTO_DIST_TEMPC_MAX_FCT 0.25 /* 0.25 */
00075 
00076 #define XCOR_YSHIFT_KS_CLIP 3      /* this should be 5 */
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                             Functions prototypes
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   //double temp=252.69284;
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   // set wavelength ranges
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   //check_nomsg(cpl_table_save(grange, NULL, NULL, "out_grange.fits", 0));
00617  
00618   /*
00619   sinfo_msg("lr20=%d",cpl_table_get_nrow(lr20));
00620   sinfo_msg("lr31=%d",cpl_table_get_nrow(lr31));
00621   sinfo_msg("lr42=%d",cpl_table_get_nrow(lr42));
00622   sinfo_msg("lr53=%d",cpl_table_get_nrow(lr53));
00623   sinfo_msg("lr64=%d",cpl_table_get_nrow(lr64));
00624   sinfo_msg("lr75=%d",cpl_table_get_nrow(lr75));
00625   sinfo_msg("lr86=%d",cpl_table_get_nrow(lr86));
00626   sinfo_msg("lr97=%d",cpl_table_get_nrow(lr97));
00627   sinfo_msg("lr00=%d",cpl_table_get_nrow(lr00));
00628 
00629   sinfo_msg("min_lrange=%f",cpl_table_get_column_min(lrange,"INDEX"));
00630   sinfo_msg("min_grange=%f",cpl_table_get_column_min(grange,"INDEX"));
00631   sinfo_msg("min_srange=%f",cpl_table_get_column_min(lambdas,"WAVE"));
00632   sinfo_msg("min_mrange=%f",cpl_table_get_column_min(mrange,"INDEX"));
00633 
00634   sinfo_msg("max_lrange=%f",cpl_table_get_column_max(lrange,"INDEX"));
00635   sinfo_msg("max_grange=%f",cpl_table_get_column_max(grange,"INDEX"));
00636   sinfo_msg("max_srange=%f",cpl_table_get_column_max(lambdas,"WAVE"));
00637   sinfo_msg("max_mrange=%f",cpl_table_get_column_max(mrange,"INDEX"));
00638   */
00639 
00640   sinfo_msg("Estimate noise");
00641   ck0(sinfo_object_estimate_noise(obj_frm,obj_noise_fit,&centre,&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   //cpl_imagelist_save(fdata,"out_fdata.fits",
00648   //                   CPL_BPP_DEFAULT,NULL,CPL_IO_DEFAULT);
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   //cpl_image_save(gpix,"out_gpix.fits",CPL_BPP_DEFAULT,NULL,CPL_IO_DEFAULT);
00656   //cpl_image_save(ratio,"out_ratio.fits",CPL_BPP_DEFAULT,NULL,CPL_IO_DEFAULT);
00657   //cpl_image_save(ima_sky,"out_ima_sky.fits",
00658   //               CPL_BPP_DEFAULT,NULL,CPL_IO_DEFAULT);
00659 
00660 
00661 
00662 
00663   // choose pixels which seems to be good sky pixels 
00664   // (95% spectral pixels are flagged) 
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   //cpl_image_save(mask,"out_mask.fits",CPL_BPP_DEFAULT,NULL,CPL_IO_DEFAULT);
00669  
00670   // threshold ratio for fraction 'minfract' of spatial pixels to be 'sky' 
00671   //sinfo_msg("To do: flag_threshold_sky_pixels");
00672 
00673 
00674   // sum spectra of flagged pixels in object and sky frames 
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   //check_nomsg(cpl_table_save(*int_obj, NULL, NULL, "out_int_obj.fits", 0));
00679   //check_nomsg(cpl_table_save(int_sky, NULL, NULL, "out_int_sky.fits", 0));
00680 
00681 
00682    // Computes thermal background
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   //check_nomsg(cpl_table_save(bkg,NULL,NULL,"out_thermal_background.fits",0));
00690 
00691   /*
00692   ck0(sinfo_pro_save_tbl(bkg,set,set,"out_thermal_background.fits",
00693              "THERMAL_BACKGROUND",NULL,cpl_func,config),
00694       "Error saving %s","THERMAL_BACKGROUND");
00695   */
00696 
00697   sinfo_msg("Subtracts thermal background from integrated OH spectrum");
00698   //sinfo_msg("nrow=%d %d",cpl_table_get_nrow(int_sky),
00699   //                       cpl_table_get_nrow(bkg));
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   //check_nomsg(cpl_table_save(int_sky, NULL, NULL, "out_int_sky_sub.fits", 0));
00705 
00706 
00707   // check for screw values at ends of spectrum
00708   sinfo_msg("Checks for screw values at ends of spectrum");
00709   //check_nomsg(cpl_table_save(*int_obj, NULL, NULL, "out_int_obj_chk.fits", 0));
00710   //check_nomsg(cpl_table_save(int_sky, NULL, NULL, "out_int_sky_chk.fits", 0));
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   //check_nomsg(cpl_table_save(*int_obj, NULL, NULL, "out_int_obj_chk.fits", 0));
00722   //check_nomsg(cpl_table_save(int_sky, NULL, NULL, "out_int_sky_chk.fits", 0));
00723   check_nomsg(wshift=sinfo_xcorr(*int_obj,int_sky,lambda,dispersion,line_hw));
00724 
00725  
00726   //wshift=-1.7164495*dispersion; 
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   //check_nomsg(cpl_table_save(*int_obj, NULL, NULL, "out_pip3_int_obj.fits", 0));
00732   //check_nomsg(cpl_table_save(int_sky, NULL, NULL, "out_pip3_int_sky.fits", 0));
00733 
00734 
00735 
00736   //DEBUG
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   //check_nomsg(cpl_table_save(*int_obj, NULL, NULL, "out_int_obj.fits", 0));
00747   //check_nomsg(cpl_table_save(int_sky, NULL, NULL, "out_int_sky.fits", 0));
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   //sinfo_free_table(&int_obj);
00786   sinfo_free_table(&int_sky);
00787 
00788   if (cpl_error_get_code() != CPL_ERROR_NONE) {
00789     //sinfo_free_imagelist(&obj_cor);
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   /* wavelength min-max J-H-K band */
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   /* wavelength boundaries between line groups corresponding 
00897      to transitions 5-2 to 9-7 */
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   /* Get Object relevant information */
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   /* defines object related wavelength ranges */
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   /* lrange */
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   /* mrange */
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   /* Get Sky relevant information */
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   /* defines sky related wavelength ranges */
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   // sum spectra of flagged spaxels 
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   //loop = where(mask > 0.5);
01263   //TO BE REMOVED: loop_tbl is not used
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   //Determines object spectrum
01272   for (i=0;i<zsize;i++) {
01273     check_nomsg(obj_slice = cpl_imagelist_get(obj,i));
01274 
01275     //pos = where(mask > 0.5 && finite(obj_slice),pos_i);
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     //int_obj[i] = mean(obj_slice[pos]);
01280         //TODO here obj_slice should be considered only on pos:
01281         //     one should do a selection/thresholding
01282         check_nomsg(cpl_table_set_double(*int_obj,"INT",i,
01283                                           cpl_image_get_mean(obj_slice)));
01284       } else {
01285         // select only poisitions where mask>0.5 and obj is finite
01286     // gslice = obj_slice[pos];
01287         //sinfo_msg("obj pos_i=%d",pos_i);
01288         
01289         check_nomsg(gslice = cpl_image_duplicate(obj_slice));
01290         check_nomsg(sinfo_image_flag_nan(&gslice));
01291     /*
01292         sinfo_msg("obj: min=%f max=%f",
01293                   cpl_image_get_min(obj_slice),
01294           cpl_image_get_max(obj_slice));
01295     */
01296         //check_nomsg(cpl_image_threshold(gslice,CPL_PIXEL_MINVAL,3.0e6,1,0));
01297         //check_nomsg(cpl_image_multiply(gslice,mask));
01298     check_nomsg(med = cpl_image_get_median(gslice));
01299     check_nomsg(sdv = cpl_image_get_stdev(gslice));
01300         //sinfo_msg("med=%f sdv=%f",med,sdv);
01301         //avg = mean(gslice[where(gslice < med+3*sdv && gslice > med-3*sdv)]);
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         //sinfo_msg("sky int=%f",avg);
01307       }
01308     }
01309 
01310     //Determines sky spectrum
01311     check_nomsg(sky_slice = cpl_imagelist_get(sky,i));
01312     //pos = where(mask > 0.5 and finite(sky_slice),pos_i);
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     //int_obj[i] = mean(obj_slice[pos]);
01317         //TODO here obj_slice should be considered only on pos:
01318         //     one should do a selection/thresholding
01319         check_nomsg(cpl_table_set_double(*int_sky,"INT",i,
01320                            cpl_image_get_mean(sky_slice)));
01321       } else {
01322         //sinfo_msg("pos_i=%d",pos_i);
01323         // select only poisitions where mask>0.5 and obj is finite
01324     // gslice = obj_slice[pos];
01325         check_nomsg(gslice = cpl_image_duplicate(sky_slice));
01326         check_nomsg(sinfo_image_flag_nan(&gslice));
01327         //check_nomsg(cpl_image_threshold(gslice,CPL_PIXEL_MINVAL,3.0e6,1,0));
01328         //check_nomsg(cpl_image_multiply(gslice,mask));
01329     check_nomsg(med = cpl_image_get_median(gslice));
01330     check_nomsg(sdv = cpl_image_get_stdev(gslice));
01331         //avg = mean(gslice[where(gslice < med+3*sdv && gslice > med-3*sdv)]);
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         if(i<100) {
01338            sinfo_msg("sky: wave=%f int=%f",
01339                       cpl_table_get_double(*int_sky,"WAVE",i,&status),avg);
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   //sinfo_msg("No bad pixels: %d",cnt);
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   //double mse=0;
01534   //double chired=0;
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   // Get Object relevant information 
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   //cpl_table_save(data_tbl, NULL, NULL, "out_data.fits", 0);
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   //sinfo_msg("value=%10.8f",value);
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     //cpl_table_save(histo, NULL, NULL, "out_pippo.fits", 0);
01598     check_nomsg(max_pos=sinfo_table_get_index_of_max(histo,"HY",CPL_TYPE_INT));
01599     //sinfo_msg("max_pos=%d",max_pos);
01600  
01601     /*
01602     check_nomsg(max_pos=sinfo_extract_table_rows(histo,"HY",
01603                                                  CPL_EQUAL_TO,max_h));
01604     sinfo_msg("size max_pos %d",cpl_table_get_nrow(max_pos));
01605     sinfo_msg("value max_pos %d",cpl_table_get_int(max_pos,"HY",0,&status));
01606     */
01607     min_x=max_pos-1;
01608     max_x=max_pos+2;
01609     //sinfo_msg("min_x=%d max_x=%d",min_x,max_x); 
01610  
01611     sinfo_free_table(&tmp_tbl1);  
01612     //sinfo_msg("x selection threshold: %f %d",
01613     //          cpl_table_get(histo,"HL",max_pos,&status),max_pos);
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     //sinfo_msg("threshold=%f",thres);
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     //cpl_table_save(min_xi, NULL, NULL, "out_min_xi.fits", 0);
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     //sinfo_msg("min_pos=%d max_pos=%d max(h)=%d min_xi_sz=%d x[maxpos[0]]=%f",
01635     //           min_pos,   max_pos,   max_h,    min_xi_sz, val);
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     //sinfo_msg("min_x=%d max_x=%d",min_x,max_x); 
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     //sinfo_msg("hmin=%f hmax=%f min_xi_sz=%d",hmin,hmax,min_xi_sz);
01648     //cpl_table_save(histo, NULL, NULL, "out_histo.fits", 0);
01649     sinfo_free_table(&histo);
01650     ck0(sinfo_histogram(data_tbl,nbins,hmin,hmax,&histo),"building histogram");
01651     //cpl_table_save(histo, NULL, NULL, "out_histo1.fits", 0);
01652     check_nomsg(cpl_table_add_scalar(histo,"HL",(hmax-hmin)/nbins/2));
01653     //cpl_table_save(histo, NULL, NULL, "out_histo2.fits", 0);
01654 
01655 
01656 
01657   }
01658   sinfo_free_table(&data_tbl);
01659   sinfo_free_table(&min_xi);
01660 
01661   //cpl_table_save(histo, NULL, NULL, "out_histo.fits", 0);
01662 
01663   check_nomsg(peak=cpl_table_get_column_max(histo,"HY"));
01664   //sinfo_msg("peak=%f",peak);
01665   sinfo_free_table(&tmp_tbl1);
01666 
01667   check_nomsg(tmp_tbl1=sinfo_extract_table_rows(histo,"HY",CPL_EQUAL_TO,peak));
01668   
01669   //cpl_table_save(tmp_tbl1, NULL, NULL, "out_tmp_tbl1.fits", 0);
01670 
01671 
01672   check_nomsg(*centre=cpl_table_get_column_mean(tmp_tbl1,"HL"));
01673   //sinfo_msg("Background level=%f",*centre);
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   //sinfo_msg("min HX %f",cpl_table_get_column_min(tmp_tbl2,"HL"));
01686   sinfo_free_table(&tmp_tbl2);
01687   //sinfo_msg("Tempc=%f",tempc);
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   //cpl_table_save(dist, NULL, NULL, "out_dist.fits", 0);
01701 
01702   //TODO
01703   //gaussfit(distx,disty,dista,nterms=3);
01704   //*noise=dista[2];
01705   *noise=tempc/2;
01706   /* THIS DOES NOT WORK */
01707   //sinfo_msg("FWHM/2=%f",*noise);
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     //sinfo_msg("Gauss fit parameters:" 
01724     //          "x0=%f sigma=%f area=%f offset=%f mse=%f chired=%f",
01725     //           x0,sigma,area,offset,mse,chired);
01726     //sinfo_msg("Background level=%f",*centre);
01727     //sinfo_msg("Noise=%f",sigma);
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   //*noise=18.7448;
01736   //*noise=20.585946;
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   /* double* phx=NULL; */
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   sinfo_msg("min=%f max=%f",
01850             cpl_table_get_column_min(dat,"DDATA"),
01851             cpl_table_get_column_max(dat,"DDATA"));
01852   */
01853   check_nomsg(cpl_table_and_selected_double(dat,"DDATA",
01854                                             CPL_NOT_GREATER_THAN,max));
01855   /*
01856   check_nomsg(cpl_table_and_selected_double(dat,"DDATA",CPL_LESS_THAN,max));
01857   */
01858   check_nomsg(tmp_tbl1=cpl_table_extract_selected(dat));
01859   /*
01860   sinfo_msg("min=%f max=%f",
01861              cpl_table_get_column_min(tmp_tbl1,"DDATA"),
01862              cpl_table_get_column_max(tmp_tbl1,"DDATA"));
01863   */
01864   /*
01865   check_nomsg(cpl_table_and_selected_double(tmp_tbl1,"DDATA",
01866                                             CPL_NOT_LESS_THAN,min));
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   sinfo_msg("min=%f max=%f",
01873              cpl_table_get_column_min(tmp_tbl2,"DDATA"),
01874              cpl_table_get_column_max(tmp_tbl2,"DDATA"));
01875   */
01876   sinfo_free_table(&tmp_tbl1);
01877   /*
01878   check_nomsg(tmp_tbl1=sinfo_extract_table_rows(dat,"DDATA",
01879                                                 CPL_NOT_GREATER_THAN,max));
01880   check_nomsg(tmp_tbl2=sinfo_extract_table_rows(tmp_tbl1,"DDATA",
01881                                                 CPL_NOT_LESS_THAN,min));
01882   */
01883 
01884   check_nomsg(ntot=cpl_table_get_nrow(tmp_tbl2));
01885   /* not necessry to sort:
01886     check_nomsg(sinfo_sort_table_1(tmp_tbl2,"DDATA",FALSE));*/
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   /* sinfo_msg("vmin=%f vmax=%f step=%f",vmin,vmax,vstp); */
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   /*check_nomsg(cpl_table_fill_column_window(*histo,"HX",0,nbins,0)); */ 
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   /*check_nomsg(phx=cpl_table_get_data_double(*histo,"HX")); */
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   //cpl_table_save(*histo, NULL, NULL, "out_histo_p5.fits", 0);
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  /* Get Object relevant information */
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   //check_nomsg(cpl_table_save(flag_tbl,NULL,NULL,"out_flag.fits",0));
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   /* Get Object relevant information */
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   /* Flag sky pixels in data cube */
02101   /* create images */
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   /* TODO */
02111   // fill image points:
02112   // g_img: mask with at least half good pixels along spectral range
02113   // r_img: mask with ratio of good pixels along spectral range
02114   // image: image with mean of spectrum over good pixels
02115 
02116   //check_nomsg(cpl_table_save(lambda, NULL, NULL, "out_lambda.fits", 0));
02117   //check_nomsg(cpl_table_save(mrange, NULL, NULL, "out_mrange.fits", 0));
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   //check_nomsg(cpl_imagelist_save(osel,"out_osel.fits", 
02128   //                               CPL_BPP_DEFAULT,NULL,CPL_IO_DEFAULT));
02129   //check_nomsg(cpl_imagelist_save(fsel,"out_fsel.fits", 
02130   //                               CPL_BPP_DEFAULT,NULL,CPL_IO_DEFAULT));
02131 
02132   for(j=0;j<ysz;j++) {
02133     for(i=0;i<xsz;i++) {
02134       // consider only planes in the proper wavelegth ranges
02135       cknull_nomsg(ospec=sinfo_slice_z(osel,i,j));
02136       cknull_nomsg(fspec=sinfo_slice_z(fsel,i,j));
02137       // consider only finite pixels
02138       check_nomsg(gpix_i=sinfo_table_extract_finite(ospec,fspec,&gpix,&gspec));
02139       //sinfo_msg("gpix_i=%d",gpix_i);
02140       if(gpix_i > 0) {
02141         // build two arrays of proper size
02142         all_pix=(double)gpix_i;
02143     /*
02144         sinfo_msg("flagspec: min=%f max=%f",
02145                   cpl_table_get_column_min(fspec,"VALUE"),
02146                   cpl_table_get_column_max(fspec,"VALUE"));
02147         sinfo_msg("good flagspec: min=%f max=%f",
02148                   cpl_table_get_column_min(gspec,"VALUE"),
02149                   cpl_table_get_column_max(gspec,"VALUE"));
02150     sinfo_msg("nfspec=%d",cpl_table_get_nrow(fspec));
02151         check_nomsg(cpl_table_save(fspec, NULL, NULL, "out_fspec.fits", 0));
02152         check_nomsg(cpl_table_save(gspec, NULL, NULL, "out_gspec.fits", 0));
02153     */
02154         //check_nomsg(flag_pix=cpl_table_and_selected_double(fspec,"VALUE",
02155         //                                              CPL_GREATER_THAN,0.5));
02156         //sinfo_msg("all_pix=%f flag_pix=%f",all_pix,flag_pix);
02157 
02158         check_nomsg(flag_pix=cpl_table_and_selected_double(gspec,"VALUE",
02159                                                         CPL_GREATER_THAN,0.5));
02160         //sinfo_msg("all_pix=%f flag_pix=%f",all_pix,flag_pix);
02161         // flag_pix = float(n_elements(where(fspec[gpix] > 0.5)));
02162         // compute the ratio between the two arrays
02163         ratio=flag_pix/all_pix;
02164         // considers only pixels with have at least half good pixels
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     //mean(ospec(gpix))
02170         check_nomsg(pimage[i+j*xsz]=cpl_table_get_column_mean(gpix,"VALUE"));
02171         //sinfo_msg("ix=%d iy=%d r=%f",i,j,ratio);
02172       }
02173       sinfo_free_table(&ospec);
02174       sinfo_free_table(&fspec);
02175       sinfo_free_table(&gpix);
02176       sinfo_free_table(&gspec);
02177 
02178     } /* end for over i */
02179   } /* end for over j */
02180   sinfo_free_imagelist(&osel);
02181   sinfo_free_imagelist(&fsel);
02182  
02183   /*
02184   cpl_image_save(*r_img,"out_r_img.fits",CPL_BPP_DEFAULT, NULL,CPL_IO_DEFAULT);
02185   cpl_image_save(*g_img,"out_g_img.fits",CPL_BPP_DEFAULT, NULL,CPL_IO_DEFAULT);
02186   cpl_image_save(*image,"out_image.fits",CPL_BPP_DEFAULT, NULL,CPL_IO_DEFAULT);
02187   */
02188   // get total(g_arr)
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   // choose pixels which seem to be sky (ie 95% of spectral pixels are flagged)
02251   check_nomsg(*mask=cpl_image_new(xsz,ysz,CPL_TYPE_DOUBLE));
02252   //r2 = where(r_img >= thres,r2i);
02253   // count good pixels: set to 0 what < thres and to 1 what > thres
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    //threshold ratio for fraction 'minfrac' of spatial pixels to be 'sky'
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      //hcum = r_img(sort(r_img));
02284      hcum = sinfo_image2table(r_img);
02285      check_nomsg(sinfo_sort_table_1(hcum,"VALUE",FALSE));
02286    
02287      //thresh = hcum[where(xcum/max(xcum) >= 1.-min_frac)];
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      //*mask[where(r_img >= thresh)] = 1;
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     //sinfo_msg("x=%g p=%g",px[i],pv[i]);
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   //TO BE FIXED: Why lrange to gat wave min and max: int_sky is sufficient
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   //Init amoeba fit parameters  
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   //check_nomsg(cpl_table_save(*bkg,NULL,NULL, "out_amoeba5.fits", 0));
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;/* There are 3 parameters */
02589     int ia[ndim];
02590 
02591   int NPOINTS=0;
02592 
02593 
02594   double temp1=0;
02595   double temp2=0;
02596 
02597   //double r0=80.306773;
02598   //double r1=450.50027;
02599   //double r2=252.17949;
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   //double ga1=0;
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   //check_nomsg(cpl_table_save(int_sky,NULL,NULL, "out_pippo.fits", 0));
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   //sinfo_msg("bkg: min=%g max=%g",bkg_min,bkg_max);
02641   //sinfo_scale_fct=sinfo_scale_fct*bkg_max;
02642   //sinfo_scale_fct=sinfo_scale_fct;
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   //To be sure one has not strange cases
02662   check_nomsg(cpl_table_and_selected_double(tmp,"INT",CPL_GREATER_THAN,-2));
02663   check_nomsg(ylr=cpl_table_extract_selected(tmp));
02664   //check_nomsg(cpl_table_save(ylr,NULL,NULL,"out_ylr_0.fits",0));
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   check_nomsg(xlr=sinfo_table_select_range(lambda,lrange,0.003));
02680   check_nomsg(ylr=sinfo_table_select_range(int_sky,lrange,0.003));
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   //ga1=p0[1];
02697   ga2=p0[2];
02698 
02699   //sinfo_msg("p= %g %g %g",p0[0],p0[1],p0[2]);
02700   check_nomsg(sinfo_table_flag_nan(&wlr,"INT"));
02701   check_nomsg(cpl_table_erase_invalid(wlr));
02702   //check_nomsg(cpl_table_save(xlr,NULL,NULL,"out_xlr.fits",0));
02703   //check_nomsg(cpl_table_save(ylr,NULL,NULL,"out_ylr.fits",0));
02704   //check_nomsg(cpl_table_save(wlr,NULL,NULL,"out_wlr.fits",0));
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   //check_nomsg(cpl_vector_fill(sy,1.));
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     sinfo_msg("before fit: a=%g %g %g",
02733               cpl_vector_get(a,0),
02734               cpl_vector_get(a,1),
02735               cpl_vector_get(a,2));
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     //sinfo_scale_fct=sinfo_scale_fct*bkg_max;
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   //ga2=252.69284;
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   // crosscorrelate obj & sky to check for lambda offset
02984 
02985   //if (mean(z[where(finite(z))]) < 0) z = z * (-1);
02986   //if sky mean is < 0 flip sky intensity
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   //check_nomsg(cpl_table_save(z,NULL,NULL,"out_z1.fits",0));
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   //check_nomsg(cpl_table_save(z,NULL,NULL,"out_z2.fits",0));
02996 
02997   //z[where(int_sky < max(int_sky[where(finite(int_sky))])/4)] = 0;
02998   // take in consideration only strong sky lines (set else to 0)
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   //flag too low values
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   //check_nomsg(cpl_table_save(z,NULL,NULL,"out_z4.fits",0));
03016   //computes gradient
03017   //z_diff = z[0:n_elements(z)-2] - z[1:n_elements(z)-1];
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   //check_nomsg(cpl_table_save(z_diff,NULL,NULL,"out_z_diff.fits",0));
03027 
03028   //identify points positions at which there is a line pick
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   //check_nomsg(cpl_table_save(z_diff,NULL,NULL,"out_z_diff.fits",0));
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   //check_nomsg(cpl_table_save(z_pos,NULL,NULL,"out_z_pos.fits",0));
03048 
03049   //Do a gaussian fit in a range of size 2*zext centered at
03050   //each line maximum position (fit the line) to get in corresponding arrays: 
03051   // 1) line lambda position of object and sky
03052   // 2) line object -sky intensity
03053   // 3) line object-sky intensity error 
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   //sinfo_msg("npos=%d z_ext=%d",npos,z_ext);
03083   //sinfo_table_column_dump(z_pos,"POS",CPL_TYPE_INT);
03084   //check_nomsg(cpl_table_save(z_pos,NULL,NULL,"out_z_pos_0.fits",0));
03085 
03086   for (jz=0;jz<npos;jz++) {
03087     check_nomsg(z1 = cpl_table_get_int(z_pos,"POS",jz,&status));
03088     //sinfo_msg("z1=%d",z1);
03089     // AMO added if check to prevent array explosion 
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         //sinfo_msg("wc_s=%f",wc_s);
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           sinfo_msg("Gauss fit parameters:");
03130           sinfo_msg("wc_s=%f sig_s=%f area_s=%f bkg_s=%f",
03131                      wc_s,sig_s,area_s,bkg_s);
03132           sinfo_msg("mse=%f chired=%f amp_s=%f",
03133                      mse,chired,amp_s);
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           sinfo_msg("Gauss fit parameters:");
03165           sinfo_msg("wc_o=%f sig_o=%f area_o=%f bkg_o=%f",
03166                      wc_o,sig_o,area_o,bkg_o);
03167           sinfo_msg("mse=%f chired=%f amp_o=%f",
03168                      mse,chired,amp_o);
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           if (lambda[z1] < 2.35 &&
03185              total(finite([l1,s1,o1])) == n_elements([l1,s1,o1])) {
03186             gs1 = float(gaussfit(l1,s1,as1,nterms=3));
03187             go1 = float(gaussfit(l1,o1,ao1,nterms=3));
03188             g_lam = [g_lam,(as1[1]+ao1[1])/2.];
03189             g_diff = [g_diff,as1[1]-ao1[1]];
03190             g_err = [g_err,sqrt(as1[2]^2+ao1[2]^2)];
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   //sinfo_table_column_dump(z_pos,"WAVEC",CPL_TYPE_DOUBLE);
03217   //sinfo_table_column_dump(z_pos,"STATUS",CPL_TYPE_INT);
03218   //sinfo_table_column_dump(z_pos,"WDIF",CPL_TYPE_DOUBLE);
03219 
03220   check_nomsg(cpl_table_and_selected_int(z_pos,"STATUS",CPL_GREATER_THAN,-1));
03221 
03222   //check_nomsg(cpl_table_save(z_pos,NULL,NULL,"out_z_pos.fits",0));
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   //check_nomsg(cpl_table_save(z_pos,NULL,NULL,"out_z_pos.fits",0));
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   //check_nomsg(cpl_table_save(z_good,NULL,NULL,"out_z_good.fits",0));  
03258   //Do a kappa-sigma clip of the differences of line positions 
03259   //as determined in the object and in the sky spectrum
03260 
03261   sinfo_msg("ks-clip1");
03262   sinfo_msg("iter mean (um) sdv (um) mean (pix) sdv (pix)"); 
03263   //sinfo_table_column_dump(z_pos,"WAVEC",CPL_TYPE_DOUBLE);
03264   //sinfo_table_column_dump(z_pos,"WDIF",CPL_TYPE_DOUBLE);
03265 
03266   for (iq = 0;iq<XCOR_YSHIFT_KS_CLIP;iq++) {
03267     //sinfo_msg("nval=%d",cpl_table_get_nrow(z_pos));
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     //z_good = where(abs(g_diff-z_mean) <= 2*z_sdv);
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     //sinfo_msg("ngood=%d",cpl_table_get_nrow(z_good));
03282     check_nomsg(cpl_table_select_all(z_pos));
03283     //z_mean = median(g_diff[z_good]);
03284     //z_sdv = stddev(g_diff[z_good]);
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   /* do a poly fit of wdif versus wave*/
03292   /*
03293   for (iq = 0; iq<3; iq++) {
03294     // sinfo_msg("%d %f %f",iq,mean(zfit),zsdv); 
03295     par1 = poly_fit(g_lam[z_good],g_diff[z_good],poly_n);
03296     z_fit = g_diff*0.;
03297     for (ii=0;ii<poly_n) z_fit = z_fit + par1[ii]*g_lam^ii;
03298     z_res = g_diff-z_fit;
03299     z_sdv = stddev(z_res[zgood]);
03300     z_good = where(abs(z_res) le 3*z_sdv);
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   //check_nomsg(cpl_table_save(z_pos,NULL,NULL,"out_z_pos2.fits",0));
03308   check_nomsg(z_good=cpl_table_duplicate(z_pos));
03309 
03310   //Do a fit of a uniform function to the residuals line position differences
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   //sinfo_table_column_dump(z_pos,"WDIF",CPL_TYPE_DOUBLE);
03315   for(iq=0;iq<XCOR_YSHIFT_KS_CLIP;iq++) {
03316     //cpl_table_dump(z_pos,0,cpl_table_get_nrow(z_pos),stdout);
03317     check_nomsg(nfit=cpl_table_get_nrow(z_good));
03318     //sinfo_msg("nfit=%d",nfit);
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     //sinfo_msg("coeff 0=%g um %g pix",zfit,zfit/dispersion);
03330 
03331     //computes residuals=difference-fit and their standard deviation
03332     //and then do a kappa-sigma clip of outliers (out of 3 sigma)
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     //cpl_table_dump(z_pos,0,cpl_table_get_nrow(z_pos),stdout);
03351     /*
03352     sinfo_msg("min=%g max=%g ndat=%d",
03353               cpl_table_get_column_min(z_pos,"WRES"),
03354               cpl_table_get_column_max(z_pos,"WRES"),
03355               cpl_table_get_nrow(z_pos));
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   //sinfo_msg(">>mean=%g",cpl_table_get_column_mean(z_good,"WDIF"));
03376 
03377   //check_nomsg(cpl_table_save(z_good,NULL,NULL,"out_z_pos3.fits",0));
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   //sinfo_msg("no=%d ns=%d nw=%d",no,ns,nw); 
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 static void
03494 sinfo_shift_sky(const int x,const int y)
03495 {
03496 
03497   //To remove compilation warnings
03498   ck0_nomsg(x);
03499   ck0_nomsg(y);
03500 
03501   // shift sky spectrum of a given amount
03502   if (max(abs(z_fit))/cdelts < 0.01) {
03503     sinfo_msg("shift <0.01 pixels will not be applied");
03504   } else {
03505     sinfo_msg("shifting sky cube by mean of %f pix wrt object",
03506              cpl_table_column_mean(z_fit,"VALUE")/cdelto);
03507     sinfo_msg("this will take a couple of minutes...");
03508     z_good = where(finite(int_sky));
03509     new_sky = spline(lambda[z_good]-z_mean,int_sky[z_good],lambda);
03510     int_sky = new_sky;
03511     sky_out = dblarr(xsize,ysize,zsize) + !values.f_nan;
03512     for (ix=0; ix<xsize;ix++) {
03513       for (iy=0;iy<ysize;iy++) {
03514     old_sky = reform(sky[ix,iy,*]);
03515     z_good = where(finite(old_sky),z_good_i);
03516     if (z_good_i > 0) {
03517           new_sky= spline(lambda[z_good]-z_fit[z_good],old_sky[z_good],lambda);
03518           new_fin= where(finite(new_sky,/infinity) ||
03519                           finite(old_sky,/nan),newfin_i);
03520       if (new_fin_i > 0) new_sky[new_fin] = !values.f_nan;
03521       sky_out[ix,iy,*] = new_sky;
03522     }
03523       }
03524     }
03525     sky = sky_out;
03526   }
03527  cleanup:
03528   return;
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; //full width in pixels of unresolved emission line
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   //Rotational parameters
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   // For each range extract proper: obj, sky, wave spectra
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     //AMO: the following 2 seems to be critical for robustness
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         // identify sky lines
03818         //sinfo_msg("finite_pix_i=%d",finite_pix_i); 
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         //sinfo_msg("sky_thresh=%f",sky_thresh);
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       //separate line and continuum regions 
03836           //by convolving with a hat region of large as a line 
03837           //sinfo_msg("xxx1_i=%d",xxx1_i); 
03838       check_nomsg(xxx2 = cpl_table_duplicate(sky_lr));
03839           //check_nomsg(cpl_table_save(sky_lr,NULL,NULL,"out_skylr.fits",0));
03840           //check_nomsg(cpl_table_save(obj_lr,NULL,NULL,"out_objlr.fits",0));
03841           //check_nomsg(cpl_table_save(xxx2,NULL,NULL,"out_xxx2_0.fits",0));
03842           ck0_nomsg(sinfo_table_threshold(&xxx2,"INT",sky_thresh,
03843                                           sky_thresh,0.,10.));
03844           //check_nomsg(cpl_table_save(xxx2,NULL,NULL,"out_xxx2_1.fits",0));
03845 
03846           
03847           /* TODO
03848       xxx2[xxx1] = 10.;
03849       */
03850           //sinfo_msg("npixw/2=%d",npixw);
03851       check_nomsg(sinfo_convolve_kernel(&xxx2,npixw/2));
03852           //check_nomsg(cpl_table_save(xxx2,NULL,NULL,"out_xxx2_2.fits",0));
03853 
03854           // get line_regions
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           // get cont_regions
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             //If we have enough line points and continuum points
03876             //sinfo_msg("line_i=%d cont_i=%d",line_i,cont_i);
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         // Fit the object profile='fline_res' of the sky line residuals
03893             // left after proper scaled sky spectrum lines (and continua) 
03894             // subtraction. Then determines median and stdev to flag outliers  
03895 
03896         //Free memory for each loop
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         //Identify obj lines and continuum, same for sky
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             check_nomsg(cpl_table_save(line_regions,NULL,NULL,
03908                         "out_line.fits",0));
03909             check_nomsg(cpl_table_save(cont_regions,NULL,NULL,
03910                         "out_cont.fits",0));
03911             check_nomsg(cpl_table_save(obj_cont,NULL,NULL,
03912                         "out_obj_cont.fits",0));
03913             check_nomsg(cpl_table_save(obj_line,NULL,NULL,
03914                         "out_obj_line.fits",0));
03915             check_nomsg(cpl_table_save(sky_line,NULL,NULL,
03916                         "out_sky_line.fits",0));
03917             check_nomsg(cpl_table_save(sky_cont,NULL,NULL,
03918                         "out_sky_cont.fits",0));
03919        
03920         */
03921 
03922 
03923 
03924 
03925 
03926 
03927             sinfo_free_table(&fline_res);
03928             //FIXME: in some cases obj_cont is empty
03929             //sinfo_msg("first line ratio determination");
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             // check if there are outliers
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               // do a k-sigma clip to select a better line region
03961               //sinfo_msg("fclip_i=%d",fclip_i);
03962               //Find again line_regions
03963               check_nomsg(line_i=cpl_table_and_selected_double(fline_res,
03964                                                                 "AINT",
03965                                                        CPL_NOT_GREATER_THAN,
03966                                                                 fmed+3*fsdv));
03967           // get new (better) line_regions
03968           sinfo_free_table(&line_regions);
03969               //sinfo_msg("line_i=%d",line_i);
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           //check_nomsg(cpl_table_save(obj_lr,NULL,NULL,"out_obj_lr.fits",0));
03978           //check_nomsg(cpl_table_save(line_regions,NULL,NULL,"out_line_regions.fits",0));
03979        
03980 
03981 
03982 
03983         // The following 2 may return an error so we do not check and
03984         // later we reset the error
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               //sinfo_msg("fline_i=%d",fline_i);
03990               if(fline_i>=3) {
03991                 // repeat the determination of the line ratio
03992                 //sinfo_msg("second line ratio determination");
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       } /* end if line_i */
04015         } /* end if xxx1_i */
04016      } /* end finite_pix_i */
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   } /* end for loop on i */
04030 
04031   sinfo_free_array(&do_hk);
04032   sinfo_free_array(&rfit);
04033 
04034   //sinfo_msg("n scale=%d",cpl_table_get_nrow(rscale0));
04035   //check_nomsg(cpl_table_save(rscale0,NULL,NULL,
04036   //                                   "out_rscale0.fits",0));
04037   check_nomsg(cpl_table_select_all(rscale0));
04038  /* TODO: here one has to implementa an interpol function 
04039   check_nomsg(range0_i=cpl_table_and_selected_double(rscale0,"RATIO",
04040                                                      CPL_NOT_EQUAL_TO,0));
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   check_nomsg(cpl_table_duplicate_column(*int_obj,"INTC",*int_obj,"INT"));
04053   check_nomsg(cpl_table_duplicate_column(*int_obj,"SKYC",*int_sky,"INTC"));
04054   check_nomsg(cpl_table_subtract_columns(*int_obj,"INTC","SKYC"));
04055   */
04056 
04057   // do simple rotational correction
04058   if (do_rot == 1) {
04059 
04060     //finite_pix = where(finite(int_sky) && finite(int_obj),finite_pix_i);
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     //sinfo_msg("min_lrange=%g max_lrange=%g",min_lrange,max_lrange);
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     //TODO: lambda invalid values need to be reset to valid (?)
04071     
04072     check_nomsg(cpl_table_erase_invalid(finite_pix));
04073  
04074  
04075     if (finite_pix_i > npixw) {
04076 
04077       //finite_pix = finite_pix[where(finite_pix > min(lrange) && 
04078       //                              finite_pix < max(lrange))];
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       //check_nomsg(cpl_table_save(rat_sky,NULL,NULL,"out_rat_sky.fits",0));
04108       //check_nomsg(cpl_table_save(finite_pix,NULL,NULL,"out_finite_pix.fits",0));
04109       //check_nomsg(cpl_table_save(sky_lr,NULL,NULL,"out_sky_lr.fits",0));
04110       //check_nomsg(cpl_table_save(obj_lr,NULL,NULL,"out_obj_lr.fits",0));
04111       //check_nomsg(cpl_table_save(wav_lr,NULL,NULL,"out_wav_lr.fits",0));
04112 
04113       //The following may fail (sky_lr may be empty) so we do not check
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         //xxx1 = where(sky_lr > median(sky_lr)+stddev(sky_lr),xxx1_i);
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         //check_nomsg(cpl_table_save(wav_lr,NULL,NULL,"out_llr_xxx1.fits",0));
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     //check_nomsg(cpl_table_dump(low_pos,0,low_pos_i,stdout));
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     //check_nomsg(cpl_table_dump(med_pos,0,med_pos_i,stdout));
04148  
04149     //TODO:
04150         //hipos = [0]
04151         //for i=0,n_elements(xxx1)-1 do begin
04152         //    x1 = where(lowpos eq i,x1_i)
04153         //    x2 = where(medpos eq i,x2_i)
04154         //    if (x1_i eq 0 and x2_i eq 0) then hipos = [hipos,i]
04155         //endfor
04156         //hipos = hipos[1:n_elements(hipos)-1]
04157     //TODO: hi_pos=sinfo_find_rot_waves(w_rot_hi,npixw,wtol,wav_lr);
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     //check_nomsg(cpl_table_dump(hi_pos,0,hi_pos_i,stdout));
04164 
04165 
04166         //xxx2[xxx1] = 10.;
04167         check_nomsg(xxx2 = cpl_table_duplicate(sky_lr));
04168         check_nomsg(nrow=cpl_table_get_nrow(sky_lr));
04169         //check_nomsg(cpl_table_save(xxx2,NULL,NULL,"out_xxx1.fits",0));
04170         //check_nomsg(cpl_table_save(xxx2,NULL,NULL,"out_xxx2_0.fits",0));
04171 
04172     // AMO: Why the following?
04173         //check_nomsg(cpl_table_fill_column_window(xxx2,"INT",0,nrow,0));
04174 
04175         //xxx2 = convol(xxx2,replicate(1,npixw),/edge_truncate,/center);
04176         //cont_regions = where(xxx2 == 0,cont_i);
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         //check_nomsg(cpl_table_save(xxx2,NULL,NULL,"out_xxx2_1.fits",0));
04181     check_nomsg(sinfo_convolve_kernel(&xxx2,npixw/2));
04182 
04183         //check_nomsg(cpl_table_save(xxx2,NULL,NULL,"out_xxx2_2.fits",0));
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         sinfo_msg("xxx1: wmin=%g wmax=%g", 
04203                    cpl_table_get_column_min(xxx1,"WAVE"),
04204                    cpl_table_get_column_max(xxx1,"WAVE"));
04205 
04206         sinfo_msg("low_pos: wmin=%g wmax=%g", 
04207                    cpl_table_get_column_min(low_pos,"WAVE"),
04208                    cpl_table_get_column_max(low_pos,"WAVE"));
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       //compute line ratio for hi_regions
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       //compute line ratio for med_regions
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       //compute line ratio for med_regions
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       } //xxx1_i > 0
04247     }//finitepix > npixw
04248   }//do_rot==1
04249   //end of new rotational bit
04250 
04251 
04252   check_nomsg(cpl_table_duplicate_column(*int_sky,"INTC",*int_sky,"INT"));
04253   //sinfo_msg("n sky=%d",cpl_table_get_nrow(*int_sky));
04254   //sinfo_msg("n scale=%d",cpl_table_get_nrow(*rscale));
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   /* Get Object relevant information */
04497   /* here one should scan the inp image constructing a wave range from it
04498      and not from another table */
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     // insert cubes with wavelengths with appropriate values only
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       pval[k++]=pim[j*sx+i];
04617       sinfo_msg("set tab %f",pim[j*sx+i]);
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   // check for screwy values at ends of spectrum
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   //gsky_pos = where(int_sky > 1.*gsky_max || int_sky < 1.*gsky_min,gskypos_i);
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   //set selected values to NAN
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   //set selected values to NAN
04686   ck0_nomsg(sinfo_table_set_column_invalid(int_sky,"INT"));
04687   //count selected values (gsky_pos table is not really necessary)
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   //check_nomsg(cpl_table_save(gsky_pos,NULL,NULL,"out_gsky_pos.fits",0));
04692   check_nomsg(gsky_pos_i=cpl_table_get_nrow(gsky_pos));
04693   //sinfo_msg("ns1=%d,ns2=%d",ns1,ns2);
04694   //sinfo_msg("gsky_pos_i=%d",gsky_pos_i);
04695 
04696   sinfo_free_table(&gsky_pos);
04697   //check_nomsg(cpl_table_save(*int_sky,NULL,NULL,"out_sky.fits",0));
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   //check_nomsg(cpl_table_save(xobj,NULL,NULL,"out_xobj.fits",0));
04705   gobj_max = (obj_max>0) ? obj_max : 0;
04706   gobj_min = (obj_min<0) ? obj_min : 0;
04707   //gobj_pos=where(int_obj > 1.*gobj_max || int_obj < 1.*gobj_min,gobj_pos_i);
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   //set selected values to NAN
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   //set selected values to NAN
04718   ck0_nomsg(sinfo_table_set_column_invalid(int_obj,"INT"));
04719   //count selected values (gsky_pos table is not really necessary)
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   //check_nomsg(cpl_table_save(gobj_pos,NULL,NULL,"out_gobj_pos.fits",0));
04724   sinfo_free_table(&tbl1);
04725   check_nomsg(gobj_pos_i=cpl_table_get_nrow(gobj_pos));
04726   //sinfo_msg("no1=%d,no2=%d",no1,no2);
04727   //sinfo_msg("gobj_pos_i=%d",gobj_pos_i);
04728   sinfo_free_table(&gobj_pos);
04729   sinfo_msg("gobj_min=%f gobj_max=%f",gobj_min,gobj_max);
04730   //check_nomsg(cpl_table_save(*int_obj,NULL,NULL,"out_obj.fits",0));
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   sinfo_msg("ninp=%d nref=%d",ninp,nref); 
04787   sinfo_msg("winp(0)=%f wref(0)=%f tol=%f diff=%f",
04788             pwin[0],pwrf[0],tol,fabs(pwin[0]-pwrf[0])); 
04789   */
04790   if(pwin[0]<=pwrf[0]) {
04791     //sinfo_msg("case 1");
04792     for(i=0;i<ninp;i++) {
04793       if(k<nref) {
04794     /*
04795         sinfo_msg("case1: %f %f %f %f %d %d",
04796                   fabs(pwin[i] - pwrf[k]),tol,pwin[i],pwrf[k],i,k);
04797     */
04798     if(fabs(pwin[i] - pwrf[k])< tol) {
04799       pcin[i]=val;
04800       k++;
04801     }
04802       }
04803     }
04804   } else {
04805 
04806     //pwin[0]>pwrf[0]
04807     //sinfo_msg("case 2");
04808     for(k=0;k<nref;k++) {
04809       if(i<ninp) {
04810     /*
04811         sinfo_msg("case2: %f %f %f %f %d %d",
04812                   fabs(pwin[i] - pwrf[k]),tol,pwin[i],pwrf[k],i,k);
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     //sinfo_msg_warning("ninp != nref");
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   sinfo_msg("nout=%d nref=%d",nout,nref); 
04892   sinfo_msg("wout(0)=%f wref(0)=%f tol=%f diff=%f",
04893             pou[0],prf[0],tol,fabs(pou[0]-prf[0])); 
04894   */
04895   if(pou[0]<=prf[0]) {
04896     //sinfo_msg("case 1");
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     //pou[0]>prf[0]
04908     //sinfo_msg("case 2");
04909     for(k=0;k<nref;k++) {
04910       if(i<nout) {
04911     /*
04912         sinfo_msg("check: %f %f %f %f",
04913                   fabs(pou[i] - prf[k]),tol,pou[i],prf[k]);
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   //check_nomsg(cpl_table_save(tmp,NULL,NULL,"out_pre0.fits",0));
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   //check_nomsg(cpl_table_save(out,NULL,NULL,"out_post0.fits",0));
04931   /* sinfo_msg("nout=%d",nout); */
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   //interpolate sky if necessary
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   //TODO
04977   cpl_table* out=NULL;
04978 
04979   //To remove compilation warnings
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  /* Get Object relevant information */
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   cknull_nomsg(*int_sky=sinfo_table_shift_column_int(int_sky_dup, 
05085                                                      "INT", zshift,&zrest));
05086   check_nomsg(cpl_table_save(*int_sky, NULL, NULL, "out_sky_shift1.fits", 0));
05087   sinfo_free_table(&(*int_sky));
05088   
05089   sinfo_msg("min=%f max=%f",min,max);
05090   check_nomsg(cpl_table_save(int_sky_dup, NULL, NULL, "out_sky_pre2.fits", 0));
05091   cknull_nomsg(*int_sky=sinfo_table_shift_column_poly(int_sky_dup, 
05092                                                       "INT", zshift,2));
05093   check_nomsg(cpl_table_save(*int_sky, NULL, NULL, "out_sky_shift2.fits", 0));
05094   */
05095   //check_nomsg(cpl_table_save(int_sky_dup, NULL, NULL, "out_sky_pre2.fits", 0));
05096   check_nomsg(*int_sky=sinfo_table_shift_simple(int_sky_dup,"INT",zshift));
05097   /*
05098     sinfo_free_table(&(*int_sky));
05099   cknull_nomsg(*int_sky=sinfo_table_shift_column_spline3(int_sky_dup, 
05100                                                          "INT", zshift));
05101   check_nomsg(cpl_table_save(*int_sky, NULL, NULL, "out_sky_shift3.fits", 0));
05102   */
05103   sinfo_free_table(&int_sky_dup);
05104   /*
05105   check_nomsg(cpl_table_select_all(*int_sky));
05106   check_nomsg(n=cpl_table_and_selected_double(*int_sky,"INT",CPL_LESS_THAN,min));
05107   ck0_nomsg(sinfo_table_set_column_invalid(int_sky,"INT"));
05108   sinfo_msg("n=%d",n);
05109   check_nomsg(cpl_table_select_all(*int_sky));
05110   check_nomsg(n=cpl_table_and_selected_double(*int_sky,"INT",CPL_GREATER_THAN,max));
05111   sinfo_msg("n=%d",n);
05112   ck0_nomsg(sinfo_table_set_column_invalid(int_sky,"INT"));
05113   check_nomsg(cpl_table_select_all(*int_sky));
05114   */
05115   //check_nomsg(cpl_table_save(*int_sky, NULL, NULL, "out_sky_shift3.fits", 0));
05116 
05117 
05118 
05119   check_nomsg(sky_shift=sinfo_cube_zshift_simple(sky_cub,(float)zshift));
05120 
05121   //check_nomsg(sky_shift=sinfo_cube_zshift(sky_cub,zshift,&zrest));
05122   //check_nomsg(cpl_imagelist_save(sky_shift,"out_sky1.fits",
05123   //                 CPL_BPP_DEFAULT,NULL,CPL_IO_DEFAULT));
05124 
05125   sinfo_free_imagelist(&sky_shift);
05126   //sinfo_free_imagelist(&sky_shift);
05127   //sinfo_msg("zrest=%f",zrest);
05128 
05129   //check_nomsg(sky_shift=sinfo_cube_zshift_poly(sky_cub,zshift,2));
05130   //check_nomsg(cpl_imagelist_save(sky_shift,"out_sky2.fits",
05131   //                               CPL_BPP_DEFAULT,NULL,CPL_IO_DEFAULT));
05132   //sinfo_free_imagelist(&sky_shift);
05133 
05134   //check_nomsg(sky_shift=sinfo_cube_zshift_spline3(sky_cub,zshift));
05135   //check_nomsg(cpl_imagelist_save(sky_shift,"out_sky3.fits",
05136   //                               CPL_BPP_DEFAULT,NULL,CPL_IO_DEFAULT));
05137   //sinfo_free_imagelist(&sky_shift);
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   /* set to 0 edges */
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     /*val*=dw; */
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   /* set to 0 edges */
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     /*val*=dw; */
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   /* set to 0 edges */
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     /*val*=dw; */
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   /* set to 0 edges */
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     /*val*=dw; */
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   int n=sizeof(x)/sizeof(double);
05494   double fmin = sinfo_fac(x[0],a[2]);
05495   double fmax = sinfo_fac(x[n-1],a[2]);
05496   sinfo_msg("n=%d",n);
05497   if(fmax < 0) sinfo_msg("fmax=%f",fmax);
05498   */
05499   //*result = a[0]+a[1]*fac/sinfo_scale_fct;
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   int n=sizeof(x)/sizeof(double);
05536   double fmin = sinfo_fac(x[0],a[2]);
05537   double fmax = sinfo_fac(x[n],a[2]);
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   //f2=f1+da*a[0];
05545   //f2=a[0]+(a[1]+da*a[1])*fac;
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   //sinfo_msg("d0=%g d1=%g d2=%g",d[0]*a[0]/f,d[1]*a[1]/f,d[2]*a[2]/f);
05551   //sinfo_msg("comp d1=%g",d[2]);
05552   //sinfo_msg("real d1=%g",(f2-f1)/(da*a[2]));
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   //return pow(x,-5.)/(exp(c/(x*fabs(t)))-1.)/sinfo_scale_fct;
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     /*sinfo_msg("check =%g thresh=%g",fabs(piw[i]-prw[k]),tol); */
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   //sinfo_msg("Object sky residuals/Sky lines ratio determination method=%d",
05869   //          method);
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   //sinfo_msg("mden=%g tden=%g",mden,tden);
05908   //sinfo_msg("mnum=%g tnum=%g",mnum,tnum);
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   // Amoeba part
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   //sinfo_msg("Before amoeba fit");
05989   //sinfo_msg("ap[0][0]=%g ap[0][1]=%g",ap[0][0],ap[1][0]);
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   //fit residual obj continuum and subtract it
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   //fit residual sky continuum and subtract it
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   //computes diff=(obj-conto)-(sky-contsky)*p[0]
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   //computes rms diff
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   //sinfo_msg("nlin=%d",nlin);
06176   check_nomsg(cpl_table_fill_column_window(out,"CONT",0,nlin,0));
06177 
06178   //do a uniform fit 
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   //sinfo_msg("coeff 0=%g",yfit);
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   //line_regions = med_regions;
06236   check_nomsg(line_regions = cpl_table_duplicate(sel_regions));
06237   //r = amoeba(1.e-5,function_name='fitsky',p0=[1.],scale=[0.1]);
06238   //Identify obj lines and continuum, same for sky
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   //fline_res = (obj_lr[line_regions]-
06248   //             interpol(obj_lr[cont_regions],llr[cont_regions],
06249   //             llr[line_regions])) -
06250   //            (sky_lr[line_regions] - 
06251   //             interpol(sky_lr[cont_regions],llr[cont_regions],
06252   //
06253   //            llr[line_regions]))*r[0];
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   //fclip = where(abs(fline_res) > fmed+3*fsdv,fclip_i);
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     //line_regions = line_regions[where(abs(fline_res) < fmed+3*fsdv)];
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      //r = amoeba(1.e-5,function_name='fitsky',p0=[1.],scale=[0.1]);
06293       ck0_nomsg(sinfo_get_line_ratio(obj_lin,obj_cnt,sky_lin,sky_cnt,meth,r));
06294     }
06295   }
06296   *r=fabs(*r);
06297   //Free memory 
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     //x = where(lambda > l_rot_low[i]-npixw*cdelto &&
06358     //          lambda < l_rot_low[i]+npixw*cdelto,x_i);
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       //sinfo_msg("i=%d x_i=%d w_min=%g w_max=%g",i,x_i,w_min,w_max);
06374       check_nomsg(cpl_table_insert(res,w_sel,r_start));
06375     }
06376     check_nomsg(cpl_table_select_all(range));
06377   }
06378 
06379   //res = range[1:cpl_table_get_nrow(res)-1];
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   //x2[x1[pos]] = 10.;
06453   //x2 = convol(x2,replicate(1,npixw),/edge_truncate,/center);
06454   //res = where(x2 > 0,hi_i);
06455   //cpl_table_save(x1, NULL, NULL, "out_x1.fits", 0);
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   //check_nomsg(cpl_table_save(low,NULL,NULL,"out_low.fits",0));
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 

Generated on Wed Jan 17 08:33:44 2007 for SINFONI Pipeline Reference Manual by  doxygen 1.4.4