utilities_scired.c

00001 #include "utilities_scired.h"
00002 /*---------------------------------------------------------------------------*/
00025 /*---------------------------------------------------------------------------*/
00026 OneCube * cube_getvig(
00027         OneCube *        cube_in,
00028         int             loleft_x,
00029         int             loleft_y,
00030         int             upright_x,
00031         int             upright_y)
00032 {
00033   const char * fctid="cube_getvig";
00034     OneCube     *        cube_out ;
00035     int                 i ;
00036     int                 outlx,
00037                                         outly ;
00038 
00039     if (cube_in==NULL) return NULL ;
00040 
00041         if ((loleft_x>upright_x) ||
00042             (loleft_y>upright_y)) {
00043             cpl_msg_error(fctid,"ill-defined slit for extraction: aborting");
00044             return NULL ;
00045         }
00046 
00047     /* Extraction coordinates include rectangular zone  */
00048     outlx = upright_x - loleft_x + 1 ;
00049     outly = upright_y - loleft_y + 1 ;
00050 
00051     cube_out = new_cube(outlx, outly, cube_in->np) ;
00052     /* Loop on all input planes */
00053     for (i=0 ; i<cube_in->np ; i++) {
00054       /*
00055         cpl_msg_warning(fctid,"extracting subimage %d %d %d", 
00056                        i, cube_in->np, 2) ;
00057       */
00058         /* Extract a slit from this plane   */
00059                 cube_out->plane[i] =
00060                         image_getvig(cube_in->plane[i],
00061                         loleft_x, loleft_y,
00062                         upright_x, upright_y) ;
00063     }
00064     return cube_out ;
00065 }
00066 /*---------------------------------------------------------------------------*/
00083 /*---------------------------------------------------------------------------*/
00084 
00085 OneImage * image_getvig(
00086         OneImage    *    image_in,
00087         int             loleft_x,
00088         int             loleft_y,
00089         int             upright_x,
00090         int             upright_y)
00091 {
00092     OneImage    *        slit_img ;
00093     int                 i, j ;
00094     register
00095     pixelvalue  *       inpt,
00096                                 *       outpt ;
00097     int         outlx, outly ;
00098     const char* fctid="image_getvig";
00099     if (image_in==NULL) return NULL ;
00100 
00101     if ((loleft_x<1) || (loleft_x>image_in->lx) ||
00102         (loleft_y<1) || (loleft_y>image_in->ly) ||
00103         (upright_x<1) || (upright_x>image_in->lx) ||
00104         (upright_y<1) || (upright_y>image_in->ly) ||
00105         (loleft_x>upright_x) || (loleft_y>upright_y)) {
00106         cpl_msg_error(fctid,"extraction zone is [%d %d] [%d %d]\n"
00107                         "cannot extract such zone: aborting slit extraction",
00108                 loleft_x, loleft_y, upright_x, upright_y) ;
00109         return NULL ;
00110     }
00111 
00112     outlx = upright_x - loleft_x + 1 ;
00113     outly = upright_y - loleft_y + 1 ;
00114     slit_img = new_image(outlx, outly) ;
00115 
00116     for (j=0 ; j<outly ; j++) {
00117         inpt = image_in->data+loleft_x-1 + (j+loleft_y-1)*image_in->lx ;
00118         outpt = slit_img->data + j*outlx ;
00119         for (i=0 ; i<outlx ; i++) {
00120             *outpt++ = *inpt++ ;
00121         }
00122     }
00123     return slit_img ;
00124 }
00125 
00126 
00127 int
00128 set_wcs_cube(cpl_imagelist* cub, char* name, double clambda, 
00129          double dis, double cpix, double cx, double cy)
00130 {
00131   const char* _id="set_wcs_cube";
00132   cpl_propertylist* plist=NULL;
00133 
00134   if ((cpl_error_code)((plist = cpl_propertylist_load(name, 0)) == NULL)) {
00135     cpl_msg_error(_id, "getting header from frame %s",name);
00136     cpl_propertylist_delete(plist) ;
00137     return -1 ;
00138   }
00139 
00140   change_plist_cube(plist, name, clambda, dis, cpix, cx, cy) ;
00141         
00142   if (cpl_imagelist_save(cub, name, CPL_BPP_DEFAULT,plist,CPL_IO_DEFAULT)!=CPL_ERROR_NONE) {
00143     cpl_msg_error(_id, "Cannot save the product %s",name);
00144     cpl_propertylist_delete(plist) ;
00145     return -1 ;
00146   }
00147   cpl_propertylist_delete(plist) ;
00148   return 0;
00149 
00150 }
00151 
00152 int
00153 set_wcs_image(cpl_image* img, char* name,  double cx, double cy)
00154 {
00155   const char* _id="set_wcs_image";
00156   cpl_propertylist* plist=NULL;
00157 
00158   if ((cpl_error_code)((plist = cpl_propertylist_load(name, 0)) == NULL)) {
00159     cpl_msg_error(_id, "getting header from frame %s",name);
00160     cpl_propertylist_delete(plist) ;
00161     return -1 ;
00162   }
00163   change_plist_image(plist, name, cx, cy) ;
00164         
00165   if (cpl_image_save(img, name, CPL_BPP_DEFAULT,plist,CPL_IO_DEFAULT)!=CPL_ERROR_NONE) {
00166     cpl_msg_error(_id, "Cannot save the product %s",name);
00167     cpl_propertylist_delete(plist) ;
00168     return -1 ;
00169   }
00170   cpl_propertylist_delete(plist) ;
00171   return 0;
00172 
00173 }
00174 
00175 
00176 
00177 int
00178 set_wcs_spectrum(cpl_image* img, char* name, double clambda, 
00179              double dis, double cpix)
00180 {
00181   const char* _id="set_wcs_image";
00182   cpl_propertylist* plist=NULL;
00183 
00184   if ((cpl_error_code)((plist = cpl_propertylist_load(name, 0)) == NULL)) {
00185     cpl_msg_error(_id, "getting header from frame %s",name);
00186     cpl_propertylist_delete(plist) ;
00187     return -1 ;
00188   }
00189   change_plist_spectrum(plist, clambda, dis,cpix) ;
00190         
00191   if (cpl_image_save(img, name, CPL_BPP_DEFAULT,plist,CPL_IO_DEFAULT)!=CPL_ERROR_NONE) {
00192     cpl_msg_error(_id, "Cannot save the product %s",name);
00193     cpl_propertylist_delete(plist) ;
00194     return -1 ;
00195   }
00196   cpl_propertylist_delete(plist) ;
00197   return 0;
00198 
00199 }
00200 
00201 void change_plist_cube (cpl_propertylist * plist, 
00202                     char * name,
00203             float cenLambda,
00204             float dispersion,
00205             int   cenpix,
00206             float center_x,
00207             float center_y )
00208 {
00209 
00210     float pixelscale ;
00211     float ra ;
00212     float dec ;
00213     float angle ;
00214     float radangle ;
00215     float cd1_1, cd1_2, cd2_1, cd2_2 ;
00216     char firsttext[2*FILENAMESZ] ;
00217  
00218     strcpy(firsttext, "spiffi_objnod -f \0") ;
00219     
00220     pixelscale = spiffi_get_pixelscale(name)/2. ;
00221     ra = spiffi_get_RA(name) ;
00222     dec = spiffi_get_DEC(name) ;
00223     angle = spiffi_get_posangle(name) ;
00224     radangle = angle * PI_NUMB / 180. ;
00225     cd1_1 = cos(radangle)  * pixelscale / 3600. ;
00226     cd1_2 = sin(radangle) * pixelscale / 3600. ;
00227     cd2_1 = -sin(radangle)  * pixelscale / 3600. ;
00228     cd2_2 = cos(radangle)  * pixelscale / 3600. ;
00229     cpl_propertylist_erase_regexp(plist, "^CTYPE1");
00230     cpl_propertylist_insert_after_string(plist, "EXPTIME", "CTYPE1", "RA---TAN" ) ;
00231     cpl_propertylist_set_comment(plist, "CTYPE1", "Projected Rectascension" ) ;
00232 
00233     cpl_propertylist_erase_regexp(plist, "^CRPIX1");
00234     cpl_propertylist_insert_after_float(plist, "CTYPE1",  "CRPIX1", center_x ) ;
00235     cpl_propertylist_set_comment(plist, "CRPIX1","Reference pixel in RA" ) ;
00236 
00237     cpl_propertylist_erase_regexp(plist, "^CRVAL1");
00238     cpl_propertylist_insert_after_float(plist, "CRPIX1",  "CRVAL1", ra ) ;
00239     cpl_propertylist_set_comment(plist, "CRVAL1","Reference RA" ) ;
00240 
00241     cpl_propertylist_erase_regexp(plist, "^CDELT1");
00242     cpl_propertylist_insert_after_float(plist, "CRVAL1",  "CDELT1", -pixelscale/3600.  ) ;
00243     cpl_propertylist_set_comment(plist, "CDELT1","pixel scale" ) ;
00244 
00245     cpl_propertylist_erase_regexp(plist, "^CUNIT1");
00246     cpl_propertylist_insert_after_string(plist, "CDELT1",  "CUNIT1", "DEGREE" ) ;
00247     cpl_propertylist_set_comment(plist, "CUNIT1","RA-UNIT" ) ;
00248 
00249     cpl_propertylist_erase_regexp(plist, "^CTYPE2");
00250     cpl_propertylist_insert_after_string(plist, "CUNIT1",  "CTYPE2", "DEC--TAN" ) ;
00251     cpl_propertylist_set_comment(plist, "CTYPE2", "Projected Declination") ;
00252 
00253     cpl_propertylist_erase_regexp(plist, "^CRPIX2");
00254     cpl_propertylist_insert_after_float(plist, "CTYPE2",  "CRPIX2", center_y ) ;
00255     cpl_propertylist_set_comment(plist, "CRPIX2", "Reference pixel in DEC") ;
00256 
00257     cpl_propertylist_erase_regexp(plist, "^CRVAL2");
00258     cpl_propertylist_insert_after_float(plist, "CRPIX2",  "CRVAL2", dec) ;
00259     cpl_propertylist_set_comment(plist, "CRVAL2", "Reference DEC") ;
00260 
00261     cpl_propertylist_erase_regexp(plist, "^CDELT2");
00262     cpl_propertylist_insert_after_float(plist, "CRVAL2",  "CDELT2", pixelscale/3600. ) ;
00263     cpl_propertylist_set_comment(plist, "CDELT2", "pixel scale") ;
00264 
00265     cpl_propertylist_erase_regexp(plist, "^CUNIT2");
00266     cpl_propertylist_insert_after_string(plist, "CDELT2",  "CUNIT2", "DEGREE" ) ;
00267     cpl_propertylist_set_comment(plist, "CUNIT2", "DEC-UNIT") ;
00268 
00269     cpl_propertylist_erase_regexp(plist, "^CTYPE3");
00270     cpl_propertylist_insert_after_string(plist, "CUNIT2",  "CTYPE3", "WAVE" ) ;
00271     cpl_propertylist_set_comment(plist, "CTYPE3", "wavelength axis in microns") ;
00272 
00273     cpl_propertylist_erase_regexp(plist, "^CRPIX3");
00274     cpl_propertylist_insert_after_int(plist, "CTYPE3",  "CRPIX3", cenpix+1 ) ;
00275     cpl_propertylist_set_comment(plist, "CRPIX3", "Reference pixel in z") ;
00276 
00277     cpl_propertylist_erase_regexp(plist, "^CRVAL3");
00278     cpl_propertylist_insert_after_float(plist, "CRPIX3",  "CRVAL3", cenLambda ) ;
00279     cpl_propertylist_set_comment(plist, "CRVAL3", "central wavelength") ;
00280 
00281     cpl_propertylist_erase_regexp(plist, "^CDELT3");
00282     cpl_propertylist_insert_after_float(plist, "CRVAL3",  "CDELT3", dispersion ) ;
00283     cpl_propertylist_set_comment(plist, "CDELT3", "microns per pixel") ;
00284 
00285     cpl_propertylist_erase_regexp(plist, "^CUNIT3");
00286     cpl_propertylist_insert_after_string(plist, "CDELT3",  "CUNIT3", "MICRON" ) ;
00287     cpl_propertylist_set_comment(plist, "CUNIT3",  "spectral unit" ) ;
00288     
00289     cpl_propertylist_erase_regexp(plist, "^CD1_1");
00290     cpl_propertylist_insert_after_float(plist, "CUNIT3",  "CD1_1", cd1_1 ) ;
00291     cpl_propertylist_set_comment(plist, "CD1_1",  "CD rotation matrix" ) ;
00292 
00293     cpl_propertylist_erase_regexp(plist, "^CD1_2");
00294     cpl_propertylist_insert_after_float(plist, "CD1_1",   "CD1_2", cd1_2 ) ;
00295     cpl_propertylist_set_comment(plist, "CD1_2",  "CD rotation matrix" ) ;
00296     
00297     cpl_propertylist_erase_regexp(plist, "^CD2_1");
00298     cpl_propertylist_insert_after_float(plist, "CD1_2",   "CD2_1", cd2_1 ) ;
00299     cpl_propertylist_set_comment(plist, "CD2_1",  "CD rotation matrix" ) ;
00300     
00301     cpl_propertylist_erase_regexp(plist, "^CD2_2");
00302     cpl_propertylist_insert_after_float(plist, "CD2_1",   "CD2_2", cd2_2 ) ;
00303     cpl_propertylist_set_comment(plist, "CD2_2",  "CD rotation matrix" ) ;
00304 
00305 
00306 }        
00307 
00308 
00309 void change_plist_spectrum (cpl_propertylist * plist, 
00310                 double cenLambda,
00311                 double dispersion,
00312                 int   cenpix)
00313 {
00314 
00315     cpl_propertylist_erase_regexp(plist, "^CTYPE3");
00316     cpl_propertylist_erase_regexp(plist, "^CRPIX3");
00317     cpl_propertylist_erase_regexp(plist, "^CRVAL3");
00318     cpl_propertylist_erase_regexp(plist, "^CDELT3");
00319     cpl_propertylist_erase_regexp(plist, "^CUNIT3");
00320 
00321     cpl_propertylist_erase_regexp(plist, "^CTYPE2");
00322     cpl_propertylist_erase_regexp(plist, "^CRPIX2");
00323     cpl_propertylist_erase_regexp(plist, "^CRVAL2");
00324     cpl_propertylist_erase_regexp(plist, "^CDELT2");
00325     cpl_propertylist_erase_regexp(plist, "^CUNIT2");
00326 
00327 
00328     cpl_propertylist_erase_regexp(plist, "^CD1_1");
00329     cpl_propertylist_erase_regexp(plist, "^CD1_2");
00330     cpl_propertylist_erase_regexp(plist, "^CD2_1");
00331     cpl_propertylist_erase_regexp(plist, "^CD2_2");
00332 
00333 
00334     cpl_propertylist_erase_regexp(plist, "^CTYPE1");
00335     cpl_propertylist_insert_after_string(plist, "EXPTIME",  "CTYPE1", "PIXEL");
00336     cpl_propertylist_set_comment(plist, "CTYPE1", "Pixel coordinate system.");
00337 
00338     cpl_propertylist_erase_regexp(plist, "^CRPIX1");
00339     cpl_propertylist_insert_after_int(plist, "CTYPE1",  "CRPIX1", 1 ) ;
00340     cpl_propertylist_set_comment(plist, "CRPIX1", "Reference pixel in x") ;
00341 
00342     cpl_propertylist_erase_regexp(plist, "^CRVAL1");
00343     cpl_propertylist_insert_after_double(plist, "CRPIX1",  "CRVAL1", 1 ) ;
00344     cpl_propertylist_set_comment(plist, "CRVAL1", "value of ref pixel.") ;
00345 
00346     cpl_propertylist_erase_regexp(plist, "^CDELT1");
00347     cpl_propertylist_insert_after_double(plist, "CRVAL1",  "CDELT1", 1 ) ;
00348     cpl_propertylist_set_comment(plist, "CDELT1", "pixel scale") ;
00349 
00350 
00351 
00352 
00353 
00354     cpl_propertylist_erase_regexp(plist, "^CUNIT1");
00355     cpl_propertylist_insert_after_string(plist, "CDELT1",  "CUNIT1", "Pixel" );
00356     cpl_propertylist_set_comment(plist, "CUNIT1",  "spectral unit" );
00357 
00358 
00359     cpl_propertylist_erase_regexp(plist, "^CTYPE2");
00360     cpl_propertylist_insert_after_string(plist, "EXPTIME","CTYPE2","WAVE" );
00361     cpl_propertylist_set_comment(plist, "CTYPE2","wavelength axis in microns");
00362 
00363 
00364     cpl_propertylist_erase_regexp(plist, "^CRPIX2");
00365     cpl_propertylist_insert_after_int(plist, "CTYPE2",  "CRPIX2",cenpix+1 ) ;
00366     cpl_propertylist_set_comment(plist, "CRPIX2", "Reference pixel in x") ;
00367 
00368     cpl_propertylist_erase_regexp(plist, "^CRVAL2");
00369     cpl_propertylist_insert_after_double(plist, "CRPIX2","CRVAL2",cenLambda ) ;
00370     cpl_propertylist_set_comment(plist, "CRVAL2", "central wavelength") ;
00371 
00372     cpl_propertylist_erase_regexp(plist, "^CDELT2");
00373     cpl_propertylist_insert_after_double(plist, "CRVAL2", "CDELT2",dispersion);
00374     cpl_propertylist_set_comment(plist, "CDELT2", "microns per pixel");
00375 
00376     cpl_propertylist_erase_regexp(plist, "^CUNIT2");
00377     cpl_propertylist_insert_after_string(plist, "CDELT2",  "CUNIT2", "MICRON");
00378     cpl_propertylist_set_comment(plist, "CUNIT2",  "spectral unit" );
00379 
00380     
00381 
00382 }   
00383 
00384 
00385 void change_plist_image (cpl_propertylist * plist, 
00386                     char * name,
00387             float center_x,
00388             float center_y )
00389 {
00390 
00391     float pixelscale ;
00392     float ra ;
00393     float dec ;
00394     float angle ;
00395     float radangle ;
00396     float cd1_1, cd1_2, cd2_1, cd2_2 ;
00397     char firsttext[2*FILENAMESZ] ;
00398  
00399     strcpy(firsttext, "spiffi_objnod -f \0") ;
00400     
00401     pixelscale = spiffi_get_pixelscale(name)/2. ;
00402     ra = spiffi_get_RA(name) ;
00403     dec = spiffi_get_DEC(name) ;
00404     angle = spiffi_get_posangle(name) ;
00405     radangle = angle * PI_NUMB / 180. ;
00406     cd1_1 = cos(radangle)  * pixelscale / 3600. ;
00407     cd1_2 = sin(radangle) * pixelscale / 3600. ;
00408     cd2_1 = -sin(radangle)  * pixelscale / 3600. ;
00409     cd2_2 = cos(radangle)  * pixelscale / 3600. ;
00410 
00411     cpl_propertylist_erase_regexp(plist, "^CTYPE1");
00412     cpl_propertylist_insert_after_string(plist, "EXPTIME", "CTYPE1", "RA---TAN" ) ;
00413     cpl_propertylist_set_comment(plist, "CTYPE1", "Projected Rectascension" ) ;
00414 
00415     cpl_propertylist_erase_regexp(plist, "^CRPIX1");
00416     cpl_propertylist_insert_after_float(plist, "CTYPE1",  "CRPIX1", center_x ) ;
00417     cpl_propertylist_set_comment(plist, "CRPIX1","Reference pixel in RA" ) ;
00418 
00419     cpl_propertylist_erase_regexp(plist, "^CRVAL1");
00420     cpl_propertylist_insert_after_float(plist, "CRPIX1",  "CRVAL1", ra ) ;
00421     cpl_propertylist_set_comment(plist, "CRVAL1","Reference RA" ) ;
00422 
00423     cpl_propertylist_erase_regexp(plist, "^CDELT1");
00424     cpl_propertylist_insert_after_float(plist, "CRVAL1",  "CDELT1", -pixelscale/3600.  ) ;
00425     cpl_propertylist_set_comment(plist, "CDELT1","pixel scale" ) ;
00426 
00427     cpl_propertylist_erase_regexp(plist, "^CUNIT1");
00428     cpl_propertylist_insert_after_string(plist, "CDELT1",  "CUNIT1", "DEGREE" ) ;
00429     cpl_propertylist_set_comment(plist, "CUNIT1","RA-UNIT" ) ;
00430 
00431     cpl_propertylist_erase_regexp(plist, "^CTYPE2");
00432     cpl_propertylist_insert_after_string(plist, "CUNIT1",  "CTYPE2", "DEC--TAN" ) ;
00433     cpl_propertylist_set_comment(plist, "CTYPE2", "Projected Declination") ;
00434 
00435     cpl_propertylist_erase_regexp(plist, "^CRPIX2");
00436     cpl_propertylist_insert_after_float(plist, "CTYPE2",  "CRPIX2", center_y ) ;
00437     cpl_propertylist_set_comment(plist, "CRPIX2", "Reference pixel in DEC") ;
00438 
00439     cpl_propertylist_erase_regexp(plist, "^CRVAL2");
00440     cpl_propertylist_insert_after_float(plist, "CRPIX2",  "CRVAL2", dec) ;
00441     cpl_propertylist_set_comment(plist, "CRVAL2", "Reference DEC") ;
00442 
00443     cpl_propertylist_erase_regexp(plist, "^CDELT2");
00444     cpl_propertylist_insert_after_float(plist, "CRVAL2",  "CDELT2", pixelscale/3600. ) ;
00445     cpl_propertylist_set_comment(plist, "CDELT2", "pixel scale") ;
00446 
00447     cpl_propertylist_erase_regexp(plist, "^CUNIT2");
00448     cpl_propertylist_insert_after_string(plist, "CDELT2",  "CUNIT2", "DEGREE" ) ;
00449     cpl_propertylist_set_comment(plist, "CUNIT2", "DEC-UNIT") ;
00450     
00451     cpl_propertylist_erase_regexp(plist, "^CD1_1");
00452     cpl_propertylist_insert_after_float(plist, "CUNIT2",  "CD1_1", cd1_1 ) ;
00453     cpl_propertylist_set_comment(plist, "CD1_1",  "CD rotation matrix" ) ;
00454 
00455     cpl_propertylist_erase_regexp(plist, "^CD1_2");
00456     cpl_propertylist_insert_after_float(plist, "CD1_1",   "CD1_2", cd1_2 ) ;
00457     cpl_propertylist_set_comment(plist, "CD1_2",  "CD rotation matrix" ) ;
00458     
00459     cpl_propertylist_erase_regexp(plist, "^CD2_1");
00460     cpl_propertylist_insert_after_float(plist, "CD1_2",   "CD2_1", cd2_1 ) ;
00461     cpl_propertylist_set_comment(plist, "CD2_1",  "CD rotation matrix" ) ;
00462     
00463     cpl_propertylist_erase_regexp(plist, "^CD2_2");
00464     cpl_propertylist_insert_after_float(plist, "CD2_1",   "CD2_2", cd2_2 ) ;
00465     cpl_propertylist_set_comment(plist, "CD2_2",  "CD rotation matrix" ) ;
00466 
00467 
00468 }        
00469 
00470 
00471 cpl_imagelist* 
00472 sinfoni_cube2imglist(OneCube* cube)
00473 {
00474   cpl_imagelist* ims=NULL;
00475   int i=0;
00476   OneImage* eimg=NULL;
00477   cpl_image* cimg=NULL;
00478   cpl_image* wimg=NULL;
00479   ims=cpl_imagelist_new(cube->lx,cube->ly,cube->np,CPL_TYPE_FLOAT);
00480  
00481   for(i=0;i<cube->np;i++) {
00482     eimg=cube->plane[i];
00483     wimg=cpl_image_wrap_float(eimg->lx,eimg->ly,eimg->data);
00484     cimg=cpl_image_duplicate(wimg);
00485     cpl_image_unwrap(wimg);
00486     cpl_imagelist_set(ims,cimg,i);
00487   }
00488 
00489   return ims;
00490 }
00491 
00492 OneCube** sinfoni_correct_median(OneCube** cubes, const int n_cubes)
00493 {
00494   const char* fctid="sinfoni_correct_median";
00495   int i=0;
00496   OneCube** cubes_cor=NULL;
00497   double median=0;
00498   pixelvalue pv=0;
00499   int z=0;
00500   if ( cubes == NULL ) {
00501     cpl_msg_error (fctid,"no cube list given!") ;
00502     return NULL ;
00503   }
00504   if ( n_cubes <= 0 ) {
00505     cpl_msg_error (fctid,"wrong number of data cubes in list!") ;
00506     return NULL ;
00507   }
00508 
00509   cubes_cor = (OneCube**) cpl_calloc (n_cubes, sizeof (OneCube*)) ;     
00510   
00511   for ( i = 0 ; i < n_cubes ; i++ ) {
00512     cubes_cor[i] = newCube(cubes[i]->lx,cubes[i]->ly,cubes[i]->np);
00513     for(z=0;z< cubes[i]->np; z++) {
00514       pv=get_median_pixelvalue_image(cubes[i]->plane[z]);
00515       median=(double)pv;
00516       if(!qfits_isnan(median)) {
00517         cubes_cor[i]->plane[z]=cst_op_image(cubes[i]->plane[z],median,'-');
00518       }
00519     }
00520   }
00521 
00522   return cubes_cor;
00523 }
00524 
00525 int sinfoni_correct_median_it(OneCube** inp)
00526 {
00527 
00528   const char* fctid="sinfoni_correct_median";
00529   double median=0;
00530   int z=0;
00531 
00532   for(z=0;z< (*inp)->np; z++) {
00533     median=median_image((*inp)->plane[z]);
00534     if(!qfits_isnan(median)) {
00535       cst_op_image_local((*inp)->plane[z],median,'-');
00536     }  else {
00537       cpl_msg_error(fctid,"median is NAN");
00538     }
00539   }
00540 
00541   return 0;
00542 }
00543 
00544 
00545 
00546 OneCube** sinfoni_correct_sky(OneCube** cubes,
00547                   const int nc,
00548                   OneCube* sky_cube)
00549 
00550 {
00551   const char* fctid="sinfoni_correct_sky";
00552   OneCube** cubes_sky=NULL;
00553   int x=0;
00554   int y=0;
00555   int z=0;
00556   int i=0;
00557   float k=0.5;
00558   int ovr=0;
00559   int ks=0;
00560   int nclip=0;
00561   double med=0;
00562   double avg=0;
00563   double sig=0;
00564   int msk_sum=0;
00565   double val_msk_sum=0;
00566   cpl_vector* val=NULL;
00567   cpl_vector* msk=NULL;
00568   
00569 
00570   if ( cubes == NULL ) {
00571     cpl_msg_error (fctid,"no cube list given!") ;
00572     return NULL ;
00573   }
00574 
00575   if ( nc <= 0 ) {
00576     cpl_msg_error (fctid,"wrong number of data cubes in list!") ;
00577     return NULL ;
00578   }
00579 
00580   cubes_sky = (OneCube**) cpl_calloc (nc, sizeof (OneCube*)) ;     
00581   for(z=0;z< cubes[0]->np; z++) {
00582     for(y=0;y< cubes[0]->ly; y++) {
00583       for(x=0;x< cubes[0]->lx; x++) {
00584         /* here we start to do a k-s clipping */
00585     msk=cpl_vector_new(nc);
00586     for (i=0;i<nc;i++) {
00587       cpl_vector_set(msk,i,1);
00588     }
00589         nclip=0;
00590         for (ks=0;ks<nc;ks++) {
00591       sig=0;
00592       med=0;
00593       ovr=0;
00594       val=cpl_vector_new(nc-nclip);
00595 
00596       for ( i = 0 ; i < nc ; i++ ) {
00597         if ((!qfits_isnan(cubes[i]->plane[z]->data[x+y*cubes[i]->lx])) &&
00598                  (cpl_vector_get(msk,i) != 0) ) {
00599           cpl_vector_set(val,ovr,(double)cubes[i]->plane[z]->data[x+y*cubes[i]->lx]); 
00600           ovr++;
00601         }
00602       }
00603 
00604       if(ovr>0) {
00605         avg=cpl_vector_get_mean(val);
00606         med=cpl_vector_get_median(val);
00607         if(ovr>1) {
00608           sig=cpl_vector_get_stdev(val);
00609         } else {
00610               sig=0;
00611         }
00612       } else {
00613         avg=cpl_vector_get(val,0);
00614         med=avg;
00615         sig=0;
00616       }
00617 
00618       cpl_vector_delete(val);
00619       for ( i = 0 ; i < nc ; i++ ) {
00620       /* Do k-s clipping at each pixel */
00621         if ((!qfits_isnan(cubes[i]->plane[z]->data[x+y*cubes[i]->lx])) && 
00622                 (cpl_vector_get(msk,i) != 0)) { 
00623           if(abs((cubes[i]->plane[z]->data[x+y*cubes[i]->lx]-med))> k*sig) {
00624         cpl_vector_set(msk,i,0);
00625                 nclip++;
00626           }
00627         }
00628       }
00629     }/* end of k-s clipping */
00630     msk_sum=0;
00631     val_msk_sum=0;
00632     for ( i = 0 ; i < nc ; i++ ) {
00633       /* computes sky at each point */
00634       if (!qfits_isnan(cubes[i]->plane[z]->data[x+y*cubes[i]->lx])) { 
00635         msk_sum+=cpl_vector_get(msk,i);
00636                 val_msk_sum+=cubes[i]->plane[z]->data[x+y*cubes[i]->lx]*
00637           cpl_vector_get(msk,i);
00638       }
00639     }
00640     sky_cube->plane[z]->data[x+y*sky_cube->lx]=val_msk_sum/msk_sum;
00641     cpl_vector_delete(msk);
00642       } /* end loop over x */
00643     } /* end loop over y */
00644   } /* end loop over z */
00645   for ( i = 0 ; i < nc ; i++ ) {
00646     cubes_sky[i]=copy_cube(cubes[i]);
00647     /* subtract the variable clean sky */
00648     op_cube(&(cubes_sky[i]),sky_cube,'-');
00649 
00650   }
00651 
00652 
00653   return cubes_sky;
00654 }
00655 
00656 
00657 OneCube** sinfoni_correct_sky2(OneCube** cubes,
00658                   const int nc,
00659                   OneCube* sky_cube,
00660                               OneCube* med_cube,
00661                               OneCube* msk_cube,
00662                               OneCube* avg_cube,
00663                               OneCube* sig_cube,
00664                               OneCube* ovr_cube)
00665 {
00666   const char* fctid="sinfoni_correct_sky";
00667   OneCube** cubes_sky=NULL;
00668   int x=0;
00669   int y=0;
00670   int z=0;
00671   int i=0;
00672   float k=0.5;
00673   int ovr=0;
00674   int ks=0;
00675   int nclip=0;
00676   double med=0;
00677   double avg=0;
00678   double sig=0;
00679   int msk_sum=0;
00680   double val_msk_sum=0;
00681   cpl_vector* val=NULL;
00682   cpl_vector* msk=NULL;
00683   
00684 
00685   if ( cubes == NULL ) {
00686     cpl_msg_error (fctid,"no cube list given!") ;
00687     return NULL ;
00688   }
00689 
00690   if ( nc <= 0 ) {
00691     cpl_msg_error (fctid,"wrong number of data cubes in list!") ;
00692     return NULL ;
00693   }
00694 
00695   cubes_sky = (OneCube**) cpl_calloc (nc, sizeof (OneCube*)) ;     
00696   for(z=0;z< cubes[0]->np; z++) {
00697     for(y=0;y< cubes[0]->ly; y++) {
00698       for(x=0;x< cubes[0]->lx; x++) {
00699         /* here we start to do a k-s clipping */
00700     msk=cpl_vector_new(nc);
00701     for (i=0;i<nc;i++) {
00702       cpl_vector_set(msk,i,1);
00703     }
00704     ovr_cube->plane[z]->data[x+y*ovr_cube->lx]=nc;
00705     msk_cube->plane[z]->data[x+y*msk_cube->lx]=nc;
00706         nclip=0;
00707         for (ks=0;ks<nc;ks++) {
00708       sig=0;
00709       med=0;
00710       ovr=0;
00711       val=cpl_vector_new(nc-nclip);
00712 
00713       for ( i = 0 ; i < nc ; i++ ) {
00714         if ((!qfits_isnan(cubes[i]->plane[z]->data[x+y*cubes[i]->lx])) &&
00715                  (cpl_vector_get(msk,i) != 0) ) {
00716           cpl_vector_set(val,ovr,(double)cubes[i]->plane[z]->data[x+y*cubes[i]->lx]); 
00717           ovr++;
00718         }
00719       }
00720 
00721       if(ovr>1) {
00722         avg=cpl_vector_get_mean(val);
00723         med=cpl_vector_get_median(val);
00724         sig=cpl_vector_get_stdev(val);
00725       } else {
00726         avg=cpl_vector_get(val,0);
00727         med=avg;
00728         sig=0;
00729       }
00730 
00731       med_cube->plane[z]->data[x+y*med_cube->lx]=med;
00732       avg_cube->plane[z]->data[x+y*avg_cube->lx]=avg;
00733       sig_cube->plane[z]->data[x+y*sig_cube->lx]=sig;
00734       cpl_vector_delete(val);
00735       for ( i = 0 ; i < nc ; i++ ) {
00736       /* Do k-s clipping at each pixel */
00737         if ((!qfits_isnan(cubes[i]->plane[z]->data[x+y*cubes[i]->lx])) && 
00738                 (cpl_vector_get(msk,i) != 0)) { 
00739           if(abs((cubes[i]->plane[z]->data[x+y*cubes[i]->lx]-med))> k*sig) {
00740         /* cubes[i]->plane[z]->data[x+y*cubes[i]->lx]=0; */
00741         msk_cube->plane[z]->data[x+y*msk_cube->lx]-=1;
00742         cpl_vector_set(msk,i,0);
00743                 nclip++;
00744           }
00745         }
00746       }
00747     }/* end of k-s clipping */
00748     msk_sum=0;
00749     val_msk_sum=0;
00750     for ( i = 0 ; i < nc ; i++ ) {
00751       /* computes sky at each point */
00752       if (!qfits_isnan(cubes[i]->plane[z]->data[x+y*cubes[i]->lx])) { 
00753         /*
00754         msk_sum+=msk_cube->plane[z]->data[x+y*msk_cube->lx];
00755                 val_msk_sum+=cubes[i]->plane[z]->data[x+y*cubes[i]->lx]*
00756           msk_cube->plane[z]->data[x+y*msk_cube->lx];
00757         */
00758         msk_sum+=cpl_vector_get(msk,i);
00759                 val_msk_sum+=cubes[i]->plane[z]->data[x+y*cubes[i]->lx]*
00760           cpl_vector_get(msk,i);
00761       }
00762     }
00763     sky_cube->plane[z]->data[x+y*sky_cube->lx]=val_msk_sum/msk_sum;
00764     cpl_vector_delete(msk);
00765       } /* end loop over x */
00766     } /* end loop over y */
00767   } /* end loop over z */
00768   for ( i = 0 ; i < nc ; i++ ) {
00769     cubes_sky[i]=copy_cube(cubes[i]);
00770     /* subtract the variable clean sky */
00771     op_cube(&(cubes_sky[i]),sky_cube,'-');
00772   }
00773 
00774 
00775   return cubes_sky;
00776 }
00777  
00778 int
00779 extract_ao_info(cpl_frameset* obj_set,cpl_frameset** set,const char* recid,cpl_parameterlist* config)
00780 {
00781 
00782   const char*  fctid=NULL;
00783   cpl_frame* frm=NULL;
00784   cpl_table* tbl=NULL;
00785   cpl_frame* product_frame=NULL;
00786   cpl_propertylist* plist=NULL;
00787   cpl_image* image=NULL;
00788 
00789   int nobj=0;
00790   int i=0;
00791   const char* name_o="out_ao_info.tfits";
00792 
00793   name_o="out_ao_info.tfits";
00794   nobj=cpl_frameset_get_size(obj_set);
00795   for(i=0;i<nobj;i++) {
00796     frm=cpl_frameset_get_frame(obj_set,i);
00797     tbl=cpl_table_load(cpl_frame_get_filename(frm),1,1);
00798     plist=cpl_propertylist_load(cpl_frame_get_filename(frm),0);
00799     if(i==0) {
00800       cpl_table_save(tbl,NULL,plist,name_o,CPL_IO_DEFAULT);
00801     } else {
00802       cpl_table_save(tbl,NULL,plist,name_o,CPL_IO_EXTEND);
00803     }
00804   }
00805   image=cpl_image_new(0,0,CPL_TYPE_FLOAT);
00806     /* Create product frame */
00807     product_frame = cpl_frame_new();
00808     cpl_frame_set_filename(product_frame, name_o) ;
00809     cpl_frame_set_tag(product_frame, PRO_AO_INFO) ;
00810     cpl_frame_set_type(product_frame, CPL_FRAME_TYPE_IMAGE) ;
00811     cpl_frame_set_group(product_frame, CPL_FRAME_GROUP_PRODUCT) ;
00812     cpl_frame_set_level(product_frame, CPL_FRAME_LEVEL_FINAL) ;
00813     if (cpl_error_get_code()) {
00814         cpl_msg_error(fctid, "Error while initialising the product frame") ;
00815         cpl_propertylist_delete(plist) ;
00816         cpl_frame_delete(product_frame) ;
00817         cpl_image_delete(image) ;
00818         return -1 ;
00819     }
00820     
00821      /* Add DataFlow keywords */
00822     if (cpl_dfs_setup_product_header(plist, product_frame, *set, config,
00823             recid, "SINFONI", "?Dictionary?") != CPL_ERROR_NONE) {
00824         cpl_msg_error(fctid, "Problem in the product DFS-compliance") ;
00825         cpl_propertylist_delete(plist) ;
00826         cpl_frame_delete(product_frame) ;
00827         cpl_image_delete(image) ;
00828         return -1 ;
00829     }
00830 
00831      /* Save the file */
00832     if (cpl_image_save(image, name_o, CPL_BPP_DEFAULT, plist,
00833                        CPL_IO_DEFAULT) != CPL_ERROR_NONE) {
00834         cpl_msg_error(fctid, "Could not save product");
00835         cpl_propertylist_delete(plist) ;
00836         cpl_frame_delete(product_frame) ;
00837         cpl_image_delete(image) ;
00838         return -1 ;
00839     }
00840      cpl_propertylist_delete(plist) ;
00841      cpl_image_delete(image) ;
00842  
00843     /* Log the saved file in the input frameset */
00844     cpl_frameset_insert(*set, product_frame) ;
00845  
00846   return 0;
00847 
00848 }
00849 
00850 int
00851 assign_offset(const int n,char* name,float* offsetx,float* offsety,
00852               const float ref_offx,const float ref_offy)
00853 {
00854 
00855   const char * _id = "assign_offset";
00856   float offx=0;
00857   float offy=0;
00858   double mjd_obs=0;
00859   /*
00860   double pixelscale=0;
00861   double angle=0;
00862   double radangle=0;
00863   double cd1_1=0;
00864   double cd1_2=0;
00865   double cd2_1=0;
00866   double cd2_2=0;
00867   double ra=0;
00868   double dec=0;
00869   */
00870 
00871   cpl_propertylist * plist=NULL;
00872   cpl_msg_info(_id,"Assign offsets");
00873   
00874   /*
00875   pixelscale = spiffi_get_pixelscale(name)/2.;
00876   angle = spiffi_get_posangle(name) ;
00877   radangle = - angle * PI_NUMB / 180. ;
00878   cd1_1 = cos(radangle) ;
00879   cd1_2 = sin(radangle) ;
00880   cd2_1 = -sin(radangle) ;
00881   cd2_2 = cos(radangle) ;
00882   
00883   cpl_msg_info(_id,"scale=%f angle=%f cd1_1=%f cd1_2=%f cd2_1=%f cd2_2=%f\n",
00884            pixelscale,angle,cd1_1,cd1_2,cd2_1,cd2_2);
00885   
00886   ra  = spiffi_get_cumoffseta(name) ;
00887   dec = spiffi_get_cumoffsetd(name) ;
00888 
00889   offx= -(ra*cd1_1 + dec*cd1_2);
00890   offy= -(ra*cd2_1 + dec*cd2_2);
00891   
00892   cpl_msg_info(_id,"RA=%f DEC=%f offx=%f offy=%f\n",ra,dec,offx,offy);
00893   offx/=pixelscale;
00894   offy/=pixelscale;
00895   cpl_msg_info(_id,"RA=%f DEC=%f offx=%f offy=%f\n",ra,dec,offx,offy);
00896   */
00897   /*
00898   offx_rot = cd1_1 * offx + cd2_1 * offy ;
00899   offy_rot = cd1_2 * offx + cd2_2 * offy ;
00900    */
00901   /* convert the coordinates to pixel units */
00902   /*
00903   offx_rot_pix = offx_rot / pixelscale ;
00904   offy_rot_pix = offy_rot / pixelscale ;
00905    */
00906   /* cpl_msg_info(_id,"offset(posangle) pix: %f,%f",offx_rot_pix,offy_rot_pix); */
00907   /*
00908   offsetx[i] = offx_rot_pix ;
00909   offsety[i] = offy_rot_pix ;
00910    */
00911 
00912   offx = spiffi_get_cumoffsetx(name) - ref_offx ;  /* was - */
00913   if (offx == FLAG) {
00914     cpl_msg_warning(_id," could not read fits header keyword cummoffsetx!") ;
00915     cpl_msg_warning(_id," Set relative offset to 0 - %f!",ref_offx) ;
00916     offx =  - ref_offx;
00917     /* return -1 ; */
00918   }
00919     
00920   offy = spiffi_get_cumoffsety(name) - ref_offy ; /* was - */
00921   if (offy == FLAG ) {
00922     cpl_msg_warning(_id," could not read fits header keyword! cumoffsety") ;
00923     cpl_msg_warning(_id," Set relative offset to 0 - %f!",ref_offx) ;
00924     offy =  - ref_offy;
00925     /* return -1 ; */
00926   }
00927   cpl_msg_info(_id,"offx=%f offy=%f",offx,offy);
00928 
00929 
00930 
00931 
00932   if ((cpl_error_code)((plist = cpl_propertylist_load(name, 0)) == NULL)) {
00933     cpl_msg_error(_id, "getting header from reference frame %s",name);
00934     cpl_propertylist_delete(plist) ;
00935     cpl_free(name);
00936     return -1 ;
00937   }
00938 
00939   if (cpl_propertylist_contains(plist, KEY_NAME_MJD_OBS)) {
00940     mjd_obs=cpl_propertylist_get_double(plist, KEY_NAME_MJD_OBS);
00941   } else {
00942     cpl_msg_error(_id,"keyword %s does not exist",KEY_NAME_MJD_OBS);
00943     cpl_propertylist_delete(plist) ;
00944     return -1;
00945   }
00946   cpl_propertylist_delete(plist) ;
00947 
00948   if (mjd_obs > 53421.58210082 ) {
00949     /* after detector's upgrade */
00950     /*
00951     array_set_value(offsetx,-offx*2,n);
00952     array_set_value(offsety,+offy*2,n);
00953     */
00954     array_set_value(offsetx,-2*offx,n);
00955     array_set_value(offsety,2*offy,n);
00956   } else {
00957     /* before detector's upgrade */
00958     /*
00959     array_set_value(offsetx,+offx*2,n);
00960     array_set_value(offsety,-offy*2,n);
00961     */
00962     array_set_value(offsetx,2*offx,n);
00963     array_set_value(offsety,-2*offy,n);
00964   }
00965 
00966   return 0;
00967 
00968 
00969 }
00970 
00971 int
00972 object_assign_offset(char* name, const int n, double* ref_offx, double* ref_offy, float** offsetx, float** offsety)
00973 {
00974 
00975   const char * _id = "object_assign_offset";
00976   float offx=0;
00977   float offy=0;
00978   double mjd_obs=0;
00979   cpl_propertylist * plist=NULL;
00980   cpl_msg_info(_id,"Assign offsets");
00981   if ( n == 0 )
00982     {
00983 
00984 
00985       *ref_offx = spiffi_get_cumoffsetx(name) ;
00986       if (*ref_offx == FLAG)
00987     {
00988       cpl_msg_error(_id," could not read fits header keyword cummoffsetx!") ;
00989       /* return -1 ; */
00990     }   
00991     
00992       *ref_offy = spiffi_get_cumoffsety(name) ;
00993       if (*ref_offy == FLAG )
00994     {
00995       cpl_msg_error(_id," could not read fits header keyword! cumoffsety") ;
00996       /* return -1 ; */
00997     }
00998       cpl_msg_info(_id,"Reference offx=%f offy=%f",*ref_offx,*ref_offy);
00999           
01000       offx = 0. ;
01001       offy = 0. ;
01002 
01003     } else {
01004        
01005       offx = spiffi_get_cumoffsetx(name) - *ref_offx ;  /* was - */
01006       if (offx == FLAG)
01007     {
01008       cpl_msg_error(_id," could not read fits header keyword cummoffsetx!") ;
01009       /* return -1 ; */
01010     }
01011     
01012       offy = spiffi_get_cumoffsety(name) - *ref_offy ; /* was - */
01013       if (offy == FLAG )
01014     {
01015       cpl_msg_error(_id," could not read fits header keyword! cumoffsety") ;
01016       /* return -1 ; */
01017     }
01018       cpl_msg_info(_id,"offx=%f offy=%f",offx,offy);
01019                
01020     }
01021            
01022     /* rotate the coordinates 
01023        offx_rot = cd1_1 * offx + cd2_1 * offy ;
01024        offy_rot = cd1_2 * offx + cd2_2 * offy ;
01025         convert the coordinates to pixel units 
01026        offx_rot_pix = offx_rot / pixelscale ;
01027        offy_rot_pix = offy_rot / pixelscale ;
01028        offsetx[i] = offx_rot_pix ;
01029        offsety[i] = offy_rot_pix ;
01030     */
01031 
01032   if ((cpl_error_code)((plist = cpl_propertylist_load(name, 0)) == NULL)) {
01033     cpl_msg_error(_id, "getting header from reference frame %s",name);
01034     cpl_propertylist_delete(plist) ;
01035     cpl_free(name);
01036     return -1 ;
01037   }
01038 
01039   if (cpl_propertylist_contains(plist, KEY_NAME_MJD_OBS)) {
01040     mjd_obs=cpl_propertylist_get_double(plist, KEY_NAME_MJD_OBS);
01041   } else {
01042     cpl_msg_error(_id,"keyword %s does not exist",KEY_NAME_MJD_OBS);
01043     cpl_propertylist_delete(plist) ;
01044     return -1;
01045   }
01046   cpl_propertylist_delete(plist) ;
01047 
01048   if (mjd_obs > 53421.58210082 ) {
01049     /* after detector's upgrade */
01050     array_set_value(*offsetx,-offx*2,n);
01051     array_set_value(*offsety,+offy*2,n);
01052   } else {
01053     /* before detector's upgrade */
01054     array_set_value(*offsetx,+offx*2,n);
01055     array_set_value(*offsety,-offy*2,n);
01056   }
01057 
01058   return 0;
01059 }
01060 
01061 
01062 
01063 
01064 OneCube*
01065 fine_tune(OneCube* cube,float* correct_dist, const char* method, const int order, const int nslits) {
01066   const char * _id= "fine_tune";
01067   int i =0;
01068   OneCube* outcube2=NULL;
01069   float* neg_dist=NULL;
01070   cpl_msg_info(_id,"Finetuning, method=%s",method);
01071   
01072   if (strcmp(method,"P")==0)
01073     {
01074       outcube2 = fineTuneCube( cube, correct_dist, order ) ;
01075       if (outcube2 == NULL)
01076     {
01077       cpl_msg_error (_id," could not fine tune the data cube\n") ;
01078       return NULL ;
01079     }
01080     }       
01081   else if (strcmp(method,"F")==0)
01082     {
01083       neg_dist=cpl_calloc(nslits,sizeof(float));
01084       for ( i = 0 ; i < nslits ; i++ )
01085     {
01086       neg_dist[i] = -correct_dist[i] ;
01087     }
01088       outcube2 = fineTuneCubeByFFT( cube, neg_dist ) ;
01089      cpl_free(neg_dist);
01090       if ( outcube2 == NULL )
01091     {
01092       cpl_msg_error (_id," could not fine tune the data cube\n") ;
01093       return NULL ;
01094     }
01095     }       
01096   else if (strcmp(method,"S")==0)
01097     {
01098       outcube2 = fineTuneCubeBySpline( cube, correct_dist ) ;
01099       if ( outcube2 == NULL )
01100     {
01101       cpl_msg_error (_id," could not fine tune the data cube\n") ;
01102       return NULL ;
01103     }
01104     }
01105   else
01106     {
01107       cpl_msg_error (_id," wrong method indicator given!") ;
01108       return NULL ;
01109     }
01110 
01111 
01112 
01113 return outcube2;
01114 
01115 }

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