00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #include "cubecreate.h"
00014 #include "sinfoni_pro_save.h"
00015 #include "cubecreate_ini_by_cpl.h"
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035 void change_plist_cubecreate (cpl_propertylist * plist,
00036 char * name,
00037 float cenLambda,
00038 float dispersion,
00039 int cenpix,
00040 float center_x,
00041 float center_y )
00042 {
00043
00044 float pixelscale ;
00045 float ra ;
00046 float dec ;
00047 float angle ;
00048 float radangle ;
00049 float cd1_1, cd1_2, cd2_1, cd2_2 ;
00050 char firsttext[2*FILENAMESZ] ;
00051 char * date ;
00052
00053 strcpy(firsttext, "spiffi_cubecreate -f \0") ;
00054
00055 pixelscale = spiffi_get_pixelscale(name) /2;
00056 ra = spiffi_get_RA(name) ;
00057 dec = spiffi_get_DEC(name) ;
00058 angle = spiffi_get_posangle(name) ;
00059 radangle = angle * PI_NUMB / 180. ;
00060 cd1_1 = -cos(radangle) * pixelscale / 3600. ;
00061 cd1_2 = sin(radangle) * pixelscale / 3600. ;
00062 cd2_1 = sin(radangle) * pixelscale / 3600. ;
00063 cd2_2 = cos(radangle) * pixelscale / 3600. ;
00064
00065
00066
00067
00068 cpl_propertylist_insert_after_string(plist, "EXPTIME", "CTYPE1", "RA---TAN" ) ;
00069 cpl_propertylist_set_comment(plist, "CTYPE1", "Projected Rectascension" ) ;
00070
00071 cpl_propertylist_insert_after_float(plist, "CTYPE1", "CRPIX1", center_x ) ;
00072 cpl_propertylist_set_comment(plist, "CRPIX1","Reference pixel in RA" ) ;
00073
00074
00075
00076 cpl_propertylist_insert_after_float(plist, "CRPIX1", "CRVAL1", ra ) ;
00077 cpl_propertylist_set_comment(plist, "CRVAL1","Reference RA" ) ;
00078
00079 cpl_propertylist_insert_after_float(plist, "CRVAL1", "CDELT1", pixelscale/3600. ) ;
00080 cpl_propertylist_set_comment(plist, "CDELT1","pixel scale" ) ;
00081
00082 cpl_propertylist_insert_after_string(plist, "CDELT1", "CUNIT1", "DEGREE" ) ;
00083 cpl_propertylist_set_comment(plist, "CUNIT1","RA-UNIT" ) ;
00084
00085 cpl_propertylist_insert_after_string(plist, "CUNIT1", "CTYPE2", "DEC--TAN" ) ;
00086 cpl_propertylist_set_comment(plist, "CTYPE2", "Projected Declination") ;
00087
00088
00089
00090
00091 cpl_propertylist_insert_after_float(plist, "CTYPE2", "CRPIX2", center_y ) ;
00092 cpl_propertylist_set_comment(plist, "CRPIX2", "Reference pixel in DEC") ;
00093
00094 cpl_propertylist_insert_after_float(plist, "CRPIX2", "CRVAL2", dec) ;
00095 cpl_propertylist_set_comment(plist, "CRVAL2", "Reference DEC") ;
00096
00097 cpl_propertylist_insert_after_float(plist, "CRVAL2", "CDELT2", pixelscale/3600. ) ;
00098 cpl_propertylist_set_comment(plist, "CDELT2", "pixel scale") ;
00099
00100 cpl_propertylist_insert_after_string(plist, "CDELT2", "CUNIT2", "DEGREE" ) ;
00101 cpl_propertylist_set_comment(plist, "CUNIT2", "DEC-UNIT") ;
00102
00103 cpl_propertylist_insert_after_string(plist, "CUNIT2", "CTYPE3", "WAVE" ) ;
00104 cpl_propertylist_set_comment(plist, "CTYPE3", "wavelength axis in microns") ;
00105
00106
00107
00108
00109 cpl_propertylist_insert_after_int(plist, "CTYPE3", "CRPIX3", cenpix+1 ) ;
00110 cpl_propertylist_set_comment(plist, "CRPIX3", "Reference pixel in z") ;
00111
00112 cpl_propertylist_insert_after_float(plist, "CRPIX3", "CRVAL3", cenLambda ) ;
00113 cpl_propertylist_set_comment(plist, "CRVAL3", "central wavelength") ;
00114
00115 cpl_propertylist_insert_after_float(plist, "CRVAL3", "CDELT3", dispersion ) ;
00116 cpl_propertylist_set_comment(plist, "CDELT3", "microns per pixel") ;
00117
00118 cpl_propertylist_insert_after_string(plist, "CDELT3", "CUNIT3", "MICRON" ) ;
00119 cpl_propertylist_set_comment(plist, "CUNIT3", "spectral unit" ) ;
00120
00121
00122
00123
00124
00125
00126 cpl_propertylist_insert_after_float(plist, "CUNIT3", "CD1_1", cd1_1 ) ;
00127 cpl_propertylist_set_comment(plist, "CD1_1", "CD rotation matrix" ) ;
00128
00129 cpl_propertylist_insert_after_float(plist, "CD1_1", "CD1_2", cd1_2 ) ;
00130 cpl_propertylist_set_comment(plist, "CD1_2", "CD rotation matrix" ) ;
00131
00132 cpl_propertylist_insert_after_float(plist, "CD1_2", "CD2_1", cd2_1 ) ;
00133 cpl_propertylist_set_comment(plist, "CD2_1", "CD rotation matrix" ) ;
00134
00135 cpl_propertylist_insert_after_float(plist, "CD2_1", "CD2_2", cd2_2 ) ;
00136 cpl_propertylist_set_comment(plist, "CD2_2", "CD rotation matrix" ) ;
00137
00138
00139 date = get_datetime_iso8601() ;
00140
00141 }
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159 int cubecreate (cpl_parameterlist* config, cpl_frameset* sof)
00160 {
00161
00162 const char* _id = "cube_create";
00163 OneImage * im=NULL ;
00164 OneImage * wavemapim=NULL ;
00165 OneImage * resampledImage=NULL ;
00166 OneImage * calim=NULL ;
00167 OneImage * halospec=NULL ;
00168 OneCube * cube=NULL ;
00169 OneCube * outcube=NULL ;
00170 OneCube * outcube2=NULL ;
00171
00172
00173 int nvalid=0;
00174 int j=0;
00175 cpl_vector* valid=NULL;
00176 cpl_imagelist* tmp_ims=NULL;
00177 cpl_image* cub_avg=NULL;
00178 double qc_ima_stdev=0;
00179 cpl_image* img_slopex=NULL;
00180 cpl_image* img_slopey=NULL;
00181 char* tmp_name=NULL;
00182
00183 fits_header * head=NULL ;
00184 double dis ;
00185 float mi ;
00186 float ma ;
00187 float fcol ;
00188 double centralLambda ;
00189 int centralpix ;
00190 float ** slit_edges ;
00191 float * correct_dist ;
00192 float * distances ;
00193 float center_x ;
00194 float center_y ;
00195 int i ;
00196
00197 cpl_imagelist* cub_ims=NULL;
00198 cpl_frameset* stk=NULL;
00199 OneImage * eimg;
00200 cpl_image* cimg=NULL;
00201
00202 int* status=NULL;
00203 char* tbl_slitpos_name=NULL;
00204 cpl_table* tbl_slitpos=NULL;
00205 char* tbl_distances_name=NULL;
00206 cpl_table* tbl_distances=NULL;
00207 char* tbl_firstcol_name=NULL;
00208 cpl_table* tbl_firstcol=NULL;
00209 float edge_x=0;
00210 float edge_y=0;
00211 float tmp_float=0;
00212 cpl_propertylist * cplist=NULL;
00213 cube_config * cfg=NULL; ;
00214
00215
00216 cfg = parse_cpl_input_cube(config,sof,&stk) ;
00217 if (cfg == NULL)
00218 {
00219 cpl_msg_error(_id," could not parse cpl input!") ;
00220 return -1 ;
00221 }
00222
00223 if(is_fits_file(cfg->inFrame) != 1) {
00224 cpl_msg_error(_id,"Input file %s is not FITS",cfg->inFrame);
00225 return -1;
00226 }
00227 if(is_fits_file(cfg->wavemap) != 1) {
00228 cpl_msg_error(_id,"Input file %s is not FITS",cfg->wavemap);
00229 return -1;
00230 }
00231
00232 im = load_image(cfg->inFrame) ;
00233 if (im == NULL)
00234 {
00235 cpl_msg_error(_id, " could not load inFrame\n" ) ;
00236 return -1 ;
00237 }
00238 wavemapim = load_image(cfg->wavemap) ;
00239 if (wavemapim == NULL )
00240 {
00241 cpl_msg_error(_id, " could not load wavemap\n" ) ;
00242 return -1 ;
00243 }
00244 head = fits_read_header(cfg->inFrame) ;
00245
00246
00247 cpl_msg_info(_id,"Resampling");
00248 resampledImage =
00249 definedResampling( im, wavemapim, cfg->ncoeffs, &cfg->nrows, &dis,
00250 &mi, &ma, ¢ralLambda, ¢ralpix) ;
00251 if ( resampledImage == NULL )
00252 {
00253 cpl_msg_error(_id, " definedResampling() failed\n" ) ;
00254 return -1 ;
00255 }
00256
00257 cpl_msg_info (_id,"dispersion %f\n", dis) ;
00258 cpl_msg_info (_id,"minimal lambda %f\n", mi ) ;
00259 cpl_msg_info (_id,"maximal lambda %f\n", ma) ;
00260 cpl_msg_info (_id,"central lambda %f\n", centralLambda) ;
00261 cpl_msg_info (_id,"central pixel %d\n", centralpix) ;
00262
00263 destroy_image(im) ;
00264 destroy_image(wavemapim) ;
00265
00266 if (cfg->halocorrectInd == 1)
00267 {
00268
00269 cpl_msg_info(_id,"Doing halocorrection");
00270 halospec = load_image(cfg->halospectrum) ;
00271 calim = multiplyImageWithSpectrum( resampledImage, halospec) ;
00272 if (calim == NULL)
00273 {
00274 cpl_msg_error(_id, " calibration failed\n" ) ;
00275 return -1 ;
00276 }
00277 destroy_image(halospec) ;
00278 destroy_image(resampledImage) ;
00279 resampledImage = calim ;
00280 }
00281
00282
00283
00284
00285
00286
00287
00288 cpl_msg_info(_id,"Cubecreation");
00289 if (NULL == (correct_dist = (float*) cpl_calloc(cfg->nslits, sizeof (float))) )
00290 {
00291 cpl_msg_error(_id," could not allocate memory!") ;
00292 return -1 ;
00293 }
00294
00295 if (cfg->northsouthInd == 0)
00296 {
00297 if (NULL == (slit_edges = (float**) cpl_calloc (cfg->nslits, sizeof(float*)) ) )
00298 {
00299 cpl_msg_error(_id," could not allocate memory!") ;
00300 return -1 ;
00301 }
00302 for ( i = 0 ; i < cfg->nslits ; i++ )
00303 {
00304 if (NULL == (slit_edges[i] = (float*) cpl_calloc (2, sizeof(float))))
00305 {
00306 cpl_msg_error(_id," could not allocate memory!") ;
00307 return -1 ;
00308 }
00309 }
00310
00311
00312
00313
00314
00315 tbl_slitpos_name = (char*) cpl_calloc(MAX_NAME_SIZE,sizeof(char*));
00316 strcpy(tbl_slitpos_name,cfg->poslist);
00317 tbl_slitpos = cpl_table_load(tbl_slitpos_name,1,0);
00318
00319 if(cpl_error_get_code() != CPL_ERROR_NONE) {
00320 cpl_msg_error(_id,"error loading tbl %s",tbl_slitpos_name);
00321 cpl_msg_error(_id,(char* ) cpl_error_get_message());
00322 return -1;
00323 }
00324 for (i =0 ; i< cfg->nslits; i++){
00325
00326 edge_x=cpl_table_get_float(tbl_slitpos,"start",i,status);
00327 edge_y=cpl_table_get_float(tbl_slitpos,"end",i,status);
00328 array2D_set_value(slit_edges,edge_x,i,0);
00329 array2D_set_value(slit_edges,edge_y,i,1);
00330 }
00331 if(cpl_error_get_code() != CPL_ERROR_NONE) {
00332 cpl_msg_error(_id,"error reading tbl %s",tbl_slitpos_name);
00333 cpl_msg_error(_id,(char* ) cpl_error_get_message());
00334 return -1;
00335 }
00336 cpl_table_delete(tbl_slitpos);
00337 cpl_free(tbl_slitpos_name);
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361
00362
00363
00364 cpl_msg_info(_id,"Make the data cube");
00365 cube = makeCubeSpi( resampledImage, slit_edges, correct_dist ) ;
00366 if (cube == NULL)
00367 {
00368 cpl_msg_error(_id," could not construct data cube!") ;
00369 return -1 ;
00370 }
00371
00372 for ( i = 0 ; i < cfg->nslits ; i++ )
00373 {
00374 cpl_free (slit_edges[i]) ;
00375 }
00376 cpl_free (slit_edges) ;
00377 }
00378 else
00379 {
00380 if ( NULL == (distances = (float*) cpl_calloc (cfg->nslits - 1, sizeof (float))) )
00381 {
00382 cpl_msg_error(_id," could allocate memory!") ;
00383 return -1 ;
00384 }
00385
00386
00387
00388 tbl_distances_name = (char*) cpl_calloc(MAX_NAME_SIZE,sizeof(char*));
00389 strcpy(tbl_distances_name,cfg->poslist);
00390 tbl_distances = cpl_table_load(tbl_distances_name,1,0);
00391 if(cpl_error_get_code() != CPL_ERROR_NONE) {
00392 cpl_msg_error(_id,"error loading tbl %s",tbl_distances_name);
00393 cpl_msg_error(_id,(char* ) cpl_error_get_message());
00394 return -1;
00395 }
00396 for (i =0 ; i< cfg->nslits-1; i++){
00397 tmp_float=cpl_table_get_float(tbl_distances,"slitlet_distance",i,status);
00398 array_set_value(distances,tmp_float,i);
00399
00400 }
00401 if(cpl_error_get_code() != CPL_ERROR_NONE) {
00402 cpl_msg_error(_id,"error reading tbl %s",tbl_distances_name);
00403 cpl_msg_error(_id,(char* ) cpl_error_get_message());
00404 return -1;
00405 }
00406
00407 cpl_table_delete(tbl_distances);
00408 cpl_free(tbl_distances_name);
00409
00410
00411
00412
00413
00414
00415
00416
00417
00418
00419
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435
00436 tbl_firstcol_name=cfg->firstCol;
00437 if ((tbl_firstcol=cpl_table_load(tbl_firstcol_name,1,0)) != NULL) {
00438 fcol=cpl_table_get_float(tbl_firstcol,"first_col",0,status);
00439 cpl_table_delete(tbl_firstcol);
00440 } else {
00441 cpl_msg_error(_id,"table %s not found! Exit!",tbl_firstcol_name);
00442 return -1 ;
00443 }
00444
00445
00446 if(cpl_error_get_code() != CPL_ERROR_NONE) {
00447 cpl_msg_error(_id,"error reading tbl %s",tbl_firstcol_name);
00448 cpl_msg_error(_id,(char* ) cpl_error_get_message());
00449 return -1;
00450 }
00451
00452
00453
00454
00455
00456
00457
00458
00459
00460
00461
00462 cpl_msg_info(_id,"Build the data cube");
00463 cube = makeCubeDist( resampledImage, fcol, distances, correct_dist ) ;
00464 if (cube == NULL)
00465 {
00466 cpl_msg_error(_id," could not construct a data cube!") ;
00467 return -1 ;
00468 }
00469 cpl_free(distances) ;
00470 }
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481 cpl_msg_info(_id,"Finetuning");
00482 if (cfg->method == 'P')
00483 {
00484 outcube = fineTuneCube( cube, correct_dist, cfg->order ) ;
00485 if (outcube == NULL)
00486 {
00487 cpl_msg_error(_id," could not fine tune the data cube") ;
00488 return -1 ;
00489 }
00490 }
00491 else if (cfg->method == 'F')
00492 {
00493 outcube = fineTuneCubeByFFT( cube, correct_dist ) ;
00494 if ( outcube == NULL )
00495 {
00496 cpl_msg_error(_id," could not fine tune the data cube") ;
00497 return -1 ;
00498 }
00499 }
00500 else if (cfg->method == 'S')
00501 {
00502 outcube = fineTuneCubeBySpline( cube, correct_dist ) ;
00503 if ( outcube == NULL )
00504 {
00505 cpl_msg_error(_id," could not fine tune the data cube") ;
00506 return -1 ;
00507 }
00508 }
00509 else
00510 {
00511 cpl_msg_error(_id," wrong method indicator given!") ;
00512 return -1 ;
00513 }
00514 outcube2=binCube(outcube,1,2,2,61,0,63);
00515 center_x = outcube2->lx / 2. + 0.5 ;
00516 center_y = outcube2->ly / 2. + 0.5 ;
00517
00518
00519
00520 cpl_msg_info(_id,"copy eclipse cube to CPL imset");
00521 cub_ims=cpl_imagelist_new(outcube2->lx,outcube2->ly,outcube->np,CPL_TYPE_FLOAT);
00522 for(i=0;i<outcube2->np;i++) {
00523 eimg=outcube2->plane[i];
00524 cimg=cpl_image_wrap_float(eimg->lx,eimg->ly,eimg->data);
00525 cpl_imagelist_set(cub_ims,cimg,i);
00526 }
00527
00528
00529
00530 nvalid=outcube2->np;
00531 valid=cpl_vector_new(nvalid);
00532 cpl_vector_fill(valid,1);
00533 for(i=0;cfg->qc_chop_zmin;i++) {
00534 cpl_vector_set(valid,i,0);
00535 }
00536 for(i=0;i<cfg->qc_chop_zmax;i++) {
00537 j=nvalid-i-1;
00538 cpl_vector_set(valid,j,0);
00539 }
00540 tmp_ims=cpl_imagelist_duplicate(cub_ims);
00541 cpl_imagelist_erase(tmp_ims,valid);
00542 cub_avg=cpl_imagelist_collapse_create(tmp_ims);
00543 qc_ima_stdev=cpl_image_get_stdev(cub_avg);
00544 img_slopex=cpl_image_collapse_create(cub_avg,0);
00545 img_slopey=cpl_image_collapse_create(cub_avg,1);
00546
00547 if(-1 == sinfoni_pro_dump_ims(cub_ims,stk,sof,cfg->outName,
00548 PRO_CUBE,_id,config)) {
00549 cpl_msg_error(_id,"cannot dump cube %s", cfg->outName);
00550 fits_header_destroy(head) ;
00551 cpl_free(correct_dist) ;
00552 destroy_image(resampledImage) ;
00553 destroy_cube (cube) ;
00554 destroy_cube (outcube) ;
00555 destroy_cube (outcube2) ;
00556 cube_free(cfg);
00557
00558 return -1;
00559 }
00560
00561
00562 cplist = cpl_propertylist_new() ;
00563 if ((cpl_error_code)((cplist = cpl_propertylist_load(cfg->outName, 0)) == NULL)) {
00564 cpl_msg_error(_id, "getting header from frame %s",cfg->outName);
00565 cpl_propertylist_delete(cplist) ;
00566 return -1 ;
00567 }
00568
00569 change_plist_cubecreate(cplist, cfg->inFrame,
00570 centralLambda, dis, centralpix, center_x, center_y) ;
00571
00572 if (cpl_imagelist_save(cub_ims, cfg->outName, CPL_BPP_DEFAULT, cplist,CPL_IO_DEFAULT)!=CPL_ERROR_NONE) {
00573 cpl_msg_error(_id, "Cannot save the product %s",cfg->outName);
00574 cpl_propertylist_delete(cplist) ;
00575 return -1 ;
00576 }
00577
00578 cpl_propertylist_delete(cplist) ;
00579
00580
00581 tmp_name = (char*) cpl_calloc(FILE_NAME_SZ, sizeof (char));
00582 strcpy(tmp_name,"avg_");
00583 strcat(tmp_name,cfg->outName);
00584 if(-1 == sinfoni_pro_dump_ima(cub_avg,stk,sof,tmp_name,
00585 PRO_CUBE_COLL,_id,config)) {
00586 cpl_msg_error(_id,"cannot dump ima %s", tmp_name);
00587 cpl_free(tmp_name);
00588 fits_header_destroy(head) ;
00589 cpl_free(correct_dist) ;
00590 destroy_image(resampledImage) ;
00591 destroy_cube (cube) ;
00592 destroy_cube (outcube) ;
00593 destroy_cube (outcube2) ;
00594 cube_free(cfg);
00595
00596 return -1;
00597 }
00598
00599 strcpy(tmp_name,"slopex_avg_");
00600 strcat(tmp_name,cfg->outName);
00601 if(-1 == sinfoni_pro_dump_ima(img_slopex,stk,sof,tmp_name,
00602 PRO_SLOPEX,_id,config)) {
00603 cpl_msg_error(_id,"cannot dump ima %s", tmp_name);
00604 cpl_free(tmp_name);
00605 fits_header_destroy(head) ;
00606 cpl_free(correct_dist) ;
00607 destroy_image(resampledImage) ;
00608 destroy_cube (cube) ;
00609 destroy_cube (outcube) ;
00610 destroy_cube (outcube2) ;
00611 cube_free(cfg);
00612
00613 return -1;
00614 }
00615
00616
00617 strcpy(tmp_name,"slopey_avg_");
00618 strcat(tmp_name,cfg->outName);
00619 if(-1 == sinfoni_pro_dump_ima(img_slopey,stk,sof,tmp_name,
00620 PRO_SLOPEY,_id,config)) {
00621 cpl_msg_error(_id,"cannot dump ima %s", tmp_name);
00622 cpl_free(tmp_name);
00623 fits_header_destroy(head) ;
00624 cpl_free(correct_dist) ;
00625 destroy_image(resampledImage) ;
00626 destroy_cube (cube) ;
00627 destroy_cube (outcube) ;
00628 destroy_cube (outcube2) ;
00629 cube_free(cfg);
00630
00631 return -1;
00632 }
00633
00634 cpl_free(tmp_name);
00635
00636
00637 fits_header_destroy(head) ;
00638 cpl_free(correct_dist) ;
00639 destroy_image(resampledImage) ;
00640 destroy_cube (cube) ;
00641 destroy_cube (outcube) ;
00642 destroy_cube (outcube2) ;
00643 cube_free(cfg);
00644
00645 return 0;
00646 }
00647
00648