sinfo_utilities_scired.c

00001 /*
00002  * This file is part of the ESO SINFONI Pipeline
00003  * Copyright (C) 2004,2005 European Southern Observatory
00004  *
00005  * This program is free software; you can redistribute it and/or modify
00006  * it under the terms of the GNU General Public License as published by
00007  * the Free Software Foundation; either version 2 of the License, or
00008  * (at your option) any later version.
00009  *
00010  * This program is distributed in the hope that it will be useful,
00011  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00012  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00013  * GNU General Public License for more details.
00014  *
00015  * You should have received a copy of the GNU General Public License
00016  * along with this program; if not, write to the Free Software
00017  * Foundation, 51 Franklin St, Fifth Floor, Boston, MA  02111-1307  USA
00018  */
00019 #ifdef HAVE_CONFIG_H
00020 #  include <config.h>
00021 #endif
00022 #include "sinfo_utilities_scired.h"
00023 #include "sinfo_functions.h"
00024 #include "sinfo_pfits.h"
00025 #include "sinfo_utils_wrappers.h"
00033 int
00034 sinfo_check_input_data(object_config* cfg)
00035 {
00036 
00037 
00038     if (cfg == NULL)
00039     {
00040         sinfo_msg_error (" could not parse cpl input!\n") ;
00041         return -1 ;
00042     }
00043 
00044 
00045     if(sinfo_is_fits_file(cfg->wavemap) != 1) {
00046       sinfo_msg_error("Input file wavemap %s is not FITS",cfg->wavemap);
00047       return -1;
00048     }
00049 
00050 
00051     if (cfg->halocorrectInd == 1)
00052     {
00053        if(sinfo_is_fits_file(cfg->halospectrum) != 1) {
00054           sinfo_msg_error("Input file %s is not FITS",cfg->halospectrum);
00055           return -1;
00056        }
00057 
00058     }
00059 
00060     if (cfg->northsouthInd == 0) {
00061       if (sinfo_is_fits_file(cfg->poslist) != 1)
00062     {
00063       sinfo_msg_error("File %s with tag %s is not FITS!",
00064                           cfg->poslist,PRO_SLIT_POS);
00065       return -1 ;
00066     }
00067     } else {
00068 
00069       if (sinfo_is_fits_file(cfg->distlist) != 1)
00070     {
00071       sinfo_msg_error("File %s with tag %s is not FITS!",
00072                            cfg->distlist,PRO_SLITLETS_DISTANCE);
00073       return -1;
00074     }
00075     }
00076 
00077 
00078     return 0;
00079      
00080 
00081 }
00082 
00083 
00084 
00096 int
00097 sinfo_auto_size_cube3(char** framelist, 
00098                     const int nframes, 
00099                     float* ref_offx, 
00100                     float* ref_offy,
00101                     int* size_x,
00102                     int* size_y)
00103 {
00104 
00105   char* name =NULL;
00106   int n=0;
00107   float offx=0;
00108   float offy=0;
00109   float min_offx=0;
00110   float max_offx=0;
00111   float min_offy=0;
00112   float max_offy=0;
00113   cpl_propertylist* plist=NULL;
00114   sinfo_msg ("Automatic computation of output cube size") ; 
00115   for ( n = 0 ; n < nframes ; n++ ) {
00116     name = framelist[n] ;
00117     plist=cpl_propertylist_load(name,0);
00118     offx = sinfo_pfits_get_cumoffsetx(plist);  /* was - */
00119     if(cpl_error_get_code() != CPL_ERROR_NONE) {
00120       sinfo_msg_warning(" could not read fits header keyword cummoffsetx!");
00121       sinfo_msg_warning(" set it to  0");
00122       offx = 0;
00123       cpl_error_reset();
00124       /* return -1 ; */
00125     }
00126     
00127     offy = sinfo_pfits_get_cumoffsety(plist); /* was - */
00128     sinfo_free_propertylist(&plist);
00129     if(cpl_error_get_code() != CPL_ERROR_NONE) {
00130       sinfo_msg_warning(" could not read fits header keyword! cumoffsety");
00131       sinfo_msg_warning(" set it to 0");
00132       offy = 0;
00133       cpl_error_reset();
00134       /* return -1 ; */
00135     }
00136     if(n==0) {
00137       min_offx=offx;
00138       min_offy=offy;
00139       max_offx=offx;
00140       max_offy=offy;
00141     } else {
00142       if(offx > max_offx) max_offx=offx;
00143       if(offy > max_offy) max_offy=offy;
00144       if(offx < min_offx) min_offx=offx;
00145       if(offy < min_offy) min_offy=offy;
00146     }
00147   }
00148 
00149   *ref_offx=(min_offx+max_offx)/2;
00150   *ref_offy=(min_offy+max_offy)/2;
00151 
00152   if(*size_x == 0) *size_x=2*floor(max_offx-min_offx+0.5)+64;
00153   if(*size_y == 0) *size_y=2*floor(max_offy-min_offy+0.5)+64;
00154   sinfo_msg("Output cube size: %d x %d",*size_x,*size_y);
00155   sinfo_msg_debug("Ref offset. x: %f y: %f",*ref_offx,*ref_offy);
00156   sinfo_msg_debug("Max offset. x: %f y: %f",max_offx,max_offy);
00157   sinfo_msg_debug("Min offset. x: %f y: %f",min_offx,min_offy);
00158   return 0;  
00159 
00160 }
00161 
00162 
00163 
00172 int
00173 sinfo_auto_size_cube2(object_config * cfg, float* ref_offx, float* ref_offy)
00174 {
00175 
00176   char* name =NULL;
00177   int n=0;
00178   float offx=0;
00179   float offy=0;
00180   float min_offx=0;
00181   float max_offx=0;
00182   float min_offy=0;
00183   float max_offy=0;
00184   cpl_propertylist * plist=NULL;
00185   sinfo_msg ("Automatic computation of output cube size") ;
00186   for ( n = 0 ; n < cfg->nframes ; n++ ) {
00187     name = cfg->framelist[n] ;    
00188     plist=cpl_propertylist_load(name,0);
00189     offx = sinfo_pfits_get_cumoffsetx(plist);  /* was - */
00190     if(cpl_error_get_code() != CPL_ERROR_NONE) {
00191       sinfo_msg_warning(" could not read fits header keyword cummoffsetx!");
00192       sinfo_msg_warning(" set it to  0");
00193       offx = 0;
00194       cpl_error_reset();
00195       /* return -1 ; */
00196     }   
00197     
00198     offy = sinfo_pfits_get_cumoffsety(plist); /* was - */
00199     sinfo_free_propertylist(&plist);
00200     if(cpl_error_get_code() != CPL_ERROR_NONE) {
00201       sinfo_msg_warning(" could not read fits header keyword! cumoffsety") ;
00202       sinfo_msg_warning(" set it to 0");
00203       offy = 0;
00204       cpl_error_reset();
00205       /* return -1 ; */
00206     }
00207     if(n==0) {
00208       min_offx=offx;
00209       min_offy=offy;
00210       max_offx=offx;
00211       max_offy=offy;
00212     } else {
00213       if(offx > max_offx) max_offx=offx;
00214       if(offy > max_offy) max_offy=offy;
00215       if(offx < min_offx) min_offx=offx;
00216       if(offy < min_offy) min_offy=offy;
00217     }
00218   }
00219   *ref_offx=(min_offx+max_offx)/2;
00220   *ref_offy=(min_offy+max_offy)/2;
00221 
00222   if(cfg->size_x == 0) cfg->size_x=2*floor(max_offx-min_offx+0.5)+64;
00223   if(cfg->size_y == 0) cfg->size_y=2*floor(max_offy-min_offy+0.5)+64;
00224   sinfo_msg("Output cube size: %d x %d",cfg->size_x,cfg->size_y);
00225   sinfo_msg("Ref offset. x: %f y: %f",*ref_offx,*ref_offy);
00226   sinfo_msg("Max offset. x: %f y: %f",max_offx,max_offy);
00227   sinfo_msg("Min offset. x: %f y: %f",min_offx,min_offy);
00228   return 0;
00229     
00230 
00231 }
00232 
00241 float*
00242 sinfo_read_distances(object_config* cfg)
00243 {
00244   int i=0;
00245   int* status=NULL;
00246   float * distances = NULL;
00247   float tmp_float=0;
00248   char tbl_distances_name[FILE_NAME_SZ];
00249   cpl_table* tbl_distances = NULL;
00250 
00251   sinfo_msg("Read distances");
00252   distances = (float*) cpl_calloc (cfg->nslits - 1, sizeof (float));
00253 
00254   if ( NULL == distances )
00255     {
00256       sinfo_msg_error ("could allocate memory!") ;
00257       return NULL ;
00258     }
00259 
00260             /*READ TFITS TABLE*/
00261   if(cpl_error_get_code() != CPL_ERROR_NONE) {
00262     sinfo_msg_error("Before loading input table");
00263     sinfo_msg_error((char* ) cpl_error_get_message());
00264     return NULL;
00265   }
00266   strcpy(tbl_distances_name,cfg->distlist);
00267   tbl_distances = cpl_table_load(tbl_distances_name,1,0);
00268   if(cpl_error_get_code() != CPL_ERROR_NONE) {
00269     sinfo_msg_error("loading input table %s",tbl_distances_name);
00270     sinfo_msg_error((char* ) cpl_error_get_message());
00271     return NULL;
00272   }
00273 
00274   for (i =0 ; i< cfg->nslits-1; i++){
00275     tmp_float=cpl_table_get_float(tbl_distances,"slitlet_distance",i,status);
00276     if(cpl_error_get_code() != CPL_ERROR_NONE) {
00277       sinfo_msg_error("reading col %s from table %s","slitlet_distance",
00278                        tbl_distances_name);
00279       sinfo_msg_error((char* ) cpl_error_get_message());
00280       return NULL;
00281     }
00282     sinfo_new_array_set_value(distances,tmp_float,i);
00283   }
00284   cpl_table_delete(tbl_distances);
00285   return distances;
00286 
00287 }
00288 
00289 
00290 
00297 float** 
00298 sinfo_read_slitlets_edges(object_config * cfg)
00299 {
00300 
00301 
00302   char tbl_slitpos_name[FILE_NAME_SZ];
00303   cpl_table* tbl_slitpos=NULL;
00304   int n=0;
00305   int i=0;
00306   int* status=NULL;
00307   float edge_x=0;
00308   float edge_y=0;
00309   float ** slit_edges = NULL;
00310 
00311   slit_edges = sinfo_new_2Dfloatarray(cfg->nslits, 2) ;
00312 
00313   strcpy(tbl_slitpos_name,cfg->poslist);
00314   tbl_slitpos = cpl_table_load(tbl_slitpos_name,1,0);
00315   if(cpl_error_get_code() != CPL_ERROR_NONE) {
00316     sinfo_msg_error("error loading tbl %s",tbl_slitpos_name);
00317     sinfo_msg_error((char* ) cpl_error_get_message());
00318     return NULL;
00319   }
00320   n = cpl_table_get_nrow(tbl_slitpos);
00321   if (n != cfg->nslits) {
00322     sinfo_msg_error("No of slitlets in table is n = %d != %d !",n,cfg->nslits);
00323     return NULL;
00324   }
00325 
00326   for (i =0 ; i< cfg->nslits; i++){
00327     edge_x=cpl_table_get_double(tbl_slitpos,"pos1",i,status);
00328     edge_y=cpl_table_get_double(tbl_slitpos,"pos2",i,status);
00329     if(cpl_error_get_code() != CPL_ERROR_NONE) {
00330       sinfo_msg_error("error reading tbl %s row %d",tbl_slitpos_name,i);
00331       sinfo_msg_error((char* ) cpl_error_get_message());
00332       return NULL;
00333     }
00334     sinfo_new_array2D_set_value(slit_edges,edge_x,i,0);
00335     sinfo_new_array2D_set_value(slit_edges,edge_y,i,1);
00336   }
00337   cpl_table_delete(tbl_slitpos);
00338   if(cpl_error_get_code() != CPL_ERROR_NONE) {
00339     sinfo_msg_error("error reading tbl %s",tbl_slitpos_name);
00340     sinfo_msg_error((char* ) cpl_error_get_message());
00341     return NULL;
00342   }
00343  
00344   return slit_edges;
00345 
00346 }
00347 
00348 
00349 /*---------------------------------------------------------------------------*/
00372 /*---------------------------------------------------------------------------*/
00373 cpl_imagelist * sinfo_new_cube_getvig(
00374         cpl_imagelist *        cube_in,
00375         int             loleft_x,
00376         int             loleft_y,
00377         int             upright_x,
00378         int             upright_y)
00379 {
00380     cpl_imagelist     *        cube_out ;
00381     int                 i ;
00382     int                 outlx,
00383                                         outly ;
00384 
00385     int ilx=0;
00386     int ily=0;
00387     int inp=0;
00388     cpl_image* i_img=NULL;
00389     cpl_image* o_img=NULL;
00390 
00391     if (cube_in==NULL) return NULL ;
00392 
00393     ilx=cpl_image_get_size_x(cpl_imagelist_get(cube_in,0));
00394     ily=cpl_image_get_size_y(cpl_imagelist_get(cube_in,0));
00395     inp=cpl_imagelist_get_size(cube_in);
00396 
00397         if ((loleft_x>upright_x) ||
00398             (loleft_y>upright_y)) {
00399             sinfo_msg_error("ill-defined slit for extraction: aborting");
00400             return NULL ;
00401         }
00402 
00403     /* Extraction coordinates include rectangular zone  */
00404     outlx = upright_x - loleft_x + 1 ;
00405     outly = upright_y - loleft_y + 1 ;
00406     /*
00407     cube_out = sinfo_new_cube(outlx, outly, cube_in->np) ;
00408     */
00409     cube_out = cpl_imagelist_new() ;
00410     /* Loop on all input planes */
00411     for (i=0 ; i<cpl_imagelist_get_size(cube_in) ; i++) {
00412 
00413       i_img=cpl_imagelist_get(cube_in,i);
00414         /* Extract a slit from this plane   */
00415                 o_img = sinfo_new_image_getvig(i_img,
00416                         loleft_x, loleft_y,
00417                         upright_x, upright_y) ;
00418         cpl_imagelist_set(cube_out,o_img,i);
00419     }
00420     return cube_out ;
00421 }
00422 /*---------------------------------------------------------------------------*/
00439 /*---------------------------------------------------------------------------*/
00440 
00441 cpl_image * sinfo_new_image_getvig(
00442         cpl_image    *    image_in,
00443         int             loleft_x,
00444         int             loleft_y,
00445         int             upright_x,
00446         int             upright_y)
00447 {
00448     cpl_image    *        slit_img ;
00449     int                 i, j ;
00450     register
00451     pixelvalue  *       inpt,
00452                                 *       outpt ;
00453     int         outlx, outly ;
00454     int ilx=0;
00455     int ily=0;
00456  
00457     float* pidata=NULL;
00458     float* podata=NULL;
00459 
00460     if (image_in==NULL) return NULL ;
00461 
00462     ilx=cpl_image_get_size_x(image_in);
00463     ily=cpl_image_get_size_y(image_in);
00464     pidata=cpl_image_get_data_float(image_in);
00465 
00466     if ((loleft_x<1) || (loleft_x>ilx) ||
00467         (loleft_y<1) || (loleft_y>ily) ||
00468         (upright_x<1) || (upright_x>ilx) ||
00469         (upright_y<1) || (upright_y>ily) ||
00470         (loleft_x>upright_x) || (loleft_y>upright_y)) {
00471         sinfo_msg_error("extraction zone is [%d %d] [%d %d]\n"
00472                         "cannot extract such zone: aborting slit extraction",
00473                 loleft_x, loleft_y, upright_x, upright_y) ;
00474         return NULL ;
00475     }
00476 
00477     outlx = upright_x - loleft_x + 1 ;
00478     outly = upright_y - loleft_y + 1 ;
00479     slit_img = cpl_image_new(outlx, outly,CPL_TYPE_FLOAT) ;
00480     podata=cpl_image_get_data_float(slit_img);
00481 
00482     for (j=0 ; j<outly ; j++) {
00483         inpt = pidata+loleft_x-1 + (j+loleft_y-1)*ilx ;
00484         outpt = podata + j*outlx ;
00485         for (i=0 ; i<outlx ; i++) {
00486             *outpt++ = *inpt++ ;
00487         }
00488     }
00489     return slit_img ;
00490 }
00491 
00503 int
00504 sinfo_new_set_wcs_cube(cpl_imagelist* cub, const char* name, double clambda, 
00505          double dis, double cpix, double cx, double cy)
00506 {
00507   cpl_propertylist* plist=NULL;
00508 
00509   if ((cpl_error_code)((plist = cpl_propertylist_load(name, 0)) == NULL)) {
00510     sinfo_msg_error( "getting header from frame %s",name);
00511     cpl_propertylist_delete(plist) ;
00512     return -1 ;
00513   }
00514 
00515   sinfo_new_change_plist_cube(plist, clambda, dis, cpix, cx, cy) ;
00516         
00517   if (cpl_imagelist_save(cub, name, CPL_BPP_DEFAULT,
00518                          plist,CPL_IO_DEFAULT)!=CPL_ERROR_NONE) {
00519     sinfo_msg_error( "Cannot save the product %s",name);
00520     cpl_propertylist_delete(plist) ;
00521     return -1 ;
00522   }
00523   cpl_propertylist_delete(plist) ;
00524   return 0;
00525 
00526 }
00535 int
00536 sinfo_new_set_wcs_image(cpl_image* img, 
00537                         const char* name,  
00538                         double cx, 
00539                         double cy)
00540 {
00541   cpl_propertylist* plist=NULL;
00542 
00543   if ((cpl_error_code)((plist = cpl_propertylist_load(name, 0)) == NULL)) {
00544     sinfo_msg_error( "getting header from frame %s",name);
00545     cpl_propertylist_delete(plist) ;
00546     return -1 ;
00547   }
00548   sinfo_new_change_plist_image(plist, cx, cy) ;
00549         
00550   if (cpl_image_save(img, name, CPL_BPP_DEFAULT,
00551                      plist,CPL_IO_DEFAULT)!=CPL_ERROR_NONE) {
00552     sinfo_msg_error( "Cannot save the product %s",name);
00553     cpl_propertylist_delete(plist) ;
00554     return -1 ;
00555   }
00556   cpl_propertylist_delete(plist) ;
00557   return 0;
00558 
00559 }
00560 
00561 
00571 int
00572 sinfo_new_set_wcs_spectrum(cpl_image* img, const char* name, double clambda, 
00573              double dis, double cpix)
00574 {
00575   cpl_propertylist* plist=NULL;
00576 
00577   if ((cpl_error_code)((plist = cpl_propertylist_load(name, 0)) == NULL)) {
00578     sinfo_msg_error( "getting header from frame %s",name);
00579     cpl_propertylist_delete(plist) ;
00580     return -1 ;
00581   }
00582   sinfo_new_change_plist_spectrum(plist, clambda, dis,cpix) ;
00583         
00584   if (cpl_image_save(img, name, CPL_BPP_DEFAULT,
00585                      plist,CPL_IO_DEFAULT)!=CPL_ERROR_NONE) {
00586     sinfo_msg_error( "Cannot save the product %s",name);
00587     cpl_propertylist_delete(plist) ;
00588     return -1 ;
00589   }
00590   cpl_propertylist_delete(plist) ;
00591   return 0;
00592 
00593 }
00605 void sinfo_new_change_plist_cube (cpl_propertylist * plist, 
00606             float cenLambda,
00607             float dispersion,
00608             int   cenpix,
00609             float center_x,
00610             float center_y )
00611 {
00612 
00613     float pixelscale ;
00614     double ra ;
00615     double dec ;
00616     double angle ;
00617     float radangle ;
00618     float cd1_1, cd1_2, cd2_1, cd2_2 ;
00619     char firsttext[2*FILE_NAME_SZ] ;
00620     strcpy(firsttext, "sinfo_rec_objnod -f \0") ;
00621    
00622     pixelscale = sinfo_pfits_get_pixscale(plist)/2. ;
00623     ra = sinfo_pfits_get_ra(plist) ;
00624     dec = sinfo_pfits_get_DEC(plist) ;
00625     angle = sinfo_pfits_get_posangle(plist) ;
00626     /* in PUPIL data there is not posangle info: we reset the error */
00627     if(cpl_error_get_code() != CPL_ERROR_NONE) {
00628        cpl_error_reset();
00629     }
00630 
00631 
00632 
00633     radangle = angle * PI_NUMB / 180. ;
00634     cd1_1 = cos(radangle)  * pixelscale / 3600. ;
00635     cd1_2 = sin(radangle) * pixelscale / 3600. ;
00636     cd2_1 = -sin(radangle)  * pixelscale / 3600. ;
00637     cd2_2 = cos(radangle)  * pixelscale / 3600. ;
00638 
00639     cpl_propertylist_erase_regexp(plist, "^CTYPE1",0);
00640     cpl_propertylist_insert_after_string(plist,"EXPTIME","CTYPE1","RA---TAN");
00641     cpl_propertylist_set_comment(plist, "CTYPE1", "Projected Rectascension");
00642     cpl_propertylist_erase_regexp(plist, "^CRPIX1",0);
00643     cpl_propertylist_insert_after_float(plist,"CTYPE1","CRPIX1", center_x ) ;
00644     cpl_propertylist_set_comment(plist, "CRPIX1","Reference pixel in RA" ) ;
00645 
00646     cpl_propertylist_erase_regexp(plist, "^CRVAL1",0);
00647     cpl_propertylist_insert_after_double(plist, "CRPIX1",  "CRVAL1", ra ) ;
00648     cpl_propertylist_set_comment(plist, "CRVAL1","Reference RA" ) ;
00649 
00650     cpl_propertylist_erase_regexp(plist, "^CDELT1",0);
00651     cpl_propertylist_insert_after_float(plist,"CRVAL1","CDELT1",
00652                                         -pixelscale/3600.  ) ;
00653     cpl_propertylist_set_comment(plist, "CDELT1","pixel scale" ) ;
00654 
00655     cpl_propertylist_erase_regexp(plist, "^CUNIT1",0);
00656     cpl_propertylist_insert_after_string(plist, "CDELT1",  "CUNIT1", "deg" ) ;
00657     cpl_propertylist_set_comment(plist, "CUNIT1","RA-UNIT" ) ;
00658 
00659     cpl_propertylist_erase_regexp(plist, "^CTYPE2",0);
00660     cpl_propertylist_insert_after_string(plist,"CUNIT1","CTYPE2","DEC--TAN");
00661     cpl_propertylist_set_comment(plist, "CTYPE2", "Projected Declination") ;
00662 
00663     cpl_propertylist_erase_regexp(plist, "^CRPIX2",0);
00664     cpl_propertylist_insert_after_float(plist, "CTYPE2",  "CRPIX2", center_y ) ;
00665     cpl_propertylist_set_comment(plist, "CRPIX2", "Reference pixel in DEC") ;
00666 
00667     cpl_propertylist_erase_regexp(plist, "^CRVAL2",0);
00668     cpl_propertylist_insert_after_double(plist, "CRPIX2",  "CRVAL2", dec) ;
00669     cpl_propertylist_set_comment(plist, "CRVAL2", "Reference DEC") ;
00670 
00671     cpl_propertylist_erase_regexp(plist, "^CDELT2",0);
00672     cpl_propertylist_insert_after_float(plist, "CRVAL2",  "CDELT2", 
00673                                         pixelscale/3600. ) ;
00674     cpl_propertylist_set_comment(plist, "CDELT2", "pixel scale") ;
00675 
00676     cpl_propertylist_erase_regexp(plist, "^CUNIT2",0);
00677     cpl_propertylist_insert_after_string(plist, "CDELT2",  "CUNIT2", "deg" ) ;
00678     cpl_propertylist_set_comment(plist, "CUNIT2", "DEC-UNIT") ;
00679 
00680     cpl_propertylist_erase_regexp(plist, "^CTYPE3",0);
00681     cpl_propertylist_insert_after_string(plist, "CUNIT2",  "CTYPE3", "WAVE" ) ;
00682     cpl_propertylist_set_comment(plist,"CTYPE3","wavelength axis in microns") ;
00683 
00684     cpl_propertylist_erase_regexp(plist, "^CRPIX3",0);
00685     cpl_propertylist_insert_after_int(plist, "CTYPE3",  "CRPIX3", cenpix ) ;
00686     cpl_propertylist_set_comment(plist, "CRPIX3", "Reference pixel in z") ;
00687 
00688     cpl_propertylist_erase_regexp(plist, "^CRVAL3",0);
00689     cpl_propertylist_insert_after_float(plist,"CRPIX3", "CRVAL3", cenLambda) ;
00690     cpl_propertylist_set_comment(plist, "CRVAL3", "central wavelength") ;
00691 
00692     cpl_propertylist_erase_regexp(plist, "^CDELT3",0);
00693 
00694     cpl_propertylist_insert_after_float(plist,"CRVAL3","CDELT3",dispersion) ;
00695     cpl_propertylist_set_comment(plist, "CDELT3", "microns per pixel") ;
00696 
00697     cpl_propertylist_erase_regexp(plist, "^CUNIT3",0);
00698     cpl_propertylist_insert_after_string(plist,"CDELT3", "CUNIT3", "mum" ) ;
00699     cpl_propertylist_set_comment(plist, "CUNIT3",  "spectral unit" ) ;
00700     
00701     cpl_propertylist_erase_regexp(plist, "^CD1_1",0);
00702     cpl_propertylist_insert_after_float(plist, "CUNIT3",  "CD1_1", cd1_1 ) ;
00703     cpl_propertylist_set_comment(plist, "CD1_1",  "CD rotation matrix" ) ;
00704 
00705     cpl_propertylist_erase_regexp(plist, "^CD1_2",0);
00706     cpl_propertylist_insert_after_float(plist, "CD1_1",   "CD1_2", cd1_2 ) ;
00707     cpl_propertylist_set_comment(plist, "CD1_2",  "CD rotation matrix" ) ;
00708     
00709     cpl_propertylist_erase_regexp(plist, "^CD2_1",0);
00710     cpl_propertylist_insert_after_float(plist, "CD1_2",   "CD2_1", cd2_1 ) ;
00711     cpl_propertylist_set_comment(plist, "CD2_1",  "CD rotation matrix" ) ;
00712     
00713     cpl_propertylist_erase_regexp(plist, "^CD2_2",0);
00714     cpl_propertylist_insert_after_float(plist, "CD2_1",   "CD2_2", cd2_2 ) ;
00715     cpl_propertylist_set_comment(plist, "CD2_2",  "CD rotation matrix" ) ;
00716 
00717 
00718 }        
00719 
00728 void sinfo_new_change_plist_spectrum (cpl_propertylist * plist, 
00729                 double cenLambda,
00730                 double dispersion,
00731                 int   cenpix)
00732 {
00733 
00734     cpl_propertylist_erase_regexp(plist, "^CTYPE3",0);
00735     cpl_propertylist_erase_regexp(plist, "^CRPIX3",0);
00736     cpl_propertylist_erase_regexp(plist, "^CRVAL3",0);
00737     cpl_propertylist_erase_regexp(plist, "^CDELT3",0);
00738     cpl_propertylist_erase_regexp(plist, "^CUNIT3",0);
00739 
00740     cpl_propertylist_erase_regexp(plist, "^CTYPE2",0);
00741     cpl_propertylist_erase_regexp(plist, "^CRPIX2",0);
00742     cpl_propertylist_erase_regexp(plist, "^CRVAL2",0);
00743     cpl_propertylist_erase_regexp(plist, "^CDELT2",0);
00744     cpl_propertylist_erase_regexp(plist, "^CUNIT2",0);
00745 
00746 
00747     cpl_propertylist_erase_regexp(plist, "^CD1_1",0);
00748     cpl_propertylist_erase_regexp(plist, "^CD1_2",0);
00749     cpl_propertylist_erase_regexp(plist, "^CD2_1",0);
00750     cpl_propertylist_erase_regexp(plist, "^CD2_2",0);
00751 
00752     cpl_propertylist_erase_regexp(plist, "^CTYPE1",0);
00753     cpl_propertylist_insert_after_string(plist, "EXPTIME",  "CTYPE1", "PIXEL");
00754     cpl_propertylist_set_comment(plist, "CTYPE1", "Pixel coordinate system.");
00755 
00756     cpl_propertylist_erase_regexp(plist, "^CRPIX1",0);
00757     cpl_propertylist_insert_after_int(plist, "CTYPE1",  "CRPIX1", 1 ) ;
00758     cpl_propertylist_set_comment(plist, "CRPIX1", "Reference pixel in x") ;
00759 
00760     cpl_propertylist_erase_regexp(plist, "^CRVAL1",0);
00761     cpl_propertylist_insert_after_double(plist, "CRPIX1",  "CRVAL1", 1 ) ;
00762     cpl_propertylist_set_comment(plist, "CRVAL1", "value of ref pixel.") ;
00763 
00764     cpl_propertylist_erase_regexp(plist, "^CDELT1",0);
00765     cpl_propertylist_insert_after_double(plist, "CRVAL1",  "CDELT1", 1 ) ;
00766     cpl_propertylist_set_comment(plist, "CDELT1", "pixel scale") ;
00767 
00768 
00769 
00770 
00771 
00772     cpl_propertylist_erase_regexp(plist, "^CUNIT1",0);
00773     cpl_propertylist_insert_after_string(plist, "CDELT1",  "CUNIT1", "Pixel" );
00774     cpl_propertylist_set_comment(plist, "CUNIT1",  "spectral unit" );
00775 
00776     cpl_propertylist_erase_regexp(plist, "^CTYPE2",0);
00777     cpl_propertylist_insert_after_string(plist, "EXPTIME","CTYPE2","WAVE" );
00778     cpl_propertylist_set_comment(plist, "CTYPE2","wavelength axis in microns");
00779 
00780     cpl_propertylist_erase_regexp(plist, "^CRPIX2",0);
00781     cpl_propertylist_insert_after_int(plist, "CTYPE2",  "CRPIX2",cenpix ) ;
00782     cpl_propertylist_set_comment(plist, "CRPIX2", "Reference pixel in x") ;
00783 
00784     cpl_propertylist_erase_regexp(plist, "^CRVAL2",0);
00785     cpl_propertylist_insert_after_double(plist, "CRPIX2","CRVAL2",cenLambda ) ;
00786     cpl_propertylist_set_comment(plist, "CRVAL2", "central wavelength") ;
00787 
00788     cpl_propertylist_erase_regexp(plist, "^CDELT2",0);
00789     cpl_propertylist_insert_after_double(plist, "CRVAL2", "CDELT2",dispersion);
00790     cpl_propertylist_set_comment(plist, "CDELT2", "microns per pixel");
00791 
00792     cpl_propertylist_erase_regexp(plist, "^CUNIT2",0);
00793     cpl_propertylist_insert_after_string(plist, "CDELT2",  "CUNIT2", "mum");
00794     cpl_propertylist_set_comment(plist, "CUNIT2",  "spectral unit" );
00795 
00796     
00797 
00798 }   
00808 void sinfo_new_change_plist_image (cpl_propertylist * plist, 
00809             float center_x,
00810             float center_y )
00811 {
00812 
00813     float pixelscale ;
00814     double ra ;
00815     double dec ;
00816     double angle ;
00817     float radangle ;
00818     float cd1_1, cd1_2, cd2_1, cd2_2 ;
00819     char firsttext[2*FILE_NAME_SZ] ;
00820  
00821     strcpy(firsttext, "sinfo_rec_objnod -f \0") ;
00822     
00823     pixelscale = sinfo_pfits_get_pixscale(plist)/2. ;
00824     ra = sinfo_pfits_get_ra(plist) ;
00825     dec = sinfo_pfits_get_DEC(plist) ;
00826     angle = sinfo_pfits_get_posangle(plist) ;
00827     /* in PUPIL data there is not posangle info: we reset the error */
00828     if(cpl_error_get_code() != CPL_ERROR_NONE) {
00829        cpl_error_reset();
00830     }
00831 
00832     radangle = angle * PI_NUMB / 180. ;
00833     cd1_1 = cos(radangle)  * pixelscale / 3600. ;
00834     cd1_2 = sin(radangle) * pixelscale / 3600. ;
00835     cd2_1 = -sin(radangle)  * pixelscale / 3600. ;
00836     cd2_2 = cos(radangle)  * pixelscale / 3600. ;
00837 
00838     cpl_propertylist_erase_regexp(plist, "^CTYPE1",0);
00839     cpl_propertylist_insert_after_string(plist,"EXPTIME","CTYPE1","RA---TAN") ;
00840     cpl_propertylist_set_comment(plist, "CTYPE1", "Projected Rectascension" ) ;
00841 
00842     cpl_propertylist_erase_regexp(plist, "^CRPIX1",0);
00843     cpl_propertylist_insert_after_float(plist, "CTYPE1",  "CRPIX1", center_x) ;
00844     cpl_propertylist_set_comment(plist, "CRPIX1","Reference pixel in RA" ) ;
00845 
00846     cpl_propertylist_erase_regexp(plist, "^CRVAL1",0);
00847     cpl_propertylist_insert_after_double(plist, "CRPIX1",  "CRVAL1", ra ) ;
00848     cpl_propertylist_set_comment(plist, "CRVAL1","Reference RA" ) ;
00849 
00850     cpl_propertylist_erase_regexp(plist, "^CDELT1",0);
00851     cpl_propertylist_insert_after_float(plist,"CRVAL1","CDELT1", 
00852                -pixelscale/3600.  ) ;
00853     cpl_propertylist_set_comment(plist, "CDELT1","pixel scale" ) ;
00854 
00855     cpl_propertylist_erase_regexp(plist, "^CUNIT1",0);
00856     cpl_propertylist_insert_after_string(plist, "CDELT1",  "CUNIT1", "deg" ) ;
00857     cpl_propertylist_set_comment(plist, "CUNIT1","RA-UNIT" ) ;
00858 
00859     cpl_propertylist_erase_regexp(plist, "^CTYPE2",0);
00860     cpl_propertylist_insert_after_string(plist,"CUNIT1","CTYPE2","DEC--TAN") ;
00861     cpl_propertylist_set_comment(plist, "CTYPE2", "Projected Declination") ;
00862 
00863     cpl_propertylist_erase_regexp(plist, "^CRPIX2",0);
00864     cpl_propertylist_insert_after_float(plist, "CTYPE2",  "CRPIX2", center_y) ;
00865     cpl_propertylist_set_comment(plist, "CRPIX2", "Reference pixel in DEC") ;
00866 
00867     cpl_propertylist_erase_regexp(plist, "^CRVAL2",0);
00868     cpl_propertylist_insert_after_double(plist, "CRPIX2",  "CRVAL2", dec) ;
00869     cpl_propertylist_set_comment(plist, "CRVAL2", "Reference DEC") ;
00870 
00871     cpl_propertylist_erase_regexp(plist, "^CDELT2",0);
00872     cpl_propertylist_insert_after_float(plist, "CRVAL2",  "CDELT2", 
00873                 pixelscale/3600. ) ;
00874     cpl_propertylist_set_comment(plist, "CDELT2", "pixel scale") ;
00875 
00876     cpl_propertylist_erase_regexp(plist, "^CUNIT2",0);
00877     cpl_propertylist_insert_after_string(plist, "CDELT2",  "CUNIT2", "deg" ) ;
00878     cpl_propertylist_set_comment(plist, "CUNIT2", "DEC-UNIT") ;
00879     
00880     cpl_propertylist_erase_regexp(plist, "^CD1_1",0);
00881     cpl_propertylist_insert_after_float(plist, "CUNIT2",  "CD1_1", cd1_1 ) ;
00882     cpl_propertylist_set_comment(plist, "CD1_1",  "CD rotation matrix" ) ;
00883 
00884     cpl_propertylist_erase_regexp(plist, "^CD1_2",0);
00885     cpl_propertylist_insert_after_float(plist, "CD1_1",   "CD1_2", cd1_2 ) ;
00886     cpl_propertylist_set_comment(plist, "CD1_2",  "CD rotation matrix" ) ;
00887     
00888     cpl_propertylist_erase_regexp(plist, "^CD2_1",0);
00889     cpl_propertylist_insert_after_float(plist, "CD1_2",   "CD2_1", cd2_1 ) ;
00890     cpl_propertylist_set_comment(plist, "CD2_1",  "CD rotation matrix" ) ;
00891     
00892     cpl_propertylist_erase_regexp(plist, "^CD2_2",0);
00893     cpl_propertylist_insert_after_float(plist, "CD2_1",   "CD2_2", cd2_2 ) ;
00894     cpl_propertylist_set_comment(plist, "CD2_2",  "CD rotation matrix" ) ;
00895 
00896 
00897 }        
00898 
00906 cpl_imagelist** 
00907 sinfo_new_sinfoni_correct_median(cpl_imagelist** cubes, const int n_cubes)
00908 {
00909   int i=0;
00910   cpl_imagelist** cubes_cor=NULL;
00911   double local_median=0;
00912   int z=0;
00913   cpl_image* i_img=NULL;
00914   cpl_image* o_img=NULL;
00915 
00916 
00917   if ( cubes == NULL ) {
00918     sinfo_msg_error ("no cube list given!") ;
00919     return NULL ;
00920   }
00921   if ( n_cubes <= 0 ) {
00922     sinfo_msg_error ("wrong number of data cubes in list!") ;
00923     return NULL ;
00924   }
00925 
00926   cubes_cor = (cpl_imagelist**) cpl_calloc (n_cubes, sizeof (cpl_imagelist*)) ;
00927   
00928   for ( i = 0 ; i < n_cubes ; i++ ) {
00929     cubes_cor[i] = cpl_imagelist_new();
00930     for(z=0;z< cpl_imagelist_get_size(cubes[i]); z++) {
00931       i_img=cpl_imagelist_get(cubes[i],z); 
00932       local_median=cpl_image_get_median(i_img);;
00933       o_img=cpl_image_duplicate(i_img);
00934       if(!isnan(local_median)) {
00935         cpl_image_subtract_scalar(o_img,local_median);
00936        }
00937       cpl_imagelist_set(cubes_cor[i],o_img,z);
00938     }
00939   }
00940 
00941   return cubes_cor;
00942 }
00943 
00950 int sinfo_new_sinfoni_correct_median_it(cpl_imagelist** inp)
00951 {
00952 
00953   double local_median=0;
00954   int z=0;
00955   cpl_image* img=NULL;
00956 
00957   for(z=0;z< cpl_imagelist_get_size((*inp)); z++) {
00958     img=cpl_imagelist_get((*inp),z);
00959     local_median=sinfo_new_my_median_image(img);
00960     if(!isnan(local_median)) {
00961       cpl_image_subtract_scalar(img,local_median);
00962     }  else {
00963       sinfo_msg_error("local_median is NAN");
00964     }
00965     cpl_imagelist_set((*inp),img,z);
00966   }
00967 
00968   return 0;
00969 }
00970 
00980 cpl_imagelist** sinfo_new_sinfoni_correct_sky(cpl_imagelist** cubes,
00981                   const int nc,
00982                   cpl_imagelist* sky_cube)
00983 
00984 {
00985   cpl_imagelist** cubes_sky=NULL;
00986   int x=0;
00987   int y=0;
00988   int z=0;
00989   int i=0;
00990   float k=0.5;
00991   int ovr=0;
00992   int ks=0;
00993   int nclip=0;
00994   double med=0;
00995   double avg=0;
00996   double sig=0;
00997   int msk_sum=0;
00998   double val_msk_sum=0;
00999   cpl_vector* val=NULL;
01000   cpl_vector* msk=NULL;
01001   int ilx=0;
01002   int ily=0;
01003   int inp=0;
01004 
01005 
01006   int sky_lx=0;
01007   int sky_ly=0;
01008 
01009   float* pidata=NULL;
01010   float* p_sky_data=NULL;
01011   cpl_image* i_img=NULL;
01012   cpl_image* sky_img=NULL;
01013   
01014 
01015   if ( cubes == NULL ) {
01016     sinfo_msg_error ("no cube list given!") ;
01017     return NULL ;
01018   }
01019   ilx=cpl_image_get_size_x(cpl_imagelist_get(cubes[0],0));
01020   ily=cpl_image_get_size_y(cpl_imagelist_get(cubes[0],0));
01021   inp=cpl_imagelist_get_size(cubes[0]);
01022 
01023 
01024   sky_lx=cpl_image_get_size_x(cpl_imagelist_get(sky_cube,0));
01025   sky_ly=cpl_image_get_size_y(cpl_imagelist_get(sky_cube,0));
01026   if ( nc <= 0 ) {
01027     sinfo_msg_error ("wrong number of data cubes in list!") ;
01028     return NULL ;
01029   }
01030 
01031   cubes_sky = (cpl_imagelist**) cpl_calloc (nc, sizeof (cpl_imagelist*)) ;     
01032   for(z=0;z< inp; z++) {
01033     sky_img=cpl_imagelist_get(sky_cube,z);
01034     p_sky_data=cpl_image_get_data_float(sky_img);
01035     for(y=0;y< ily; y++) {
01036       for(x=0;x< ilx; x++) {
01037         /* here we start to do a k-s clipping */
01038     msk=cpl_vector_new(nc);
01039     for (i=0;i<nc;i++) {
01040       cpl_vector_set(msk,i,1);
01041     }
01042         nclip=0;
01043         for (ks=0;ks<nc;ks++) {
01044       sig=0;
01045       med=0;
01046       ovr=0;
01047       val=cpl_vector_new(nc-nclip);
01048 
01049       for ( i = 0 ; i < nc ; i++ ) {
01050             i_img=cpl_imagelist_get(cubes[i],z);
01051             pidata=cpl_image_get_data_float(i_img);
01052         if ((!isnan(pidata[x+y*ilx])) &&
01053                  (cpl_vector_get(msk,i) != 0) ) {
01054           cpl_vector_set(val,ovr,(double)pidata[x+y*ilx]); 
01055           ovr++;
01056         }
01057       }
01058 
01059       if(ovr>0) {
01060         avg=cpl_vector_get_mean(val);
01061         med=cpl_vector_get_median(val);
01062         if(ovr>1) {
01063           sig=cpl_vector_get_stdev(val);
01064         } else {
01065               sig=0;
01066         }
01067       } else {
01068         avg=cpl_vector_get(val,0);
01069         med=avg;
01070         sig=0;
01071       }
01072 
01073       cpl_vector_delete(val);
01074       for ( i = 0 ; i < nc ; i++ ) {
01075             i_img=cpl_imagelist_get(cubes[i],z);
01076             pidata=cpl_image_get_data_float(i_img);
01077       /* Do k-s clipping at each pixel */
01078         if ((!isnan(pidata[x+y*ilx])) && 
01079                 (cpl_vector_get(msk,i) != 0)) { 
01080           if(abs((pidata[x+y*ilx]-med))> k*sig) {
01081         cpl_vector_set(msk,i,0);
01082                 nclip++;
01083           }
01084         }
01085       }
01086     }/* end of k-s clipping */
01087     msk_sum=0;
01088     val_msk_sum=0;
01089     for ( i = 0 ; i < nc ; i++ ) {
01090       /* computes sky at each point */
01091       if (!isnan(pidata[x+y*ilx])) { 
01092         msk_sum+=cpl_vector_get(msk,i);
01093                 val_msk_sum+=pidata[x+y*ilx]*
01094           cpl_vector_get(msk,i);
01095       }
01096     }
01097     p_sky_data[x+y*sky_lx]=val_msk_sum/msk_sum;
01098     cpl_vector_delete(msk);
01099       } /* end loop over x */
01100     } /* end loop over y */
01101   } /* end loop over z */
01102   for ( i = 0 ; i < nc ; i++ ) {
01103     cubes_sky[i]=cpl_imagelist_duplicate(cubes[i]);
01104     /* subtract the variable clean sky */
01105     cpl_imagelist_subtract(cubes_sky[i],sky_cube);
01106 
01107   }
01108 
01109 
01110   return cubes_sky;
01111 }
01112 
01127 cpl_imagelist** sinfo_new_sinfoni_correct_sky2(cpl_imagelist** cubes,
01128                   const int nc,
01129                   cpl_imagelist* sky_cube,
01130                               cpl_imagelist* med_cube,
01131                               cpl_imagelist* msk_cube,
01132                               cpl_imagelist* avg_cube,
01133                               cpl_imagelist* sig_cube,
01134                               cpl_imagelist* ovr_cube)
01135 {
01136   cpl_imagelist** cubes_sky=NULL;
01137   int x=0;
01138   int y=0;
01139   int z=0;
01140   int i=0;
01141   float k=0.5;
01142   int ovr=0;
01143   int ks=0;
01144   int nclip=0;
01145   double med=0;
01146   double avg=0;
01147   double sig=0;
01148   int msk_sum=0;
01149   double val_msk_sum=0;
01150   cpl_vector* val=NULL;
01151   cpl_vector* msk=NULL;
01152   
01153   int ilx=0;
01154   int ily=0;
01155   int inp=0;
01156 
01157   int ovr_lx=0;
01158   int msk_lx=0;
01159   int avg_lx=0;
01160   int sig_lx=0;
01161   int sky_lx=0;
01162   int med_lx=0;
01163 
01164 
01165   float* p_ovr_data=NULL;
01166   float* p_msk_data=NULL;
01167   float* p_avg_data=NULL;
01168   float* p_sig_data=NULL;
01169   float* p_sky_data=NULL;
01170   float* p_med_data=NULL;
01171 
01172   float* pidata=NULL;
01173 
01174   if ( cubes == NULL ) {
01175     sinfo_msg_error ("no cube list given!") ;
01176     return NULL ;
01177   }
01178 
01179   ilx=cpl_image_get_size_x(cpl_imagelist_get(cubes[0],0));
01180   ily=cpl_image_get_size_y(cpl_imagelist_get(cubes[0],0));
01181   inp=cpl_imagelist_get_size(cubes[0]);
01182 
01183   if ( nc <= 0 ) {
01184     sinfo_msg_error ("wrong number of data cubes in list!") ;
01185     return NULL ;
01186   }
01187 
01188   cubes_sky = (cpl_imagelist**) cpl_calloc (nc, sizeof (cpl_imagelist*)) ;     
01189 
01190   ovr_lx=ilx;
01191   msk_lx=ilx;
01192   avg_lx=ilx;
01193   sig_lx=ilx;
01194   sky_lx=ilx;
01195   med_lx=ilx;
01196 
01197   for(z=0;z< inp; z++) {
01198     p_ovr_data=cpl_image_get_data_float(cpl_imagelist_get(ovr_cube,z));
01199     p_msk_data=cpl_image_get_data_float(cpl_imagelist_get(msk_cube,z));
01200     p_avg_data=cpl_image_get_data_float(cpl_imagelist_get(avg_cube,z));
01201     p_sig_data=cpl_image_get_data_float(cpl_imagelist_get(sig_cube,z));
01202     p_sky_data=cpl_image_get_data_float(cpl_imagelist_get(sky_cube,z));
01203     p_med_data=cpl_image_get_data_float(cpl_imagelist_get(med_cube,z));
01204 
01205 
01206     for(y=0;y< ily; y++) {
01207       for(x=0;x< ilx; x++) {
01208         /* here we start to do a k-s clipping */
01209     msk=cpl_vector_new(nc);
01210     for (i=0;i<nc;i++) {
01211       cpl_vector_set(msk,i,1);
01212     }
01213     p_ovr_data[x+y*ovr_lx]=nc;
01214     p_msk_data[x+y*msk_lx]=nc;
01215         nclip=0;
01216         for (ks=0;ks<nc;ks++) {
01217       sig=0;
01218       med=0;
01219       ovr=0;
01220       val=cpl_vector_new(nc-nclip);
01221 
01222       for ( i = 0 ; i < nc ; i++ ) {
01223             pidata=cpl_image_get_data_float(cpl_imagelist_get(cubes[i],z));
01224         if ((!isnan(pidata[x+y*ilx])) &&
01225                  (cpl_vector_get(msk,i) != 0) ) {
01226           cpl_vector_set(val,ovr,(double)pidata[x+y*ilx]); 
01227           ovr++;
01228         }
01229       }
01230 
01231       if(ovr>1) {
01232         avg=cpl_vector_get_mean(val);
01233         med=cpl_vector_get_median(val);
01234         sig=cpl_vector_get_stdev(val);
01235       } else {
01236         avg=cpl_vector_get(val,0);
01237         med=avg;
01238         sig=0;
01239       }
01240 
01241       p_med_data[x+y*med_lx]=med;
01242       p_avg_data[x+y*avg_lx]=avg;
01243       p_sig_data[x+y*sig_lx]=sig;
01244       cpl_vector_delete(val);
01245       for ( i = 0 ; i < nc ; i++ ) {
01246            pidata=cpl_image_get_data_float(cpl_imagelist_get(cubes[i],z));
01247       /* Do k-s clipping at each pixel */
01248         if ((!isnan(pidata[x+y*ilx])) && 
01249                 (cpl_vector_get(msk,i) != 0)) { 
01250           if(abs((pidata[x+y*ilx]-med))> k*sig) {
01251         /* pidata[x+y*ilx]=0; */
01252         p_msk_data[x+y*msk_lx]-=1;
01253         cpl_vector_set(msk,i,0);
01254                 nclip++;
01255           }
01256         }
01257       }
01258     }/* end of k-s clipping */
01259     msk_sum=0;
01260     val_msk_sum=0;
01261     for ( i = 0 ; i < nc ; i++ ) {
01262            pidata=cpl_image_get_data_float(cpl_imagelist_get(cubes[i],z));
01263       /* computes sky at each point */
01264       if (!isnan(pidata[x+y*ilx])) { 
01265         /*
01266         msk_sum+=p_msk_data[x+y*msk_lx];
01267                 val_msk_sum+=pidata[x+y*cubes[i]->lx]*
01268           p_msk_data[x+y*msk_lx];
01269         */
01270         msk_sum+=cpl_vector_get(msk,i);
01271                 val_msk_sum+=pidata[x+y*ilx]*
01272           cpl_vector_get(msk,i);
01273       }
01274     }
01275     p_sky_data[x+y*sky_lx]=val_msk_sum/msk_sum;
01276     cpl_vector_delete(msk);
01277       } /* end loop over x */
01278     } /* end loop over y */
01279   } /* end loop over z */
01280   for ( i = 0 ; i < nc ; i++ ) {
01281     cubes_sky[i]=cpl_imagelist_duplicate(cubes[i]);
01282     /* subtract the variable clean sky */
01283     cpl_imagelist_subtract(cubes_sky[i],sky_cube);
01284   }
01285 
01286 
01287   return cubes_sky;
01288 }
01298 int
01299 sinfo_new_extract_ao_info(cpl_frameset* obj_set,
01300                           cpl_frameset** set,
01301                           const char* recid,
01302                           cpl_parameterlist* config)
01303 {
01304 
01305   cpl_frame* frm=NULL;
01306   cpl_table* tbl=NULL;
01307   cpl_frame* product_frame=NULL;
01308   cpl_propertylist* plist=NULL;
01309   cpl_image* image=NULL;
01310 
01311   int nobj=0;
01312   int i=0;
01313   const char* name_o="out_ao_info.fits";
01314 
01315   name_o="out_ao_info.fits";
01316   nobj=cpl_frameset_get_size(obj_set);
01317   for(i=0;i<nobj;i++) {
01318     frm=cpl_frameset_get_frame(obj_set,i);
01319     tbl=cpl_table_load(cpl_frame_get_filename(frm),1,1);
01320     plist=cpl_propertylist_load(cpl_frame_get_filename(frm),0);
01321     if(i==0) {
01322       cpl_table_save(tbl,NULL,plist,name_o,CPL_IO_DEFAULT);
01323     } else {
01324       cpl_table_save(tbl,NULL,plist,name_o,CPL_IO_EXTEND);
01325     }
01326   }
01327   image=cpl_image_new(0,0,CPL_TYPE_FLOAT);
01328     /* Create product frame */
01329     product_frame = cpl_frame_new();
01330     cpl_frame_set_filename(product_frame, name_o) ;
01331     cpl_frame_set_tag(product_frame, PRO_AO_INFO) ;
01332     cpl_frame_set_type(product_frame, CPL_FRAME_TYPE_IMAGE) ;
01333     cpl_frame_set_group(product_frame, CPL_FRAME_GROUP_PRODUCT) ;
01334     cpl_frame_set_level(product_frame, CPL_FRAME_LEVEL_FINAL) ;
01335     if (cpl_error_get_code()) {
01336         sinfo_msg_error( "Error while initialising the product frame") ;
01337         cpl_propertylist_delete(plist) ;
01338         cpl_frame_delete(product_frame) ;
01339         cpl_image_delete(image) ;
01340         return -1 ;
01341     }
01342     
01343      /* Add DataFlow keywords */
01344     if (cpl_dfs_setup_product_header(plist, product_frame, *set, config,
01345             recid, "SINFONI", "?Dictionary?") != CPL_ERROR_NONE) {
01346         sinfo_msg_error( "Problem in the product DFS-compliance") ;
01347         cpl_propertylist_delete(plist) ;
01348         cpl_frame_delete(product_frame) ;
01349         cpl_image_delete(image) ;
01350         return -1 ;
01351     }
01352 
01353      /* Save the file */
01354     if (cpl_image_save(image, name_o, CPL_BPP_DEFAULT, plist,
01355                        CPL_IO_DEFAULT) != CPL_ERROR_NONE) {
01356         sinfo_msg_error( "Could not save product");
01357         cpl_propertylist_delete(plist) ;
01358         cpl_frame_delete(product_frame) ;
01359         cpl_image_delete(image) ;
01360         return -1 ;
01361     }
01362      cpl_propertylist_delete(plist) ;
01363      cpl_image_delete(image) ;
01364  
01365     /* Log the saved file in the input frameset */
01366     cpl_frameset_insert(*set, product_frame) ;
01367  
01368   return 0;
01369 
01370 }
01383 int
01384 sinfo_new_assign_offset(const int n,
01385                         const char* name,
01386                         float* offsetx,
01387                         float* offsety,
01388                         const float ref_offx,
01389                         const float ref_offy)
01390 {
01391 
01392   float offx=0;
01393   float offy=0;
01394   double mjd_obs=0;
01395   /*
01396   double pixelscale=0;
01397   double angle=0;
01398   double radangle=0;
01399   double cd1_1=0;
01400   double cd1_2=0;
01401   double cd2_1=0;
01402   double cd2_2=0;
01403   double ra=0;
01404   double dec=0;
01405   */
01406 
01407   cpl_propertylist * plist=NULL;
01408   sinfo_msg_debug("Assign offsets");
01409   
01410   if ((cpl_error_code)((plist = cpl_propertylist_load(name, 0)) == NULL)) {
01411     sinfo_msg_error( "getting header from reference frame %s",name);
01412     cpl_propertylist_delete(plist) ;
01413     return -1 ;
01414   }
01415 
01416   offx = sinfo_pfits_get_cumoffsetx(plist) - ref_offx ;  /* was - */
01417   if(cpl_error_get_code() != CPL_ERROR_NONE) {
01418     sinfo_msg_warning(" could not read fits header keyword cummoffsetx!") ;
01419     sinfo_msg_warning(" Set relative offset to 0 - %f!",ref_offx) ;
01420     offx =  - ref_offx;
01421     cpl_error_reset();
01422     /* return -1 ; */
01423   }
01424     
01425   offy = sinfo_pfits_get_cumoffsety(plist) - ref_offy ; /* was - */
01426   if(cpl_error_get_code() != CPL_ERROR_NONE) {
01427     sinfo_msg_warning(" could not read fits header keyword! cumoffsety") ;
01428     sinfo_msg_warning(" Set relative offset to 0 - %f!",ref_offx) ;
01429     offy =  - ref_offy;
01430     cpl_error_reset();
01431     /* return -1 ; */
01432   }
01433   sinfo_msg_debug("offx=%f offy=%f",offx,offy);
01434 
01435   if (cpl_propertylist_contains(plist, KEY_NAME_MJD_OBS)) {
01436     mjd_obs=cpl_propertylist_get_double(plist, KEY_NAME_MJD_OBS);
01437   } else {
01438     sinfo_msg_error("keyword %s does not exist",KEY_NAME_MJD_OBS);
01439     cpl_propertylist_delete(plist) ;
01440     return -1;
01441   }
01442   cpl_propertylist_delete(plist) ;
01443 
01444   if (mjd_obs > 53825. ) {
01445     /* April 1st 2006 */
01446     sinfo_msg("New cumoffset setting convention");
01447     sinfo_new_array_set_value(offsetx,2*offx,n);
01448     sinfo_new_array_set_value(offsety,2*offy,n);
01449   } else if ((mjd_obs > 53421.58210082 ) && (mjd_obs <= 53825.)){
01450     /* after detector's upgrade */
01451     /*
01452     sinfo_new_array_set_value(offsetx,-offx*2,n);
01453     sinfo_new_array_set_value(offsety,+offy*2,n);
01454     */
01455     sinfo_new_array_set_value(offsetx,-2*offx,n);
01456     sinfo_new_array_set_value(offsety,2*offy,n);
01457   } else {
01458     /* before detector's upgrade */
01459     /*
01460     sinfo_new_array_set_value(offsetx,+offx*2,n);
01461     sinfo_new_array_set_value(offsety,-offy*2,n);
01462     */
01463     sinfo_new_array_set_value(offsetx,2*offx,n);
01464     sinfo_new_array_set_value(offsety,-2*offy,n);
01465   }
01466 
01467   return 0;
01468 
01469 
01470 }
01483 int
01484 sinfo_new_object_assign_offset(const char* name, 
01485                                const int n, 
01486                                double* ref_offx, 
01487                                double* ref_offy, 
01488                                float** offsetx, 
01489                                float** offsety)
01490 {
01491 
01492   float offx=0;
01493   float offy=0;
01494   double mjd_obs=0;
01495   cpl_propertylist * plist=NULL;
01496   sinfo_msg_debug("Assign offsets");
01497 
01498   if ((cpl_error_code)((plist = cpl_propertylist_load(name, 0)) == NULL)) {
01499     sinfo_msg_error( "getting header from reference frame %s",name);
01500     cpl_propertylist_delete(plist) ;
01501     return -1 ;
01502   }
01503   if ( n == 0 ) {
01504 
01505     *ref_offx = sinfo_pfits_get_cumoffsetx(plist) ;
01506     if(cpl_error_get_code() != CPL_ERROR_NONE) {
01507       sinfo_msg_error(" could not read fits header keyword cummoffsetx!") ;
01508       /* return -1 ; */
01509       cpl_error_reset();
01510     }
01511     
01512     *ref_offy = sinfo_pfits_get_cumoffsety(plist) ;
01513     if(cpl_error_get_code() != CPL_ERROR_NONE) {
01514       sinfo_msg_error(" could not read fits header keyword! cumoffsety") ;
01515       cpl_error_reset();
01516        /* return -1 ; */
01517     }
01518     sinfo_msg_debug("Reference offx=%f offy=%f",*ref_offx,*ref_offy);
01519           
01520     offx = 0. ;
01521     offy = 0. ;
01522 
01523   } else {
01524        
01525     offx = sinfo_pfits_get_cumoffsetx(plist) - *ref_offx ;  /* was - */
01526     if(cpl_error_get_code() != CPL_ERROR_NONE) {
01527       sinfo_msg_error(" could not read fits header keyword cummoffsetx!") ;
01528       cpl_error_reset();
01529       /* return -1 ; */
01530     }
01531     
01532     offy = sinfo_pfits_get_cumoffsety(plist) - *ref_offy ; /* was - */
01533     if(cpl_error_get_code() != CPL_ERROR_NONE) {
01534       sinfo_msg_error(" could not read fits header keyword! cumoffsety") ;
01535       /* return -1 ; */
01536       cpl_error_reset();
01537     }
01538     sinfo_msg_debug("offx=%f offy=%f",offx,offy);               
01539   }
01540            
01541 
01542     /* rotate the coordinates 
01543        offx_rot = cd1_1 * offx + cd2_1 * offy ;
01544        offy_rot = cd1_2 * offx + cd2_2 * offy ;
01545         convert the coordinates to pixel units 
01546        offx_rot_pix = offx_rot / pixelscale ;
01547        offy_rot_pix = offy_rot / pixelscale ;
01548        offsetx[i] = offx_rot_pix ;
01549        offsety[i] = offy_rot_pix ;
01550     */
01551 
01552 
01553 
01554   if (cpl_propertylist_contains(plist, KEY_NAME_MJD_OBS)) {
01555     mjd_obs=cpl_propertylist_get_double(plist, KEY_NAME_MJD_OBS);
01556   } else {
01557     sinfo_msg_error("keyword %s does not exist",KEY_NAME_MJD_OBS);
01558     cpl_propertylist_delete(plist) ;
01559     return -1;
01560   }
01561   cpl_propertylist_delete(plist) ;
01562 
01563   if (mjd_obs > 53825. ) {
01564     /* April 1st 2006 */
01565     sinfo_msg("New cumoffset setting convention");
01566     sinfo_new_array_set_value(*offsetx,2*offx,n);
01567     sinfo_new_array_set_value(*offsety,2*offy,n);
01568   } else if ((mjd_obs > 53421.58210082 ) && (mjd_obs <= 53825.)){
01569     /* after detector's upgrade */
01570     sinfo_new_array_set_value(*offsetx,-offx*2,n);
01571     sinfo_new_array_set_value(*offsety,+offy*2,n);
01572   } else {
01573     /* before detector's upgrade */
01574     sinfo_new_array_set_value(*offsetx,+offx*2,n);
01575     sinfo_new_array_set_value(*offsety,-offy*2,n);
01576   }
01577 
01578   return 0;
01579 }
01580 
01592 cpl_imagelist*
01593 sinfo_new_fine_tune(cpl_imagelist* cube,
01594               float* correct_dist, 
01595               const char* method, 
01596               const int order, 
01597               const int nslits) {
01598   int i =0;
01599   cpl_imagelist* outcube2=NULL;
01600   float* neg_dist=NULL;
01601   sinfo_msg("Finetuning, method=%s",method);
01602   
01603   if (strcmp(method,"P")==0)
01604     {
01605       outcube2 = sinfo_new_fine_tune_cube( cube, correct_dist, order ) ;
01606       if (outcube2 == NULL)
01607     {
01608       sinfo_msg_error (" could not fine tune the data cube\n") ;
01609       return NULL ;
01610     }
01611     }       
01612   else if (strcmp(method,"F")==0)
01613     {
01614       neg_dist=cpl_calloc(nslits,sizeof(float));
01615       for ( i = 0 ; i < nslits ; i++ )
01616     {
01617       neg_dist[i] = -correct_dist[i] ;
01618     }
01619       outcube2 = sinfo_new_fine_tune_cube_by_FFT( cube, neg_dist ) ;
01620      cpl_free(neg_dist);
01621       if ( outcube2 == NULL )
01622     {
01623       sinfo_msg_error (" could not fine tune the data cube\n") ;
01624       return NULL ;
01625     }
01626     }       
01627   else if (strcmp(method,"S")==0)
01628     {
01629       outcube2 = sinfo_new_fine_tune_cube_by_spline( cube, correct_dist ) ;
01630       if ( outcube2 == NULL )
01631     {
01632       sinfo_msg_error (" could not fine tune the data cube\n") ;
01633       return NULL ;
01634     }
01635     }
01636   else
01637     {
01638       sinfo_msg_error (" wrong method indicator given!") ;
01639       return NULL ;
01640     }
01641 
01642 
01643 
01644 return outcube2;
01645 
01646 }

Generated on Wed Jan 17 08:33:44 2007 for SINFONI Pipeline Reference Manual by  doxygen 1.4.4