00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
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
00025
00026
00027 #define STREHL_M1 8.0
00028 #define STREHL_M2 1.1
00029 #define STREHL_BOX_SIZE 64
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
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
00137
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
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
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
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
00238 nsample=(int)((wend-wstart-wstep)/wstep);
00239
00240
00241
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
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
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
00279
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
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
00301
00302
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
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
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