sinfo_new_psf.c

00001 /*
00002  * This file is part of the ESO SINFONI Pipeline
00003  * Copyright (C) 2004,2005 European Southern Observatory
00004  *
00005  * This program is free software; you can redistribute it and/or modify
00006  * it under the terms of the GNU General Public License as published by
00007  * the Free Software Foundation; either version 2 of the License, or
00008  * (at your option) any later version.
00009  *
00010  * This program is distributed in the hope that it will be useful,
00011  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00012  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00013  * GNU General Public License for more details.
00014  *
00015  * You should have received a copy of the GNU General Public License
00016  * along with this program; if not, write to the Free Software
00017  * Foundation, 51 Franklin St, Fifth Floor, Boston, MA  02111-1307  USA
00018  */
00019 /*----------------------------------------------------------------------------
00020  
00021      File name    :     sinfo_new_psf.c
00022    Author       :   A. Modigliani
00023    Created on   :   Sep 17, 2003
00024    Description  : 
00025 
00026  sinfo_new_psf.py does the image reconstruction of a set of 
00027  sky-subtracted, flatfielded, 
00028  bad pixel corrected and slope of the spectra aligned exposures of a bright 
00029  star with continuum spectrum. The resulting image can be used to determine 
00030  the PSF
00031  ---------------------------------------------------------------------------*/
00032 #ifdef HAVE_CONFIG_H
00033 #  include <config.h>
00034 #endif
00035 
00036 /*----------------------------------------------------------------------------
00037                                 Includes
00038  ---------------------------------------------------------------------------*/
00039 #include <irplib_utils.h>
00040 #include <irplib_error.h>
00041 #include <assert.h>
00042 #include "sinfo_new_psf.h"
00043 #include "sinfo_pro_save.h"
00044 #include "sinfo_key_names.h"
00045 #include "sinfo_psf_ini.h"
00046 #include "sinfo_psf_ini_by_cpl.h"
00047 #include "sinfo_utilities_scired.h"
00048 #include "sinfo_hidden.h"
00049 #include "sinfo_pfits.h"
00050 #include "sinfo_functions.h"
00051 #include "sinfo_error.h"
00052 #include "sinfo_utils_wrappers.h"
00053 #include "sinfo_globals.h"
00054 
00055 /*----------------------------------------------------------------------------
00056                                 Defines
00057  ---------------------------------------------------------------------------*/
00058 
00059 #define SINFO_STREHL_M1                       8.0
00060 #define SINFO_STREHL_M2                       1.1
00061 #define SINFO_STREHL_BOX_SIZE                 64
00062 #define SINFO_STREHL_WINDOW                   6
00063 #define SINFO_PSF_SZ                          3
00064 #define SINFO_PI                   3.1415926535897932384626433
00065 #define SINFO_STREHL_ERROR_COEFFICIENT    SINFO_PI * 0.007 / 0.0271
00066 /*----------------------------------------------------------------------------
00067                              Function Definitions
00068  ---------------------------------------------------------------------------*/
00069 static void 
00070 sinfo_check_borders(int* val,const int max,const int thresh);
00071 
00072 static void 
00073 sinfo_get_safe_box(int* llx, 
00074                    int* lly, 
00075                    int* urx, 
00076                    int* ury, 
00077                    const int xpos,
00078                    const int ypos,
00079                    const int box,
00080                    const int szx, 
00081                    const int szy);
00082 
00083 static int 
00084 sinfo_get_strehl_from_slice(cpl_imagelist* cube,
00085                             double disp, 
00086                             double cWave,
00087                             double ws,
00088                             double we,
00089                             double pscale,
00090                             double strehl_star_radius,
00091                             double strehl_bg_r1,
00092                             double strehl_bg_r2,
00093                             double* strehl,
00094                             double* strehl_err);
00095 
00096 
00097 static cpl_table* 
00098 sinfo_get_encircled_energy(cpl_frameset* sof,
00099                            cpl_image* img, 
00100                            double* fwhm_x, 
00101                double* fwhm_y,
00102                            cpl_table** qclog);
00103 
00104 static double 
00105 sinfo_get_strehl_from_ima(cpl_image* ima,
00106                           char* name, 
00107                           cpl_frame* frame);
00108 
00109 static int 
00110 sinfo_get_strehl_from_image(cpl_image* img,
00111                             double ws,
00112                             double we,
00113                             double pscale,
00114                             double strehl_star_radius,
00115                             double strehl_bg_r1,
00116                             double strehl_bg_r2,
00117                             double* strehl,
00118                 double* strehl_err);
00119 
00120 static cpl_table* 
00121 sinfo_get_strehl_from_cube(cpl_imagelist* cube,
00122                            char* name, 
00123                            cpl_frame* frame);
00124 
00125 static int 
00126 sinfo_get_frm12(cpl_frameset* sof,cpl_frame** frm1,cpl_frame** frm2);
00127 
00128 static cpl_table* 
00129 sinfo_get_strehl_from_2cubes(cpl_imagelist* cube1, 
00130                              cpl_imagelist* cube2,
00131                              const char* fname1, 
00132                              cpl_frame* frm1,
00133                              cpl_frame* frm2);
00134 
00135 static int sinfo_strehl_compute(
00136         const cpl_image *   im1,
00137         const cpl_image *   im2,
00138         double              m1,
00139         double              m2,
00140         double              lam,
00141         double              dlam,
00142         double              pscale1,
00143         double              pscale2,
00144         int                 size,
00145         int                 xpos1,
00146         int                 ypos1,
00147         int                 xpos2,
00148         int                 ypos2,
00149         double              r1,
00150         double              r2,
00151         double              r3,
00152         int                 noise_box_sz,
00153         int                 noise_nsamples,
00154         double          *   strehl,
00155         double          *   strehl_err,
00156         double          *   star_bg,
00157         double          *   star_peak,
00158         double          *   star_flux,
00159         double          *   psf_peak,
00160         double          *   psf_flux,
00161         double          *   bg_noise);
00162 
00171 /*----------------------------------------------------------------------------
00172    Function     :       sinfo_new_psf()
00173    In           :       ini_file: file name of according .ini file
00174    Out          :       integer (0 if it worked, -1 if it doesn't) 
00175    Job          : 
00176 
00177  sinfo_new_psf.py does the image reconstruction of a set of sky-subtracted, 
00178  flatfielded, 
00179  bad pixel corrected and slope of the spectra aligned exposures of a bright 
00180  star with continuum spectrum. The resulting image can be used to determine 
00181  the PSF
00182 
00183  ---------------------------------------------------------------------------*/
00184 
00185 int 
00186 sinfo_new_psf (const char* plugin_id,
00187                cpl_parameterlist* config,
00188                cpl_frameset* sof)
00189 {
00190 
00191   cpl_imagelist* cube1=NULL;
00192   cpl_imagelist* cube2=NULL;
00193   cpl_image * med_img1=NULL ;
00194   cpl_image * med_img2=NULL ;
00195 
00196   cpl_table* ao_performance=NULL;
00197   cpl_table* enc_energy=NULL;
00198 
00199   cpl_frame* frm1=NULL;
00200   cpl_frame* frm2=NULL;
00201 
00202   cpl_table* qclog_tbl=NULL;
00203   cpl_frameset* stk=NULL;
00204   cpl_propertylist* plist =NULL;
00205 
00206   psf_config * cfg =NULL;
00207 
00208   int nsample=0;
00209   int i = 0;
00210   int status=0;
00211 
00212  
00213 
00214   int strehl_sw=0;
00215   int ilx1=0;
00216   int ily1=0;
00217   int ilx2=0;
00218   int ily2=0;
00219 
00220   float cx1=0;
00221   float cy1=0;
00222   float cx2=0;
00223   float cy2=0;
00224  
00225   double fwhm_x=0;
00226   double fwhm_y=0;
00227   double lam=0;
00228   double strehl=0;
00229  
00230   char fname1[MAX_NAME_SIZE];
00231   char fname2[MAX_NAME_SIZE];
00232 
00233   char key_name[MAX_NAME_SIZE];
00234 
00235   char obs_name1[MAX_NAME_SIZE];
00236   char hlamp_st='F';
00237   char shut2_st='F';
00238 
00239 
00240 
00241   /* 
00242        -----------------------------------------------------------------
00243        1) parse the file names and parameters to the psf_config data 
00244           structure cfg
00245        -----------------------------------------------------------------
00246   */
00247 
00248   sinfo_msg("Parsing cpl input");
00249   check_nomsg(stk=cpl_frameset_new());
00250 
00251   cknull(cfg = sinfo_parse_cpl_input_psf(sof,&stk),
00252       "error parsing cpl input");
00253 
00254   /* TODO the following generate a small leak of 72 bytes */
00255   strehl_sw=sinfo_get_strehl_type(sof);
00256 
00257   if(strehl_sw==0) {
00258     sinfo_msg("One target Strehl computation");
00259     if(sinfo_is_fits_file(cfg->inFrame) != 1) {
00260       sinfo_msg_error("Input file %s is not FITS",cfg->inFrame);
00261       goto cleanup;
00262     } else {
00263       strcpy(fname1,cfg->inFrame);
00264     }
00265 
00266      if(NULL != cpl_frameset_find(sof,PRO_COADD_PSF)) {
00267       frm1 = cpl_frameset_find(sof,PRO_COADD_PSF);
00268     } else if(NULL != cpl_frameset_find(sof,PRO_OBS_PSF)) {
00269       frm1 = cpl_frameset_find(sof,PRO_OBS_PSF);
00270     } else if(NULL != cpl_frameset_find(sof,PRO_COADD_STD)) {
00271       frm1 = cpl_frameset_find(sof,PRO_COADD_STD);
00272     } else if(NULL != cpl_frameset_find(sof,PRO_OBS_STD)) {
00273       frm1 = cpl_frameset_find(sof,PRO_OBS_STD);
00274     } else if(NULL != cpl_frameset_find(sof,PRO_COADD_OBJ)) {
00275       frm1 = cpl_frameset_find(sof,PRO_COADD_OBJ);
00276     } else if(NULL != cpl_frameset_find(sof,PRO_OBS_OBJ)) {
00277       frm1 = cpl_frameset_find(sof,PRO_OBS_OBJ);
00278     } else {
00279       sinfo_msg_error("Frame %s  or %s or %s  or %s or %s  or %s not found!",
00280               PRO_COADD_PSF,PRO_OBS_PSF,
00281               PRO_COADD_STD,PRO_OBS_STD,
00282               PRO_COADD_OBJ,PRO_OBS_OBJ);
00283       goto cleanup;
00284     }
00285 
00286     sinfo_get_obsname(frm1,obs_name1);
00287     check_nomsg(hlamp_st=sinfo_get_keyvalue_bool(frm1,KEY_NAME_LAMP_HALO));
00288     check_nomsg(shut2_st=sinfo_get_keyvalue_bool(frm1,KEY_NAME_SHUT2_ST));
00289  
00290     check_nomsg(cube1 = cpl_imagelist_load(fname1,CPL_TYPE_FLOAT,0));
00291     cknull(med_img1=sinfo_new_median_cube(cube1),
00292       " could not do sinfo_medianCube()");
00293 
00294     check_nomsg(ilx1=cpl_image_get_size_x(med_img1));
00295     check_nomsg(ily1=cpl_image_get_size_y(med_img1));
00296 
00297     cx1 = ilx1 / 2. + 0.5;
00298     cy1 = ily1 / 2. + 0.5;
00299     cknull(ao_performance=sinfo_get_strehl_from_cube(cube1,fname1,frm1),
00300        "error computing strehl");
00301     sinfo_free_imagelist(&cube1);
00302   } else {
00303     sinfo_msg("Two target Strehl computation");
00304     sinfo_get_frm12(sof,&frm1,&frm2);
00305     strcpy(fname1,cpl_frame_get_filename(frm1));
00306     strcpy(fname2,cpl_frame_get_filename(frm2));
00307 
00308     check_nomsg(cube1 = cpl_imagelist_load(fname1,CPL_TYPE_FLOAT,0));
00309     check_nomsg(cube2 = cpl_imagelist_load(fname2,CPL_TYPE_FLOAT,0));
00310     cknull(med_img1=sinfo_new_median_cube(cube1),"Computing median on cube");
00311     cknull(med_img2=sinfo_new_median_cube(cube2),"Computing median on cube");
00312 
00313 
00314     check_nomsg(ilx1=cpl_image_get_size_x(med_img1));
00315     check_nomsg(ily1=cpl_image_get_size_y(med_img1));
00316     check_nomsg(ilx2=cpl_image_get_size_x(med_img2));
00317     check_nomsg(ily2=cpl_image_get_size_y(med_img2));
00318 
00319     cx1 = ilx1 / 2. + 0.5;
00320     cy1 = ily1 / 2. + 0.5;
00321     cx2 = ilx2 / 2. + 0.5;
00322     cy2 = ily2 / 2. + 0.5;
00323     cknull(ao_performance=sinfo_get_strehl_from_2cubes(cube1,
00324                                cube2,
00325                                fname1,
00326                                frm1,
00327                                frm2),
00328        "Computing strehl");
00329 
00330     sinfo_free_imagelist(&cube1);
00331     sinfo_free_imagelist(&cube2);
00332 
00333   }
00334 
00335   /* STREHL computation */
00336 
00337   check_nomsg(nsample=cpl_table_get_nrow(ao_performance));
00338   cknull_nomsg(qclog_tbl = sinfo_qclog_init());
00339  
00340   /* QC log */
00341   cknull(plist = cpl_propertylist_load(fname1, 0),
00342      "getting header from reference ima frame %s",fname1);
00343 
00344   if (cpl_propertylist_contains(plist, KEY_NAME_LOOP_STATE)) {
00345     sinfo_qclog_add_string(qclog_tbl,KEY_NAME_LOOP_STATE,
00346                       cpl_propertylist_get_string(plist,KEY_NAME_LOOP_STATE),
00347                       KEY_HELP_LOOP_STATE,"%s");
00348   }
00349 
00350 
00351   if (cpl_propertylist_contains(plist, KEY_NAME_LOOP_LGS)) {
00352     sinfo_qclog_add_int(qclog_tbl,KEY_NAME_LOOP_LGS,
00353                       cpl_propertylist_get_int(plist,KEY_NAME_LOOP_LGS),
00354                       KEY_HELP_LOOP_LGS,"%d");
00355   }
00356 
00357   if (cpl_propertylist_contains(plist, KEY_NAME_INS1_MODE)) {
00358     sinfo_qclog_add_string(qclog_tbl,KEY_NAME_INS1_MODE,
00359                       cpl_propertylist_get_string(plist,KEY_NAME_INS1_MODE),
00360                       KEY_HELP_INS1_MODE,"%s");
00361   }
00362   sinfo_free_propertylist(&plist);
00363 
00364   check_nomsg(strehl=cpl_table_get_column_median(ao_performance,"strehl"));
00365 
00366   ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,"QC STREHL MED",strehl,
00367             "STREHL MEDIAN","%f"));
00368   strehl=sinfo_get_strehl_from_ima(med_img1,fname1,frm1);
00369 
00370   ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,"QC STREHL AVG",strehl,
00371             "STREHL AVERAGE","%f"));
00372 
00373   for(i=1;i<nsample;i++) {
00374 
00375     check_nomsg(strehl=cpl_table_get_double(ao_performance,"strehl",
00376                                             i,&status));
00377     snprintf(key_name,MAX_NAME_SIZE-1,"%s%d","QC STREHL",i);
00378     ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,key_name,strehl,"STREHL","%f"));
00379     
00380     check_nomsg(lam=cpl_table_get_double(ao_performance,"wavelength",
00381                                          i,&status));
00382     snprintf(key_name,MAX_NAME_SIZE-1,"%s%d","QC LAMBDA",i);
00383     ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,key_name,lam,
00384                                      "WAVELENGTH","%f"));
00385         
00386   }
00387 
00388   check_nomsg(strehl=cpl_table_get_column_median(ao_performance,
00389                                                  "strehl_error"));
00390   ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,"QC STREHL MEDERR",strehl,
00391                                    "STREHL ERROR MEDIAN","%f"));
00392   ck0_nomsg(sinfo_qclog_add_string(qclog_tbl,"OBS NAME",obs_name1,
00393                                    "OB name","%s"));
00394   ck0_nomsg(sinfo_qclog_add_bool(qclog_tbl,PAF_NAME_LAMP_HALO,hlamp_st,
00395                                   KEY_NAME_LAMP_HALO,"%d"));
00396   ck0_nomsg(sinfo_qclog_add_bool(qclog_tbl,PAF_NAME_SHUT2_ST,shut2_st,
00397                                   KEY_NAME_SHUT2_ST,"%d"));
00398 
00399   ck0(sinfo_pro_save_tbl(ao_performance,stk,sof,
00400                   PSF_AO_PERFORMANCE_OUT_FILENAME,
00401              PRO_AO_PERFORMANCE,qclog_tbl,plugin_id,config),
00402     "cannot save tbl %s", PSF_AO_PERFORMANCE_OUT_FILENAME);
00403 
00404   sinfo_free_table(&qclog_tbl);
00405   sinfo_free_table(&ao_performance);
00406 
00407   /* Encircled energy & FWHM computation */
00408   cknull_nomsg(qclog_tbl=sinfo_qclog_init());
00409   cknull(enc_energy=sinfo_get_encircled_energy(sof,
00410                            med_img1,
00411                            &fwhm_x,
00412                            &fwhm_y,
00413                            &qclog_tbl),
00414            "Computing encircled energy");
00415 
00416   ck0(sinfo_pro_save_tbl(enc_energy,stk,sof,PSF_ENC_ENERGY_OUT_FILENAME,
00417              PRO_ENC_ENERGY,qclog_tbl,plugin_id,config),
00418       "cannot save tbl %s", PSF_ENC_ENERGY_OUT_FILENAME);
00419 
00420   sinfo_free_table(&qclog_tbl);
00421   sinfo_free_table(&enc_energy);
00422 
00423   /* QC log */
00424   cknull_nomsg(qclog_tbl = sinfo_qclog_init());
00425   ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,"QC FWHMX",fwhm_x,
00426                                    "QC FWHM X","%f"));
00427   ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,"QC FWHMY",fwhm_y,
00428                                    "QC FWHM Y","%f"));
00429   ck0_nomsg(sinfo_qclog_add_bool(qclog_tbl,PAF_NAME_LAMP_HALO,
00430                                  hlamp_st,KEY_NAME_LAMP_HALO,"%d"));
00431   ck0_nomsg(sinfo_qclog_add_bool(qclog_tbl,PAF_NAME_SHUT2_ST,shut2_st,
00432                                  KEY_NAME_SHUT2_ST,"%d"));
00433 
00434   ck0(sinfo_pro_save_ima(med_img1,stk,sof,cfg->outName,PRO_PSF,
00435              qclog_tbl,plugin_id,config),
00436       "cannot save ima %s", cfg->outName);
00437 
00438   sinfo_free_table(&qclog_tbl);
00439   sinfo_new_set_wcs_image(med_img1,cfg->outName,cx1, cy1);
00440   sinfo_free_image(&med_img1);
00441   sinfo_free_frameset(&stk);
00442   sinfo_free_psf(&cfg);
00443   return 0;
00444 
00445  cleanup:
00446 
00447   sinfo_free_table(&qclog_tbl);
00448   sinfo_free_imagelist(&cube2);
00449   sinfo_free_imagelist(&cube1);
00450   sinfo_free_table(&enc_energy);
00451   sinfo_free_image(&med_img1);
00452   sinfo_free_table(&ao_performance);
00453   sinfo_free_propertylist(&plist) ;
00454   sinfo_free_psf(&cfg);
00455   sinfo_free_frameset(&stk);
00456   return -1 ;
00457 
00458 }
00459 
00460 static int 
00461 sinfo_get_strehl_from_image(cpl_image* img,
00462                             double ws,
00463                             double we,
00464                             double pscale,
00465                             double strehl_star_radius,
00466                             double strehl_bg_r1,
00467                             double strehl_bg_r2,
00468                             double* strehl,
00469                             double* strehl_err)
00470 {
00471 
00472   cpl_image* img_dup=NULL;
00473 
00474   double dlam=0.;
00475   double lam=0.;
00476   double norm=0.;
00477   double xc=0.;
00478   double yc=0.;
00479   double sx=0.;
00480   double sy=0.;
00481   double fwhm_x=0.;
00482   double fwhm_y=0.;
00483   double max_ima_cx=0.;
00484   double max_ima_cy=0.;
00485   double bkg=0.;
00486   double psf_peak=0.;
00487   double psf_flux=0.;
00488   double bkg_noise=0.;
00489   double star_bkg=0.;
00490   double star_peak=0.;
00491   double star_flux=0.;
00492 
00493   int max_ima_x=0;
00494   int max_ima_y=0;
00495   int wllx=0;
00496   int wlly=0;
00497   int wurx=0;
00498   int wury=0;
00499   int ima_szx=0;
00500   int ima_szy=0;
00501 
00502 
00503   lam = (double)0.5*(ws+we);
00504   dlam=we-ws;
00505   sinfo_msg_warning("ws=%f we=%f dl=%f",ws,we,dlam);
00506   check_nomsg(img_dup=cpl_image_duplicate(img));
00507   sinfo_clean_nan(&img_dup);
00508   check_nomsg(cpl_image_get_maxpos(img_dup,&max_ima_x,&max_ima_y));
00509   sinfo_free_image(&img_dup);
00510 
00511   check_nomsg(ima_szx=cpl_image_get_size_x(img));
00512   check_nomsg(ima_szy=cpl_image_get_size_y(img));
00513   sinfo_check_borders(&max_ima_x,ima_szx,SINFO_STREHL_WINDOW);
00514   sinfo_check_borders(&max_ima_y,ima_szy,SINFO_STREHL_WINDOW);
00515   sinfo_get_safe_box(&wllx,&wlly,&wurx,&wury,max_ima_x,max_ima_y,SINFO_PSF_SZ,
00516                                ima_szx,ima_szy);
00517 
00518   /* 
00519   cpl_image_get_maxpos_window(img,wllx,wlly,wurx,wury,&max_ima_x,&max_ima_y); 
00520    */
00521   check_nomsg(max_ima_cx=cpl_image_get_centroid_x_window(img,wllx,wlly,
00522                                                          wurx,wury));
00523   check_nomsg(max_ima_cy=cpl_image_get_centroid_y_window(img,wllx,wlly,
00524                                                          wurx,wury)); 
00525 
00526   if(CPL_ERROR_NONE != cpl_image_fit_gaussian(img,max_ima_x,max_ima_y,
00527                                               SINFO_PSF_SZ,&norm,&xc,&yc,
00528                                               &sx,&sy,&fwhm_x,&fwhm_y)) {
00529                                            
00530     sinfo_msg_warning("Gaussian fit failed");
00531     cpl_error_reset();
00532 
00533   }
00534 
00535 
00536   check_nomsg(bkg=irplib_strehl_ring_background(img,max_ima_x,max_ima_y,
00537                     19,21,IRPLIB_BG_METHOD_AVER_REJ)) ;
00538 
00539 
00540   if(CPL_ERROR_NONE != irplib_strehl_compute(img,
00541                                              SINFO_STREHL_M1,
00542                                              SINFO_STREHL_M2,
00543                                              lam,
00544                                              dlam,
00545                                              pscale,
00546                                  SINFO_STREHL_BOX_SIZE,
00547                                              max_ima_x,
00548                                              max_ima_y,
00549                                  strehl_star_radius,
00550                                              strehl_bg_r1,
00551                                              strehl_bg_r2,
00552                                  NOISE_HSIZE,
00553                                              NOISE_NSAMPLES,
00554                                  strehl,
00555                                              strehl_err,
00556                                              &star_bkg,
00557                                              &star_peak,
00558                                              &star_flux,
00559                                  &psf_peak,
00560                                              &psf_flux,
00561                                              &bkg_noise)) {
00562 
00563 
00564      sinfo_msg_warning("Problem computing strehl");
00565      *strehl=-1;
00566      *strehl_err=0;
00567 
00568      irplib_error_reset();
00569      cpl_error_reset();
00570   }
00571 
00572   return 0;
00573 
00574  cleanup:
00575 
00576   return -1;
00577 
00578 }
00579 
00580 
00581 
00582 
00583 static int 
00584 sinfo_get_strehl_from_slice(cpl_imagelist* cube,
00585                             double disp, 
00586                             double cWave,
00587                             double ws,
00588                             double we,
00589                             double pscale,
00590                             double strehl_star_radius,
00591                             double strehl_bg_r1,
00592                             double strehl_bg_r2,
00593                             double* strehl,
00594                             double* strehl_err)
00595 {
00596   cpl_image* img_dup=NULL;
00597   cpl_image* img=NULL;
00598 
00599   double dlam=0.;
00600   double lam=0.;
00601   double norm=0.;
00602   double xc=0.;
00603   double yc=0.;
00604   double sx=0.;
00605   double sy=0.;
00606   double fwhm_x=0.;
00607   double fwhm_y=0.;
00608   double max_ima_cx=0.;
00609   double max_ima_cy=0.;
00610   double bkg=0.;
00611   double psf_peak=0.;
00612   double psf_flux=0.;
00613   double bkg_noise=0.;
00614   double star_bkg=0.;
00615   double star_peak=0.;
00616   double star_flux=0.;
00617 
00618   int max_ima_x=0;
00619   int max_ima_y=0;
00620   int wllx=0;
00621   int wlly=0;
00622   int wurx=0;
00623   int wury=0;
00624   int ima_szx=0;
00625   int ima_szy=0;
00626 
00627 
00628   lam = (double)0.5*(ws+we);
00629   dlam=we-ws;
00630 
00631   img=sinfo_new_average_cube_to_image_between_waves(cube,disp,cWave,ws,we);
00632   check_nomsg(img_dup=cpl_image_duplicate(img));
00633   sinfo_clean_nan(&img_dup);
00634   check_nomsg(cpl_image_get_maxpos(img_dup,&max_ima_x,&max_ima_y));
00635   check_nomsg(cpl_image_delete(img_dup));
00636 
00637   check_nomsg(ima_szx=cpl_image_get_size_x(img));
00638   check_nomsg(ima_szy=cpl_image_get_size_y(img));
00639   sinfo_check_borders(&max_ima_x,ima_szx,SINFO_STREHL_WINDOW);
00640   sinfo_check_borders(&max_ima_y,ima_szy,SINFO_STREHL_WINDOW);
00641   sinfo_get_safe_box(&wllx,&wlly,&wurx,&wury,max_ima_x,max_ima_y,SINFO_PSF_SZ,
00642                                ima_szx,ima_szy);
00643 
00644   /* 
00645   cpl_image_get_maxpos_window(img,wllx,wlly,wurx,wury,&max_ima_x,&max_ima_y); 
00646    */
00647   check_nomsg(max_ima_cx=cpl_image_get_centroid_x_window(img,wllx,wlly,
00648                                                          wurx,wury));
00649   check_nomsg(max_ima_cy=cpl_image_get_centroid_y_window(img,wllx,wlly,
00650                                                          wurx,wury)); 
00651 
00652   if(CPL_ERROR_NONE != cpl_image_fit_gaussian(img,max_ima_x,max_ima_y,
00653                                               SINFO_PSF_SZ,&norm,&xc,&yc,
00654                                               &sx,&sy,&fwhm_x,&fwhm_y)) {
00655                                            
00656     sinfo_msg_warning("Gaussian fit failed");
00657     cpl_error_reset();
00658 
00659   }
00660 
00661 
00662   check_nomsg(bkg=irplib_strehl_ring_background(img,max_ima_x,max_ima_y,
00663                     19,21,IRPLIB_BG_METHOD_AVER_REJ)) ;
00664 
00665 
00666   if(CPL_ERROR_NONE != irplib_strehl_compute(img,
00667                                              SINFO_STREHL_M1,
00668                                              SINFO_STREHL_M2,
00669                                              lam,
00670                                              dlam,
00671                                              pscale,
00672                                  SINFO_STREHL_BOX_SIZE,
00673                                              max_ima_x,
00674                                              max_ima_y,
00675                                  strehl_star_radius,
00676                                              strehl_bg_r1,
00677                                              strehl_bg_r2,
00678                                  NOISE_HSIZE,
00679                                              NOISE_NSAMPLES,
00680                                  strehl,
00681                                              strehl_err,
00682                                              &star_bkg,
00683                                              &star_peak,
00684                                              &star_flux,
00685                                  &psf_peak,
00686                                              &psf_flux,
00687                                              &bkg_noise)) {
00688 
00689 
00690      sinfo_msg_warning("Problem computing strehl");
00691      *strehl=-1;
00692      *strehl_err=0;
00693 
00694      irplib_error_reset();
00695      cpl_error_reset();
00696   }
00697 
00698 
00699   sinfo_free_image(&img);
00700 
00701 
00702   return 0;
00703 
00704  cleanup:
00705   return -1;
00706 
00707 }
00708 
00709 
00710 
00711 
00712 
00713 static int sinfo_get_strehl_from_2slices(cpl_imagelist* cube1,
00714                      cpl_imagelist* cube2,
00715                                 double disp, 
00716                                 double cWave,
00717                                 double ws,
00718                                 double we,
00719                                 double pscale1,
00720                                 double pscale2,
00721                                 double strehl_star_rad1,
00722                                 double strehl_bg_rmin1,
00723                                 double strehl_bg_rmax1,
00724                                 double* strehl,
00725                                 double* strehl_err)
00726 {
00727   cpl_image* img1=NULL;
00728   cpl_image* img2=NULL;
00729 
00730   cpl_image* img_dup=NULL;
00731 
00732   double dlam=0.;
00733   double lam=0.;
00734   double norm=0.;
00735   double xc=0.;
00736   double yc=0.;
00737   double sx=0.;
00738   double sy=0.;
00739   double fwhm_x=0.;
00740   double fwhm_y=0.;
00741   double max_ima1_cx=0.;
00742   double max_ima1_cy=0.;
00743 
00744   double bkg1=0.;
00745   double psf_peak=0.;
00746   double psf_flux=0.;
00747   double bkg_noise=0.;
00748   double star_bkg=0.;
00749   double star_peak=0.;
00750   double star_flux=0.;
00751 
00752   int max_ima1_x=0;
00753   int max_ima1_y=0;
00754 
00755   int max_ima2_x=0;
00756   int max_ima2_y=0;
00757 
00758   int wllx=0;
00759   int wlly=0;
00760   int wurx=0;
00761   int wury=0;
00762   int ima_szx=0;
00763   int ima_szy=0;
00764 
00765 
00766   lam = (double)0.5*(ws+we);
00767   dlam=we-ws;
00768 
00769   img1=sinfo_new_average_cube_to_image_between_waves(cube1,disp,cWave,ws,we);
00770   img2=sinfo_new_average_cube_to_image_between_waves(cube2,disp,cWave,ws,we);
00771 
00772   check_nomsg(img_dup=cpl_image_duplicate(img1));
00773   sinfo_clean_nan(&img_dup);
00774   check_nomsg(cpl_image_get_maxpos(img_dup,&max_ima1_x,&max_ima1_y));
00775   sinfo_free_image(&img_dup);
00776 
00777   check_nomsg(img_dup=cpl_image_duplicate(img2));
00778   sinfo_clean_nan(&img_dup);
00779   check_nomsg(cpl_image_get_maxpos(img_dup,&max_ima2_x,&max_ima2_y));
00780   sinfo_free_image(&img_dup);
00781 
00782   check_nomsg(ima_szx=cpl_image_get_size_x(img1));
00783   check_nomsg(ima_szy=cpl_image_get_size_y(img1));
00784   sinfo_check_borders(&max_ima1_x,ima_szx,SINFO_STREHL_WINDOW);
00785   sinfo_check_borders(&max_ima1_y,ima_szy,SINFO_STREHL_WINDOW);
00786   sinfo_get_safe_box(&wllx,&wlly,&wurx,&wury,max_ima1_x,max_ima1_y,
00787                          SINFO_PSF_SZ,ima_szx,ima_szy);
00788 
00789   /* 
00790   cpl_image_get_maxpos_window(img,wllx,wlly,wurx,wury,&max_ima_x,&max_ima_y); 
00791    */
00792   check_nomsg(max_ima1_cx=cpl_image_get_centroid_x_window(img1,wllx,wlly,
00793                                                           wurx,wury));
00794   check_nomsg(max_ima1_cy=cpl_image_get_centroid_y_window(img1,wllx,wlly,
00795                                                           wurx,wury)); 
00796 
00797   if(CPL_ERROR_NONE != cpl_image_fit_gaussian(img1,max_ima1_x,max_ima1_y,
00798                                               SINFO_PSF_SZ,
00799                                   &norm,&xc,&yc,&sx,&sy,
00800                                               &fwhm_x,&fwhm_y)) {
00801 
00802     sinfo_msg_warning("problem in Gaussian fit");
00803     cpl_error_reset();
00804   }
00805 
00806 
00807   check_nomsg(bkg1=irplib_strehl_ring_background(img1,max_ima1_x,max_ima1_y,
00808                                 19,21,
00809                                                 IRPLIB_BG_METHOD_AVER_REJ));
00810 
00811 
00812   if(CPL_ERROR_NONE != sinfo_strehl_compute(img1,img2,
00813                         SINFO_STREHL_M1,SINFO_STREHL_M2,lam,dlam,
00814                         pscale1,pscale2,
00815             SINFO_STREHL_BOX_SIZE,
00816                         max_ima1_x,max_ima1_y,max_ima2_x,max_ima2_y,
00817             strehl_star_rad1,strehl_bg_rmin1,strehl_bg_rmax1,
00818             NOISE_HSIZE,NOISE_NSAMPLES,
00819             strehl,strehl_err,&star_bkg,&star_peak,&star_flux,
00820                         &psf_peak,&psf_flux,&bkg_noise))
00821     {
00822       sinfo_msg_warning("Problem computing strehl, set it to -1");
00823       *strehl=-1;
00824       *strehl_err=0;
00825       irplib_error_reset();
00826       cpl_error_reset();
00827 
00828 
00829     }
00830 
00831   sinfo_free_image(&img1);
00832 
00833   return 0;
00834  cleanup:
00835 
00836   sinfo_free_image(&img_dup);
00837   return -1;
00838 
00839 
00840 }
00841 
00842 
00843 
00844 
00845 cpl_table* sinfo_get_encircled_energy(cpl_frameset* sof,
00846                                       cpl_image* img, 
00847                                       double* fwhm_x, 
00848                                       double* fwhm_y,
00849                                       cpl_table** qclog_tbl)
00850 {
00851 
00852   cpl_image* img_dup=NULL;
00853   int max_ima_x=0;
00854   int max_ima_y=0;
00855   int wllx=0;
00856   int wlly=0;
00857   int wurx=0;
00858   int wury=0;
00859   const double d_mirror = 8.;
00860   const double factor = 180/PI_NUMB*3600.;
00861   double max_ima_cx=0;
00862   double max_ima_cy=0;
00863   double norm=0.;
00864   double xc=0.;
00865   double yc=0.;
00866   double sx=0.;
00867   double sy=0.;
00868   double flux=0;
00869   double flux_max=0;
00870   double pix_scale=0;
00871   double lam=0.;
00872   double pscale=0.;
00873   int dr_difr=0;
00874 
00875   double r=0.;
00876   double bkg=0.;
00877   int i=0;
00878   int ni=0;
00879   int ir_difr=0;
00880   int dr=0;
00881   int rmin=0;
00882 
00883   char band[MAX_NAME_SIZE];
00884   char spat_res[MAX_NAME_SIZE];
00885 
00886   cpl_table* enc_energy=NULL;
00887   cpl_frame* frame=NULL;
00888 
00889   int ima_szx=0;
00890   int ima_szy=0;
00891 
00892 
00893 
00894   if(NULL != cpl_frameset_find(sof,PRO_COADD_PSF)) {
00895     frame = cpl_frameset_find(sof,PRO_COADD_PSF);
00896   } else if(NULL != cpl_frameset_find(sof,PRO_OBS_PSF)) {
00897     frame = cpl_frameset_find(sof,PRO_OBS_PSF);
00898   } else if(NULL != cpl_frameset_find(sof,PRO_COADD_STD)) {
00899     frame = cpl_frameset_find(sof,PRO_COADD_STD);
00900   } else if(NULL != cpl_frameset_find(sof,PRO_OBS_STD)) {
00901     frame = cpl_frameset_find(sof,PRO_OBS_STD);
00902   } else if(NULL != cpl_frameset_find(sof,PRO_COADD_OBJ)) {
00903     frame = cpl_frameset_find(sof,PRO_COADD_OBJ);
00904   } else if(NULL != cpl_frameset_find(sof,PRO_OBS_OBJ)) {
00905     frame = cpl_frameset_find(sof,PRO_OBS_OBJ);
00906   } else {
00907     sinfo_msg_error("Frame %s  or %s or  %s  or %s or %s  or %s not found!", 
00908             PRO_COADD_PSF,PRO_OBS_PSF,
00909             PRO_COADD_STD, PRO_OBS_STD,
00910             PRO_COADD_OBJ, PRO_OBS_OBJ);
00911     return NULL;
00912   }
00913 
00914   sinfo_get_spatial_res(frame,spat_res);
00915   sinfo_get_band(frame,band);
00916   pix_scale=atof(spat_res);
00917   lam=sinfo_get_wave_cent(band);
00918   /* factor 2 due to change of detector to 2K */
00919   pscale=0.5*pix_scale;
00920 
00921 
00922 
00923   dr_difr=factor*1.22*lam*1.e-6/d_mirror/pscale;
00924   ir_difr=floor(dr_difr+0.5);
00925   if (pix_scale==0.025) {
00926     ni=10;
00927     rmin=ir_difr;
00928     dr=rmin;
00929   } else {
00930     ni=15;
00931     sinfo_msg_warning("Reset diffraction limit");
00932     ir_difr=10;
00933     rmin=1;
00934     dr=2;
00935   }
00936 
00937   sinfo_msg("Diffraction limit: %d",ir_difr);
00938 
00939   check_nomsg(img_dup=cpl_image_duplicate(img));
00940   sinfo_clean_nan(&img_dup);
00941   check_nomsg(cpl_image_get_maxpos(img_dup,&max_ima_x,&max_ima_y));
00942   sinfo_free_image(&img_dup);
00943 
00944 
00945 
00946   check_nomsg(ima_szx=cpl_image_get_size_x(img));
00947   check_nomsg(ima_szy=cpl_image_get_size_y(img));
00948   sinfo_check_borders(&max_ima_x,ima_szx,SINFO_STREHL_WINDOW);
00949   sinfo_check_borders(&max_ima_y,ima_szy,SINFO_STREHL_WINDOW);
00950   sinfo_get_safe_box(&wllx,&wlly,&wurx,&wury,max_ima_x,max_ima_y,SINFO_PSF_SZ,
00951                                ima_szx,ima_szy);
00952 
00953   check_nomsg(max_ima_cx=cpl_image_get_centroid_x_window(img,wllx,wlly,
00954                                                          wurx,wury));
00955   check_nomsg(max_ima_cy=cpl_image_get_centroid_y_window(img,wllx,wlly,
00956                                                          wurx,wury));
00957 
00958   if(CPL_ERROR_NONE != cpl_image_fit_gaussian(img,max_ima_x,max_ima_y,
00959                                               SINFO_PSF_SZ,
00960                                               &norm,&xc,&yc,&sx,&sy,
00961                                               fwhm_x,fwhm_y)) {
00962 
00963 
00964     sinfo_msg_warning("problem in Gaussian fit");
00965     cpl_error_reset();
00966   }
00967  check_nomsg(enc_energy = cpl_table_new(ni));
00968  check_nomsg(cpl_table_new_column(enc_energy,"r_pix", CPL_TYPE_INT));
00969  check_nomsg(cpl_table_new_column(enc_energy,"r_mas", CPL_TYPE_DOUBLE));
00970  check_nomsg(cpl_table_new_column(enc_energy,"r_dif", CPL_TYPE_DOUBLE));
00971  check_nomsg(cpl_table_new_column(enc_energy,"abs_energy" , CPL_TYPE_DOUBLE));
00972  check_nomsg(cpl_table_new_column(enc_energy,"rel_energy" , CPL_TYPE_DOUBLE));
00973  /* encircled energy computation */
00974  check_nomsg(bkg=irplib_strehl_ring_background(img,max_ima_x,max_ima_y,
00975                     25,30,IRPLIB_BG_METHOD_AVER_REJ)) ;
00976  r=rmin+(ni-1)*dr;
00977  check_nomsg(flux_max=irplib_strehl_disk_flux(img,max_ima_x,max_ima_y,r,bkg));
00978  r=rmin;
00979 
00980  for(i=0; i<ni; i++)
00981    {
00982      check_nomsg(flux=irplib_strehl_disk_flux(img,max_ima_x,max_ima_y,r,bkg));
00983      check_nomsg(cpl_table_set_int(enc_energy,"r_pix",i,r));
00984      check_nomsg(cpl_table_set_double(enc_energy,"r_mas",i,r*pscale));
00985      check_nomsg(cpl_table_set_double(enc_energy,"r_dif",i,r/ir_difr));
00986      check_nomsg(cpl_table_set_double(enc_energy,"abs_energy",i,flux));
00987      check_nomsg(cpl_table_set_double(enc_energy,"rel_energy",i,flux/flux_max));
00988      r+=dr;
00989 
00990    }
00991  
00992  sinfo_msg("max ima=%d %d\n",max_ima_x,max_ima_y);
00993  sinfo_msg("centroid ima=%f %f\n",max_ima_cx,max_ima_cy);
00994  sinfo_msg("gauss info=%f %f %f %f %f %f %f\n",
00995                           norm,xc,yc,sx,sy,*fwhm_x,*fwhm_y);
00996 
00997  check_nomsg(flux=irplib_strehl_disk_flux(img,max_ima_x,max_ima_y,
00998                                           ir_difr,bkg));
00999  ck0_nomsg(sinfo_qclog_add_double(*qclog_tbl,"QC ENC CORE",
01000                                   flux/flux_max,
01001                                   "Encircled energy within PSF core","%f"));
01002 
01003  return enc_energy;
01004 
01005  cleanup:
01006   sinfo_free_image(&img_dup);
01007 
01008   return NULL;
01009 }
01010 
01011 
01012 static cpl_table* sinfo_get_strehl_from_cube(cpl_imagelist* cube,
01013                                             char* name, 
01014                                             cpl_frame* frame)
01015 {
01016   cpl_table* strehl_tbl=NULL;
01017 
01018   double dispersion=0.;
01019   double centralWave=0.;
01020   double wrange=0;
01021   double wstart=0;
01022   double wstep=0;
01023   double wend=0;
01024   double ws=0;
01025   double we=0;
01026   double pix_scale=0;
01027   double lam=0;
01028   double dlam=0;
01029   double pscale = 0;
01030 
01031   double strehl_star_radius=0;
01032   double strehl_bg_r1=0;
01033   double strehl_bg_r2=0;
01034   double strehl=0;
01035   double strehl_err=0;
01036   char spat_res[MAX_NAME_SIZE];
01037   cpl_propertylist* plist=NULL;
01038 
01039   int naxis3=0;
01040   int nsample=0;
01041   int i=0;
01042 
01043 
01044   sinfo_get_spatial_res(frame,spat_res);
01045   pix_scale=atof(spat_res);
01046   sinfo_msg("pix scale=%f",pix_scale);
01047   /* factor 2 due to change of detector to 2K */
01048   pscale=0.5*pix_scale;
01049 
01050   strehl_star_radius=25*pscale;
01051   strehl_bg_r1=25*pscale;
01052   strehl_bg_r2=30*pscale;
01053 
01054   plist=cpl_propertylist_load(name,0);
01055   dispersion=sinfo_pfits_get_cdelt3(plist);
01056   centralWave=sinfo_pfits_get_crval3(plist);
01057   naxis3=sinfo_pfits_get_naxis3(plist);
01058   sinfo_free_propertylist(&plist);
01059   wrange=dispersion*naxis3;
01060 
01061   wstart = centralWave - (float) (cpl_imagelist_get_size(cube) / 2)*
01062                                  dispersion+dispersion;
01063   wend  =wstart + dispersion * cpl_imagelist_get_size(cube);
01064   wstep=0.025;
01065  /* 
01066    note: 
01067     -wstep as we do not hit the borders where the 
01068     sinfo_gaussian fit has a problem 
01069   */
01070   nsample=(int)((wend-wstart-wstep)/wstep);
01071   check_nomsg(strehl_tbl = cpl_table_new(nsample));
01072   check_nomsg(cpl_table_new_column(strehl_tbl,"wavelength",CPL_TYPE_DOUBLE));
01073   check_nomsg(cpl_table_new_column(strehl_tbl,"strehl",CPL_TYPE_DOUBLE));
01074   check_nomsg(cpl_table_new_column(strehl_tbl,"strehl_error",CPL_TYPE_DOUBLE));
01075 
01076 
01077   for(i=1;i<nsample;i++) {
01078 
01079     ws=wstart+wstep*i;
01080     we=ws+wstep;
01081 
01082     lam = (double)0.5*(ws+we);
01083     dlam=wstep;
01084 
01085     check(sinfo_get_strehl_from_slice(cube,
01086                                 dispersion, 
01087                                 centralWave,
01088                                 ws,
01089                                 we,
01090                                 pscale,
01091                                 strehl_star_radius,
01092                                 strehl_bg_r1,
01093                                 strehl_bg_r2,
01094                                 &strehl,
01095                       &strehl_err),"Error computing strehl");
01096 
01097        if((isnan(lam) ==0) &&
01098           (isnan(lam) ==0) && 
01099           (isnan(lam) ==0)) {
01100      check_nomsg(cpl_table_set_double(strehl_tbl,"wavelength",i,lam));
01101      check_nomsg(cpl_table_set_double(strehl_tbl,"strehl",i,strehl));
01102      check_nomsg(cpl_table_set_double(strehl_tbl,"strehl_error",i,
01103                      strehl_err));
01104     
01105        }
01106   }
01107 
01108   return strehl_tbl;
01109 
01110  cleanup:
01111   return NULL;
01112 
01113 
01114 }
01115 
01116 static double sinfo_get_strehl_from_ima(cpl_image* ima,
01117                                         char* name, 
01118                                         cpl_frame* frame)
01119 {
01120 
01121   double dispersion=0.;
01122   double centralWave=0.;
01123   double wrange=0;
01124   double wstart=0;
01125   double wend=0;
01126   double pix_scale=0;
01127   double pscale = 0;
01128 
01129   double strehl_star_radius=0;
01130   double strehl_bg_r1=0;
01131   double strehl_bg_r2=0;
01132   double strehl=0;
01133   double strehl_err=0;
01134   char spat_res[MAX_NAME_SIZE];
01135 
01136   int naxis3=0;
01137   cpl_propertylist* plist=NULL;
01138 
01139 
01140   sinfo_get_spatial_res(frame,spat_res);
01141   pix_scale=atof(spat_res); 
01142   /* factor 2 due to change of detector to 2K */
01143   pscale=0.5*pix_scale;
01144 
01145 
01146 
01147   strehl_star_radius=25*pscale;
01148   strehl_bg_r1=25*pscale;
01149   strehl_bg_r2=30*pscale;
01150 
01151   plist=cpl_propertylist_load(name,0);
01152   dispersion=sinfo_pfits_get_cdelt3(plist);
01153   centralWave=sinfo_pfits_get_crval3(plist);
01154   naxis3=sinfo_pfits_get_naxis3(plist);
01155   sinfo_free_propertylist(&plist);
01156   wrange=dispersion*naxis3;
01157 
01158   wstart = centralWave - (wrange / 2);
01159   wend   = centralWave + (wrange / 2);
01160 
01161   check(sinfo_get_strehl_from_image(ima,
01162                                 wstart,
01163                                 wend,
01164                                 pscale,
01165                                 strehl_star_radius,
01166                                 strehl_bg_r1,
01167                                 strehl_bg_r2,
01168                                 &strehl,
01169                 &strehl_err),"Computing Strehl");
01170 
01171 
01172 
01173   
01174 
01175  cleanup:
01176   return strehl;
01177 
01178 
01179 }
01180 
01181 static int 
01182 sinfo_get_frm12(cpl_frameset* sof,cpl_frame** frm1,cpl_frame** frm2){
01183 
01184   cpl_frameset* obs=NULL;
01185   int nobs=0;
01186   float eps=0.0001;
01187   float* pix_scale=NULL;
01188   int i=0;
01189   cpl_frame* frame=NULL;
01190 
01191   obs = cpl_frameset_new();
01192   sinfo_contains_frames_kind(sof,obs,PRO_OBS_PSF);
01193   nobs=cpl_frameset_get_size(obs);
01194   if (nobs < 1) {
01195      sinfo_contains_frames_kind(sof,obs,PRO_OBS_STD);
01196      nobs=cpl_frameset_get_size(obs);
01197   }
01198 
01199   nobs=cpl_frameset_get_size(obs);
01200   if (nobs < 1) {
01201     return -1;
01202   } else {
01203     pix_scale=cpl_calloc(nobs,sizeof(float));
01204     for(i=0;i<nobs;i++) {
01205       frame=cpl_frameset_get_frame(obs,i);
01206       pix_scale[i]=sinfo_pfits_get_pixelscale(
01207                            (char*)cpl_frame_get_filename(frame));
01208       if(fabs(pix_scale[i]-0.025)< eps) {
01209         *frm1=cpl_frame_duplicate(frame);
01210       } else if (fabs(pix_scale[i]-0.1) <eps) {
01211         *frm2=cpl_frame_duplicate(frame);        
01212       } else {
01213         sinfo_msg_error("No proper frame found for strehl computation");
01214     return -1;
01215       }
01216     }
01217   }
01218   cpl_free(pix_scale);
01219   cpl_frameset_delete(obs);
01220 
01221   return 0;
01222 
01223 }
01224 
01225 static cpl_table* 
01226 sinfo_get_strehl_from_2cubes(cpl_imagelist* cube1, 
01227                              cpl_imagelist* cube2,
01228                              const char* fname1, 
01229                              cpl_frame* frm1,
01230                              cpl_frame* frm2)
01231 {
01232 
01233   cpl_table* strehl_tbl=NULL;
01234 
01235 
01236   double dispersion=0.;
01237   double centralWave=0.;
01238   double wrange=0;
01239   double wstart=0;
01240   double wstep=0;
01241   double wend=0;
01242   double ws=0;
01243   double we=0;
01244   double pix_scale1=0;
01245   double pix_scale2=0;
01246   double lam=0;
01247   double dlam=0;
01248   double pscale1 = 0;
01249   double pscale2 = 0;
01250 
01251   double strehl_star_rad1=0;
01252   double strehl_star_rad2=0;
01253   double strehl_bg_rmin1=0;
01254   double strehl_bg_rmin2=0;
01255   double strehl_bg_rmax1=0;
01256   double strehl_bg_rmax2=0;
01257   double strehl=0;
01258   double strehl_err=0;
01259   char res1[MAX_NAME_SIZE];
01260   char res2[MAX_NAME_SIZE];
01261 
01262   int naxis3=0;
01263   int nsample=0;
01264   int i=0;
01265   cpl_propertylist* plist=NULL;
01266 
01267   sinfo_get_spatial_res(frm1,res1);
01268   sinfo_get_spatial_res(frm2,res2);
01269   pix_scale1=atof(res1);
01270   pix_scale2=atof(res2);
01271   /* factor 2 due to change of detector to 2K */
01272   pscale1=0.5*pix_scale1;
01273   pscale2=0.5*pix_scale2;
01274 
01275 
01276 
01277   strehl_star_rad1=25*pscale1;
01278   strehl_bg_rmin1=25*pscale1;
01279   strehl_bg_rmax1=30*pscale1;
01280 
01281   strehl_star_rad2=25*pscale2;
01282   strehl_bg_rmin2=25*pscale2;
01283   strehl_bg_rmax2=30*pscale2;
01284   plist=cpl_propertylist_load(fname1,0);
01285   dispersion=sinfo_pfits_get_cdelt3(plist);
01286   centralWave=sinfo_pfits_get_crval3(plist);
01287   naxis3=sinfo_pfits_get_naxis3(plist);
01288   sinfo_free_propertylist(&plist);
01289   wrange=dispersion*naxis3;
01290 
01291   wstart = centralWave - (float) (cpl_imagelist_get_size(cube1) / 2)*
01292                                   dispersion+dispersion;
01293   wend   = wstart + dispersion * cpl_imagelist_get_size(cube1);
01294   wstep  = 0.025;
01295  /* 
01296    note: 
01297     -wstep as we do not hit the borders where the 
01298     sinfo_gaussian fit has a problem 
01299   */
01300   nsample=(int)((wend-wstart-wstep)/wstep);
01301 
01302   check_nomsg(strehl_tbl = cpl_table_new(nsample));
01303   check_nomsg(cpl_table_new_column(strehl_tbl,"wavelength",CPL_TYPE_DOUBLE));
01304   check_nomsg(cpl_table_new_column(strehl_tbl,"strehl",CPL_TYPE_DOUBLE));
01305   check_nomsg(cpl_table_new_column(strehl_tbl,"strehl_error",CPL_TYPE_DOUBLE));
01306 
01307 
01308   for(i=1;i<nsample;i++) {
01309 
01310     ws=wstart+wstep*i;
01311     we=ws+wstep;
01312 
01313     lam = (double)0.5*(ws+we);
01314     dlam=wstep;
01315 
01316     sinfo_get_strehl_from_2slices(cube1,
01317                                   cube2,
01318                                 dispersion, 
01319                                 centralWave,
01320                                 ws,
01321                                 we,
01322                                 pscale1,
01323                                 pscale2,
01324                                 strehl_star_rad1,
01325                                 strehl_bg_rmin1,
01326                                 strehl_bg_rmax1,
01327                                 &strehl,
01328                                 &strehl_err);
01329 
01330 
01331        if((isnan(lam) ==0) &&
01332           (isnan(lam) ==0) && 
01333           (isnan(lam) ==0)) {
01334      check_nomsg(cpl_table_set_double(strehl_tbl,"wavelength",i,lam));
01335      check_nomsg(cpl_table_set_double(strehl_tbl,"strehl",i,strehl));
01336      check_nomsg(cpl_table_set_double(strehl_tbl,"strehl_error",
01337                                           i,strehl_err));
01338     
01339        }
01340 
01341   }
01342 
01343   return strehl_tbl;
01344  cleanup:
01345   return NULL;
01346 }
01347 
01348 
01349 
01350 /*---------------------------------------------------------------------------*/
01383 /*---------------------------------------------------------------------------*/
01384 int sinfo_strehl_compute(
01385         const cpl_image *   im1,
01386         const cpl_image *   im2,
01387         double              m1,
01388         double              m2,
01389         double              lam,
01390         double              dlam,
01391         double              pscale1,
01392         double              pscale2,
01393         int                 size,
01394         int                 xpos1,
01395         int                 ypos1,
01396         int                 xpos2,
01397         int                 ypos2,
01398         double              r1,
01399         double              r2,
01400         double              r3,
01401         int                 noise_box_sz,
01402         int                 noise_nsamples,
01403         double          *   strehl,
01404         double          *   strehl_err,
01405         double          *   star_bg,
01406         double          *   star_peak,
01407         double          *   star_flux,
01408         double          *   psf_peak,
01409         double          *   psf_flux,
01410         double          *   bg_noise)
01411 {
01412     cpl_image * psf ;
01413     const int   window_size = 5;
01414     int         star_radius ;
01415     int         ring[4] ;
01416     double peak2 =0;
01417 
01418     /* Test inputs */
01419     irplib_assure_code(im1 != NULL,         CPL_ERROR_NULL_INPUT);
01420     irplib_assure_code(im2 != NULL,         CPL_ERROR_NULL_INPUT);
01421     irplib_assure_code(strehl != NULL,     CPL_ERROR_NULL_INPUT);
01422     irplib_assure_code(strehl_err != NULL, CPL_ERROR_NULL_INPUT);
01423     irplib_assure_code(star_bg != NULL,    CPL_ERROR_NULL_INPUT);
01424     irplib_assure_code(star_peak != NULL,  CPL_ERROR_NULL_INPUT);
01425     irplib_assure_code(star_flux != NULL,  CPL_ERROR_NULL_INPUT);
01426     irplib_assure_code(psf_peak != NULL,   CPL_ERROR_NULL_INPUT);
01427     irplib_assure_code(psf_flux != NULL,   CPL_ERROR_NULL_INPUT);
01428 
01429     irplib_assure_code(pscale1 > 0.0,      CPL_ERROR_ILLEGAL_INPUT);
01430     irplib_assure_code(pscale2 > 0.0,      CPL_ERROR_ILLEGAL_INPUT);
01431 
01432     irplib_assure_code(xpos1-window_size > 0, CPL_ERROR_ACCESS_OUT_OF_RANGE);
01433     irplib_assure_code(ypos1-window_size > 0, CPL_ERROR_ACCESS_OUT_OF_RANGE);
01434     irplib_assure_code(xpos2-window_size > 0, CPL_ERROR_ACCESS_OUT_OF_RANGE);
01435     irplib_assure_code(ypos2-window_size > 0, CPL_ERROR_ACCESS_OUT_OF_RANGE);
01436 
01437     irplib_assure_code(xpos1+window_size <= cpl_image_get_size_x(im1),
01438                        CPL_ERROR_ACCESS_OUT_OF_RANGE);
01439     irplib_assure_code(ypos1+window_size <= cpl_image_get_size_y(im1),
01440                        CPL_ERROR_ACCESS_OUT_OF_RANGE);
01441 
01442     irplib_assure_code(xpos2+window_size <= cpl_image_get_size_x(im2),
01443                        CPL_ERROR_ACCESS_OUT_OF_RANGE);
01444     irplib_assure_code(ypos2+window_size <= cpl_image_get_size_y(im2),
01445                        CPL_ERROR_ACCESS_OUT_OF_RANGE);
01446 
01447     irplib_assure_code(r1 > 0.0,      CPL_ERROR_ILLEGAL_INPUT);
01448     irplib_assure_code(r2 > 0.0,      CPL_ERROR_ILLEGAL_INPUT);
01449     irplib_assure_code(r3 > r2,       CPL_ERROR_ILLEGAL_INPUT);
01450 
01451     /* Computing a Strehl ratio is a story between an ideal PSF */
01452     /* and a candidate image supposed to approximate this ideal PSF. */
01453 
01454     /* Generate first appropriate PSF to find max peak: same pscale as
01455        the one of the image where we compute the flux */
01456     psf = irplib_strehl_generate_psf(m1, m2, lam, dlam, pscale2, size);
01457     irplib_assure_code(psf != NULL,      CPL_ERROR_ILLEGAL_OUTPUT);
01458 
01459     /* Compute flux in PSF and find max peak */
01460     *psf_peak = cpl_image_get_max(psf) ;
01461     /*
01462     sinfo_msg_warning("Total flux in PSF is: %g\n",cpl_image_get_flux(psf));
01463     */
01464     cpl_image_delete(psf) ;
01465     assert( *psf_peak > 0.0); /* The ideal PSF has a positive maximum */
01466     *psf_flux = 1.0; /* The psf flux, cpl_image_get_flux(psf), is always 1 */
01467 
01468     ring[0] = xpos2 ;
01469     ring[1] = ypos2 ;
01470     ring[2] = (int)(r2/pscale2);
01471     ring[3] = (int)(r3/pscale2);
01472 
01473     /* Measure the background in the candidate image */
01474     *star_bg = irplib_strehl_ring_background(im2, xpos2, ypos2, 
01475                                                   ring[2], ring[3],
01476                                              IRPLIB_BG_METHOD_AVER_REJ);
01477 
01478     /* Compute star_radius in pixels */
01479     star_radius = (int)(r1/pscale2);
01480 
01481     /* Measure the flux on the candidate image */
01482     *star_flux = irplib_strehl_disk_flux(im2, xpos2, ypos2, star_radius, 
01483                                          *star_bg);
01484 
01485     irplib_assure_code(*star_flux > 0.0,      CPL_ERROR_ILLEGAL_OUTPUT);
01486 
01487     /* Measure the peak value on the candidate image */
01488     /* FIXME: Arbitrary choice of image border */
01489     /* here we should not subtract the background */
01490     *star_peak = cpl_image_get_max_window(im1, 
01491                                           xpos1-window_size,
01492                                           ypos1-window_size, 
01493                                           xpos1+window_size, 
01494                                           ypos1+window_size) - *star_bg ;
01495     peak2 = cpl_image_get_max_window(im2, 
01496                                      xpos2-window_size,
01497                      ypos2-window_size, 
01498                                      xpos2+window_size, 
01499                                      ypos2+window_size) - *star_bg ;
01500     if ( fabs(peak2-*star_peak)/peak2 > 0.25) {
01501       sinfo_msg_warning("rel diff: %g",fabs(peak2-*star_peak)/peak2);
01502     }
01503     irplib_assure_code(*star_peak > 0.0,      CPL_ERROR_ILLEGAL_OUTPUT);
01504 
01505     /* Compute Strehl */
01506     *strehl = (*star_peak / *star_flux) / (*psf_peak / *psf_flux) ;
01507 
01508     if (*strehl > 1)
01509         cpl_msg_warning(cpl_func, "Extreme Strehl-ratio=%g, star_peak=%g, "
01510                         "star_flux=%g, psf_peak=%g, psf_flux=%g", *strehl,
01511                         *star_peak, *star_flux, *psf_peak, *psf_flux);
01512 
01513     /* Compute Strehl error */
01514     if (cpl_flux_get_noise_ring(im2, ring, noise_box_sz, noise_nsamples, 
01515                                 bg_noise, NULL) == CPL_ERROR_NONE) {
01516         *strehl_err = SINFO_STREHL_ERROR_COEFFICIENT * (*bg_noise) * pscale2 * 
01517             star_radius * star_radius / *star_flux;
01518         irplib_assure_code(*strehl_err >= 0.0,       CPL_ERROR_ILLEGAL_OUTPUT);
01519     }
01520 
01521     return CPL_ERROR_NONE;
01522 }
01523 
01524 static void 
01525 sinfo_check_borders(int* val,const int max,const int thresh)
01526 {
01527    
01528   *val = ((*val-thresh) > 0) ? *val : thresh;
01529   *val = ((*val+thresh) < max) ? *val : max-thresh-1;
01530   return;
01531 }
01532 
01533 static void 
01534 sinfo_get_safe_box(int* llx, 
01535                    int* lly, 
01536                    int* urx, 
01537                    int* ury, 
01538                    const int xpos,
01539                    const int ypos,
01540                    const int box,
01541                    const int szx, 
01542                    const int szy)
01543 
01544 {
01545   *llx= ((xpos-box)>0)   ? (xpos-box) : 1;
01546   *lly= ((ypos-box)>0)   ? (ypos-box) : 1;
01547   *urx= ((xpos+box)<szx) ? (xpos+box) : szx-1 ;
01548   *ury= ((ypos+box)<szy) ? (ypos+box) : szy-1 ;
01549   return;
01550 }

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