sinfo_new_stdstar.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_stdstar.c
00022    Author       :   J. Schreiber
00023    Created on   :   December 3, 2003
00024    Description  :   this routine doess the optimal extraction of a spectrum 
00025                         of an already reduced data cube of a standard star 
00026                         observation. Additionally, a conversion factor from 
00027                         mag to counts/sec can be determined if the magnitude 
00028                         of the standard star is known. 
00029                         This is done for a number of jittered data cubes and 
00030                         the results are averaged by rejecting the extreme 
00031                         values
00032 
00033  ---------------------------------------------------------------------------*/
00034 
00035 #ifdef HAVE_CONFIG_H
00036 #  include <config.h>
00037 #endif
00038 /*----------------------------------------------------------------------------
00039                                 Includes
00040  ---------------------------------------------------------------------------*/
00041 #include <irplib_std.h>
00042 #include <math.h>
00043 #include "sinfo_new_stdstar.h"
00044 #include "sinfo_standstar_ini_by_cpl.h"
00045 #include "sinfo_pro_save.h"
00046 #include "sinfo_pfits.h"
00047 #include "sinfo_utilities_scired.h"
00048 #include "sinfo_spectrum_ops.h"
00049 #include "sinfo_hidden.h"
00050 #include "sinfo_functions.h"
00051 #include "sinfo_error.h"
00052 #include "sinfo_utils_wrappers.h"
00053 /*----------------------------------------------------------------------------
00054                                 Defines
00055  ---------------------------------------------------------------------------*/
00056 
00057 
00058 
00059 /*----------------------------------------------------------------------------
00060                              Function Definitions
00061  ---------------------------------------------------------------------------*/
00062 static int  sinfo_compute_efficiency(const char* name,
00063                                      standstar_config ** cfg,
00064                      cpl_imagelist  * cube, 
00065                                      cpl_table** tbl_spectrum);
00066 
00074 /*----------------------------------------------------------------------------
00075    Function     :       sinfo_stdstar()
00076    In           :       ini_file: file name of according .ini file
00077    Out          :       integer (0 if it worked, -1 if it doesn't) 
00078    Job          :     this routine carries through the data cube creation of an 
00079                         object science observation using object-sky nodding
00080                         and jittering. This script expects jittered frames that
00081                 were already sky-subtracted
00082                         averaged, flatfielded, spectral tilt corrected and 
00083             interleaved if necessary
00084  ---------------------------------------------------------------------------*/
00085 int 
00086 sinfo_new_stdstar(const char* plugin_id,
00087                   cpl_parameterlist* config, 
00088                   cpl_frameset* sof)
00089 {
00090 
00091     standstar_config * cfg=NULL ;
00092     cpl_imagelist  * cube=NULL ;
00093 
00094     cpl_imagelist  * list_object=NULL ;
00095     cpl_image ** spectrum=NULL ;
00096     cpl_image * img_spct=NULL ;
00097 
00098 
00099   
00100     cpl_frameset* raw=NULL;
00101     cpl_frame* frame=NULL;
00102     cpl_image* std_med_ima=NULL;
00103     cpl_image* std_med_dup=NULL;
00104 
00105     cpl_table* qclog_tbl=NULL;
00106     cpl_table* tbl_spectrum=NULL;
00107     cpl_propertylist* plist=NULL;
00108 
00109     char * name=NULL ;
00110     int fra=0;
00111     float exptime=0 ;
00112     double convfactor=0;
00113     double cleanfactor=0;
00114     float* factor=NULL;
00115 
00116     char std_med_filename[MAX_NAME_SIZE];
00117     char std_cub_filename[MAX_NAME_SIZE];
00118 
00119     double max_ima_cx=0;
00120     double max_ima_cy=0;
00121     int max_ima_x=0;
00122     int max_ima_y=0;
00123     double norm=0;
00124     double xcen=0;
00125     double ycen=0;
00126     double sig_x=0;
00127     double sig_y=0;
00128     double fwhm_x=0;
00129     double fwhm_y=0;
00130     double disp=0;
00131     double dispersion=0;
00132     int i=0;
00133     int wllx=0;
00134     int wlly=0;
00135     int wurx=0;
00136     int wury=0;
00137     int psf_sz=20;
00138     int qc_info=0;
00139     int ima_szx=0;
00140     int ima_szy=0;
00141     int check1=0;
00142     int check2=0;
00143     int check3=0;
00144     int check4=0;
00145     double xshift=0;
00146     double yshift=0;
00147 
00148     float cenpix = 0;
00149     float cenLambda = 0;
00150 
00151     int no=0;
00152     double lo_cut=0.;
00153     double hi_cut=0.;
00154   
00155 
00156 
00157 
00158     double fpar[7];
00159     double dpar[7];
00160     int mpar[7];
00161 
00162     /*
00163        parse the file names and parameters to the cube_config 
00164        data structure cfg
00165      */
00166  
00167     /* sinfo_msg("Parse cpl input"); */
00168      check_nomsg(raw=cpl_frameset_new());
00169  
00170     cknull(cfg=sinfo_parse_cpl_input_standstar(config,sof,&raw),
00171        "could not parse cpl input!") ;
00172 
00173     cknull(list_object = cpl_imagelist_new (),
00174        "could not allocate memory");
00175 
00176 
00177     if (cfg->convInd == 1) {
00178       factor = sinfo_new_floatarray(cfg->nframes);
00179     }
00180 
00181   
00182     if(NULL != cpl_frameset_find(sof,PRO_COADD_STD)) {
00183       frame = cpl_frameset_find(sof,PRO_COADD_STD);
00184       strcpy(std_cub_filename,cpl_frame_get_filename(frame));
00185     } else if(NULL != cpl_frameset_find(sof,PRO_COADD_PSF)) {
00186       frame = cpl_frameset_find(sof,PRO_COADD_PSF);
00187       strcpy(std_cub_filename,cpl_frame_get_filename(frame));
00188     } else if(NULL != cpl_frameset_find(sof,PRO_COADD_OBJ)) {
00189       frame = cpl_frameset_find(sof,PRO_COADD_OBJ);
00190       strcpy(std_cub_filename,cpl_frame_get_filename(frame));
00191     } else if(NULL != cpl_frameset_find(sof,PRO_COADD_PUPIL)) {
00192       frame = cpl_frameset_find(sof,PRO_COADD_PUPIL);
00193       strcpy(std_cub_filename,cpl_frame_get_filename(frame));
00194     } else {
00195       sinfo_msg_error("Frame %s, %s, %s or %s not found! Exit!", 
00196                        PRO_COADD_STD,PRO_COADD_PSF,
00197                        PRO_COADD_OBJ,PRO_COADD_PUPIL );
00198       goto cleanup;
00199     }
00200 
00201     cknull(plist = cpl_propertylist_load(std_cub_filename, 0),
00202       "getting header from reference ima frame %s",std_cub_filename);
00203 
00204     cenpix = sinfo_pfits_get_crpix3(plist);
00205     cenLambda = sinfo_pfits_get_crval3(plist);
00206     dispersion = sinfo_pfits_get_cdelt3(plist);
00207  
00208 
00209     if (cpl_propertylist_contains(plist, KEY_NAME_CDELT3)) {
00210       disp=cpl_propertylist_get_double(plist, KEY_NAME_CDELT3);
00211     } else {
00212       sinfo_msg_warning("Keyword %s not found.",KEY_NAME_CDELT3);
00213     }
00214     sinfo_free_propertylist(&plist) ;
00215     
00216     /* we find automatiocally extraction parameters */
00217 
00218     if(NULL != cpl_frameset_find(sof,PRO_MED_COADD_STD)) {
00219       frame = cpl_frameset_find(sof,PRO_MED_COADD_STD);
00220       strcpy(std_med_filename,cpl_frame_get_filename(frame));
00221       check_nomsg(std_med_ima=cpl_image_load(std_med_filename,
00222                                              CPL_TYPE_FLOAT,0,0));    
00223     } else if(NULL != cpl_frameset_find(sof,PRO_OBS_STD)) {
00224       check_nomsg(frame = cpl_frameset_find(sof,PRO_OBS_STD));
00225       strcpy(std_cub_filename,cpl_frame_get_filename(frame));
00226       check_nomsg(cube = cpl_imagelist_load(std_cub_filename,
00227                                              CPL_TYPE_FLOAT,0));
00228       strcpy(std_med_filename,STDSTAR_OUT_MED_CUBE);
00229       check_nomsg(std_med_ima=cpl_imagelist_collapse_median_create(cube));
00230       sinfo_free_imagelist(&cube);
00231 
00232       ck0(sinfo_pro_save_ima(std_med_ima,sof,sof,STDSTAR_OUT_MED_CUBE,
00233                  PRO_MED_OBS_PSF,NULL,plugin_id,config),
00234       "Error saving image %s tag %s",STDSTAR_OUT_MED_CUBE,PRO_MED_OBS_PSF);
00235 
00236     } else if(NULL != cpl_frameset_find(sof,PRO_MED_COADD_PSF)) {
00237       check_nomsg(frame = cpl_frameset_find(sof,PRO_MED_COADD_PSF));
00238       strcpy(std_med_filename,cpl_frame_get_filename(frame));
00239       check_nomsg(std_med_ima=cpl_image_load(std_med_filename,
00240                                              CPL_TYPE_FLOAT,0,0));    
00241     } else if(NULL != cpl_frameset_find(sof,PRO_OBS_PSF)) {
00242       check_nomsg(frame = cpl_frameset_find(sof,PRO_OBS_PSF));
00243       strcpy(std_cub_filename,cpl_frame_get_filename(frame));
00244       check_nomsg(cube=cpl_imagelist_load(std_cub_filename,CPL_TYPE_FLOAT,0));
00245       strcpy(std_med_filename,STDSTAR_OUT_MED_CUBE);
00246       check_nomsg(std_med_ima=cpl_imagelist_collapse_median_create(cube));
00247       sinfo_free_imagelist(&cube);
00248 
00249       ck0(sinfo_pro_save_ima(std_med_ima,sof,sof,STDSTAR_OUT_MED_CUBE,
00250                  PRO_MED_OBS_PSF,NULL,plugin_id,config),
00251       "Error saving image %s tag %s",STDSTAR_OUT_MED_CUBE,PRO_MED_OBS_PSF);
00252 
00253     } else if(NULL != cpl_frameset_find(sof,PRO_MED_COADD_OBJ)) {
00254       check_nomsg(frame = cpl_frameset_find(sof,PRO_MED_COADD_OBJ));
00255       strcpy(std_med_filename,cpl_frame_get_filename(frame));
00256       check_nomsg(std_med_ima=cpl_image_load(std_med_filename,
00257                                              CPL_TYPE_FLOAT,0,0));    
00258     } else if(NULL != cpl_frameset_find(sof,PRO_OBS_OBJ)) {
00259       check_nomsg(frame = cpl_frameset_find(sof,PRO_OBS_OBJ));
00260       strcpy(std_cub_filename,cpl_frame_get_filename(frame));
00261       check_nomsg(cube = cpl_imagelist_load(std_cub_filename,
00262                                              CPL_TYPE_FLOAT,0));
00263       strcpy(std_med_filename,STDSTAR_OUT_MED_CUBE);
00264       check_nomsg(std_med_ima=cpl_imagelist_collapse_median_create(cube));
00265       sinfo_free_imagelist(&cube);
00266 
00267       ck0(sinfo_pro_save_ima(std_med_ima,sof,sof,STDSTAR_OUT_MED_CUBE,
00268                  PRO_MED_OBS_OBJ,NULL,plugin_id,config),
00269       "Error saving image %s tag %s",STDSTAR_OUT_MED_CUBE,PRO_MED_OBS_OBJ);
00270     } else {
00271       sinfo_msg_error("Frame %s %s %s %s %s %s not found! Exit!",
00272               PRO_MED_COADD_STD, PRO_OBS_STD, 
00273                       PRO_MED_COADD_PSF, PRO_OBS_PSF,
00274                       PRO_MED_COADD_OBJ, PRO_OBS_OBJ);
00275       goto cleanup;
00276     }
00277 
00278     check_nomsg(std_med_dup=cpl_image_duplicate(std_med_ima));
00279     sinfo_clean_nan(&std_med_dup);
00280     check_nomsg(cpl_image_get_maxpos(std_med_dup,&max_ima_x,&max_ima_y));
00281     sinfo_free_image(&std_med_dup);
00282 
00283     ima_szx=cpl_image_get_size_x(std_med_ima);
00284     ima_szy=cpl_image_get_size_y(std_med_ima);
00285     wllx= ((max_ima_x-psf_sz)>0)       ? (max_ima_x-psf_sz) : 1;
00286     wlly= ((max_ima_y-psf_sz)>0)       ? (max_ima_y-psf_sz) : 1;
00287     wurx= ((max_ima_x+psf_sz)<ima_szx) ? (max_ima_x+psf_sz) : ima_szx ;
00288     wury= ((max_ima_y+psf_sz)<ima_szy) ? (max_ima_y+psf_sz) : ima_szy ;
00289     /*
00290     sinfo_msg("wllx=%d wlly=%d wurx=%d wury=%d\n",wllx,wlly,wurx,wury);
00291     cpl_image_get_maxpos_window(std_med_ima,wllx,wlly,wurx,wury,
00292                                 &max_ima_x,&max_ima_y); 
00293     */
00294     check_nomsg(qclog_tbl = sinfo_qclog_init());
00295     check_nomsg(max_ima_cx=cpl_image_get_centroid_x_window(std_med_ima,wllx,
00296                                                            wlly,wurx,wury));
00297     check_nomsg(max_ima_cy=cpl_image_get_centroid_y_window(std_med_ima,wllx,
00298                                                            wlly,wurx,wury));
00299 
00300 
00301       xshift=max_ima_cx-ima_szx/2;
00302       yshift=max_ima_cy-ima_szy/2;
00303 
00304       ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,"QC SHIFTX",xshift,
00305                                        "X shift centroid - center image","%f"));
00306 
00307       ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,"QC SHIFTY",yshift,
00308                                        "Y shift centroid - center image","%f"));
00309 
00310 
00311     if(
00312           ((max_ima_x-psf_sz) < 1) ||
00313           ((max_ima_y-psf_sz) < 1) ||
00314           ((max_ima_x+psf_sz) > ima_szx) ||
00315           ((max_ima_x+psf_sz) > ima_szy) 
00316       ) 
00317       {
00318         psf_sz = (psf_sz < (max_ima_x-1))     ? psf_sz : (max_ima_x-1);
00319         psf_sz = (psf_sz < (max_ima_y-1))     ? psf_sz : (max_ima_y-1);
00320         psf_sz = (psf_sz < ima_szx-max_ima_x) ? psf_sz : (ima_szx-max_ima_x);
00321         psf_sz = (psf_sz < ima_szy-max_ima_y) ? psf_sz : (ima_szy-max_ima_y);
00322       }
00323 
00324 
00325       ck0_nomsg(sinfo_qclog_add_int(qclog_tbl,"QC FWHM LLX",cfg->llx,
00326                                     "STD star FWHM LLX","%d"));
00327       ck0_nomsg(sinfo_qclog_add_int(qclog_tbl,"QC FWHM LLY",cfg->lly,
00328                                     "STD star FWHM LLY","%d"));
00329       ck0_nomsg(sinfo_qclog_add_int(qclog_tbl,"QC FWHM HBX",cfg->halfbox_x,
00330                                     "STD star FWHM HBX","%d"));
00331       ck0_nomsg(sinfo_qclog_add_int(qclog_tbl,"QC FWHM HBX",cfg->halfbox_y,
00332                                     "STD star FWHM HBY","%d"));
00333 
00334 
00335     /* call the 2D-Gaussian fit routine */
00336     for ( i = 0 ; i < 7 ; i++ )
00337     {
00338         mpar[i] = 1 ;
00339     }
00340 
00341      if(-1 == sinfo_new_fit_2d_gaussian(std_med_ima, 
00342                                   fpar, 
00343                                   dpar, 
00344                                   mpar, 
00345                                   cfg->llx, 
00346                                   cfg->lly, 
00347                                   cfg->halfbox_x, 
00348                                   cfg->halfbox_y, 
00349                                   &check4 ) ) {
00350 
00351       sinfo_msg_warning("2d Gaussian fit failed");
00352       cpl_error_reset();
00353       irplib_error_reset();
00354       /* return 0; */
00355 
00356      } else {
00357 
00358       ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,"QC FWHM MAJ",fpar[4],
00359                                        "STD star FWHM on major axis","%f"));
00360       ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,"QC FWHM MIN",fpar[5],
00361                                        "STD star FWHM on minor axis","%f"));
00362       ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,"QC THETA",fpar[6],
00363                                        "STD star ellipsis angle theta","%f"));
00364 
00365 
00366      }
00367 
00368     /*
00369     sinfo_msg("Gauss fit params: xc,yc,amp,bkg,fwhm_x,fwhm_y,angle\n");
00370     for ( i = 0 ; i < 7 ; i++ )
00371       {
00372         sinfo_msg("fpar[%d]=%f dpar[%d]=%f\n",i,fpar[i],i,dpar[i]);
00373       }
00374     */
00375 
00376     if(CPL_ERROR_NONE == cpl_image_fit_gaussian(std_med_ima,
00377                                                 max_ima_x,
00378                                                 max_ima_y,
00379                                                 psf_sz,
00380                                     &norm,
00381                                                 &xcen,
00382                                                 &ycen,
00383                                                 &sig_x,
00384                                                 &sig_y,
00385                                                 &fwhm_x,
00386                                                 &fwhm_y)) {
00387 
00388       ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,"QC FWHMX",fwhm_x,
00389                                        "STD star FWHM on X","%f"));
00390       ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,"QC FWHMY",fwhm_y,
00391                                        "STD star FWHM on Y","%f"));
00392       cfg -> halfbox_x =  (floor)(0.5*(fwhm_x+fwhm_y)*cfg->fwhm_factor+0.5);
00393       cfg -> halfbox_y =  (floor)(0.5*(fwhm_x+fwhm_y)*cfg->fwhm_factor+0.5);
00394 
00395     } else {
00396       sinfo_msg_warning("Problem fitting Gaussian");
00397       irplib_error_reset();
00398 
00399     }
00400     sinfo_free_image(&std_med_ima);
00401 
00402       /*
00403       sinfo_msg("max ima=%d %d psf_sz=%d\n",max_ima_x,max_ima_y,psf_sz);
00404       sinfo_msg("centroid ima=%f %f\n",max_ima_cx,max_ima_cy);
00405       sinfo_msg("gauss=norm=%f xcen=%f ycen=%f sig_x=%f "
00406                 "sig_y=%f fwhm_x=%f fwhm_y=%f\n",
00407                           norm,xcen,ycen,sig_x,sig_y,fwhm_x,fwhm_y);
00408       */
00409 
00410       cfg -> llx = (int)(xcen-cfg->halfbox_x);
00411       cfg -> llx = (cfg -> llx  > 0 ) ? cfg -> llx  : 1;
00412 
00413       if((cfg->llx+2*cfg->halfbox_x) >= ima_szx) {
00414     cfg -> halfbox_x=(int) ((ima_szx-cfg->llx-1)/2);
00415         check1++;
00416       } 
00417 
00418       cfg -> lly = (int)(ycen-cfg->halfbox_y);
00419       cfg -> lly = (cfg -> lly  > 0 ) ? cfg -> lly  : 1;
00420       if((cfg->lly+2*cfg->halfbox_y) >= ima_szy) {
00421     cfg -> halfbox_y=(int) ((ima_szy-cfg->lly-1)/2);
00422         check1++;
00423       } 
00424      ck0_nomsg(sinfo_qclog_add_int(qclog_tbl,"QC CHECK1",check1,
00425                                    "Check on evaluation box","%d"));
00426      /*     
00427       sinfo_msg("llx= %d lly= %d\n",cfg->llx, cfg->lly);
00428       sinfo_msg("halfbox_x=%d halfbox_y=%d\n",cfg->halfbox_x,cfg->halfbox_y);
00429      */
00430  
00431     /*
00432 #----------------------------------------------------------------------
00433 #---------------------------EXTRACTION---------------------------------
00434 #----------------------------------------------------------------------
00435     */
00436 
00437       sinfo_msg("Extraction");
00438       cknull(spectrum = (cpl_image**) cpl_calloc (cfg -> nframes,
00439                                                   sizeof(cpl_image*)),
00440                             "Could not allocate memory for spectrum image");
00441 
00442       for (fra=0; fra < cfg->nframes; fra++) {
00443          name = cfg->inFrameList[fra];
00444          if(sinfo_is_fits_file(name) != 1) {
00445             sinfo_msg_error("Input file %s is not FITS",name);
00446             goto cleanup;
00447      }
00448          cknull(cube = cpl_imagelist_load(name,CPL_TYPE_FLOAT,0),
00449         "could not load data cube" );
00450        
00451          if (exptime == FLAG) {
00452              sinfo_msg_error("could not find exposure time in the fits header" );
00453              return -1;
00454      }
00455          exptime = sinfo_pfits_get_ditndit(name);
00456 
00457          sinfo_msg("cfg->gain %f",cfg->gain);
00458          check_nomsg(tbl_spectrum=cpl_table_new(cpl_imagelist_get_size(cube)));
00459          if(NULL==(spectrum[fra]=sinfo_new_optimal_extraction_from_cube( cube, 
00460                                       cfg->llx,
00461                                       cfg->lly,
00462                                       cfg->halfbox_x,
00463                                       cfg->halfbox_y,
00464                                       cfg->fwhm_factor,
00465                                       BKG_VARIANCE,
00466                                       SKY_FLUX,
00467                                       cfg->gain,
00468                                       exptime,
00469                                       name,
00470                                       &tbl_spectrum,
00471                                       qc_info,
00472                                       &check2))){
00473        sinfo_msg_warning("could not do sinfo_optimalExtractionFromCube" );
00474        cpl_error_reset();
00475      } else {
00476      check_nomsg(cpl_imagelist_set(list_object,
00477                      cpl_image_duplicate(spectrum[fra]), fra)); 
00478 
00479      }
00480             
00481 
00482 
00483 
00484      ck0_nomsg(sinfo_qclog_add_int(qclog_tbl,"QC CHECK2",
00485                                        check2,"Check on evaluation box","%d"));
00486 
00487      if(-1 == sinfo_compute_efficiency(name,&cfg,cube,&tbl_spectrum)) {
00488            sinfo_msg_warning("Problems computing efficiency");
00489            cpl_error_reset();
00490          } 
00491 
00492         ck0(sinfo_pro_save_tbl(tbl_spectrum,raw,sof,
00493                                       (char*)STDSTAR_OUT_TABLE,
00494                       PRO_STD_STAR_SPECTRA,qclog_tbl,
00495                 plugin_id,config),
00496          "cannot dump ima %s", "out_std_star_spectrum.fits");
00497 
00498      
00499 
00500      
00501      sinfo_free_table(&qclog_tbl);
00502 
00503 
00504      /*
00505          if(spectrum[fra] != NULL ) {
00506        sinfo_free_image(&(spectrum)[fra]);
00507          }
00508      */
00509      /*----determine the intensity conversion factor if wished--------*/
00510      if (cfg->convInd == 1) {
00511        sinfo_msg("Determines convertion factor");
00512        
00513            convfactor = sinfo_new_determine_conversion_factor( cube, 
00514                                                 cfg->mag, 
00515                                                 exptime, 
00516                                                 cfg->llx, 
00517                                                 cfg->lly, 
00518                                                 cfg->halfbox_x, 
00519                                                 cfg->halfbox_y, 
00520                                                 &check3 );
00521         
00522  
00523             if (convfactor < -100000.) {
00524            sinfo_msg_warning("could not do sinfo_determineConversionFactor!" );
00525            /* goto cleanup; */
00526         } else {
00527              sinfo_new_array_set_value(factor, convfactor, fra);
00528         }
00529     }
00530         sinfo_free_imagelist(&cube);
00531       } /* end loop over fra */
00532 
00533       sinfo_free_table(&tbl_spectrum);
00534       sinfo_free_image_array(&spectrum,cfg->nframes);
00535       if (cfg->convInd == 1) {
00536     sinfo_msg("Determines clean factor");
00537         cleanfactor = sinfo_new_clean_mean(factor, 
00538                                  cfg->nframes, 
00539                                  cfg->lo_reject*100., 
00540                  cfg->hi_reject*100.);
00541       }
00542       if (cleanfactor > 100000. || cleanfactor == FLAG) {
00543     sinfo_msg_error("could not do sinfo_clean_mean!" );
00544         goto cleanup;
00545       }
00546 
00547 
00548   /*---read the fits header to change the gain and noise parameter-----*/
00549       sinfo_msg("average with rejection");
00550 
00551       no=cpl_imagelist_get_size(list_object);
00552       lo_cut=(floor)(cfg->lo_reject*no+0.5);
00553       hi_cut=(floor)(cfg->hi_reject*no+0.5);
00554       if(no > 0) {
00555          cknull(img_spct=cpl_imagelist_collapse_minmax_create(list_object,
00556                                                               lo_cut,hi_cut),
00557                           "sinfo_average_with_rejection failed" );
00558       }
00559 
00560       sinfo_free_imagelist(&list_object);
00561       if(no > 0) {
00562     check_nomsg(qclog_tbl = sinfo_qclog_init());
00563  
00564         ck0_nomsg(sinfo_qclog_add_double(qclog_tbl,"QC CONVFCT",cleanfactor,
00565                                          "Conversion factor","%g"));
00566 
00567         ck0_nomsg(sinfo_qclog_add_int(qclog_tbl,"QC CHECK3",check3,
00568                                       "Check evaluation box","%d"));
00569   
00570 
00571         ck0(sinfo_pro_save_ima(img_spct,raw,sof,cfg->outName,
00572                     PRO_STD_STAR_SPECTRUM,qclog_tbl,
00573                  plugin_id,config),
00574       "cannot dump ima %s", cfg->outName);
00575       
00576         sinfo_new_set_wcs_spectrum(img_spct,cfg->outName,cenLambda,disp,cenpix);
00577         sinfo_free_table(&qclog_tbl);
00578       }
00579       /*#---free memory---*/
00580       if(factor != NULL) sinfo_new_destroy_array(&factor);
00581       sinfo_free_image(&img_spct);
00582       sinfo_stdstar_free(&cfg);
00583       sinfo_free_frameset(&raw);
00584 
00585       return 0;
00586 
00587  cleanup:
00588       sinfo_free_table(&tbl_spectrum);
00589       sinfo_free_table(&qclog_tbl);
00590       sinfo_free_imagelist(&list_object);
00591       if(spectrum != NULL) sinfo_free_image_array(&spectrum,cfg->nframes);
00592       sinfo_free_image(&std_med_ima);
00593       sinfo_free_image(&std_med_dup);
00594       sinfo_free_imagelist(&cube);
00595       sinfo_free_propertylist(&plist) ;
00596       if(factor != NULL) sinfo_new_destroy_array(&factor);
00597       sinfo_free_image(&img_spct);
00598       sinfo_stdstar_free (&cfg);
00599       sinfo_free_frameset(&raw);
00600     return -1;
00601  
00602 }
00603 
00604 static int  sinfo_compute_efficiency(const char* name,
00605                                      standstar_config ** cfg,
00606                      cpl_imagelist  * cube, 
00607                                      cpl_table** tbl_spectrum)
00608 {
00609 
00610 
00611   int dit=0;
00612   int i=0;
00613   int status=0;
00614 
00615 
00616   double dispersion=0;
00617   double h=0;
00618   double c=0;
00619   double k=0;
00620   double temp=0;
00621   double surface=0;
00622   double gain=0;
00623   double num=0;
00624   double den=0;
00625   double fctr=0;
00626   double wc=0;
00627   double wave=0;
00628   double std_ra=0;
00629   double std_dec=0;
00630   double mag=0;
00631 
00632 
00633   double int_spct=0;
00634   double bb_flux_norm=0;
00635   double efficiency=0;
00636   double f0=0;
00637   char band[FILE_NAME_SZ];
00638   irplib_band std_band=0;
00639   cpl_propertylist* plist=NULL;
00640 
00641   plist=cpl_propertylist_load(name,0);
00642   strcpy(band,sinfo_pfits_get_band(plist));
00643   std_ra =(double) sinfo_pfits_get_ra(plist);
00644   std_dec=(double) sinfo_pfits_get_dec(plist);
00645   sinfo_free_propertylist(&plist);
00646 
00647   if (strcmp(band,"H+K") == 0) {
00648     f0=8.00e-11;   /* erg/s/cm/cm/A */
00649     wc=1.950;
00650         std_band=BAND_K;
00651   } else if (strcmp(band,"K") == 0) {
00652     f0=4.10e-11; 
00653     wc=2.175;
00654         std_band=BAND_K;
00655   } else if (strcmp(band,"J") == 0) {
00656     f0=3.11e-10;   
00657     wc=1.225;
00658         std_band=BAND_J;
00659   } else if (strcmp(band,"H") == 0) {
00660     f0=1.15e-10;   
00661     wc=1.675;
00662         std_band=BAND_H;
00663   }
00664 
00665 
00666   irplib_std_get_mag(std_ra, std_dec, std_band, &mag) ;
00667   (*cfg)->mag=mag; 
00668   sinfo_msg("Catalogue magnitude=%g",mag);
00669 
00670   plist=cpl_propertylist_load(name,0);
00671   dispersion = sinfo_pfits_get_cdelt3(plist);
00672   dit = sinfo_pfits_get_dit(plist); // s 
00673   sinfo_free_propertylist(&plist);
00674   dispersion = 1e4 * dispersion;   // A/pix 
00675   h=PLANCK;                        // J s 
00676   c=SPEED_OF_LIGHT;                 // m / s 
00677   k= BOLTZMANN;                     // J / K 
00678   temp=1e4;                         // K, star temperature 
00679   surface=1e4 * TELESCOPE_SURFACE;  // cm * cm         
00680   gain = 2.42;                      // e / ADU 
00681 
00682   num= gain * pow(10,(mag/2.5)) * 1e7 * h * c; // e / ADU erg s  m / s 
00683   den= dit * surface * dispersion * wc * 1e-6; // s cm cm  A / pix m  
00684 
00685   fctr = num/den;                              // e pix/ADU erg/s/cm/cm/A
00686   sinfo_msg("Dispersion: %f A/pix",dispersion);
00687         
00688      // erg/s/cm/cm/A 
00689      // int_spct * fctr  = e erg/s/cm/cm/A 
00690 
00691   check_nomsg(cpl_table_new_column(*tbl_spectrum,"efficiency",CPL_TYPE_FLOAT));
00692 
00693 
00694   for (i=0;i<cpl_imagelist_get_size(cube);i++) {
00695   
00696     check_nomsg(wave=cpl_table_get_float(*tbl_spectrum,"wavelength",i,&status));
00697     check_nomsg(int_spct=cpl_table_get_float(*tbl_spectrum,"counts_bkg",
00698                                              i,&status));
00699 
00700     bb_flux_norm = pow(wc/wave,5) * (exp(h*c/(wc * 1e-6*k*temp))-1)/ 
00701                                        (exp(h*c/(wave*1e-6*k*temp))-1);
00702 
00703 
00704     efficiency = int_spct*fctr/(bb_flux_norm *f0);
00705     check_nomsg(cpl_table_set_float(*tbl_spectrum,"efficiency",i,efficiency));
00706 
00707 
00708   } 
00709 
00710 return 0;
00711 
00712  cleanup:
00713   sinfo_free_propertylist(&plist);
00714  return -1;
00715 
00716 }
00717 
00718 
00719 

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