00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014 #include "objnod.h"
00015 #include "sinfoni_pro_save.h"
00016 #include "objnod_ini_by_cpl.h"
00017 #include "sinfoni_functions.h"
00018 #include "utilities_scired.h"
00019 #include <sinfoni_memory.h>
00020
00021
00022
00023 #define PI_NUMB (3.1415926535897932384626433832795)
00024
00025
00026
00027
00028
00029
00030
00031 int
00032 check_input_data(object_config* cfg)
00033 {
00034
00035 const char* _id = "check_input_data";
00036
00037 if (cfg == NULL)
00038 {
00039 cpl_msg_error (_id," could not parse cpl input!\n") ;
00040 return -1 ;
00041 }
00042
00043
00044 if(is_fits_file(cfg->wavemap) != 1) {
00045 cpl_msg_error(_id,"Input file wavemap %s is not FITS",cfg->wavemap);
00046 return -1;
00047 }
00048
00049
00050 if (cfg->halocorrectInd == 1)
00051 {
00052 if(is_fits_file(cfg->halospectrum) != 1) {
00053 cpl_msg_error(_id,"Input file %s is not FITS",cfg->halospectrum);
00054 return -1;
00055 }
00056
00057 }
00058
00059 if (cfg->northsouthInd == 0) {
00060 if (is_fits_file(cfg->poslist) != 1)
00061 {
00062 cpl_msg_error(_id,"File %s with tag %s is not FITS!",cfg->poslist,PRO_SLIT_POS);
00063 return -1 ;
00064 }
00065 } else {
00066
00067 if (is_fits_file(cfg->distlist) != 1)
00068 {
00069 cpl_msg_error(_id,"File %s with tag %s is not FITS!",cfg->distlist,PRO_SLITLETS_DISTANCE);
00070 return -1;
00071 }
00072 }
00073
00074
00075 return 0;
00076
00077
00078 }
00079
00080
00081 int
00082 auto_size_cube(object_config * cfg)
00083 {
00084
00085 const char* _id = "auto_size_cube";
00086 char* name =NULL;
00087 int n=0;
00088 float offx=0;
00089 float offy=0;
00090 float ref_offx=0;
00091 float ref_offy=0;
00092 float min_offx=0;
00093 float max_offx=0;
00094 float min_offy=0;
00095 float max_offy=0;
00096 float extra_offx=0;
00097 float extra_offy=0;
00098
00099 cpl_msg_info (_id,"Automatic computation of output cube size") ;
00100 for ( n = 0 ; n < cfg->nframes ; n++ )
00101 {
00102 name = cfg->framelist[n] ;
00103 if ( n == 0 )
00104 {
00105
00106
00107 ref_offx = spiffi_get_cumoffsetx(name) ;
00108 if (ref_offx == FLAG)
00109 {
00110 cpl_msg_warning(_id," could not read fits header keyword cummoffsetx! Set it to 0") ;
00111 ref_offx = 0;
00112
00113 }
00114
00115 ref_offy = spiffi_get_cumoffsety(name) ;
00116 if (ref_offy == FLAG )
00117 {
00118 cpl_msg_warning(_id," could not read fits header keyword cumoffsety! Set it to 0") ;
00119 ref_offy = 0;
00120
00121 }
00122
00123
00124 offx = 0. ;
00125 offy = 0. ;
00126 min_offx=0;
00127 max_offx=0;
00128 min_offy=0;
00129 max_offy=0;
00130
00131
00132 }
00133 else
00134 {
00135
00136 offx = spiffi_get_cumoffsetx(name) - ref_offx ;
00137 if (offx == FLAG)
00138 {
00139 cpl_msg_warning(_id," could not read fits header keyword cummoffsetx!") ;
00140 cpl_msg_warning(_id," set it to - %f",ref_offx);
00141 offx = - ref_offx;
00142
00143 }
00144
00145 offy = spiffi_get_cumoffsety(name) - ref_offy ;
00146 if (offy == FLAG )
00147 {
00148 cpl_msg_warning(_id," could not read fits header keyword! cumoffsety") ;
00149 cpl_msg_warning(_id," set it to - %f",ref_offy);
00150 offy = - ref_offy;
00151
00152 }
00153
00154 if(fabs(offx) > extra_offx) extra_offx=fabs(offx);
00155 if(fabs(offy) > extra_offy) extra_offy=fabs(offy);
00156
00157 if(offx < min_offx) min_offx=offx;
00158 if(offy < min_offy) min_offy=offy;
00159 if(offx > max_offx) max_offx=offx;
00160 if(offy > max_offy) max_offy=offy;
00161
00162
00163 }
00164
00165
00166 }
00167
00168
00169
00170
00171
00172 if(cfg->size_x == 0) cfg->size_x=SLITLENGTH+4*floor(extra_offx+0.5);
00173 if(cfg->size_y == 0) cfg->size_y=SLITLENGTH+4*floor(extra_offy+0.5);
00174 cpl_msg_info(_id,"Output cube size: %d x %d",cfg->size_x,cfg->size_y);
00175 return 0;
00176
00177
00178 }
00179
00180
00181
00182
00183 int
00184 auto_size_cube2(object_config * cfg, float* ref_offx, float* ref_offy)
00185 {
00186
00187 const char* _id = "auto_size_cube";
00188 char* name =NULL;
00189 int n=0;
00190 float offx=0;
00191 float offy=0;
00192 float min_offx=0;
00193 float max_offx=0;
00194 float min_offy=0;
00195 float max_offy=0;
00196
00197 cpl_msg_info (_id,"Automatic computation of output cube size") ;
00198 for ( n = 0 ; n < cfg->nframes ; n++ )
00199 {
00200 name = cfg->framelist[n] ;
00201
00202 offx = spiffi_get_cumoffsetx(name);
00203 if (offx == FLAG)
00204 {
00205 cpl_msg_warning(_id," could not read fits header keyword cummoffsetx!") ;
00206 cpl_msg_warning(_id," set it to 0");
00207 offx = 0;
00208
00209 }
00210
00211 offy = spiffi_get_cumoffsety(name);
00212 if (offy == FLAG )
00213 {
00214 cpl_msg_warning(_id," could not read fits header keyword! cumoffsety") ;
00215 cpl_msg_warning(_id," set it to 0");
00216 offy = 0;
00217
00218 }
00219 if(n==0) {
00220 min_offx=offx;
00221 min_offy=offy;
00222 max_offx=offx;
00223 max_offy=offy;
00224 } else {
00225 if(offx > max_offx) max_offx=offx;
00226 if(offy > max_offy) max_offy=offy;
00227 if(offx < min_offx) min_offx=offx;
00228 if(offy < min_offy) min_offy=offy;
00229 }
00230 }
00231 *ref_offx=(min_offx+max_offx)/2;
00232 *ref_offy=(min_offy+max_offy)/2;
00233
00234 if(cfg->size_x == 0) cfg->size_x=2*floor(max_offx-min_offx+0.5)+64;
00235 if(cfg->size_y == 0) cfg->size_y=2*floor(max_offy-min_offy+0.5)+64;
00236 cpl_msg_info(_id,"Output cube size: %d x %d",cfg->size_x,cfg->size_y);
00237 cpl_msg_info(_id,"Ref offset. x: %f y: %f",*ref_offx,*ref_offy);
00238 cpl_msg_info(_id,"Max offset. x: %f y: %f",max_offx,max_offy);
00239 cpl_msg_info(_id,"Min offset. x: %f y: %f",min_offx,min_offy);
00240 return 0;
00241
00242
00243 }
00244
00245
00246 float*
00247 read_distances(object_config* cfg)
00248 {
00249 const char * _id = "read_distances";
00250 int i=0;
00251 int* status=NULL;
00252 float * distances = NULL;
00253 float tmp_float=0;
00254 char tbl_distances_name[FILE_NAME_SZ];
00255 cpl_table* tbl_distances = NULL;
00256
00257 cpl_msg_info(_id,"Read distances");
00258 distances = (float*) cpl_calloc (cfg->nslits - 1, sizeof (float));
00259
00260 if ( NULL == distances )
00261 {
00262 cpl_msg_error (_id,"could allocate memory!") ;
00263 return NULL ;
00264 }
00265
00266
00267 if(cpl_error_get_code() != CPL_ERROR_NONE) {
00268 cpl_msg_error(_id,"Before loading input table");
00269 cpl_msg_error(_id,(char* ) cpl_error_get_message());
00270 return NULL;
00271 }
00272 strcpy(tbl_distances_name,cfg->distlist);
00273 tbl_distances = cpl_table_load(tbl_distances_name,1,0);
00274 if(cpl_error_get_code() != CPL_ERROR_NONE) {
00275 cpl_msg_error(_id,"loading input table %s",tbl_distances_name);
00276 cpl_msg_error(_id,(char* ) cpl_error_get_message());
00277 return NULL;
00278 }
00279
00280 for (i =0 ; i< cfg->nslits-1; i++){
00281 tmp_float=cpl_table_get_float(tbl_distances,"slitlet_distance",i,status);
00282 if(cpl_error_get_code() != CPL_ERROR_NONE) {
00283 cpl_msg_error(_id,"reading col %s from table %s","slitlet_distance",tbl_distances_name);
00284 cpl_msg_error(_id,(char* ) cpl_error_get_message());
00285 return NULL;
00286 }
00287 array_set_value(distances,tmp_float,i);
00288 }
00289 cpl_table_delete(tbl_distances);
00290 return distances;
00291
00292 }
00293
00294
00295
00296
00297 float**
00298 read_slitlets_edges(object_config * cfg)
00299 {
00300
00301 const char * _id="read_slitlets_edge";
00302
00303 char tbl_slitpos_name[FILE_NAME_SZ];
00304 cpl_table* tbl_slitpos=NULL;
00305 int n=0;
00306 int i=0;
00307 int* status=NULL;
00308 float edge_x=0;
00309 float edge_y=0;
00310 float ** slit_edges = NULL;
00311
00312 slit_edges = new_2Dfloat_array(cfg->nslits, 2) ;
00313
00314 strcpy(tbl_slitpos_name,cfg->poslist);
00315 tbl_slitpos = cpl_table_load(tbl_slitpos_name,1,0);
00316 if(cpl_error_get_code() != CPL_ERROR_NONE) {
00317 cpl_msg_error(_id,"error loading tbl %s",tbl_slitpos_name);
00318 cpl_msg_error(_id,(char* ) cpl_error_get_message());
00319 return NULL;
00320 }
00321 n = cpl_table_get_nrow(tbl_slitpos);
00322 if (n != cfg->nslits) {
00323 cpl_msg_error (_id,"No of slitlets in table is n = %d != %d !",n,cfg->nslits);
00324 return NULL;
00325 }
00326
00327 for (i =0 ; i< cfg->nslits; i++){
00328 edge_x=cpl_table_get_double(tbl_slitpos,"pos1",i,status);
00329 edge_y=cpl_table_get_double(tbl_slitpos,"pos2",i,status);
00330 if(cpl_error_get_code() != CPL_ERROR_NONE) {
00331 cpl_msg_error(_id,"error reading tbl %s row %d",tbl_slitpos_name,i);
00332 cpl_msg_error(_id,(char* ) cpl_error_get_message());
00333 return NULL;
00334 }
00335 array2D_set_value(slit_edges,edge_x,i,0);
00336 array2D_set_value(slit_edges,edge_y,i,1);
00337 }
00338 cpl_table_delete(tbl_slitpos);
00339 if(cpl_error_get_code() != CPL_ERROR_NONE) {
00340 cpl_msg_error(_id,"error reading tbl %s",tbl_slitpos_name);
00341 cpl_msg_error(_id,(char* ) cpl_error_get_message());
00342 return NULL;
00343 }
00344
00345 return slit_edges;
00346
00347 }
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358
00359
00360
00361 int objnod (const char* plugin_id,cpl_parameterlist* config,
00362 cpl_frameset* sof, const char* procatg)
00363 {
00364
00365 const char* _id = "objnod";
00366 object_config * cfg=NULL ;
00367 OneImage * im=NULL ;
00368 OneImage * wavemapim=NULL ;
00369 OneImage * resampledImage=NULL ;
00370 OneImage * calim=NULL ;
00371 OneImage * halospec=NULL ;
00372 OneImage * sky_im=NULL;
00373 OneImage* res_flat=NULL;
00374 OneImage* res_sky=NULL;
00375 OneImage* flat_im=NULL;
00376 OneImage* jitter_image=NULL;
00377 OneImage* eima_avg=NULL;
00378 OneImage* eima_med=NULL;
00379
00380
00381 OneCube * cube=NULL ;
00382 OneCube * outcube=NULL ;
00383 OneCube * outcube2=NULL ;
00384 OneCube ** cubeobject ;
00385 OneCube ** cube_tmp ;
00386 OneCube * jittercube=NULL ;
00387
00388 OneCube * maskcube=NULL ;
00389 OneCube* cflat=NULL;
00390 OneCube* cflat2=NULL;
00391 OneCube* csky=NULL;
00392 OneCube* csky2=NULL;
00393
00394
00395 int i, n ;
00396 int partind = 0 ;
00397 int centralpix ;
00398 int z_siz=0;
00399 int z_min=0;
00400 int z_max=0;
00401 int z=0;
00402 int z_stp=100;
00403 int* status=NULL;
00404 int scales_sky=0;
00405 int ks_clip=0;
00406 double kappa=2.0;
00407
00408 float ref_offx=0;
00409 float ref_offy=0;
00410 float mi ;
00411 float ma ;
00412 float fcol ;
00413 float center_x, newcenter_x ;
00414 float center_y, newcenter_y ;
00415 float cd1_1, cd1_2, cd2_1, cd2_2 ;
00416 float pixelscale, angle, radangle ;
00417 float exptime ;
00418
00419 float * correct_dist ;
00420 float * distances ;
00421 float * times ;
00422 float * offsetx, * offsety ;
00423 float ** slit_edges ;
00424
00425 double dis ;
00426 double centralLambda ;
00427
00428 char name_jitter[FILENAMESZ] ;
00429 char pro_mjit[MAX_NAME_SIZE];
00430 char pro_obs[MAX_NAME_SIZE];
00431 char pro_med[MAX_NAME_SIZE];
00432
00433
00434 char * name ;
00435 char * partname, * partname2 ;
00436 char tbl_firstcol_name[FILE_NAME_SZ];
00437 char file_name[FILE_NAME_SZ];
00438 int vllx=0;
00439 int vlly=0;
00440 int vurx=0;
00441 int vury=0;
00442
00443
00444 cpl_table* tbl_firstcol=NULL;
00445
00446
00447 cpl_imagelist* cub_ims=NULL;
00448 cpl_imagelist* msk_ims=NULL;
00449 cpl_imagelist* jit_ims=NULL;
00450
00451
00452 cpl_image* img=NULL;
00453 cpl_image* wimg=NULL;
00454 cpl_image* ima_med=NULL;
00455 cpl_image* ima_avg=NULL;
00456
00457 cpl_frameset* stk=NULL;
00458 cpl_parameter* p=NULL;
00459
00460
00461
00462 if (strcmp(procatg,PRO_COADD_STD) == 0) {
00463 strcpy(pro_mjit,PRO_MASK_COADD_STD);
00464 strcpy(pro_obs,PRO_OBS_STD);
00465 strcpy(pro_med,PRO_MED_COADD_STD);
00466
00467 } else if (strcmp(procatg,PRO_COADD_PSF) == 0) {
00468 strcpy(pro_mjit,PRO_MASK_COADD_PSF);
00469 strcpy(pro_obs,PRO_OBS_PSF);
00470 strcpy(pro_med,PRO_MED_COADD_PSF);
00471 } else {
00472 strcpy(pro_mjit,PRO_MASK_COADD_OBJ);
00473 strcpy(pro_obs,PRO_OBS_OBJ);
00474 strcpy(pro_med,PRO_MED_COADD_OBJ);
00475 }
00476
00477
00478
00479 stk = cpl_frameset_new();
00480 cfg = parse_cpl_input_objnod(config,sof,&stk) ;
00481
00482
00483
00484 if(-1 == check_input_data(cfg)) {
00485 cpl_msg_error(_id,"error checking input");
00486 cpl_frameset_delete(stk);
00487 return -1;
00488 }
00489
00490 if ( cfg->jitterind == 1 )
00491 {
00492 times = (float*) cpl_calloc (cfg->nframes, sizeof (float)) ;
00493 if ( times == NULL )
00494 {
00495 cpl_msg_error(_id," could not allocate memory!") ;
00496 objnod_free(cfg);
00497 cpl_frameset_delete(stk);
00498 return -1 ;
00499 }
00500 offsetx = (float*) cpl_calloc (cfg->nframes, sizeof(float)) ;
00501 if ( offsetx == NULL )
00502 {
00503 cpl_msg_error (_id," could not allocate memory!") ;
00504 cpl_free(times);
00505 objnod_free(cfg);
00506 cpl_frameset_delete(stk);
00507 return -1 ;
00508 }
00509 offsety = (float*) cpl_calloc (cfg->nframes, sizeof(float)) ;
00510 if ( offsety == NULL )
00511 {
00512 cpl_msg_error (_id," could not allocate memory!") ;
00513 cpl_free(times);
00514 cpl_free(offsetx);
00515 objnod_free(cfg);
00516 cpl_frameset_delete(stk);
00517 return -1 ;
00518 }
00519 }
00520
00521 if (cfg->jitterind == 0)
00522 {
00523 if ( NULL != (partname = strtok(cfg->outName, ".")))
00524 {
00525 partname2 = strtok (NULL, ".") ;
00526 partind = 1 ;
00527 }
00528 }
00529
00530
00531 if( -1 == auto_size_cube2(cfg,&ref_offx,&ref_offy) ) {
00532 cpl_msg_error(_id,"Error resizing cube");
00533 cpl_free(times);
00534 cpl_free(offsetx);
00535 cpl_free(offsety);
00536 objnod_free(cfg);
00537 cpl_frameset_delete(stk);
00538 return -1;
00539 }
00540
00541 if(cfg->northsouthInd !=0 ) {
00542 cpl_msg_info(_id,"Read table %s",PRO_FIRST_COL);
00543 strcpy(tbl_firstcol_name,cfg->firstCol);
00544 tbl_firstcol = cpl_table_load(tbl_firstcol_name,1,0);
00545 if(cpl_error_get_code() != CPL_ERROR_NONE) {
00546 cpl_msg_error(_id,"loading input table %s",tbl_firstcol_name);
00547 cpl_msg_error(_id,(char* ) cpl_error_get_message());
00548 return -1;
00549 }
00550
00551 fcol=cpl_table_get_float(tbl_firstcol,"first_col",0,status);
00552 cpl_table_delete(tbl_firstcol);
00553
00554 }
00555
00556
00557
00558
00559 for ( n = 0 ; n < cfg->nframes ; n++ )
00560 {
00561
00562 cpl_msg_info(_id,"Read FITS information");
00563 name = cfg->framelist[n] ;
00564 if (n == 0)
00565 {
00566 strcpy (name_jitter, name) ;
00567 }
00568 if( is_fits_file(name) != 1) {
00569 cpl_msg_error(_id,"Input file %s is not FITS",name);
00570 cpl_free(times);
00571 cpl_free(offsetx);
00572 cpl_free(offsety);
00573 objnod_free(cfg);
00574 cpl_frameset_delete(stk);
00575 return -1;
00576 }
00577
00578
00579
00580
00581
00582 pixelscale = spiffi_get_pixelscale(name) /2;
00583 angle = spiffi_get_posangle(name) ;
00584 radangle = angle * PI_NUMB / 180. ;
00585 cd1_1 = cos(radangle) ;
00586 cd1_2 = sin(radangle) ;
00587 cd2_1 = -sin(radangle) ;
00588 cd2_2 = cos(radangle) ;
00589
00590 cpl_msg_info(_id,"frame no.: %d, name: %s\n", n, name) ;
00591 im = load_image(name) ;
00592
00593 if (im == NullImage)
00594 {
00595 cpl_msg_error(_id," could not load frame %s!",name) ;
00596 return -1 ;
00597 }
00598
00599 if (cfg->jitterind == 1)
00600 {
00601 exptime = spiffi_get_ditndit(name) ;
00602
00603 if (exptime == FLAG)
00604 {
00605 cpl_msg_error(_id,"could not read fits header keyword exptime!");
00606 return -1 ;
00607 }
00608 times[n] = exptime ;
00609
00610 if(-1==assign_offset(n,name,offsetx,offsety,ref_offx,ref_offy)) {
00611 cpl_msg_error(_id,"Error assigning offsets");
00612 return -1;
00613 }
00614
00615 }
00616
00617
00618
00619
00620
00621
00622
00623
00624 cpl_msg_info(_id,"Resampling object");
00625 wavemapim = load_image(cfg->wavemap);
00626
00627 if(wavemapim == NULL) {
00628 cpl_msg_error(_id,"could not load wavemap\n");
00629 return -1;
00630 }
00631 resampledImage = definedResampling( im,
00632 wavemapim,
00633 cfg->ncoeffs,
00634 &cfg->nrows,
00635 &dis,
00636 &mi,
00637 &ma,
00638 ¢ralLambda,
00639 ¢ralpix) ;
00640
00641
00642 if ( resampledImage == NULL )
00643 {
00644 cpl_msg_error ( _id," definedResampling() failed\n" ) ;
00645
00646
00647 destroy_image(wavemapim);
00648 destroy_image(im);
00649 cpl_free(times);
00650 cpl_free(offsetx);
00651 cpl_free(offsety);
00652 objnod_free(cfg);
00653 cpl_frameset_delete(stk);
00654
00655 return -1 ;
00656 }
00657
00658
00659 if(n ==0) {
00660 if(strcmp(cfg->mflat_dist,"not_found") != 0) {
00661 cpl_msg_info(_id,"Resampling master flat");
00662 flat_im = load_image(cfg->mflat_dist);
00663 if( flat_im == NULL) {
00664 cpl_msg_error(_id,"Distorted master flat field not found");
00665 cpl_msg_error(_id,"You may have set --stack-flat_ind=FALSE");
00666 cpl_msg_error(_id,"Flat field resampling skipped");
00667 destroy_image(wavemapim);
00668 destroy_image(im);
00669 objnod_free(cfg);
00670 cpl_frameset_delete(stk);
00671 return -1;
00672 }
00673
00674 res_flat = definedResampling(flat_im,
00675 wavemapim,
00676 cfg->ncoeffs,
00677 &cfg->nrows,
00678 &dis,
00679 &mi,
00680 &ma,
00681 ¢ralLambda,
00682 ¢ralpix) ;
00683
00684
00685
00686 if ( res_flat == NULL )
00687 {
00688 cpl_msg_error ( _id," definedResampling() failed\n" ) ;
00689
00690
00691 destroy_image(resampledImage);
00692 destroy_image(wavemapim);
00693 destroy_image(im);
00694 cpl_free(times);
00695 cpl_free(offsetx);
00696 cpl_free(offsety);
00697 objnod_free(cfg);
00698 cpl_frameset_delete(stk);
00699 destroy_image(flat_im) ;
00700
00701 return -1 ;
00702 }
00703 destroy_image(flat_im) ;
00704 }
00705
00706
00707 if(strcmp(cfg->sky_dist,"no_sky")!=0) {
00708 cpl_msg_info(_id,"Resampling sky");
00709 sky_im = load_image(cfg->sky_dist);
00710 res_sky = definedResampling(sky_im,
00711 wavemapim,
00712 cfg->ncoeffs,
00713 &cfg->nrows,
00714 &dis,
00715 &mi,
00716 &ma,
00717 ¢ralLambda,
00718 ¢ralpix) ;
00719
00720
00721
00722
00723 if ( res_sky == NULL ) {
00724 cpl_msg_error ( _id," definedResampling() failed\n" ) ;
00725
00726
00727 destroy_image(resampledImage);
00728 destroy_image(wavemapim);
00729 destroy_image(im);
00730 cpl_free(times);
00731 cpl_free(offsetx);
00732 cpl_free(offsety);
00733 objnod_free(cfg);
00734 cpl_frameset_delete(stk);
00735 destroy_image(flat_im) ;
00736 destroy_image(sky_im) ;
00737
00738 return -1 ;
00739 }
00740 destroy_image(sky_im) ;
00741
00742 }
00743 }
00744
00745
00746 cpl_msg_info (_id,"dispersion %f\n", dis) ;
00747 cpl_msg_info (_id,"lambda min %f max %f cent %f\n", mi,ma,centralLambda ) ;
00748 cpl_msg_info (_id,"central pixel %d\n", centralpix) ;
00749
00750 destroy_image(im) ;
00751 destroy_image(wavemapim) ;
00752
00753
00754
00755
00756
00757
00758
00759 if (cfg->halocorrectInd == 1)
00760 {
00761 cpl_msg_info(_id,"Calibration");
00762 halospec = load_image(cfg->halospectrum) ;
00763 calim = multiplyImageWithSpectrum( resampledImage, halospec) ;
00764 if (calim == NULL)
00765 {
00766 cpl_msg_error ( _id," multiplyImageWithSpectrum() failed\n" ) ;
00767 return -1 ;
00768 }
00769 destroy_image(halospec) ;
00770 destroy_image(resampledImage) ;
00771 resampledImage = copy_image(calim) ;
00772 destroy_image(calim);
00773 }
00774
00775
00776
00777
00778
00779
00780
00781
00782 cpl_msg_info(_id,"Cube creation");
00783
00784 if (cfg->northsouthInd == 0) {
00785 cpl_msg_info(_id,"cfg->northsouthInd == 0");
00786 if( NULL == (slit_edges = read_slitlets_edges(cfg)) ) {
00787 cpl_msg_error(_id,"error reading slitlets edges");
00788 return -1;
00789 }
00790 } else {
00791 cpl_msg_info(_id,"cfg->northsouthInd != 0");
00792
00793 if( NULL == (distances = read_distances(cfg))) {
00794 cpl_msg_error(_id,"error reading distances");
00795 destroy_image(res_flat);
00796 destroy_image(res_sky);
00797 destroy_image(resampledImage);
00798 cpl_free(times);
00799 cpl_free(offsetx);
00800 cpl_free(offsety);
00801 objnod_free(cfg);
00802 cpl_frameset_delete(stk);
00803 return -1;
00804 }
00805
00806 }
00807
00808 correct_dist = (float*) cpl_calloc(cfg->nslits, sizeof (float));
00809 if (NULL == correct_dist )
00810 {
00811 cpl_msg_error(_id," could not allocate memory!\n") ;
00812
00813
00814 cpl_free(distances);
00815 destroy_image(res_flat);
00816 destroy_image(resampledImage);
00817 cpl_free(times);
00818 cpl_free(offsetx);
00819 cpl_free(offsety);
00820 objnod_free(cfg);
00821 cpl_frameset_delete(stk);
00822
00823 return -1 ;
00824 }
00825
00826
00827
00828 cpl_msg_info(_id,"Create cube object");
00829 if (cfg->northsouthInd ==0 )
00830 {
00831
00832 cube = makeCubeSpi( resampledImage, slit_edges, correct_dist ) ;
00833
00834 if (cube == NULL)
00835 {
00836 cpl_msg_error (_id," could not construct data cube!\n") ;
00837
00838
00839 cpl_free(correct_dist);
00840 cpl_free(distances);
00841 destroy_image(res_flat);
00842 destroy_image(resampledImage);
00843 cpl_free(times);
00844 cpl_free(offsetx);
00845 cpl_free(offsety);
00846 objnod_free(cfg);
00847 cpl_frameset_delete(stk);
00848
00849 return -1 ;
00850 }
00851 } else {
00852 cube = makeCubeDist(resampledImage,fcol,distances,correct_dist);
00853
00854 if (cube == NULL)
00855 {
00856 cpl_msg_error (_id," could not construct a data cube!\n") ;
00857
00858
00859 cpl_free(correct_dist);
00860 cpl_free(distances);
00861 destroy_image(res_flat);
00862 destroy_image(resampledImage);
00863 cpl_free(times);
00864 cpl_free(offsetx);
00865 cpl_free(offsety);
00866 objnod_free(cfg);
00867 cpl_frameset_delete(stk);
00868
00869 return -1 ;
00870 }
00871
00872 }
00873 destroy_image(resampledImage);
00874
00875
00876 if(n==0) {
00877 if(strcmp(cfg->mflat_dist,"not_found")!=0) {
00878 cpl_msg_info(_id,"Create cube master flat");
00879 if (cfg->northsouthInd ==0 ) {
00880 cflat = makeCubeSpi( res_flat, slit_edges, correct_dist ) ;
00881 if (cflat == NULL) {
00882 cpl_msg_error (_id," could not construct data cube!\n") ;
00883 destroy_image(res_sky);
00884 destroy_image(res_flat);
00885 cpl_free(correct_dist);
00886 destroy_cube(cube);
00887 destroy_2Darray(slit_edges,cfg->nslits);
00888 cpl_free(times);
00889 cpl_free(offsetx);
00890 cpl_free(offsety);
00891 cpl_frameset_delete(stk);
00892 objnod_free(cfg);
00893
00894 return -1 ;
00895 }
00896 } else {
00897 cflat = makeCubeDist( res_flat, fcol, distances, correct_dist ) ;
00898
00899 if (cflat == NULL) {
00900 cpl_msg_error (_id," could not construct a data cube!\n") ;
00901 return -1 ;
00902 }
00903 }
00904
00905 destroy_image(res_flat);
00906 }
00907 if(strcmp(cfg->sky_dist,"no_sky")!=0) {
00908
00909 cpl_msg_info(_id,"Create cube sky");
00910 if (cfg->northsouthInd ==0 ) {
00911 csky = makeCubeSpi( res_sky, slit_edges, correct_dist ) ;
00912 if (csky == NULL) {
00913 cpl_msg_error (_id," could not construct data cube!\n") ;
00914 return -1 ;
00915 }
00916 } else {
00917 csky = makeCubeDist( res_sky, fcol, distances, correct_dist ) ;
00918 if (csky == NULL) {
00919 cpl_msg_error (_id," could not construct a data cube!\n") ;
00920 return -1 ;
00921 }
00922 }
00923 destroy_image(res_sky);
00924 }
00925 }
00926
00927
00928 if (cfg->northsouthInd ==0 )
00929 {
00930 destroy_2Darray(slit_edges,cfg->nslits);
00931 }
00932 else
00933 {
00934 destroy_array(distances);
00935 }
00936
00937
00938
00939
00940
00941
00942
00943
00944
00945
00946
00947
00948
00949 if(n==0) {
00950
00951
00952 if(strcmp(cfg->sky_dist,"no_sky")!=0) {
00953
00954 if(NULL == (csky2=fine_tune(csky,correct_dist,cfg->method,cfg->order,cfg->nslits))){
00955 cpl_msg_error (_id," could not fine tune the data cube\n") ;
00956 return -1 ;
00957 }
00958 destroy_cube(csky);
00959 cpl_msg_info(_id,"Stretch output cube along Y direction");
00960
00961 csky = binCube(csky2,1,2,0,63,0,63);
00962 destroy_cube(csky2);
00963
00964
00965 cub_ims=sinfoni_cube2imglist(csky);
00966 if(-1 == sinfoni_pro_dump_ims(cub_ims,sof,sof,"out_sky_cube.fits",
00967 "OBS_SKY",plugin_id,config))
00968 {
00969 cpl_msg_error(_id,"cannot dump cube %s", "out_sky_cube.fits");
00970 cpl_free(correct_dist);
00971 destroy_cube(csky) ;
00972 destroy_cube(cflat) ;
00973 destroy_cube(cube);
00974
00975 cpl_imagelist_delete(cub_ims);
00976 cpl_free(times);
00977 cpl_free(offsetx);
00978 cpl_free(offsety);
00979 objnod_free(cfg);
00980 cpl_frameset_delete(stk);
00981 return -1;
00982 }
00983 eima_med=medianCube(csky);
00984 center_x = eima_med->lx / 2. + 0.5 ;
00985 center_y = eima_med->ly / 2. + 0.5 ;
00986
00987 set_wcs_cube(cub_ims, "out_sky_cube.fits", centralLambda,
00988 dis, centralpix, center_x, center_y);
00989
00990 cpl_imagelist_delete(cub_ims);
00991 destroy_cube(csky) ;
00992
00993
00994 ima_med=cpl_image_wrap_float(eima_med->lx,eima_med->ly,eima_med->data);
00995 if (ima_med == NULL) {
00996 cpl_msg_error(_id,"Creating an average image");
00997 }
00998 if(-1 == sinfoni_pro_dump_ima(ima_med,sof,sof,"out_sky_med.fits",
00999 "SKY_MED",plugin_id,config))
01000 {
01001 cpl_msg_error(_id,"cannot dump ima %s", "out_sky_med.fits");
01002 cpl_free(correct_dist);
01003 cpl_image_unwrap(ima_med);
01004 destroy_image(eima_med);
01005 destroy_cube(cflat) ;
01006 destroy_cube(cube);
01007 cpl_free(times);
01008 cpl_free(offsetx);
01009 cpl_free(offsety);
01010 objnod_free(cfg);
01011 cpl_frameset_delete(stk);
01012
01013 return -1;
01014 }
01015 set_wcs_image(ima_med, "out_sky_med.fits", center_x, center_y);
01016 cpl_image_unwrap(ima_med);
01017 destroy_image(eima_med);
01018 }
01019
01020
01021
01022 if(strcmp(cfg->mflat_dist,"not_found")!=0) {
01023
01024 if(NULL == (cflat2=fine_tune(cflat,correct_dist,cfg->method,cfg->order,cfg->nslits))){
01025 cpl_msg_error (_id," could not fine tune the data cube\n") ;
01026 return -1 ;
01027 }
01028 destroy_cube(cflat);
01029 cpl_msg_info(_id,"Stretch output cube along Y direction");
01030
01031 cflat = binCube(cflat2,1,2,0,63,0,63);
01032 destroy_cube(cflat2);
01033
01034 cub_ims=sinfoni_cube2imglist(cflat);
01035 if(-1 == sinfoni_pro_dump_ims(cub_ims,sof,sof,"out_mflat_cube.fits",
01036 "MFLAT_CUBE",plugin_id,config))
01037 {
01038 cpl_msg_error(_id,"cannot dump cube %s", "out_mflat_cube.fits");
01039 cpl_free(correct_dist);
01040 cpl_imagelist_delete(cub_ims);
01041 destroy_cube(cflat) ;
01042 destroy_cube(cube);
01043 cpl_free(times);
01044 cpl_free(offsetx);
01045 cpl_free(offsety);
01046 objnod_free(cfg);
01047 cpl_frameset_delete(stk);
01048
01049 return -1;
01050 }
01051 eima_avg=averageCubeToImage(cflat);
01052 ima_avg=cpl_image_wrap_float(eima_avg->lx,eima_avg->ly,eima_avg->data);
01053 if (ima_avg == NULL) {
01054 cpl_msg_error(_id,"Creating an average image");
01055 }
01056 if(-1 == sinfoni_pro_dump_ima(ima_avg,sof,sof,"out_mflat_avg.fits",
01057 "MFLAT_AVG",plugin_id,config))
01058 {
01059 cpl_msg_error(_id,"cannot dump ima %s", "out_mflat_avg.fits");
01060 cpl_free(correct_dist);
01061 cpl_imagelist_delete(cub_ims);
01062 destroy_cube(cflat) ;
01063 destroy_cube(cube);
01064 cpl_free(times);
01065 cpl_free(offsetx);
01066 cpl_free(offsety);
01067 cpl_image_unwrap(ima_avg);
01068 destroy_image(eima_avg);
01069 objnod_free(cfg);
01070 cpl_frameset_delete(stk);
01071 return -1;
01072
01073 }
01074 cpl_image_unwrap(ima_avg);
01075 destroy_image(eima_avg);
01076 cpl_imagelist_delete(cub_ims);
01077
01078 eima_med=medianCube(cflat);
01079 ima_med=cpl_image_wrap_float(eima_med->lx,eima_med->ly,eima_med->data);
01080 if(-1 == sinfoni_pro_dump_ima(ima_med,sof,sof,"out_mflat_med.fits",
01081 "MFLAT_MED",plugin_id,config))
01082 {
01083 cpl_msg_error(_id,"cannot dump ima %s", "out_mflat_med.fits");
01084 cpl_free(correct_dist);
01085 cpl_image_unwrap(ima_med);
01086 destroy_image(eima_med);
01087 destroy_cube(cflat) ;
01088 destroy_cube(cube);
01089 cpl_free(times);
01090 cpl_free(offsetx);
01091 cpl_free(offsety);
01092 objnod_free(cfg);
01093 cpl_frameset_delete(stk);
01094
01095 return -1;
01096 }
01097 destroy_cube(cflat);
01098 cpl_image_unwrap(ima_med);
01099 destroy_image(eima_med);
01100 }
01101 }
01102
01103
01104 if(NULL == (outcube2=fine_tune(cube,correct_dist,cfg->method,cfg->order,cfg->nslits))){
01105 cpl_msg_error (_id," could not fine tune the data cube\n") ;
01106 return -1 ;
01107 }
01108 cpl_msg_info(_id,"Stretch output cube along Y direction");
01109 outcube = binCube(outcube2,1,2,0,63,0,63);
01110 destroy_cube(cube);
01111
01112 sprintf(file_name,"%s%2.2d%s","out_cube_obj",n,".fits");
01113 cub_ims=sinfoni_cube2imglist(outcube);
01114 if(-1 == sinfoni_pro_dump_ims(cub_ims,sof,sof,file_name,
01115 pro_obs,plugin_id,config))
01116 {
01117 cpl_msg_error(_id,"cannot dump cube %s", file_name);
01118 cpl_free(correct_dist);
01119 destroy_cube(outcube2) ;
01120 destroy_cube(outcube) ;
01121 cpl_free(times);
01122 cpl_free(offsetx);
01123 cpl_free(offsety);
01124
01125 cpl_imagelist_delete(cub_ims);
01126 objnod_free(cfg);
01127 cpl_frameset_delete(stk);
01128 return -1;
01129
01130 }
01131 center_x = outcube->lx / 2. + 0.5 ;
01132 center_y = outcube->ly / 2. + 0.5 ;
01133
01134 set_wcs_cube(cub_ims, file_name, centralLambda, dis,
01135 centralpix, center_x, center_y);
01136 cpl_imagelist_delete(cub_ims);
01137
01138
01139
01140
01141 destroy_cube(outcube2);
01142 destroy_cube(outcube) ;
01143 cpl_free(correct_dist) ;
01144
01145 }
01146
01147
01148
01149 if(cfg->jitterind == 0) {
01150 goto exit;
01151 }
01152
01153
01154 if( (cfg->size_x*cfg->size_y*cfg->nframes) > (130*130*21) ) {
01155 cpl_msg_warning(_id,"Coadd cube size:%d,%d. N frames: %d",
01156 cfg->size_x,cfg->size_y,cfg->nframes);
01157 cpl_msg_warning(_id,"Max allowed should be such that sixeX*sixeY*Nframes < 130*130*21");
01158 goto exit;
01159 }
01160
01161
01162
01163 if ( cfg->jitterind == 1 )
01164 {
01165 p = cpl_parameterlist_find(config, "sinfoni.objnod.vllx");
01166 vllx = cpl_parameter_get_int(p);
01167 p = cpl_parameterlist_find(config, "sinfoni.objnod.vlly");
01168 vlly = cpl_parameter_get_int(p);
01169 p = cpl_parameterlist_find(config, "sinfoni.objnod.vurx");
01170 vurx = cpl_parameter_get_int(p);
01171 p = cpl_parameterlist_find(config, "sinfoni.objnod.vury");
01172 vury = cpl_parameter_get_int(p);
01173
01174
01175 cube_tmp = (OneCube**) cpl_calloc (cfg->nframes, sizeof (OneCube*));
01176 cubeobject = (OneCube**) cpl_calloc (cfg->nframes, sizeof (OneCube*));
01177 if ( cubeobject == NULL ) {
01178 cpl_msg_error(_id," could not allocate memory!") ;
01179 return -1 ;
01180 }
01181 if ( cube_tmp == NULL ) {
01182 cpl_msg_error(_id," could not allocate memory!") ;
01183 return -1 ;
01184 }
01185
01186 for ( n = 0 ; n < cfg->nframes ; n++ ) {
01187 sprintf(file_name,"%s%2.2d%s","out_cube_obj",n,".fits");
01188 cube_tmp[n] = load_cube(file_name);
01189 cubeobject[n] = cube_getvig(cube_tmp[n],1+vllx,1+vlly,
01190 64-vurx,64-vury);
01191 destroy_cube(cube_tmp[n]);
01192 }
01193 cpl_free(cube_tmp);
01194 }
01195
01196
01197
01198
01199
01200
01201
01202
01203
01204 if (cfg->jitterind == 1)
01205 {
01206 cpl_msg_info(_id,"Jittering...");
01207
01208 cpl_msg_info(_id,"Coadded cube size. x: %d y: %d",
01209 cfg->size_x,cfg->size_y);
01210
01211 jittercube = newCube(cfg->size_x, cfg->size_y, cubeobject[0]->np) ;
01212
01213
01214
01215
01216
01217
01218
01219
01220
01221 p = cpl_parameterlist_find(config,"sinfoni.objnod.scales_sky");
01222 scales_sky = cpl_parameter_get_bool(p);
01223 p = cpl_parameterlist_find(config,"sinfoni.objnod.ks_clip");
01224 ks_clip = cpl_parameter_get_bool(p);
01225 p = cpl_parameterlist_find(config,"sinfoni.objnod.kappa");
01226 kappa = cpl_parameter_get_double(p);
01227
01228
01229 if(scales_sky == 1) {
01230 cpl_msg_info(_id,"Subtract spatial median to each cube plane");
01231 for(n=0;n<cfg->nframes;n++) {
01232 cpl_msg_info(_id,"process cube %d\n",n);
01233 sinfoni_correct_median_it(&(cubeobject[n]));
01234 }
01235 }
01236 if (NullCube==(maskcube=newCube(jittercube->lx, jittercube->ly, jittercube->np)))
01237 {
01238 cpl_msg_error ("combineJitteredCubes:","could not allocate cube!");
01239 return -1;
01240 }
01241
01242
01243 cpl_msg_info(_id,"Combine jittered cubes");
01244
01245
01246 if(ks_clip == 1) {
01247 cpl_msg_info(_id,"Cube coaddition with kappa-sigma");
01248 }
01249
01250 for(z=0;z<cubeobject[0]->np;z+=z_stp) {
01251 z_siz=(z_stp < (cubeobject[0]->np-z)) ? z_stp : (cubeobject[0]->np-z);
01252 z_min=z;
01253 z_max=z_min+z_siz;
01254 cpl_msg_info(_id,"Coadding cubes: range [%4.4d,%4.4d) of 0-%d\n",
01255 z_min,z_max,cubeobject[0]->np);
01256
01257 if(ks_clip == 1) {
01258 combineJitteredCubesThomasRange(cubeobject,jittercube,maskcube,
01259 cfg->nframes,offsetx,offsety,
01260 times,cfg->kernel_type,
01261 z_min,z_max,kappa);
01262
01263 } else {
01264 combineJitteredCubesRange(cubeobject, jittercube, maskcube,
01265 cfg->nframes, offsetx,
01266 offsety, times,
01267 cfg->kernel_type,z_min,z_max) ;
01268 }
01269 }
01270 convert0ToZEROForCubes(jittercube) ;
01271
01272 if (jittercube == NULL)
01273 {
01274 cpl_msg_error(_id," could not allocate new data cube!") ;
01275 return -1 ;
01276 }
01277
01278 if (maskcube == NULL)
01279 {
01280 cpl_msg_error(_id," could not merge the jittered data cubes\n") ;
01281 cpl_free(cubeobject);
01282 return -1 ;
01283 }
01284
01285 for ( i = 0 ; i <cfg->nframes ; i++ ) {
01286 destroy_cube(cubeobject[i]);
01287 }
01288 cpl_free(cubeobject);
01289
01290 newcenter_x = cfg->size_x / 2. + 0.5 ;
01291 newcenter_y = cfg->size_y / 2. + 0.5 ;
01292
01293 jit_ims=sinfoni_cube2imglist(jittercube);
01294
01295 if(-1 == sinfoni_pro_dump_ims(jit_ims,sof,sof,cfg->outName,
01296 procatg,plugin_id,config)) {
01297 cpl_msg_error(_id,"cannot dump cube %s", cfg->outName);
01298 if (cfg->northsouthInd == 0)
01299 {
01300 for ( i = 0 ; i < cfg-> nslits ; i++ )
01301 {
01302 cpl_free(slit_edges[i]) ;
01303 }
01304 cpl_free(slit_edges) ;
01305 }
01306 destroy_cube(jittercube) ;
01307 destroy_cube(maskcube) ;
01308 cpl_imagelist_delete(jit_ims);
01309 cpl_free(times);
01310 cpl_free(offsetx);
01311 cpl_free(offsety);
01312 objnod_free(cfg);
01313 cpl_frameset_delete(stk);
01314
01315 return -1;
01316 }
01317 set_wcs_cube(jit_ims, cfg->outName, centralLambda, dis, centralpix, center_x, center_y);
01318
01319 cpl_imagelist_delete(jit_ims);
01320
01321 jitter_image = medianCube(jittercube);
01322 if (jitter_image == NULL) {
01323 cpl_msg_error(_id," could not do medianCube()\n");
01324 return -1;
01325 }
01326 wimg=cpl_image_wrap_float(jitter_image->lx,jitter_image->ly,jitter_image->data);
01327 img=cpl_image_duplicate(wimg);
01328 cpl_image_unwrap(wimg);
01329
01330 if(-1 == sinfoni_pro_dump_ima(img,sof,sof,cfg->med_cube_name,
01331 pro_med,plugin_id,config)) {
01332 cpl_msg_error(_id,"cannot dump ima %s", cfg->outName);
01333 cpl_image_delete(img);
01334 destroy_image(jitter_image);
01335 destroy_cube(maskcube) ;
01336 destroy_cube(jittercube);
01337 cpl_free(offsetx);
01338 cpl_free(offsety);
01339 cpl_free(times);
01340 objnod_free(cfg);
01341 cpl_frameset_delete(stk);
01342 return -1;
01343 }
01344 set_wcs_image(img, cfg->med_cube_name, center_x, center_y);
01345
01346 cpl_image_delete(img);
01347 destroy_image(jitter_image);
01348
01349
01350 msk_ims=sinfoni_cube2imglist(maskcube);
01351 if(-1 == sinfoni_pro_dump_ims(msk_ims,sof,sof,cfg->maskname,
01352 pro_mjit,plugin_id,config)) {
01353 cpl_msg_error(_id,"cannot dump cube %s", cfg->maskname);
01354 if (cfg->northsouthInd == 0)
01355 {
01356 for ( i = 0 ; i < cfg-> nslits ; i++ )
01357 {
01358 cpl_free(slit_edges[i]) ;
01359 }
01360 cpl_free(slit_edges) ;
01361 }
01362 destroy_cube(jittercube) ;
01363 destroy_cube(maskcube) ;
01364 cpl_imagelist_delete(msk_ims);
01365 cpl_free(times);
01366 cpl_free(offsetx);
01367 cpl_free(offsety);
01368 objnod_free(cfg);
01369 cpl_frameset_delete(stk);
01370 return -1;
01371 }
01372 set_wcs_cube(msk_ims, cfg->maskname, centralLambda, dis, centralpix, center_x, center_y);
01373
01374 cpl_free(times) ;
01375 cpl_free(offsetx) ;
01376 cpl_free(offsety) ;
01377 destroy_cube(maskcube) ;
01378 destroy_cube(jittercube) ;
01379 cpl_imagelist_delete(msk_ims);
01380
01381 }
01382
01383 exit:
01384
01385
01386 objnod_free(cfg);
01387 cpl_frameset_delete(stk);
01388 sinfoni_memory_status();
01389 return 0;
01390
01391
01392 for ( i = 0 ; i < cfg->nframes ; i++ ) {
01393 destroy_cube(cubeobject[i]) ;
01394 }
01395 cpl_free(cubeobject);
01396
01397 objnod_free(cfg);
01398 cpl_frameset_delete(stk);
01399 return 0;
01400
01401
01402 if (cfg->northsouthInd == 0)
01403 {
01404 for ( i = 0 ; i < cfg-> nslits ; i++ )
01405 {
01406 cpl_free(slit_edges[i]) ;
01407 }
01408 cpl_free(slit_edges) ;
01409 }
01410 objnod_free(cfg);
01411
01412 cpl_frameset_delete(stk);
01413 return 0;
01414
01415 }
01416
01417