objspider.c

00001 /*----------------------------------------------------------------------------
00002  
00003      File name    :     objspider.c
00004    Author       :   A. Modigliani
00005    Created on   :   Sep 17, 2003
00006    Description  : 
00007 
00008 This recipe carries through the data cube creation of an object science 
00009 observation using the sky spider and jittering. This script expects jittered 
00010 frames that were already dark-subtracted averaged, flatfielded, spectral tilt 
00011 corrected and interleaved if necessary
00012 
00013 
00014  ---------------------------------------------------------------------------*/
00015 
00016 /*----------------------------------------------------------------------------
00017                                 Includes
00018  ---------------------------------------------------------------------------*/
00019 #include "objspider.h"
00020 #include "sinfoni_pro_save.h"
00021 #include "objspider_ini_by_cpl.h"
00022 
00023 /*----------------------------------------------------------------------------
00024                                 Defines
00025  ---------------------------------------------------------------------------*/
00026 
00027 
00028 /*----------------------------------------------------------------------------
00029                              Function Definitions
00030  ---------------------------------------------------------------------------*/
00031 
00032 /*----------------------------------------------------------------------------
00033    Function     :       objspider()
00034    In           :       ini_file: file name of according .ini file
00035    Out          :       integer (0 if it worked, -1 if it doesn't) 
00036    Job          : 
00037 This recipe carries through the data cube creation of an object science 
00038 observation using the sky spider and jittering. This script expects jittered 
00039 frames that were already dark-subtracted averaged, flatfielded, spectral tilt 
00040 corrected and interleaved if necessary
00041 
00042 
00043  ---------------------------------------------------------------------------*/
00044 
00045 
00046 
00047 void 
00048 change_plist_objspider(cpl_propertylist * plist, char * name,
00049                         char * outName, 
00050                         double cenLambda, double dispersion, int cenpix, 
00051                         double center_x, double center_y) 
00052 {
00053 
00054       float pixelscale ;
00055       float ra ;
00056       float dec ;
00057       float angle ;
00058       float radangle ;
00059       float cd1_1, cd1_2, cd2_1, cd2_2 ;
00060       char firsttext[2*FILENAMESZ] ;
00061       char * date ;
00062 
00063       pixelscale = spiffi_get_pixelscale(name)/2;
00064       ra = spiffi_get_RA(name);
00065       dec = spiffi_get_DEC(name);
00066       angle = spiffi_get_posangle(name);
00067       radangle = angle * PI_NUMB / 180.;
00068       cd1_1 = -cos(radangle) * pixelscale / 3600.;
00069       cd1_2 = sin(radangle) * pixelscale / 3600.;
00070       cd2_1 = sin(radangle) * pixelscale / 3600.;
00071       cd2_2 = cos(radangle) * pixelscale / 3600.;
00072 
00073 
00074 
00075       cpl_propertylist_insert_after_string(plist, "EXPTIME", "CTYPE1", "RA---TAN" ) ;
00076       cpl_propertylist_set_comment(plist, "CTYPE1", "Projected Rectascension" ) ;
00077 
00078       cpl_propertylist_insert_after_float(plist, "CTYPE1",  "CRPIX1", center_x ) ;
00079       cpl_propertylist_set_comment(plist, "CRPIX1","Reference pixel in RA" ) ;
00080 
00081  
00082 
00083       cpl_propertylist_insert_after_float(plist, "CRPIX1",  "CRVAL1", ra ) ;
00084       cpl_propertylist_set_comment(plist, "CRVAL1","Reference RA" ) ;
00085 
00086       cpl_propertylist_insert_after_float(plist, "CRVAL1",  "CDELT1", -pixelscale/3600.  ) ;
00087       cpl_propertylist_set_comment(plist, "CDELT1","pixel scale" ) ;
00088 
00089       cpl_propertylist_insert_after_string(plist, "CDELT1",  "CUNIT1", "DEGREE" ) ;
00090       cpl_propertylist_set_comment(plist, "CUNIT1","RA-UNIT" ) ;
00091 
00092       cpl_propertylist_insert_after_string(plist, "CUNIT1",  "CTYPE2", "DEC--TAN" ) ;
00093       cpl_propertylist_set_comment(plist, "CTYPE2", "Projected Declination") ;
00094 
00095 
00096 
00097 
00098       cpl_propertylist_insert_after_float(plist, "CTYPE2",  "CRPIX2", center_y ) ;
00099       cpl_propertylist_set_comment(plist, "CRPIX2", "Reference pixel in DEC") ;
00100 
00101       cpl_propertylist_insert_after_float(plist, "CRPIX2",  "CRVAL2", dec) ;
00102       cpl_propertylist_set_comment(plist, "CRVAL2", "Reference DEC") ;
00103 
00104       cpl_propertylist_insert_after_float(plist, "CRVAL2",  "CDELT2", pixelscale/3600. ) ;
00105       cpl_propertylist_set_comment(plist, "CDELT2", "pixel scale") ;
00106 
00107       cpl_propertylist_insert_after_string(plist, "CDELT2",  "CUNIT2", "DEGREE" ) ;
00108       cpl_propertylist_set_comment(plist, "CUNIT2", "DEC-UNIT") ;
00109 
00110       cpl_propertylist_insert_after_string(plist, "CUNIT2",  "CTYPE3", "WAVE" ) ;
00111       cpl_propertylist_set_comment(plist, "CTYPE3", "wavelength axis in microns") ;
00112 
00113 
00114 
00115 
00116       cpl_propertylist_insert_after_int(plist, "CTYPE3",  "CRPIX3", cenpix+1 ) ;
00117       cpl_propertylist_set_comment(plist, "CRPIX3", "Reference pixel in z") ;
00118 
00119       cpl_propertylist_insert_after_float(plist, "CRPIX3",  "CRVAL3", cenLambda ) ;
00120       cpl_propertylist_set_comment(plist, "CRVAL3", "central wavelength") ;
00121 
00122       cpl_propertylist_insert_after_float(plist, "CRVAL3",  "CDELT3", dispersion ) ;
00123       cpl_propertylist_set_comment(plist, "CDELT3", "microns per pixel") ;
00124 
00125       cpl_propertylist_insert_after_string(plist, "CDELT3",  "CUNIT3", "MICRON" ) ;
00126       cpl_propertylist_set_comment(plist, "CUNIT3",  "spectral unit" ) ;
00127 
00128 
00129 
00130       cpl_propertylist_insert_after_float(plist, "CUNIT3",  "CD1_1", cd1_1 ) ;
00131       cpl_propertylist_set_comment(plist, "CD1_1",  "CD rotation matrix" ) ;
00132 
00133       cpl_propertylist_insert_after_float(plist, "CD1_1",   "CD1_2", cd1_2 ) ;
00134       cpl_propertylist_set_comment(plist, "CD1_2",  "CD rotation matrix" ) ;
00135 
00136       cpl_propertylist_insert_after_float(plist, "CD1_2",   "CD2_1", cd2_1 ) ;
00137       cpl_propertylist_set_comment(plist, "CD2_1",  "CD rotation matrix" ) ;
00138 
00139       cpl_propertylist_insert_after_float(plist, "CD2_1",   "CD2_2", cd2_2 ) ;
00140       cpl_propertylist_set_comment(plist, "CD2_2",  "CD rotation matrix" ) ;
00141 
00142 
00143 
00144       date     = get_datetime_iso8601() ;
00145 
00146       strcpy(firsttext, "spiffi_objspider -f \0") ;
00147       
00148       if (strlen(firsttext) > 30) {
00149         cpl_propertylist_set_string(plist,KEY_NAME_PIPEFILE,  outName ) ;
00150         cpl_propertylist_set_string(plist,KEY_NAME_HPRO_TYPE,  "REDUCED" ) ;
00151         cpl_propertylist_set_string(plist,KEY_NAME_HPRO_CATG,  "Science Observation" ) ;
00152         cpl_propertylist_set_string(plist,KEY_NAME_HPRO_STATUS, "OK" ) ;
00153         cpl_propertylist_set_string(plist,KEY_NAME_HPRO_DATE, date ) ;
00154         cpl_propertylist_set_string(plist,KEY_NAME_HPRO_RECID, firsttext ) ;
00155         cpl_propertylist_set_string(plist,KEY_NAME_HPRO_DRSID, PACKAGE_VERSION ) ;  
00156 
00157       }
00158 
00159 }
00160 
00161 
00162 
00163 
00164 
00165 int objspider (cpl_parameterlist* config, cpl_frameset* sof)
00166 {
00167       const char* _id = "objspider";
00168     object_config * cfg ;
00169 
00170     OneImage * im ;
00171     OneImage * im_wave_map ;
00172     OneImage * im_resampled ;
00173     OneImage * im_halo_spec ;
00174     OneImage * im_cal ;
00175 
00176 
00177     OneCube * cube;
00178     OneCube * cube_out;
00179     OneCube ** cube_sub;
00180     OneCube * cube_jitter;
00181     OneCube * cube_mask;
00182     OneCube ** cube_object;
00183     char blank[FILENAMESZ] ;
00184     fits_header ** head;
00185 
00186     Vector * vec_sky_spec;
00187 
00188 
00189     int i = 0;
00190     int fra=0;
00191     int cenpix = 0;
00192     int lx=0;
00193     int ly=0;
00194 
00195     int* pos;
00196     int* posindex;
00197     int* centralpix;
00198 
00199     char* name;
00200   
00201     char* name_jitter;
00202   
00203     float exptime=0;
00204     float offx=0;
00205     float offy=0;
00206     float offx_rot=0;
00207     float offy_rot=0;
00208     float offx_rot_pix=0;
00209     float offy_rot_pix=0;
00210     float pixelscale ;
00211     float angle ;
00212     float radangle ;
00213     float cd1_1, cd1_2, cd2_1, cd2_2 ;
00214     float min = 0;
00215     float max = 0;
00216     float center_x=0;
00217     float center_y=0;
00218     float fc=0;
00219     float new_center_x=0;
00220     float new_center_y=0;
00221       
00222     float* distances;
00223     float* times;
00224     float* offsetx;
00225     float* offsety;
00226     float* correct_dist;
00227     float* mi;
00228     float* ma;
00229 
00230     float** slit_edges;
00231 
00232     double disp = 0;
00233     double cenLam = 0;
00234 
00235     double* dis;
00236     double* centralLambda;
00237 
00238 
00239     float edge_x=0;
00240     float edge_y=0;
00241     float tmp_float=0;
00242     char* tbl_slitpos_name=NULL;
00243     cpl_table* tbl_slitpos=NULL;
00244     char* tbl_distances_name=NULL;
00245     cpl_table* tbl_distances=NULL;
00246     char* tbl_firstcol_name=NULL;
00247     cpl_table* tbl_firstcol=NULL;
00248     int* status=NULL;
00249 
00250     cpl_imagelist* cub_ims=NULL;
00251     cpl_imagelist* msk_ims=NULL;
00252     cpl_imagelist* jit_ims=NULL;
00253 
00254     OneImage * eimg=NULL;
00255     cpl_image* cimg=NULL;
00256     cpl_frameset* stk=NULL;
00257     char* key_value=NULL;
00258     cpl_propertylist * cplist=NULL;
00259 
00260     
00261     /* 
00262        -----------------------------------------------------------------
00263        1) parse the file names and parameters to the object_config data 
00264           structure cfg
00265        -----------------------------------------------------------------
00266      */
00267     cpl_msg_info(_id,"Parse cpl input");
00268     cfg = parse_cpl_input_objspider(config,sof,&stk);
00269     if(cpl_error_get_code() != CPL_ERROR_NONE) {
00270       cpl_msg_error(_id,"error parsing cpl input");
00271       cpl_msg_error(_id,(char* ) cpl_error_get_message());
00272       return -1;
00273     }
00274     if (cfg == NULL) {
00275       cpl_msg_error(_id," could not parse cpl input!\n");
00276       return -1;
00277     }
00278 
00279     if(is_fits_file(cfg->wavemap) != 1) {
00280       cpl_msg_error(_id,"Input file %s is not FITS",cfg->wavemap);
00281       return -1;
00282     }
00283 
00284     if (cfg->halocorrectInd == 1) {
00285        if(is_fits_file(cfg->halospectrum) != 1) {
00286          cpl_msg_error(_id,"Input file %s is not FITS",cfg->halospectrum);
00287          return -1;
00288        }
00289     }
00290 
00291     if (is_fits_file(cfg->poslist) != 1) {
00292       cpl_msg_error(_id,"poslist not given!");
00293       return -1;
00294     }
00295 
00296     im_wave_map = load_image(cfg->wavemap);
00297     if (im_wave_map == NULL){
00298       cpl_msg_error(_id, "could not load wavemap\n" );
00299       return -1;
00300     }
00301 
00302     /* read either the distances list or the slitlet edge position list */
00303     if (cfg->northsouthInd == 0) {
00304       cpl_msg_info(_id,"Read slitlet edges positions");
00305       slit_edges = new_2Dfloat_array(cfg->nslits, 2);
00306       /* OLD WAY
00307       if ( NULL == (slit_list = fopen(cfg->poslist,  "r" )) ) {
00308     
00309     cpl_msg_error(_id," could not open poslist!\n") ;
00310     return -1 ;
00311       }
00312       */
00313 
00314        /*READ TFITS TABLE*/
00315         if(cpl_error_get_code() != CPL_ERROR_NONE) {
00316           cpl_msg_error(_id,"error before loading cpl input");
00317           cpl_msg_error(_id,(char* ) cpl_error_get_message());
00318           return -1;
00319         }
00320         tbl_slitpos_name = (char*) cpl_calloc(MAX_NAME_SIZE,sizeof(char*));
00321         strcpy(tbl_slitpos_name,cfg->poslist);
00322         tbl_slitpos = cpl_table_load(tbl_slitpos_name,1,0);
00323 
00324         if(cpl_error_get_code() != CPL_ERROR_NONE) {
00325            cpl_msg_error(_id,"error loading cpl table %s",tbl_slitpos_name);
00326           return -1;
00327         }
00328         for (i =0 ; i< cfg->nslits; i++){
00329            edge_x=cpl_table_get_float(tbl_slitpos,"start",i,status);
00330 
00331            edge_y=cpl_table_get_float(tbl_slitpos,"end",i,status);
00332            array2D_set_value(slit_edges,edge_x,i,0);
00333            array2D_set_value(slit_edges,edge_y,i,1);
00334         }        
00335         if(cpl_error_get_code() != CPL_ERROR_NONE) {
00336            cpl_msg_error(_id,"error reading cpl table %s",tbl_slitpos_name);
00337           return -1;
00338         }
00339 
00340         cpl_table_delete(tbl_slitpos);
00341         cpl_free(tbl_slitpos_name);
00342     }
00343     else {
00344       /* open the ASCII list of the determined distances */
00345       cpl_msg_info(_id,"Read slitlets distances");
00346         if ( NULL == (distances = (float*) 
00347                       cpl_calloc (cfg->nslits - 1, sizeof (float))) ) {
00348     
00349         cpl_msg_error (_id,"could allocate memory!\n") ;
00350             return -1 ;
00351     }
00352 
00353       /*READ TFITS TABLE*/
00354         tbl_distances_name = (char*) cpl_calloc(MAX_NAME_SIZE,sizeof(char*));
00355         strcpy(tbl_distances_name,cfg->poslist);
00356         tbl_distances = cpl_table_load(tbl_distances_name,1,0);
00357         if(cpl_error_get_code() != CPL_ERROR_NONE) {
00358            cpl_msg_error(_id,"loading cpl table %s",tbl_distances_name);
00359           return -1;
00360         }
00361         for (i =0 ; i< cfg->nslits-1; i++){
00362            tmp_float=cpl_table_get_float(tbl_distances,"slitlet_distance",i,status);
00363         if(cpl_error_get_code() != CPL_ERROR_NONE) {
00364            cpl_msg_error(_id,"reading table %s raw %d val %f",tbl_distances_name,i,tmp_float);
00365           return -1;
00366         }
00367            array_set_value(distances,tmp_float,i);
00368         }
00369         cpl_table_delete(tbl_distances);
00370         cpl_free(tbl_distances_name);
00371         if(cpl_error_get_code() != CPL_ERROR_NONE) {
00372            cpl_msg_error(_id,"reading cpl table %s",tbl_distances_name);
00373           return -1;
00374         }
00375 
00376       /* open the ASCII list of the determined distances */
00377         if ( NULL ==  (distances = (float*) 
00378                cpl_calloc (cfg->nslits - 1, sizeof (float))) ) {
00379         cpl_msg_error (_id,"could allocate memory!\n") ;
00380             return -1 ;
00381     } 
00382     }
00383     if (cfg->jitterind == 1) {
00384     cube_object = new_cube_list(cfg->nframes);
00385     if (cube_object == NULL) {
00386       cpl_msg_error(_id," could not allocate memory\n");
00387       return -1;
00388     }
00389     times   = new_float_array(cfg->nframes);
00390     offsetx = new_float_array(cfg->nframes);
00391     offsety = new_float_array(cfg->nframes);
00392     }
00393 
00394     head = (fits_header**) cpl_calloc(cfg -> nframes, sizeof(fits_header*));
00395     cube_sub = (OneCube**) cpl_calloc(cfg -> nframes, sizeof(OneCube*));
00396 
00397     for (fra=0; fra< cfg->nframes; fra++) {
00398       cpl_msg_info(_id,"Read FITS information");
00399     name = cfg->framelist[fra];
00400     if (fra == 0) {
00401        name_jitter = name;
00402     }
00403         if(is_fits_file(name) != 1) {
00404            cpl_msg_error(_id,"Input file %s is not FITS",name);
00405            return -1;
00406     }
00407     /* get some header values and compute the CD-matrix */
00408     pixelscale = spiffi_get_pixelscale(name)/2;
00409     angle = spiffi_get_posangle(name);
00410     radangle = angle * PI_NUMB / 180.;
00411     cd1_1 = cos(radangle);
00412     cd1_2 = -sin(radangle);
00413     cd2_1 = sin(radangle);
00414     cd2_2 = cos(radangle);
00415 
00416     cpl_msg_info(_id,"fra=%d name=%s\n", fra, name);
00417 
00418     im = load_image(name);
00419     if (im == NULL) {
00420       cpl_msg_error(_id, " could not load frame\n" );
00421       return -1;
00422     }
00423     head[fra]=fits_read_header(name);
00424     if (cfg->jitterind == 1) {
00425       /* exptime = spiffi_get_exptime(name); */
00426       exptime = spiffi_get_ditndit(name);
00427 
00428       if (exptime == FLAG) {
00429         cpl_msg_error(_id,"could not read fits header keyword!\n");
00430         return -1;
00431       }
00432       array_set_value(times, exptime, fra);
00433 
00434       if (fra == 0) {
00435         offx = 0.;
00436         offy = 0.;
00437       }
00438       else {
00439         offx    = spiffi_get_cumoffsetx(name);
00440         if (offx == FLAG) {
00441           cpl_msg_error(_id,"could not read fits header keyword!\n");
00442           return -1;
00443         }
00444         offy    = spiffi_get_cumoffsety(name);
00445         if (offy == FLAG) {
00446           cpl_msg_error(_id,"could not read fits header keyword!\n");
00447           return -1;
00448         }
00449 
00450         /* rotate the coordinates */
00451         offx_rot = cd1_1 * offx + cd2_1 * offy;
00452         offy_rot = cd1_2 * offx + cd2_2 * offy;
00453         /* convert the coordinates to pixel units */
00454         offx_rot_pix = offx_rot / pixelscale;
00455         offy_rot_pix = offy_rot / pixelscale;
00456         array_set_value(offsetx, offx_rot_pix, fra);
00457         array_set_value(offsety, offy_rot_pix, fra);
00458 
00459   
00460       }
00461 
00462     }
00463 
00464 
00465     /*
00466      ------------------------------------------------------------------------
00467      --------------------------RESAMPLING------------------------------------
00468      ------------------------------------------------------------------------
00469     */
00470       cpl_msg_info(_id,"Resampling");
00471     dis = new_double_array(1);
00472     mi  = new_float_array(1);
00473     ma  = new_float_array(1);
00474     centralLambda  = new_double_array(1);
00475     centralpix = new_int_array(1);
00476 
00477     doublearray_set_value(dis, 0., 0);
00478     array_set_value(mi, 0., 0);
00479     array_set_value(ma, 0., 0);
00480     doublearray_set_value(centralLambda, 0., 0);
00481     intarray_set_value(centralpix, 0, 0);
00482 
00483         im_resampled = definedResampling( im, im_wave_map, 
00484                                           cfg->ncoeffs, &cfg->nrows, 
00485                                         dis, mi, ma, centralLambda, 
00486                                         centralpix );
00487 
00488     if (im_resampled == NULL) {
00489       cpl_msg_error(_id, "definedResampling failed\n" );
00490       return -1;
00491     }
00492     disp = doublearray_get_value(dis, 0);
00493     min = array_get_value(mi, 0);
00494     max = array_get_value(ma, 0);
00495     cenLam = doublearray_get_value(centralLambda, 0);
00496     cenpix = intarray_get_value(centralpix, 0);
00497 
00498     cpl_msg_info(_id,"dispersion %f\n", disp);
00499     cpl_msg_info(_id,"minimum wavelength %f\n", min);
00500     cpl_msg_info(_id,"maximum wavelength %f\n", max);
00501     cpl_msg_info(_id,"central wavelength %f\n", cenLam);
00502     cpl_msg_info(_id,"central pixel %d\n", cenpix);
00503     destroy_doublearray(dis);
00504     destroy_array(mi);
00505     destroy_array(ma);
00506     destroy_doublearray(centralLambda);
00507     destroy_intarray(centralpix);
00508     destroy_image(im);
00509 
00510     /*
00511       ------------------------------------------------------------------------
00512       --------------------------CALIBRATION-----------------------------------
00513       ------------------------------------------------------------------------
00514       ---Multiply with calibrated halogen lamp spectrum-----------------------
00515     */
00516 
00517     if (cfg->halocorrectInd == 1) {
00518           cpl_msg_info(_id,"Calibration");
00519       im_halo_spec = load_image(cfg->halospectrum);
00520       im_cal    = multiplyImageWithSpectrum( im_resampled, im_halo_spec);
00521       if (im_cal == NULL) {
00522         cpl_msg_error(_id, " calibration failed\n" );
00523         return -1;
00524       }
00525       destroy_image(im_halo_spec);
00526       destroy_image(im_resampled);
00527       im_resampled = im_cal;
00528     }
00529 
00530 
00531       /*
00532       -------------------------------------------------------------------------
00533       -----------------------------CUBECREATION--------------------------------
00534       -------------------------------------------------------------------------
00535       ---select north-south-test or fitting of slitlet edges---
00536       */
00537 
00538           cpl_msg_info(_id,"Cube creation");
00539     correct_dist = new_float_array(cfg->nslits) ;
00540     if (cfg->northsouthInd == 0) {
00541       /* build the data cube */
00542       cube = makeCubeSpi( im_resampled, slit_edges, correct_dist );
00543 
00544       if (cube == NULL) {
00545         cpl_msg_error(_id,"could not construct data cube!\n");
00546         return -1;
00547       }
00548     }
00549     else {
00550            /*READ TFITS TABLE*/
00551             tbl_firstcol_name = (char*) cpl_calloc(MAX_NAME_SIZE,sizeof(char*));
00552             strcpy(tbl_firstcol_name,cfg->firstCol);
00553             tbl_firstcol = cpl_table_load(tbl_firstcol_name,1,0);
00554             if(cpl_error_get_code() != CPL_ERROR_NONE) {
00555                cpl_msg_error(_id,"error loading cpl table %s",tbl_firstcol_name);
00556                return -1;
00557         }
00558 
00559             fc=cpl_table_get_float(tbl_firstcol,"first_col",0,status);
00560             if(cpl_error_get_code() != CPL_ERROR_NONE) {
00561                cpl_msg_error(_id,"error reading cpl table %s",tbl_firstcol_name);
00562                return -1;
00563         }
00564 
00565             cpl_table_delete(tbl_firstcol);
00566             cpl_free(tbl_firstcol_name);
00567 
00568       cube = makeCubeDist( im_resampled, fc, distances, correct_dist );
00569 
00570       if (cube == NULL) {
00571         cpl_msg_error(_id," could not construct a data cube!\n");
00572         return -1;
00573       }
00574     }
00575 
00576 
00577       /*
00578     --------------------------------------------------------------------------
00579     ------------------------FINETUNING----------------------------------------
00580     --------------------------------------------------------------------------
00581     shift the rows of the reconstructed images of the data cube to the 
00582     correct sub pixel position select the shift method: polynomial 
00583     interpolation, FFT or cubic spline interpolation
00584       */
00585 
00586          cpl_msg_info(_id,"Fine tuning, method: %s", cfg->method);
00587     if (strcmp(cfg->method, "P")==0) {
00588       cube_out = fineTuneCube( cube, correct_dist, cfg->order );
00589 
00590       if (cube_out == NULL) {
00591         cpl_msg_error(_id," could not fine tune the data cube\n");
00592         return -1;
00593       }
00594     }
00595     else if (strcmp(cfg->method, "F")==0) {
00596       for (i=0; i< cfg->nslits; i++) {
00597         array_set_value(correct_dist, -array_get_value(correct_dist,i), i);
00598       }
00599       cube_out = fineTuneCubeByFFT( cube, correct_dist );
00600       if (cube_out == NULL) {
00601         cpl_msg_error(_id," could not fine tune the data cube\n");
00602         return -1;
00603       }
00604     }
00605     else if (strcmp(cfg->method, "S")==0) {
00606       cube_out = fineTuneCubeBySpline( cube, correct_dist );
00607       if (cube_out == NULL) {
00608           cpl_msg_error(_id," could not fine tune the data cube\n");
00609           return -1;
00610       }
00611     }
00612     else {
00613         cpl_msg_error(_id," wrong method indicator given!");
00614         return -1; 
00615     }
00616     
00617     /* ---free memory---- */
00618     destroy_image(im_resampled);
00619     destroy_cube (cube);
00620     destroy_array(correct_dist);
00621 
00622 
00623       /*
00624     ---------------------------------------------------------------------------
00625     -----------------------SKY EXTRACTION AND SUBTRACTION----------------------
00626     ---------------------------------------------------------------------------
00627       */
00628         cpl_msg_info(_id,"Sky extraction and subtraction");
00629     posindex = new_int_array(1);
00630     pos = spiffi_get_skyspiderposindex(name, posindex);
00631     if (pos == NULL) {
00632       cpl_msg_error(_id,"could not find sky spider position index in fits header\n");
00633       /* return -1; */
00634     }
00635 
00636         vec_sky_spec = extractSkyFromCube(cube_out, cfg->loReject, 
00637                                    cfg->hiReject, pos, 
00638                                    cfg->tolerance, 
00639                                   intarray_get_value(posindex, 0));
00640 
00641     if (vec_sky_spec == NULL) {
00642       cpl_msg_error(_id,"could not extract sky spectrum from data cube\n");
00643       return -1;
00644     }
00645     destroy_intarray(posindex);
00646     destroy_intarray(pos);
00647 
00648     cube_sub[fra]=subSpectrumFromCube(cube_out, vec_sky_spec);
00649 
00650     if (cube_sub[fra] == NULL) {
00651       cpl_msg_error(_id," could not subtract sky spectrum from data cube\n");
00652       return -1;
00653     }
00654 
00655     destroyVector(vec_sky_spec);
00656     destroy_cube(cube_out);
00657 
00658       /*
00659         --------------------------------------------------------------------
00660         -----------------JITTERING------------------------------------------
00661         --------------------------------------------------------------------
00662       */
00663         cpl_msg_info(_id,"Jittering");
00664     if (cfg->jitterind == 1) {
00665       insert_cube(cube_object, cube_sub[fra],  fra);
00666     }
00667     if (cfg->jitterind == 0) {
00668 
00669         if (strstr(cfg->outName, "." ) != NULL ) {
00670           sprintf(name_jitter, "%s%d%s",
00671               get_rootname(cfg->outName),fra,strstr(cfg->outName,"."));
00672 
00673         }
00674             else {
00675           sprintf(name_jitter,"%s%d",cfg->outName,fra);
00676         }
00677 
00678         lx = cube_out->lx;
00679         ly = cube_out->ly;
00680         center_x = lx / 2. + 0.5;
00681         center_y = ly / 2. + 0.5;
00682 
00683             cpl_msg_info(_id,"copy eclipse cube to CPL imset");
00684             cub_ims=cpl_imagelist_new(cube_sub[fra]->lx,cube_sub[fra]->ly,
00685                                   cube_sub[fra]->np,CPL_TYPE_FLOAT);
00686             for(i=0;i<cube_sub[fra]->np;i++) {
00687               eimg=cube_sub[fra]->plane[i];
00688               cimg=cpl_image_wrap_float(eimg->lx,eimg->ly,eimg->data);
00689               cpl_imagelist_set(cub_ims,cimg,i);
00690             }
00691 
00692 
00693             if(-1 == sinfoni_pro_dump_ims(cub_ims,stk,sof,name_jitter,
00694                     PRO_CUBE,_id,config)) {
00695                 cpl_msg_error(_id,"cannot dump cube %s", name_jitter);
00696                 return -1;
00697             }
00698 
00699             cplist = cpl_propertylist_new() ;
00700             if ((cpl_error_code)((cplist = cpl_propertylist_load(name_jitter, 0)) == NULL)) {
00701                cpl_msg_error(_id, "getting header from frame %s",name_jitter);
00702                cpl_propertylist_delete(cplist) ;
00703                return -1 ;
00704         }
00705 
00706         change_plist_objspider(cplist, name, name_jitter, cenLam, disp, cenpix,
00707                                    center_x, center_y );
00708 
00709         if (fra != 0) {
00710           sprintf(blank,"%f",offx_rot_pix);
00711           cpl_propertylist_set_string(cplist,  "HIERARCH ESO SEQ CUMOFFSETX", blank);
00712 
00713           sprintf(blank,"%f",offy_rot_pix);
00714           cpl_propertylist_set_string(cplist,  "HIERARCH ESO SEQ CUMOFFSETY", blank);
00715         }
00716 
00717             if (cpl_imagelist_save(cub_ims, name_jitter, CPL_BPP_DEFAULT, cplist,CPL_IO_DEFAULT)!=CPL_ERROR_NONE) {
00718                cpl_msg_error(_id, "Cannot save the product %s",name_jitter);
00719                cpl_propertylist_delete(cplist) ;
00720                return -1 ;
00721         }
00722             cpl_propertylist_delete(cplist) ;
00723 
00724 
00725             cpl_free(key_value);
00726 
00727     }
00728     
00729     }
00730 
00731     if (cfg->jitterind == 1) {
00732       /* Check removed
00733          if (cfg->size_x < cube_sub[0]->lx ||
00734              cfg->size_x > 2*(cube_sub[0]->lx) ||
00735              cfg->size_y < cube_sub[0]->ly ||
00736          cfg->size_y > 2*(cube_sub[0]->ly)) {
00737            cpl_msg_error(_id,"given final combined cube sizes are wrong!\n");
00738            return -1;
00739      }
00740       */
00741          cube_jitter = newCube(cfg->size_x, cfg->size_y, 
00742                                cube_sub[0]->np);       
00743      if (cube_jitter == NULL) {
00744        cpl_msg_error(_id,"could allocate new data cube\n");
00745        return -1;
00746      }
00747          cube_mask = combineJitteredCubes(cube_object, cube_jitter, 
00748                                           cfg->nframes, 
00749                           offsetx, offsety, 
00750                                           times, cfg->kernel_type);
00751 
00752      if (cube_mask == NULL) {
00753        cpl_msg_error(_id,"could not merge the jittered data cubes\n");
00754        return -1;
00755      }
00756 
00757      new_center_x = cfg->size_x / 2. + 0.5;
00758      new_center_y = cfg->size_y / 2. + 0.5;
00759 
00760 
00761 
00762 
00763             cpl_msg_info(_id,"copy eclipse cube to CPL imset");
00764             jit_ims=cpl_imagelist_new(cube_jitter->lx,cube_jitter->ly,cube_jitter->np,CPL_TYPE_FLOAT);
00765             for(i=0;i<cube_jitter->np;i++) {
00766               eimg=cube_jitter->plane[i];
00767               cimg=cpl_image_wrap_float(eimg->lx,eimg->ly,eimg->data);
00768               cpl_imagelist_set(cub_ims,cimg,i);
00769             }
00770 
00771             if(-1 == sinfoni_pro_dump_ims(jit_ims,stk,sof,cfg->outName,
00772                     PRO_CUBE,_id,config)) {
00773                 cpl_msg_error(_id,"cannot dump cube %s", cfg->outName);
00774                 return -1;
00775             }
00776 
00777 
00778             cplist = cpl_propertylist_new() ;
00779             if ((cpl_error_code)((cplist = cpl_propertylist_load(cfg->outName, 0)) == NULL)) {
00780                cpl_msg_error(_id, "getting header from frame %s",cfg->outName);
00781                cpl_propertylist_delete(cplist) ;
00782                return -1 ;
00783         }
00784 
00785     
00786             change_plist_objspider(cplist, name_jitter, cfg->outName, cenLam, disp, 
00787              cenpix, new_center_x, new_center_y);
00788 
00789             sprintf(blank,"%d",cfg->nframes);
00790             cpl_propertylist_set_string(cplist,  "HIERARCH ESO PRO DATANCOM",blank);
00791 
00792             if (cpl_imagelist_save(jit_ims, cfg->outName, CPL_BPP_DEFAULT, cplist,CPL_IO_DEFAULT)!=CPL_ERROR_NONE) {
00793                cpl_msg_error(_id, "Cannot save the product %s",cfg->outName);
00794                cpl_propertylist_delete(cplist) ;
00795                return -1 ;
00796         }
00797             cpl_propertylist_delete(cplist) ;
00798 
00799             cpl_free(key_value);
00800 
00801             cpl_msg_info(_id,"copy eclipse cube to CPL imset");
00802             msk_ims=cpl_imagelist_new(cube_mask->lx,cube_mask->ly,cube_mask->np,CPL_TYPE_FLOAT);
00803             for(i=0;i<cube_mask->np;i++) {
00804               eimg=cube_mask->plane[i];
00805               cimg=cpl_image_wrap_float(eimg->lx,eimg->ly,eimg->data);
00806               cpl_imagelist_set(msk_ims,cimg,i);
00807             }
00808 
00809 
00810 
00811            if(-1 == sinfoni_pro_dump_ims(msk_ims,stk,sof,cfg->maskname,
00812                     PRO_MASK_CUBE,_id,config)) {
00813                 cpl_msg_error(_id,"cannot dump cube %s", cfg->maskname);
00814                 return -1;
00815             }
00816 
00817 
00818             cplist = cpl_propertylist_new() ;
00819             if ((cpl_error_code)((cplist = cpl_propertylist_load(cfg->maskname, 0)) == NULL)) {
00820                cpl_msg_error(_id, "getting header from frame %s",cfg->maskname);
00821                cpl_propertylist_delete(cplist) ;
00822                return -1 ;
00823         }
00824 
00825     
00826             change_plist_objspider(cplist, name_jitter, cfg->maskname, cenLam, disp, 
00827              cenpix, new_center_x, new_center_y);
00828 
00829             sprintf(blank,"%d",cfg->nframes);
00830             cpl_propertylist_set_string(cplist,  "HIERARCH ESO PRO DATANCOM",blank);
00831 
00832             if (cpl_imagelist_save(msk_ims, cfg->maskname, CPL_BPP_DEFAULT, cplist,CPL_IO_DEFAULT)!=CPL_ERROR_NONE) {
00833                cpl_msg_error(_id, "Cannot save the product %s",cfg->maskname);
00834                cpl_propertylist_delete(cplist) ;
00835                return -1 ;
00836         }
00837             cpl_propertylist_delete(cplist) ;
00838 
00839 
00840            cpl_free(key_value);
00841       
00842             /* Get the reference file  */
00843 
00844      destroy_array(times);
00845      destroy_array(offsetx);
00846      destroy_array(offsety);
00847      destroy_cube(cube_mask);
00848      destroy_cube(cube_jitter);
00849      destroy_cube_list(cube_object);
00850     }
00851 
00852     /*---free memory----*/
00853     destroy_image(im_wave_map);
00854     for (i=0; i< cfg->nframes; i++) {
00855     destroy_cube(cube_sub[i]);
00856     fits_header_destroy(head[i]);
00857     }
00858     cpl_free(head);
00859     cpl_free(cube_sub);
00860 
00861     if (cfg->northsouthInd == 0) {
00862     destroy_2Darray(slit_edges, cfg->nslits);
00863     }
00864     else {
00865       cpl_free(distances);
00866     }
00867 
00868     objspider_free(cfg);
00869     return 0;
00870 }

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