00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020 #include "stdstar.h"
00021 #include "standstar_ini_by_cpl.h"
00022 #include "sinfoni_pro_save.h"
00023 #include "utilities_scired.h"
00024 #include <irplib_std.h>
00025 #include <math.h>
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035 int
00036 sinfoni_clean_nan(cpl_image** im);
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048 int stdstar (const char* plugin_id,cpl_parameterlist* config, cpl_frameset* sof)
00049 {
00050
00051 const char* _id = "stdstar";
00052 standstar_config * cfg ;
00053 OneCube * cube ;
00054 OneCube * cube_object ;
00055
00056 OneImage ** list_object ;
00057 OneImage ** spectrum ;
00058 OneImage * averaged_spectrum ;
00059
00060 FILE * factorfile ;
00061 char * name ;
00062 int fra=0;
00063 float exptime=0 ;
00064 double convfactor=0;
00065 double cleanfactor=0;
00066 char band[FILE_NAME_SZ];
00067 float* factor;
00068 cpl_parameter* p;
00069
00070 cpl_image* asp_img=NULL;
00071 cpl_image* asp_imgw=NULL;
00072
00073 cpl_frameset* raw=NULL;
00074
00075 cpl_frame* frame=NULL;
00076 cpl_image* std_med_ima=NULL;
00077 cpl_image* std_med_dup=NULL;
00078
00079 char* std_med_filename=NULL;
00080 char* std_cub_filename=NULL;
00081 cpl_table* qclog_tbl=NULL;
00082 char* key_value=NULL;
00083 double max_ima_cx=0;
00084 double max_ima_cy=0;
00085 int max_ima_x=0;
00086 int max_ima_y=0;
00087 double norm=0;
00088 double xcen=0;
00089 double ycen=0;
00090 double sig_x=0;
00091 double sig_y=0;
00092 double fwhm_x=0;
00093 double fwhm_y=0;
00094 cpl_table* tbl_spectrum=NULL;
00095 double disp=0;
00096 cpl_propertylist* plist=NULL;
00097 double h=0;
00098 double c=0;
00099 double f0=0;
00100 double wc=0;
00101 double k= 0;
00102 double temp=0;
00103 double dispersion=0;
00104 double mag=0;
00105 double surface=0;
00106 int i=0;
00107 float int_spct=0;
00108 float bb_flux_norm=0;
00109 float wave=0;
00110 float efficiency=0;
00111 float dit=0;
00112 float fctr=0;
00113 int status=0;
00114 float gain=2.42;
00115 float den=0.;
00116 float num=0.;
00117 int wllx=0;
00118 int wlly=0;
00119 int wurx=0;
00120 int wury=0;
00121 int psf_sz=20;
00122 int qc_info=0;
00123 double std_ra=0;
00124 double std_dec=0;
00125 irplib_band std_band;
00126 int ima_szx=0;
00127 int ima_szy=0;
00128 int check1=0;
00129 int check2=0;
00130 int check3=0;
00131
00132 float cenpix = 0;
00133 float cenLambda = 0;
00134
00135
00136 if(cpl_error_get_code() != CPL_ERROR_NONE) {
00137 cpl_msg_error(_id,"error checking code");
00138 cpl_msg_error(_id,(char* ) cpl_error_get_message());
00139 return -1;
00140 }
00141
00142 raw=cpl_frameset_new();
00143 cfg = parse_cpl_input_standstar(config,sof,&raw) ;
00144
00145 if (cfg == NULL)
00146 {
00147 cpl_msg_error(_id,"could not parse cpl input!\n") ;
00148 cpl_frameset_delete(raw);
00149 return -1 ;
00150 }
00151
00152 list_object = new_image_list (cfg->nframes);
00153 if (list_object == NULL) {
00154 cpl_msg_error(_id,"could not allocate memory\n");
00155 stdstar_free (cfg);
00156 cpl_frameset_delete(raw);
00157 return -1;
00158 }
00159
00160 if (cfg->convInd == 1) {
00161 factor = new_float_array(cfg->nframes);
00162 }
00163
00164 std_cub_filename = (char*) cpl_calloc(MAX_NAME_SIZE,sizeof(char));
00165 if(NULL != cpl_frameset_find(sof,PRO_COADD_STD)) {
00166 frame = cpl_frameset_find(sof,PRO_COADD_STD);
00167 strcpy(std_cub_filename,cpl_frame_get_filename(frame));
00168 } else {
00169 cpl_msg_error(_id,"Frame %s not found! Exit!", PRO_COADD_STD);
00170 destroy_array(factor);
00171 cpl_free(std_cub_filename);
00172 destroy_list(list_object);
00173 stdstar_free (cfg);
00174 cpl_frameset_delete(raw);
00175 return -1;
00176 }
00177
00178
00179
00180
00181
00182 if ((cpl_error_code)((plist = cpl_propertylist_load(std_cub_filename, 0)) == NULL)) {
00183 cpl_msg_error(_id, "getting header from reference ima frame %s",std_cub_filename);
00184 cpl_propertylist_delete(plist) ;
00185 destroy_array(factor);
00186 destroy_list(list_object);
00187 stdstar_free (cfg);
00188 cpl_frameset_delete(raw);
00189 return -1 ;
00190 }
00191
00192
00193 cenpix = spiffi_get_cenpix(std_cub_filename);
00194 cenLambda = spiffi_get_cenLambda(std_cub_filename);
00195 dispersion = spiffi_get_dispersion(std_cub_filename);
00196
00197 cpl_free(std_cub_filename);
00198
00199
00200
00201
00202
00203 if (cpl_propertylist_contains(plist, KEY_NAME_CDELT3)) {
00204 disp=cpl_propertylist_get_double(plist, KEY_NAME_CDELT3);
00205 } else {
00206 cpl_msg_warning(_id,"Keyword %s not found.",KEY_NAME_CDELT3);
00207 }
00208 cpl_propertylist_delete(plist) ;
00209
00210
00211 std_med_filename = (char*) cpl_calloc(MAX_NAME_SIZE,sizeof(char));
00212
00213 if(NULL != cpl_frameset_find(sof,PRO_MED_COADD_STD)) {
00214 frame = cpl_frameset_find(sof,PRO_MED_COADD_STD);
00215 strcpy(std_med_filename,cpl_frame_get_filename(frame));
00216 } else {
00217 cpl_msg_error(_id,"Frame %s not found! Exit!", PRO_MED_COADD_STD);
00218 cpl_free( std_med_filename);
00219 destroy_array(factor);
00220 destroy_list(list_object);
00221 stdstar_free (cfg);
00222 cpl_frameset_delete(raw);
00223 return -1;
00224 }
00225
00226
00227
00228 std_med_ima=cpl_image_load(std_med_filename,CPL_TYPE_FLOAT,0,0);
00229 std_med_dup=cpl_image_duplicate(std_med_ima);
00230 sinfoni_clean_nan(&std_med_dup);
00231 cpl_image_get_maxpos(std_med_dup,&max_ima_x,&max_ima_y);
00232 cpl_image_delete(std_med_dup);
00233 ima_szx=cpl_image_get_size_x(std_med_ima);
00234 ima_szy=cpl_image_get_size_y(std_med_ima);
00235 wllx= ((max_ima_x-psf_sz)>0) ? (max_ima_x-psf_sz) : 1;
00236 wlly= ((max_ima_y-psf_sz)>0) ? (max_ima_y-psf_sz) : 1;
00237 wurx= ((max_ima_x+psf_sz)<ima_szx) ? (max_ima_x+psf_sz) : ima_szx ;
00238 wury= ((max_ima_y+psf_sz)<ima_szy) ? (max_ima_y+psf_sz) : ima_szy ;
00239
00240 cpl_msg_info(_id,"wllx=%d wlly=%d wurx=%d wury=%d\n",wllx,wlly,wurx,wury);
00241
00242 max_ima_cx=cpl_image_get_centroid_x_window(std_med_ima,wllx,wlly,wurx,wury);
00243 max_ima_cy=cpl_image_get_centroid_y_window(std_med_ima,wllx,wlly,wurx,wury);
00244
00245
00246
00247 if(
00248 ((max_ima_x-psf_sz) < 1) ||
00249 ((max_ima_y-psf_sz) < 1) ||
00250 ((max_ima_x+psf_sz) > ima_szx) ||
00251 ((max_ima_x+psf_sz) > ima_szy)
00252 )
00253 {
00254 psf_sz = (psf_sz < (max_ima_x-1)) ? psf_sz : (max_ima_x-1);
00255 psf_sz = (psf_sz < (max_ima_y-1)) ? psf_sz : (max_ima_y-1);
00256 psf_sz = (psf_sz < ima_szx-max_ima_x) ? psf_sz : (ima_szx-max_ima_x);
00257 psf_sz = (psf_sz < ima_szy-max_ima_y) ? psf_sz : (ima_szy-max_ima_y);
00258 }
00259
00260
00261 cpl_image_fit_gaussian(std_med_ima,max_ima_x,max_ima_y,psf_sz,
00262 &norm,&xcen,&ycen,&sig_x,&sig_y,&fwhm_x,&fwhm_y);
00263 cpl_image_delete(std_med_ima);
00264
00265
00266
00267
00268
00269 qclog_tbl = sinfoni_qclog_init(4);
00270 key_value = cpl_calloc(FILE_NAME_SZ,sizeof(char));
00271 sprintf(key_value,"%f",fwhm_x);
00272
00273 sinfoni_qclog_add(qclog_tbl,0,"QC FWHMX","CPL_TYPE_DOUBLE",key_value,"STD star FWHM on X");
00274 sprintf(key_value,"%f",fwhm_y);
00275 sinfoni_qclog_add(qclog_tbl,1,"QC FWHMY","CPL_TYPE_DOUBLE",key_value,"STD star FWHM on Y");
00276
00277
00278 cpl_msg_info(_id,"max ima=%d %d psf_sz=%d\n",max_ima_x,max_ima_y,psf_sz);
00279 cpl_msg_info(_id,"centroid ima=%f %f\n",max_ima_cx,max_ima_cy);
00280 cpl_msg_info(_id,"gauss info=%f %f %f %f %f %f %f\n",
00281 norm,xcen,ycen,sig_x,sig_y,fwhm_x,fwhm_y);
00282
00283
00284 cfg -> halfbox_x = (floor)(0.5*(fwhm_x+fwhm_y)*cfg->fwhm_factor+0.5);
00285 cfg -> llx = (int)(xcen-cfg->halfbox_x);
00286 cfg -> llx = (cfg -> llx > 0 ) ? cfg -> llx : 1;
00287
00288 if((cfg->llx+2*cfg->halfbox_x) > ima_szx) {
00289 cfg -> halfbox_x=(int) ((ima_szx-cfg->llx-1)/2);
00290 check1++;
00291 } else {
00292 }
00293
00294 cfg -> halfbox_y = (floor)(0.5*(fwhm_x+fwhm_y)*cfg->fwhm_factor+0.5);
00295 cfg -> lly = (int)(ycen-cfg->halfbox_y);
00296 cfg -> lly = (cfg -> lly > 0 ) ? cfg -> lly : 1;
00297 if((cfg->lly+2*cfg->halfbox_y) > ima_szy) {
00298 cfg -> halfbox_y=(int) ((ima_szy-cfg->lly-1)/2);
00299 check1++;
00300 } else {
00301 }
00302
00303 sinfoni_qclog_add(qclog_tbl,2,"QC CHECK1","CPL_TYPE_INT",key_value,"Check evaluation box");
00304
00305 cpl_msg_info(_id,"llx= %d lly= %d\n",cfg->llx, cfg->lly);
00306 cpl_msg_info(_id,"halfbox_x=%d halfbox_y=%d\n",cfg->halfbox_x,cfg->halfbox_y);
00307 cpl_free(std_med_filename);
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321 cpl_msg_info(_id,"Extraction");
00322 spectrum = (OneImage**) cpl_calloc (cfg -> nframes, sizeof(OneImage*)) ;
00323
00324
00325
00326 sinfoni_get_band(frame,band);
00327 name = (char*) cpl_frame_get_filename(frame);
00328
00329
00330 std_ra =(double) spiffi_get_RA(name);
00331 std_dec=(double) spiffi_get_DEC(name);
00332
00333 if (strcmp(band,"H+K") == 0) {
00334 f0=8.00e-11;
00335 wc=1.950;
00336 std_band=BAND_K;
00337 } else if (strcmp(band,"K") == 0) {
00338 f0=4.10e-11;
00339 wc=2.175;
00340 std_band=BAND_K;
00341 } else if (strcmp(band,"J") == 0) {
00342 f0=3.11e-10;
00343 wc=1.225;
00344 std_band=BAND_J;
00345 } else if (strcmp(band,"H") == 0) {
00346 f0=1.15e-10;
00347 wc=1.675;
00348 std_band=BAND_H;
00349 }
00350
00351
00352 irplib_std_get_mag(std_ra, std_dec, std_band, &mag) ;
00353 cfg->mag=mag;
00354 cpl_msg_info(_id,"Catalogue magnitude=%g",cfg->mag);
00355
00356
00357
00358
00359
00360
00361
00362 for (fra=0; fra < cfg->nframes; fra++) {
00363 name = cfg->inFrameList[fra];
00364 if(is_fits_file(name) != 1) {
00365 cpl_msg_error(_id,"Input file %s is not FITS",name);
00366 return -1;
00367 }
00368
00369
00370 cube = load_cube(name);
00371 if (cube == NULL) {
00372 cpl_msg_error (_id, "could not load data cube" );
00373 return -1;
00374 }
00375
00376 if (exptime == FLAG) {
00377 cpl_msg_error(_id,"could not find exposure time in the fits header" );
00378 return -1;
00379 }
00380
00381
00382 exptime = spiffi_get_ditndit(name);
00383
00384 cpl_msg_info(_id,"cfg->gain %f",cfg->gain);
00385 tbl_spectrum = cpl_table_new(cube->np);
00386
00387
00388 spectrum[fra]=optimalExtractionFromCube( cube,
00389 cfg->llx,
00390 cfg->lly,
00391 cfg->halfbox_x,
00392 cfg->halfbox_y,
00393 cfg->fwhm_factor,
00394 BKG_VARIANCE,
00395 SKY_FLUX,
00396 cfg->gain,
00397 exptime,
00398 name,
00399 &tbl_spectrum, qc_info, &check2);
00400
00401
00402
00403
00404 sprintf(key_value,"%d",check2);
00405 sinfoni_qclog_add(qclog_tbl,3,"QC CHECK2","CPL_TYPE_INT",key_value,"Check evaluation box");
00406
00407
00408 if (spectrum[fra] == NULL){
00409 cpl_msg_error(_id,"could not do optimalExtractionFromCube" );
00410 stdstar_free (cfg);
00411 destroy_cube(cube);
00412 cpl_table_delete(tbl_spectrum);
00413 destroy_list(list_object);
00414 cpl_frameset_delete(raw);
00415 cpl_table_delete(qclog_tbl);
00416 cpl_free(key_value);
00417 destroy_array(factor);
00418 cpl_free(spectrum);
00419
00420 return 0;
00421 }
00422 insert_image(list_object, spectrum[fra], fra);
00423
00424
00425 dispersion = spiffi_get_dispersion(name);
00426 dispersion = 1e4 * dispersion;
00427 h=PLANCK;
00428 c=SPEED_OF_LIGHT;
00429 k= BOLTZMANN;
00430 temp=1e4;
00431 surface=1e4 * TELESCOPE_SURFACE;
00432 dit = spiffi_get_dit(name);
00433 mag=cfg->mag;
00434
00435
00436
00437 gain = 2.42;
00438
00439 num= gain * pow(10,(mag/2.5)) * 1e7 * h * c;
00440 den= dit * surface * dispersion * wc * 1e-6;
00441
00442 fctr = num/den;
00443 cpl_msg_info(_id,"Dispersion: %f A/pix",disp);
00444
00445
00446
00447
00448 cpl_table_new_column(tbl_spectrum,"bb_flux_norm", CPL_TYPE_FLOAT);
00449 cpl_table_new_column(tbl_spectrum,"efficiency" , CPL_TYPE_FLOAT);
00450
00451
00452 for (i=0;i<cube->np;i++) {
00453
00454 wave = cpl_table_get_float(tbl_spectrum,"wavelength",i,&status);
00455 int_spct = cpl_table_get_float(tbl_spectrum,"counts_bkg",i,&status);
00456
00457 bb_flux_norm = pow(wc/wave,5) * (exp(h*c/(wc * 1e-6*k*temp))-1)/
00458 (exp(h*c/(wave*1e-6*k*temp))-1);
00459
00460
00461 efficiency = int_spct*fctr/(bb_flux_norm *f0);
00462 cpl_table_set_float(tbl_spectrum,"bb_flux_norm",i,bb_flux_norm);
00463 cpl_table_set_float(tbl_spectrum,"efficiency",i,efficiency);
00464
00465
00466 }
00467
00468 if(-1 == sinfoni_pro_save_tbl(tbl_spectrum,raw,sof,
00469 (char*)STDSTAR_OUT_TABLE,
00470 PRO_STD_STAR_SPECTRA,qclog_tbl,
00471 plugin_id,config)) {
00472 cpl_msg_error(_id,"cannot dump ima %s", "out_std_star_spectrum.tfits");
00473 }
00474 cpl_table_delete(qclog_tbl);
00475 cpl_free(key_value);
00476
00477
00478
00479 if(spectrum) {
00480 cpl_free(spectrum);
00481 }
00482
00483
00484
00485 if (cfg->convInd == 1) {
00486 cpl_msg_info(_id,"Determines convertion factor");
00487
00488 convfactor = determineConversionFactor( cube,
00489 cfg->mag,
00490 exptime,
00491 cfg->llx,
00492 cfg->lly,
00493 cfg->halfbox_x,
00494 cfg->halfbox_y, &check3 );
00495
00496
00497 if (convfactor < -100000.) {
00498 cpl_msg_error(_id,"could not do determineConversionFactor!" );
00499 stdstar_free (cfg);
00500 destroy_cube(cube);
00501 cpl_table_delete(tbl_spectrum);
00502 destroy_list(list_object);
00503 cpl_frameset_delete(raw);
00504 destroy_array(factor);
00505 return 0;
00506 }
00507 array_set_value(factor, convfactor, fra);
00508
00509 }
00510 destroy_cube(cube);
00511 }
00512 cpl_table_delete(tbl_spectrum);
00513
00514
00515 if (cfg->convInd == 1) {
00516 cpl_msg_info(_id,"Clean mean");
00517 cleanfactor = clean_mean(factor,
00518 cfg->nframes,
00519 cfg->lo_reject*100.,
00520 cfg->hi_reject*100.);
00521 }
00522 if (cleanfactor > 100000. || cleanfactor == FLAG) {
00523 cpl_msg_error(_id,"could not do clean_mean!" );
00524 return -1;
00525 }
00526
00527 factorfile = fopen(cfg->convName,"w");
00528 fprintf(factorfile,"%14.12g",cleanfactor);
00529 fclose(factorfile);
00530 if (cfg->convInd == 1) {
00531 destroy_array(factor);
00532 }
00533
00534
00535 cpl_msg_info(_id,"average with rejection");
00536 cube_object = list_make_cube(list_object, cfg->nframes);
00537 averaged_spectrum = average_with_rejection(cube_object,
00538 cfg->lo_reject,
00539 cfg->hi_reject );
00540
00541
00542
00543
00544
00545 destroy_cube(cube_object);
00546
00547 if (averaged_spectrum == NULL) {
00548 cpl_msg_error(_id, "average_with_rejection failed" );
00549 destroy_list(list_object);
00550 stdstar_free (cfg);
00551 cpl_frameset_delete(raw);
00552 return -1;
00553 }
00554 destroy_list(list_object);
00555
00556 qclog_tbl = sinfoni_qclog_init(2);
00557 key_value = cpl_calloc(FILE_NAME_SZ,sizeof(char));
00558
00559 asp_imgw=cpl_image_wrap_float(averaged_spectrum->lx,averaged_spectrum->ly,
00560 averaged_spectrum->data);
00561 asp_img=cpl_image_duplicate(asp_imgw);
00562 cpl_image_unwrap(asp_imgw);
00563
00564 cpl_msg_info(_id,"lx=%d ly=%d\n",averaged_spectrum->lx,averaged_spectrum->ly);
00565
00566 sprintf(key_value,"%g",cleanfactor);
00567 sinfoni_qclog_add(qclog_tbl,0,"QC CONVFCT","CPL_TYPE_DOUBLE",key_value,"Conversion factor");
00568 sprintf(key_value,"%d",check3);
00569 sinfoni_qclog_add(qclog_tbl,1,"QC CHECK3","CPL_TYPE_INT",key_value,"Check evaluation box");
00570
00571
00572
00573 if(-1 == sinfoni_pro_save_ima(asp_img,raw,sof,cfg->outName,
00574 PRO_STD_STAR_SPECTRUM,qclog_tbl,
00575 plugin_id,config)) {
00576 cpl_msg_error(_id,"cannot dump ima %s", cfg->outName);
00577 }
00578
00579
00580 set_wcs_spectrum (asp_img, cfg->outName,cenLambda, disp, cenpix);
00581 cpl_image_delete(asp_img);
00582 cpl_table_delete(qclog_tbl);
00583 cpl_free(key_value);
00584
00585
00586
00587
00588 destroy_image(averaged_spectrum);
00589 stdstar_free (cfg);
00590 cpl_frameset_delete(raw);
00591 return 0;
00592
00593 }
00594
00595
00596
00597
00598