psf.c

00001 /*----------------------------------------------------------------------------
00002  
00003      File name    :     psf.c
00004    Author       :   A. Modigliani
00005    Created on   :   Sep 17, 2003
00006    Description  : 
00007 
00008  psf.py does the image reconstruction of a set of sky-subtracted, flatfielded, 
00009  bad pixel corrected and slope of the spectra aligned exposures of a bright 
00010  star with continuum spectrum. The resulting image can be used to determine 
00011  the PSF
00012  ---------------------------------------------------------------------------*/
00013 
00014 /*----------------------------------------------------------------------------
00015                                 Includes
00016  ---------------------------------------------------------------------------*/
00017 #include "psf.h"
00018 #include "sinfoni_pro_save.h"
00019 #include "psf_ini.h"
00020 #include "psf_ini_by_cpl.h"
00021 #include "fitshead.h"
00022 #include "utilities_scired.h"
00023 /*----------------------------------------------------------------------------
00024                                 Defines
00025  ---------------------------------------------------------------------------*/
00026 
00027 #define STREHL_M1                       8.0
00028 #define STREHL_M2                       1.1
00029 #define STREHL_BOX_SIZE                 64
00030 
00031 /*----------------------------------------------------------------------------
00032                              Function Definitions
00033  ---------------------------------------------------------------------------*/
00034 
00035 /*----------------------------------------------------------------------------
00036    Function     :       psf()
00037    In           :       ini_file: file name of according .ini file
00038    Out          :       integer (0 if it worked, -1 if it doesn't) 
00039    Job          : 
00040 
00041 psf.py does the image reconstruction of a set of sky-subtracted, flatfielded, 
00042  bad pixel corrected and slope of the spectra aligned exposures of a bright 
00043  star with continuum spectrum. The resulting image can be used to determine 
00044  the PSF
00045 
00046  ---------------------------------------------------------------------------*/
00047 
00048 int psf (const char* plugin_id,cpl_parameterlist* config,cpl_frameset* sof)
00049 {
00050     const char* _id = "psf";
00051     psf_config * cfg ;
00052 
00053     OneCube* coadd_cube=NULL;
00054     OneImage * tmp_img=NULL;
00055 
00056     OneImage * out_image=NULL ;
00057 
00058     int naxis3=0;
00059     int nsample=0;
00060     int i = 0;
00061     int lx=0;
00062     int ly=0;
00063 
00064     int ni=0;
00065     int wllx=0;
00066     int wlly=0;
00067     int wurx=0;
00068     int wury=0;
00069     int psf_sz=3;
00070     int max_ima_x=0;
00071     int max_ima_y=0;
00072     int r=0;
00073     int dr=0;
00074     int rmin=0;
00075     int rdifr=0;
00076 
00077     float center_x=0;
00078     float center_y=0;
00079 
00080     double dispersion=0.;
00081     double centralWave=0.;
00082     double wrange=0;
00083     double wstart=0;
00084     double wstep=0;
00085     double wend=0;
00086     double ws=0;
00087     double we=0;
00088     double pix_scale=0;
00089     double max_ima_cx=0;
00090     double max_ima_cy=0;
00091     double norm=0;
00092     double xcen=0;
00093     double ycen=0;
00094     double sig_x=0;
00095     double sig_y=0;
00096     double fwhm_x=0;
00097     double fwhm_y=0;
00098     double flux=0.;
00099     double flux_max=0.;
00100     double lam=0;
00101     double dlam=0;
00102     double pscale=0.;
00103     double strehl=0;
00104     double strehl_err=0;
00105     double star_bg=0;
00106     double star_peak=0;
00107     double star_flux=0;
00108     double psf_peak=0;
00109     double psf_flux=0;
00110     double bg_noise=0;
00111     double bkg=0;
00112 
00113     double strehl_star_radius=0;
00114     double strehl_bg_r1=0;
00115     double strehl_bg_r2=0;
00116 
00117 
00118 
00119     char psf_cub_filename[FILE_NAME_SZ];
00120     char band[FILE_NAME_SZ];
00121     char key_value[FILE_NAME_SZ];
00122     char spat_res[FILE_NAME_SZ];
00123 
00124 
00125     cpl_table* ao_performance=NULL;
00126     cpl_table* enc_energy=NULL;
00127     cpl_frame* frame=NULL;
00128     cpl_table* qclog_tbl=NULL;
00129     cpl_image* img=NULL;
00130     cpl_image* imgw=NULL;
00131     cpl_image* img_dup=NULL;
00132     cpl_frameset* stk=NULL;
00133 
00134     /* 
00135        -----------------------------------------------------------------
00136        1) parse the file names and parameters to the psf_config data 
00137           structure cfg
00138        -----------------------------------------------------------------
00139      */
00140     cpl_msg_info(_id,"Parsing cpl input");
00141     stk=cpl_frameset_new();
00142     cfg = parse_cpl_input_psf(config,sof,&stk) ;
00143     if(cpl_error_get_code() != CPL_ERROR_NONE) {
00144       cpl_msg_error(_id,"error parsing cpl input");
00145       cpl_msg_error(_id,(char* ) cpl_error_get_message());
00146       return -1;
00147     }
00148 
00149     if (cfg == NULL)
00150     {
00151         cpl_msg_error(_id,"could not parse cpl input!\n") ;
00152         cpl_frameset_delete(stk);
00153         return -1 ;
00154     }
00155     
00156    if(is_fits_file(cfg->inFrame) != 1) {
00157        cpl_msg_error(_id,"Input file %s is not FITS",cfg->inFrame);
00158        cpl_frameset_delete(stk);
00159        return -1;
00160    } else {
00161       strcpy(psf_cub_filename,cfg->inFrame);
00162    }
00163     if(NULL != cpl_frameset_find(sof,PRO_COADD_PSF)) {
00164       frame = cpl_frameset_find(sof,PRO_COADD_PSF);
00165     } else if(NULL != cpl_frameset_find(sof,PRO_OBS_PSF)) {
00166       frame = cpl_frameset_find(sof,PRO_OBS_PSF);
00167     } else {
00168       cpl_msg_error(_id,"Frame %s  of %s not found!", PRO_COADD_PSF, PRO_OBS_PSF);
00169       return -1;
00170     }
00171 
00172 
00173  sinfoni_get_spatial_res(frame,spat_res);
00174  sinfoni_get_band(frame,band);
00175  pix_scale=atof(spat_res); 
00176  /* factor 2 due to change of detector to 2K */
00177  pscale=0.5*pix_scale;
00178 
00179  if (strcmp(band,"H+K") == 0) {
00180    lam=1.950;
00181  } else if (strcmp(band,"K") == 0) {
00182    lam=2.175;
00183  } else if (strcmp(band,"J") == 0) {
00184    lam=1.225;
00185  } else if (strcmp(band,"H") == 0) { 
00186    lam=1.675;
00187  }
00188 
00189  dlam=spiffi_get_dispersion(psf_cub_filename);
00190 
00191  /*   printf("diffraction limit=%f\n",3600*180/PI_NUMB*lam*1.e-6*1.22/8/pscale); */
00192    rdifr=floor(3600*180/PI_NUMB*lam*1.e-6*1.22/8/pscale+0.5);
00193  if (pix_scale==0.025) {
00194    ni=10;
00195    rmin=rdifr;
00196    dr=rmin;
00197  } else {
00198    ni=15;
00199    cpl_msg_warning(_id,"Reset diffraction limit");
00200    rdifr=10;
00201    rmin=1;
00202    dr=2;
00203  }
00204 
00205  cpl_msg_info(_id,"Diffraction limit: %d",rdifr);
00206 
00207  /* rdifr=floor(25*pscale+0.5); */
00208 
00209  strehl_star_radius=25*pscale;
00210  strehl_bg_r1=25*pscale;
00211  strehl_bg_r2=30*pscale;
00212 
00213      
00214  coadd_cube = load_cube(psf_cub_filename);
00215 
00216 
00217  out_image=medianCube(coadd_cube);
00218  if (out_image == NULL) {
00219     cpl_msg_error(_id," could not do medianCube()\n");
00220     return -1;
00221  }
00222  lx = out_image->lx;
00223  ly = out_image->ly;
00224 
00225  center_x = lx / 2. + 0.5;
00226  center_y = ly / 2. + 0.5;
00227 
00228 
00229  dispersion=spiffi_get_dispersion(psf_cub_filename);
00230  centralWave=spiffi_get_cenLambda(psf_cub_filename);
00231  naxis3=spiffi_get_naxis3(psf_cub_filename);
00232  wrange=dispersion*naxis3;
00233 
00234  wstart = centralWave - (float) (coadd_cube->np / 2)*dispersion+dispersion;
00235  wend  =wstart + dispersion * coadd_cube->np;
00236  wstep=0.025;
00237  /* note: -wstep as we do not hit the borders where the gaussian fit has a problem */
00238  nsample=(int)((wend-wstart-wstep)/wstep);
00239  /*
00240  printf("dis=%f cw=%f n3=%d wstart=%f  wend=%f nsample=%d pscale=%f\n",
00241          dispersion,centralWave,naxis3,wstart,wend,nsample,pscale);
00242  */
00243  ao_performance = cpl_table_new(nsample);
00244  cpl_table_new_column(ao_performance,"wavelength", CPL_TYPE_DOUBLE);
00245  cpl_table_new_column(ao_performance,"strehl" , CPL_TYPE_DOUBLE);
00246  cpl_table_new_column(ao_performance,"strehl_error" , CPL_TYPE_DOUBLE);
00247 
00248  
00249  for(i=1;i<nsample;i++) {
00250 
00251       ws=wstart+wstep*i;
00252       we=ws+wstep;
00253       /* printf("ws=%f we=%f\n",ws,we); */
00254       lam = (double)0.5*(ws+we);
00255       dlam=wstep;
00256       tmp_img=averageCubeToImageBetweenWaves(coadd_cube,dispersion,centralWave,ws,we);
00257       imgw=cpl_image_wrap_float(tmp_img->lx,tmp_img->ly,tmp_img->data);
00258    
00259       img=cpl_image_duplicate(imgw);
00260       cpl_image_unwrap(imgw);
00261       img_dup=cpl_image_duplicate(img);
00262       sinfoni_clean_nan(&img_dup);
00263       cpl_image_get_maxpos(img_dup,&max_ima_x,&max_ima_y);
00264       cpl_image_delete(img_dup);
00265 
00266       wllx=max_ima_x-psf_sz;
00267       wlly=max_ima_y-psf_sz;
00268       wurx=max_ima_x+psf_sz;
00269       wury=max_ima_y+psf_sz;
00270       
00271       /* cpl_image_get_maxpos_window(img,wllx,wlly,wurx,wury,&max_ima_x,&max_ima_y); */
00272       max_ima_cx=cpl_image_get_centroid_x_window(img,wllx,wlly,wurx,wury);
00273       max_ima_cy=cpl_image_get_centroid_y_window(img,wllx,wlly,wurx,wury); 
00274 
00275       cpl_image_fit_gaussian(img,max_ima_x,max_ima_y,psf_sz,
00276                             &norm,&xcen,&ycen,&sig_x,&sig_y,&fwhm_x,&fwhm_y);
00277       /*
00278       printf("norm %f,xcen %f,ycen %f,sig_x %f,sig_y, %f fwhm_x, %f fwhm_y %f\n",
00279               norm,xcen,ycen,sig_x,sig_y,fwhm_x,fwhm_y);
00280       */
00281 
00282       bkg=irplib_strehl_ring_background(img,max_ima_x,max_ima_y,
00283                    19,21,IRPLIB_BG_METHOD_AVER_REJ) ;
00284 
00285       /*
00286       printf("lam %f dlam %f pscale=%f size%d max_x %d max_y %d star_radius %f bg_r1 %f bg_r2 %f hsz %d nsmpl %d\n",lam,dlam,pscale,STREHL_BOX_SIZE,max_ima_x,max_ima_y,strehl_star_radius,strehl_bg_r1,strehl_bg_r2,NOISE_HSIZE,NOISE_NSAMPLES);
00287       */
00288       irplib_strehl_compute(img,STREHL_M1,STREHL_M2,lam,dlam,pscale,STREHL_BOX_SIZE,
00289                       max_ima_x,max_ima_y,strehl_star_radius,strehl_bg_r1,
00290                                       strehl_bg_r2,NOISE_HSIZE,NOISE_NSAMPLES,
00291                          &strehl,&strehl_err,&star_bg,&star_peak,&star_flux,
00292                 &psf_peak,&psf_flux,&bg_noise);
00293 
00294 
00295 
00296       cpl_image_delete(img);
00297       destroy_image(tmp_img);
00298 
00299       /*
00300       cpl_msg_info(_id,"strehl %f strehl_err %f star_bg %f star_peak %f star_flux %f psf_peak %f psf_flux %f bg_noise %f",strehl,strehl_err,star_bg,star_peak,star_flux,psf_peak,psf_flux,bg_noise);
00301       */
00302       /* printf("strehl=%f\n",strehl); */
00303       if((qfits_isnan(lam) ==0) &&(qfits_isnan(lam) ==0) && (qfits_isnan(lam) ==0)) {
00304     cpl_table_set_double(ao_performance,"wavelength",i,lam);
00305     cpl_table_set_double(ao_performance,"strehl",i,strehl);
00306     cpl_table_set_double(ao_performance,"strehl_error",i,strehl_err);
00307       }
00308 
00309  }
00310  destroy_cube(coadd_cube);
00311 
00312 
00313  if(-1 == sinfoni_pro_dump_tbl(ao_performance,stk,sof,
00314                                PSF_AO_PERFORMANCE_OUT_FILENAME,
00315                                PRO_AO_PERFORMANCE,plugin_id,config)) {
00316      cpl_msg_error(_id,"cannot dump tbl %s", PSF_AO_PERFORMANCE_OUT_FILENAME);
00317 
00318      destroy_image(out_image);
00319      cpl_table_delete(ao_performance);
00320      free_psf(cfg);
00321      cpl_frameset_delete(stk);
00322      return -1;
00323  }
00324  cpl_table_delete(ao_performance);
00325  imgw=cpl_image_wrap_float(out_image->lx,out_image->ly,out_image->data);
00326  img=cpl_image_duplicate(imgw);
00327  cpl_image_unwrap(imgw);
00328  /* FWHM computation */
00329 
00330 
00331  img_dup=cpl_image_duplicate(img);
00332  sinfoni_clean_nan(&img_dup);
00333  cpl_image_get_maxpos(img_dup,&max_ima_x,&max_ima_y);
00334  cpl_image_delete(img_dup);
00335 
00336  wllx=max_ima_x-psf_sz;
00337  wlly=max_ima_y-psf_sz;
00338  wurx=max_ima_x+psf_sz;
00339  wury=max_ima_y+psf_sz;
00340 
00341  max_ima_cx=cpl_image_get_centroid_x_window(img,wllx,wlly,wurx,wury);
00342  max_ima_cy=cpl_image_get_centroid_y_window(img,wllx,wlly,wurx,wury);
00343 
00344  cpl_image_fit_gaussian(img,max_ima_x,max_ima_y,psf_sz,
00345                             &norm,&xcen,&ycen,&sig_x,&sig_y,&fwhm_x,&fwhm_y);
00346 
00347 
00348  enc_energy = cpl_table_new(ni);
00349  cpl_table_new_column(enc_energy,"r_pix", CPL_TYPE_INT);
00350  cpl_table_new_column(enc_energy,"r_mas", CPL_TYPE_DOUBLE);
00351  cpl_table_new_column(enc_energy,"r_dif", CPL_TYPE_DOUBLE);
00352  cpl_table_new_column(enc_energy,"abs_energy" , CPL_TYPE_DOUBLE);
00353  cpl_table_new_column(enc_energy,"rel_energy" , CPL_TYPE_DOUBLE);
00354  /* encircled energy computation */
00355  bkg=irplib_strehl_ring_background(img,max_ima_x,max_ima_y,
00356                     25,30,IRPLIB_BG_METHOD_AVER_REJ) ;
00357  r=rmin+(ni-1)*dr;
00358  flux_max=irplib_strehl_disk_flux(img,max_ima_x,max_ima_y,r,bkg);
00359  r=rmin;
00360 
00361  for(i=0; i<ni; i++)
00362    {
00363      flux=irplib_strehl_disk_flux(img,max_ima_x,max_ima_y,r,bkg);
00364      cpl_table_set_int(enc_energy,"r_pix",i,r);
00365      cpl_table_set_double(enc_energy,"r_mas",i,r*pscale);
00366      cpl_table_set_double(enc_energy,"r_dif",i,r/rdifr);
00367      cpl_table_set_double(enc_energy,"abs_energy",i,flux);
00368      cpl_table_set_double(enc_energy,"rel_energy",i,flux/flux_max);
00369      r+=dr;
00370 
00371    }
00372  
00373  cpl_msg_info(_id,"max ima=%d %d\n",max_ima_x,max_ima_y);
00374  cpl_msg_info(_id,"centroid ima=%f %f\n",max_ima_cx,max_ima_cy);
00375  cpl_msg_info(_id,"gauss info=%f %f %f %f %f %f %f\n",
00376                           norm,xcen,ycen,sig_x,sig_y,fwhm_x,fwhm_y);
00377  
00378  if(-1 == sinfoni_pro_dump_tbl(enc_energy,stk,sof,PSF_ENC_ENERGY_OUT_FILENAME,
00379                                PRO_ENC_ENERGY,plugin_id,config)) {
00380      cpl_msg_error(_id,"cannot dump ima %s", PSF_ENC_ENERGY_OUT_FILENAME);
00381      destroy_image(out_image);
00382      cpl_image_delete(img);
00383      cpl_table_delete(enc_energy);
00384      free_psf(cfg);
00385      cpl_frameset_delete(stk);
00386      return -1;
00387  }
00388  cpl_table_delete(enc_energy);
00389  
00390  qclog_tbl = sinfoni_qclog_init(2);
00391  sprintf(key_value,"%f",fwhm_x);
00392  sinfoni_qclog_add(qclog_tbl,0,"QC FWHMX","CPL_TYPE_DOUBLE",key_value,"QC FWHM X");
00393  sprintf(key_value,"%f",fwhm_y);
00394  sinfoni_qclog_add(qclog_tbl,1,"QC FWHMY","CPL_TYPE_DOUBLE",key_value,"QC FWHM Y");
00395 
00396  if(-1 == sinfoni_pro_save_ima(img,stk,sof,cfg->outName,PRO_PSF,
00397                                qclog_tbl,plugin_id,config)) {
00398      cpl_msg_error(_id,"cannot dump ima %s", cfg->outName);
00399      destroy_image(out_image);
00400      cpl_image_delete(img);
00401      cpl_table_delete(qclog_tbl);
00402      free_psf(cfg);
00403      cpl_frameset_delete(stk);
00404      return -1;
00405  }
00406  cpl_table_delete(qclog_tbl);
00407  set_wcs_image(img,cfg->outName,center_x, center_y);
00408  cpl_image_delete(img);
00409  destroy_image(out_image);
00410  cpl_frameset_delete(stk);
00411  free_psf(cfg);
00412  return 0;
00413 
00414 }
00415 

Generated on Wed Oct 26 13:08:54 2005 for SINFONI Pipeline Reference Manual by doxygen1.2.13.1 written by Dimitri van Heesch, © 1997-2001