cube_create_from_resampled.c

00001 /*----------------------------------------------------------------------------
00002    
00003    File name    :       cube_create_from_resampled.c
00004    Author       :   J. Schreiber
00005    Created on   :   May 7, 2003
00006    Description  :   Creates a data cube and can determine the edge positions
00007                         of the slitlets
00008  ---------------------------------------------------------------------------*/
00009 
00010 /*----------------------------------------------------------------------------
00011                                 Includes
00012  ---------------------------------------------------------------------------*/
00013 #include "cube_create_from_resampled.h"
00014 
00015 /*----------------------------------------------------------------------------
00016                                 Defines
00017  ---------------------------------------------------------------------------*/
00018 
00019 
00020 /*----------------------------------------------------------------------------
00021                              Function Definitions
00022  ---------------------------------------------------------------------------*/
00023 
00024 /*----------------------------------------------------------------------------
00025    Function     :       change_plist_cube_create_from_resampled()
00026    In           :       many needed inputs
00027    Out          :       nothing 
00028    Job          :       this routine changes the fits header for the data cube
00029                         file using appropriate values.
00030  ---------------------------------------------------------------------------*/
00031 
00032 void change_plist_cube_create_from_resampled (cpl_propertylist * plist, 
00033                     char * ini_name,
00034                     char * name,
00035             char * outName,
00036             float cenLambda,
00037             float dispersion,
00038             int   cenpix,
00039             float center_x,
00040             float center_y )
00041 {
00042 
00043     float pixelscale ;
00044     float ra ;
00045     float dec ;
00046     float angle ;
00047     float radangle ;
00048     float cd1_1, cd1_2, cd2_1, cd2_2 ;
00049     char firsttext[2*FILENAMESZ] ;
00050     char * date ;
00051 
00052     strcpy(firsttext, "spiffi_cube_create_from_resampled -f \0") ;
00053     
00054     pixelscale = spiffi_get_pixelscale(name) ;
00055     ra = spiffi_get_RA(name) ;
00056     dec = spiffi_get_DEC(name) ;
00057     angle = spiffi_get_posangle(name) ;
00058     radangle = angle * PI_NUMB / 180. ;
00059     cd1_1 = cos(radangle)  * pixelscale / 3600. ;
00060     cd1_2 = -sin(radangle) * pixelscale / 3600. ;
00061     cd2_1 = sin(radangle)  * pixelscale / 3600. ;
00062     cd2_2 = cos(radangle)  * pixelscale / 3600. ;
00063 
00064 
00065 
00066     cpl_propertylist_insert_after_string(plist, "EXPTIME", "CTYPE1", "RA---TAN" ) ;
00067     cpl_propertylist_set_comment(plist, "CTYPE1", "Projected Rectascension" ) ;
00068 
00069     cpl_propertylist_insert_after_float(plist, "CTYPE1",  "CRPIX1", center_x ) ;
00070     cpl_propertylist_set_comment(plist, "CRPIX1","Reference pixel in RA" ) ;
00071 
00072 
00073  
00074     cpl_propertylist_insert_after_float(plist, "CRPIX1",  "CRVAL1", ra ) ;
00075     cpl_propertylist_set_comment(plist, "CRVAL1","Reference RA" ) ;
00076 
00077     cpl_propertylist_insert_after_float(plist, "CRVAL1",  "CDELT1", -pixelscale/3600.  ) ;
00078     cpl_propertylist_set_comment(plist, "CDELT1","pixel scale" ) ;
00079 
00080     cpl_propertylist_insert_after_string(plist, "CDELT1",  "CUNIT1", "DEGREE" ) ;
00081     cpl_propertylist_set_comment(plist, "CUNIT1","RA-UNIT" ) ;
00082 
00083     cpl_propertylist_insert_after_string(plist, "CUNIT1",  "CTYPE2", "DEC--TAN" ) ;
00084     cpl_propertylist_set_comment(plist, "CTYPE2", "Projected Declination") ;
00085 
00086 
00087 
00088 
00089     cpl_propertylist_insert_after_float(plist, "CTYPE2",  "CRPIX2", center_y ) ;
00090     cpl_propertylist_set_comment(plist, "CRPIX2", "Reference pixel in DEC") ;
00091 
00092     cpl_propertylist_insert_after_float(plist, "CRPIX2",  "CRVAL2", dec) ;
00093     cpl_propertylist_set_comment(plist, "CRVAL2", "Reference DEC") ;
00094 
00095     cpl_propertylist_insert_after_float(plist, "CRVAL2",  "CDELT2", pixelscale/3600. ) ;
00096     cpl_propertylist_set_comment(plist, "CDELT2", "pixel scale") ;
00097 
00098     cpl_propertylist_insert_after_string(plist, "CDELT2",  "CUNIT2", "DEGREE" ) ;
00099     cpl_propertylist_set_comment(plist, "CUNIT2", "DEC-UNIT") ;
00100 
00101     cpl_propertylist_insert_after_string(plist, "CUNIT2",  "CTYPE3", "WAVE" ) ;
00102     cpl_propertylist_set_comment(plist, "CTYPE3", "wavelength axis in microns") ;
00103 
00104 
00105     cpl_propertylist_insert_after_int(plist, "CTYPE3",  "CRPIX3", cenpix+1 ) ;
00106     cpl_propertylist_set_comment(plist, "CRPIX3", "Reference pixel in z") ;
00107 
00108     cpl_propertylist_insert_after_float(plist, "CRPIX3",  "CRVAL3", cenLambda ) ;
00109     cpl_propertylist_set_comment(plist, "CRVAL3", "central wavelength") ;
00110 
00111     cpl_propertylist_insert_after_float(plist, "CRVAL3",  "CDELT3", dispersion ) ;
00112     cpl_propertylist_set_comment(plist, "CDELT3", "microns per pixel") ;
00113 
00114     cpl_propertylist_insert_after_string(plist, "CDELT3",  "CUNIT3", "MICRON" ) ;
00115     cpl_propertylist_set_comment(plist, "CUNIT3",  "spectral unit" ) ;
00116 
00117 
00118 
00119     cpl_propertylist_insert_after_float(plist, "CUNIT3",  "CD1_1", cd1_1 ) ;
00120     cpl_propertylist_set_comment(plist, "CD1_1",  "CD rotation matrix" ) ;
00121 
00122     cpl_propertylist_insert_after_float(plist, "CD1_1",   "CD1_2", cd1_2 ) ;
00123     cpl_propertylist_set_comment(plist, "CD1_2",  "CD rotation matrix" ) ;
00124 
00125     cpl_propertylist_insert_after_float(plist, "CD1_2",   "CD2_1", cd2_1 ) ;
00126     cpl_propertylist_set_comment(plist, "CD2_1",  "CD rotation matrix" ) ;
00127 
00128     cpl_propertylist_insert_after_float(plist, "CD2_1",   "CD2_2", cd2_2 ) ;
00129     cpl_propertylist_set_comment(plist, "CD2_2",  "CD rotation matrix" ) ;
00130 
00131     date     = get_datetime_iso8601() ;
00132     strcat (firsttext , ini_name) ;
00133 
00134 
00135     cpl_propertylist_set_string(plist,KEY_NAME_PIPEFILE,  outName ) ;
00136     cpl_propertylist_set_string(plist,KEY_NAME_HPRO_TYPE,  "REDUCED" ) ;
00137     cpl_propertylist_set_string(plist,KEY_NAME_HPRO_CATG,  "Science Observation" ) ;
00138     cpl_propertylist_set_string(plist,KEY_NAME_HPRO_STATUS, "OK" ) ;
00139     cpl_propertylist_set_string(plist,KEY_NAME_HPRO_DATE, date ) ;
00140     cpl_propertylist_set_string(plist,KEY_NAME_HPRO_RECID, firsttext ) ;
00141     cpl_propertylist_set_string(plist,KEY_NAME_HPRO_DRSID, PACKAGE_VERSION ) ;  
00142 
00143 
00144 }           
00145 
00146 
00147 /*----------------------------------------------------------------------------
00148    Function     :       cube_create_from_resampled()
00149    In           :       ini_file: file name of according .ini file
00150    Out          :       integer (0 if it worked, -1 if it doesn't) 
00151    Job          :       this routine does the resampling of an offset-corrected,
00152                         flatfielded, bad pixel corrected and
00153                         eventually interleaved data frame. Additionally, an intensity 
00154             calibration is carried through by using
00155                         a standard star or a black body measurement. 
00156             The spectral features of the flatfield halogen lamp are corrected.
00157                         Afterwards a data cube is created out of the resampled image.
00158             It is the users choice to use either
00159                         the fitted edge positions of the slitlets or the distances
00160             of the slitlets gained from a north-south-test. 
00161  ---------------------------------------------------------------------------*/
00162 int cube_create_from_resampled (cpl_parameterlist* config, cpl_frameset* sof)
00163 {
00164 
00165   const char* _id = "cube_create_from_resampled";
00166 
00167     cube_config * cfg ;
00168     OneImage * im ;
00169     OneCube  * cube ;
00170     OneCube  * outcube ;
00171     fits_header * head ;
00172     float fcol ;
00173     FILE * slitlist ;
00174     FILE * list ;
00175     FILE * fc_list ;
00176     float ** slit_edges ;
00177     float *  correct_dist ;
00178     float * distances ;
00179     int i ;
00180     cpl_frameset* raw=NULL;
00181     
00182     /*----parse the file names and parameters to the cube_config data structure cfg---*/
00183     cfg = parse_cpl_input_cube(config,sof,&raw) ;
00184 
00185 
00186     if (cfg == NULL)
00187     {
00188         cpl_msg_error(_id," could not parse .ini file!") ;
00189         return -1 ;
00190     }
00191     
00192     if(is_fits_file(cfg->inFrame) != 1 ){
00193       cpl_msg_error(_id,"Input file %s is not FITS",cfg->inFrame);
00194       return -1;
00195     }
00196 
00197     if (cfg->halocorrectInd == 1)
00198     {
00199       if(is_fits_file(cfg->halospectrum) != 1) {
00200       cpl_msg_error(_id,"Input file %s is not FITS",cfg->halospectrum);
00201       return -1;
00202       }
00203     }
00204     if (is_fits_file(cfg->poslist) != 1)
00205     {
00206          cpl_msg_error(_id,"Input file %s is not FITS",cfg->poslist) ;
00207      return -1 ;
00208     } 
00209     
00210     im = load_image(cfg->inFrame) ;
00211     if (im == NULL)
00212     {
00213         cpl_msg_error(_id, " could not load inFrame" ) ;
00214         return -1 ;
00215     }
00216     head = fits_read_header(cfg->inFrame) ;
00217     
00218     /*----------------------------------------------------------------------
00219      *---------------------CUBECREATION-------------------------------------
00220      *----------------------------------------------------------------------
00221      *---select north-south-test or fitting of slitlet edges---
00222      */
00223     if (NULL == (correct_dist = (float*) cpl_calloc(cfg->nslits, sizeof (float))) )
00224     {
00225         cpl_msg_error(_id," could not allocate memory!") ;
00226     return -1 ;   
00227     }
00228     
00229     if (cfg->northsouthInd == 0)
00230     {
00231     if (NULL == (slit_edges = (float**) cpl_calloc (cfg->nslits, sizeof(float*)) ) )
00232     {
00233         cpl_msg_error(_id," could not allocate memory!") ;
00234         return -1 ;   
00235     }
00236     for ( i = 0 ; i < cfg->nslits ; i++ )
00237     {
00238         if (NULL == (slit_edges[i] = (float*) cpl_calloc (2, sizeof(float))))
00239         {
00240             cpl_msg_error(_id," could not allocate memory!") ;
00241             return -1 ;   
00242         }
00243     }
00244    
00245         if ( NULL == (slitlist = fopen(cfg->poslist,  "r" )) )
00246     {
00247         cpl_msg_error(_id," could not open poslist!") ;
00248         return -1 ;
00249     }
00250         
00251     i = 0 ;
00252     while (EOF != fscanf(slitlist, "%f %f\n", &slit_edges[i][0] , &slit_edges[i][1]) )
00253     {
00254         i++ ;
00255     } 
00256     if ( i != cfg->nslits)
00257     {
00258         cpl_msg_error(_id," wrong file format!") ;
00259         fclose(slitlist) ;
00260         return -1 ;
00261     }
00262         fclose(slitlist) ;
00263 
00264         /* build the data cube */
00265         cube = makeCubeSpi( im, slit_edges, correct_dist ) ;
00266         if (cube == NULL)
00267     {
00268             cpl_msg_error(_id," could not construct data cube!") ;
00269             return -1 ;
00270     }
00271         
00272     for ( i = 0 ; i < cfg->nslits ; i++ )
00273     {
00274         cpl_free (slit_edges[i]) ;
00275     }
00276         cpl_free (slit_edges) ;
00277     }
00278     else
00279     {
00280         if ( NULL == (distances = (float*) cpl_calloc (cfg->nslits - 1, sizeof (float))) )
00281     {
00282         cpl_msg_error(_id," could allocate memory!") ;
00283             return -1 ;
00284     } 
00285 
00286         if ( NULL == (list = fopen(cfg->poslist,  "r" ) ) ) 
00287     {
00288         cpl_msg_error(_id," could not open poslist!") ;
00289         return -1 ;
00290     }
00291         
00292     i = 0 ;
00293     while (EOF != fscanf(list, "%f\n", &distances[i] ))
00294     {
00295         i++ ;
00296     } 
00297         
00298     if ( i != cfg->nslits - 1)
00299     {
00300         cpl_msg_error(_id," wrong file format!") ;
00301         fclose(list) ;
00302         return -1 ;
00303     }
00304         
00305     fclose(list) ;
00306 
00307         if ( NULL == (fc_list = fopen(cfg->firstCol,  "r" ) ) ) 
00308     {
00309         cpl_msg_error(_id," could not open firstCol list!") ;
00310         return -1 ;
00311     }
00312         fscanf(fc_list, "%f\n", &fcol) ; 
00313         fclose(fc_list) ;
00314     
00315     /*build the data cube */
00316         cube = makeCubeDist( im, fcol, distances, correct_dist ) ;
00317         if (cube == NULL)
00318     {
00319             cpl_msg_error(_id," could not construct a data cube!") ;
00320             return -1 ;
00321         }
00322         cpl_free(distances) ;
00323     }
00324     
00325     /*-------------------------------------------------------------------------------------
00326      *------------------------FINETUNING------------------------------------------------------
00327      *-------------------------------------------------------------------------------------
00328      * shift the rows of the reconstructed images of the data cube to the correct sub pixel position
00329      * select the shift method: polynomial interpolation, FFT or cubic spline interpolation
00330      */
00331      
00332     if (cfg->method == 'P')
00333     {
00334         outcube = fineTuneCube( cube, correct_dist, cfg->order ) ;
00335         if (outcube == NULL)
00336     {
00337             cpl_msg_error(_id," could not fine tune the data cube") ;
00338             return -1 ;
00339         }
00340     }       
00341     else if (cfg->method == 'F')
00342     {
00343         for ( i = 0 ; i < cfg->nslits ; i++ )
00344     {
00345          correct_dist[i] = -correct_dist[i] ;
00346     }
00347         outcube = fineTuneCubeByFFT( cube, correct_dist ) ;
00348         if ( outcube == NULL )
00349         {
00350         cpl_msg_error(_id," could not fine tune the data cube") ;
00351             return -1 ;
00352         }
00353     }       
00354     else if (cfg->method == 'S')
00355     {
00356         outcube = fineTuneCubeBySpline( cube, correct_dist ) ;
00357         if ( outcube == NULL )
00358     {
00359             cpl_msg_error(_id," could not fine tune the data cube") ;
00360             return -1 ;
00361         }
00362     }
00363     else
00364     {
00365         cpl_msg_error(_id," wrong method indicator given!") ;
00366         return -1 ;
00367     }
00368     
00369     save_cube_to_fits_hdr_dump( outcube, cfg -> outName, head ) ;
00370 
00371     /*---free memory----*/
00372     fits_header_destroy(head) ;
00373     cpl_free(correct_dist) ;
00374     destroy_image(im) ;
00375     destroy_cube (cube) ;
00376     destroy_cube (outcube) ;
00377     cube_cfg_destroy (cfg) ;
00378     return 0 ;
00379 }
00380 
00381 

Generated on Wed Oct 26 13:08:51 2005 for SINFONI Pipeline Reference Manual by doxygen1.2.13.1 written by Dimitri van Heesch, © 1997-2001