objnod_rel.c

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

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