spectrum_ops.c

00001 /*******************************************************************************
00002 * E.S.O. - VLT project
00003 *
00004 *
00005 *
00006 * who       when      what
00007 * --------  --------  ----------------------------------------------
00008 * schreib  25/05/00  created
00009 */
00010 
00011 /************************************************************************
00012 *   NAME
00013 *        spectrum_ops.c -
00014 *        some vector procedures to operate on spectra
00015 *
00016 *   SYNOPSIS
00017 *   #include "spectrum_ops.h"
00018 *
00019 *   1) Vector * newVector( ulong32 n_elements )
00020 *   2) void * destroyVector( Vector *vector )
00021 *   3) OneImage * vectorToImage( Vector * spectrum )
00022 *   4) Vector * imageToVector( OneImage * spectrum )
00023 *   5) OneImage * extractSpectrumFromResampledFlat( OneImage * resflat,
00024 *                                                   float      loreject,
00025 *                                                   float      hireject ) 
00026 *   6) OneImage * multiplyImageWithSpectrum( OneImage * image, OneImage * spectrum )
00027 *   7) OneImage * optimalExtractionFromCube( OneCube * cube, 
00028 *                                            int       halfbox_x, 
00029 *                                            int       halfbox_y,
00030 *                                            float     fwhm_factor,
00031 *                                            float     backvariance,
00032 *                                            float     sky,
00033 *                                            float     gain,
00034 *                                            float     exptime)
00035 *   8) Vector * extractSkyFromCube( OneCube * cube,
00036 *                                   float     loReject,
00037 *                                   float     hiReject,
00038 *                                   int     * position,
00039 *                                   int       tolerance,
00040 *                                   int       posindicator )
00041 *    9) Vector * sumRectangleOfCubeSpectra( OneCube * cube,
00042 *                                     int llx,
00043 *                                     int lly,
00044 *                                     int urx,
00045 *                                     int ury )
00046 *   10) Vector * sumCircleOfCubeSpectra( OneCube * cube,
00047 *                                        int       centerx,
00048 *                                        int       centery,
00049 *                                        int       radius )
00050 *   11) Vector * meanRectangleOfCubeSpectra( OneCube * cube,
00051 *                                            int llx,
00052 *                                            int lly,
00053 *                                            int urx,
00054 *                                            int ury )
00055 *   12) Vector * meanCircleOfCubeSpectra( OneCube * cube,
00056 *                                         int       centerx,
00057 *                                         int       centery,
00058 *                                         int       radius )
00059 *   13) Vector * blackbodySpectrum( char * templateSpec, double temp )
00060 *   14) Vector * medianRectangleOfCubeSpectra( OneCube * cube,
00061 *                                               int llx,
00062 *                                               int lly,
00063 *                                               int urx,
00064 *                                               int ury )
00065 *   15) Vector * medianCircleOfCubeSpectra( OneCube * cube,
00066 *                                           int       centerx,
00067 *                                           int       centery,
00068 *                                           int       radius )
00069 *   16) Vector * cleanmeanRectangleOfCubeSpectra( OneCube * cube,
00070 *                                                 int llx,
00071 *                                                 int lly,
00072 *                                                 int urx,
00073 *                                                 int ury,
00074 *                                                 float lo_reject,
00075 *                                                 float hi_reject )
00076 *   17) Vector * cleanmeanCircleOfCubeSpectra( OneCube * cube,
00077 *                                              int       centerx,
00078 *                                              int       centery,
00079 *                                              int       radius,
00080 *                                              float     lo_reject,
00081 *                                              float     hi_reject )
00082 *   18) float * shiftArray ( float * input, int n_elements, float shift, double * ker ) 
00083 *
00084 *   DESCRIPTION
00085 *   1) allocates memory for a new vector
00086 *   2) frees memory of a vector
00087 *   3) converts a spectral vector to a fits image
00088 *      remark: vector object spectrum is destroyed
00089 *   4) converts a fits image to a spectral vector
00090 *      remark: input image is destroyed
00091 *   5) builds one spectrum in a fits image out of a resampled
00092 *      flatfield frame by taking a clean mean along the spatial pixels
00093 *   6) multiplys a resampled image with a resampled spectrum
00094 *      (calibrated halogen lamp spectrum) in the same spectral range
00095 *      that means all image columns are multiplied with the same spectrum
00096 *   7) does the optimal extraction of a standard star spectrum
00097 *      according to the equation:
00098 *       S = sum { (P^2 / V) * (I - B) / P } / sum{ P^2 / V } 
00099 *       S: spectral flux at a particular wavelength
00100 *       P: normalized PSF (determined by a 2D-Gaussian fit)
00101 *       I: pixel value
00102 *       B: background pixel value determined by the background parameter of the 2D-Gaussian fit
00103 *       V: estimated variance of a pixel: V = [R^2 + D + sky + I,c/exptime]/gain
00104 *          where R is the read noise, and D the dark current variance.
00105 *          backvariance is R^2 + D in counts/sec. I,c is the source intensity in counts
00106 *          Remember: sigma,e[e-] = gain[e/count] * sigma,c [counts] = sqrt(I,e) = sqrt(gain*I,c)
00107 *          => V,c = sigma,c^2 = sigma,e^2/gain^2 
00108 *          => sigma,c = sqrt(I,c/gain) => V,c = I,c/gain
00109 *   8) extracts a sky spectrum from a reduced sky spider observation, that
00110 *      means from a data cube. Therefore, the position of the sky within the
00111 *      field of view must be first read from the fits header.  A pixel tolerance is subtracted.
00112 *      The found sky spectra are averaged by rejecting the extreme high and low values.
00113 *   9) summing routine for a reduced data to get a better spectral S/N
00114 *      only for a rectangular aperture.
00115 *   10) summing routine for a reduced data to get a better spectral S/N
00116 *       only for a circular aperture.
00117 *   11) averaging routine for a reduced data to get a better spectral S/N
00118 *       only for a rectangular aperture.
00119 *   12) averaging routine for a reduced data to get a better spectral S/N
00120 *       only for a circular aperture.
00121 *   13) computes a blackbody spectral intensity distribution
00122 *       (W/(m^2 lambda ster)) 
00123 *   14) median routine for a reduced data to get a better spectral S/N
00124 *       only for a rectangular aperture.
00125 *   15) median routine for a reduced data to get a better spectral S/N
00126 *       only for a circular aperture.
00127 *   16) clean averaging routine for a reduced data to get a better spectral S/N
00128 *       only for a rectangular aperture.
00129 *   17) clean averaging routine for a reduced data to get a better spectral S/N
00130 *       only for a circular aperture.
00131 *   18) shifts an array by a sub-pixel shift value using a tanh
00132 *       interpolation kernel
00133 *
00134 *   FILES
00135 *
00136 *   ENVIRONMENT
00137 *
00138 *   RETURN VALUES
00139 *
00140 *   CAUTIONS
00141 *
00142 *   EXAMPLES
00143 *
00144 *   SEE ALSO
00145 *
00146 *   BUGS
00147 *
00148 *------------------------------------------------------------------------
00149 */
00150 
00151 #define POSIX_SOURCE 1
00152 #include "vltPort.h"
00153 
00154 /*
00155  * System Headers
00156  */
00157 
00158 /*
00159  * Local Headers
00160  */
00161 
00162 #include "spectrum_ops.h"
00163 
00164 /*----------------------------------------------------------------------------
00165  *                          Function codes
00166  *--------------------------------------------------------------------------*/
00167 
00168 /*----------------------------------------------------------------------------
00169    Function     :       newVector()
00170    In           :       number of vector elements
00171    Out          :       Vector 
00172    Job          :       allocates memory for a new vector
00173  ---------------------------------------------------------------------------*/
00174 
00175 Vector * newVector( ulong32 n_elements )
00176 {
00177     Vector * new_vector ;
00178 
00179     if ( n_elements <= 0 )
00180     {
00181         cpl_msg_error ("newVector:"," wrong number of elements\n") ;
00182         return NullVector ;
00183     }
00184     
00185     /* allocate memory for a vector with the given number of elements */
00186     new_vector = (Vector *) cpl_malloc (sizeof (Vector)) ; 
00187     new_vector -> n_elements = n_elements ;
00188     new_vector -> data = (pixelvalue *) cpl_calloc (n_elements, sizeof (pixelvalue)) ;
00189 
00190     return new_vector ;
00191 }
00192 
00193 
00194 /*----------------------------------------------------------------------------
00195    Function     :       destroyVector()
00196    In           :       vector to destroy
00197    Out          :       nothing
00198    Job          :       frees memory of a vector
00199  ---------------------------------------------------------------------------*/
00200 
00201 void destroyVector( Vector *vector )
00202 {
00203     if ( vector == NULL )   
00204     {
00205         cpl_msg_error("destroyVector:"," NULL Vector given!\n") ;
00206         return ;
00207     }    
00208     
00209     cpl_free ( vector -> data ) ;
00210     cpl_free ( vector ) ;
00211 }
00212 
00213 /*----------------------------------------------------------------------------
00214    Function     :       vectorToImage()
00215    In           :       spectral vector that should be converted to a fits image
00216    Out          :       image with lx = 1 and ly = length of vector
00217    Job          :       converts a spectral vector to a fits image
00218                         remark: vector object spectrum is destroyed
00219  ---------------------------------------------------------------------------*/
00220 
00221 OneImage * vectorToImage( Vector * spectrum )
00222 {
00223     OneImage * returnIm ;
00224     int i ;
00225 
00226     if ( spectrum == NULL )
00227     {
00228         cpl_msg_error("vectorToImage:"," no spectrum given!\n") ;
00229         return NullImage ;
00230     }
00231 
00232     /* allocate memory */
00233     if ( NullImage == (returnIm = new_image(1, spectrum->n_elements)) )
00234     {
00235         cpl_msg_error("vectorToImage:"," no spectrum given!\n") ;
00236         destroyVector(spectrum) ;
00237         return NullImage ;
00238     }
00239 
00240     for ( i = 0 ; i < spectrum->n_elements ; i++ )
00241     {
00242         returnIm -> data[i] = spectrum -> data[i] ;
00243     }
00244 
00245     destroyVector (spectrum) ;
00246     return returnIm ;
00247 }
00248 
00249 /*----------------------------------------------------------------------------
00250    Function     :       imageToVector()
00251    In           :       1-D Fits image that should be converted to a spectral vector
00252    Out          :       spectral vector with length ly 
00253    Job          :       converts a fits image to a spectral vector
00254                         remark: input image is destroyed
00255  ---------------------------------------------------------------------------*/
00256 
00257 Vector * imageToVector( OneImage * spectrum )
00258 {
00259     Vector * returnVector ;
00260     int i ;
00261 
00262     if ( spectrum == NULL )
00263     {
00264         cpl_msg_error("imageToVector:"," no spectrum given!\n") ;
00265         return NULL ;
00266     }
00267 
00268     /* allocate memory */
00269     if ( NULL == (returnVector = newVector(spectrum->nbpix)) )
00270     {
00271         cpl_msg_error("imageToVector:"," cannot allocate memory!\n") ;
00272         destroy_image(spectrum) ;
00273         return NULL ;
00274     }
00275 
00276     for ( i = 0 ; i < (int) spectrum->nbpix ; i++ )
00277     {
00278         returnVector -> data[i] = spectrum -> data[i] ;
00279     }
00280 
00281     destroy_image (spectrum) ;
00282     return returnVector ;
00283 }
00284 
00285 /*----------------------------------------------------------------------------
00286    Function     :       extractSpectrumFromResampledFlat()
00287    In           :       resflat: resampled halogen lamp frame, bad pixel corrected
00288                         loreject, hireject: percentage of extreme low and high intensity values
00289                                             to be rejected from averaging
00290    Out          :       fits image that contains the final halogen lamp spectrum
00291    Job          :       builds one spectrum in a fits image out of a resampled
00292                         flatfield frame by taking a clean mean along the spatial pixels
00293  ---------------------------------------------------------------------------*/
00294 
00295 OneImage * extractSpectrumFromResampledFlat( OneImage * resflat,
00296                                              float      loreject,
00297                                              float      hireject ) 
00298 {
00299     OneImage * retIm ;
00300     int col, row ;
00301     int n ;
00302     float array[resflat->ly] ;
00303     float cleanMean ;
00304     Vector * spectrum ;
00305 
00306     if ( resflat == NullImage )
00307     {
00308         cpl_msg_error("extractSpectrumFromResampledFlat:"," no flatfield given!\n") ;
00309         return NullImage ;
00310     }
00311 
00312     /* allocate memory */
00313     if ( NullVector == (spectrum = newVector(resflat->ly) ) )
00314     {
00315         cpl_msg_error("extractSpectrumFromResampledFlat:"," could not allocate memory!\n") ;
00316         return NullImage ;
00317     }
00318 
00319 
00320     for ( row = 0 ; row < resflat -> ly ; row++ )
00321     {
00322         n = 0 ;
00323         for ( col = 0 ; col < resflat -> lx ; col++ )
00324         {
00325             if ( !qfits_isnan(resflat -> data[col + row*resflat->lx]) )
00326             {
00327                 array[n] = resflat->data[col+row*resflat->lx] ;
00328                 n++ ;
00329             }
00330         }
00331         if ( n == 0 )
00332         {
00333             cpl_msg_warning("extractSpectrumFromResampledFlat:"," only bad pixels in row: %d!\n", row) ;
00334             cleanMean = ZERO ;
00335         }
00336         else
00337         {
00338             if ( FLT_MAX == (cleanMean = clean_mean(array, n, loreject, hireject)) )
00339             {
00340                 cpl_msg_error("extractSpectrumFromResampledFlat:"," could not do clean_mean!\n") ;
00341                 destroyVector(spectrum) ;
00342                 return NullImage ;
00343             }
00344         }
00345         spectrum->data[row] = cleanMean ; 
00346     }
00347     if ( NullImage == ( retIm = vectorToImage( spectrum ) ) )
00348     {
00349         cpl_msg_error("extractSpectrumFromResampledFlat:"," could not do vectorToImage!\n") ;
00350         destroyVector(spectrum) ;
00351         return NullImage ;
00352     }
00353     return retIm ;
00354 }
00355 
00356 /*----------------------------------------------------------------------------
00357    Function     :       multiplyImageWithSpectrum()
00358    In           :       image: resampled image
00359                         spectrum: resampled spectrum in image format
00360    Out          :       resulting image
00361    Job          :       multiplys a resampled image with a resampled spectrum
00362                         (calibrated halogen lamp spectrum) in the same spectral range
00363                         that means all image columns are multiplied with the same spectrum
00364  ---------------------------------------------------------------------------*/
00365 
00366 OneImage * multiplyImageWithSpectrum( OneImage * image, OneImage * spectrum )
00367 {
00368     int col, row ;
00369     OneImage * retImage ;
00370 
00371     if ( image == NullImage )
00372     {
00373         cpl_msg_error("multiplyImageWithSpectrum:"," no image given!\n") ;
00374         return NullImage ;
00375     }
00376 
00377     if ( spectrum == NullImage )
00378     {
00379         cpl_msg_error("multiplyImageWithSpectrum:"," no spectrum image given!\n") ;
00380         return NullImage ;
00381     }
00382     if ( spectrum -> ly != image -> ly )
00383     {
00384         cpl_msg_error("multiplyImageWithSpectrum:"," images are not compatible in pixel length!\n") ;
00385         return NullImage ;
00386     }
00387 
00388     if ( NullImage == (retImage = copy_image(image)) )
00389     {
00390         cpl_msg_error("multiplyImageWithSpectrum:"," could not copy original image!\n") ;
00391         return NullImage ;
00392     }
00393 
00394     for ( col = 0 ; col < image -> lx ; col++ )
00395     {
00396         for ( row = 0 ; row < image -> ly ; row++ )
00397         {
00398             if ( !qfits_isnan(image->data[col+row*image->lx]) &&
00399                  !qfits_isnan(spectrum->data[col+row*image->lx]))
00400             {
00401                 retImage->data[col+row*image->lx] = 
00402                     image->data[col+row*image->lx] * spectrum->data[row] ;
00403                                                    
00404             }
00405         }
00406     }
00407     return retImage ;
00408 }
00409 
00410 /*----------------------------------------------------------------------------
00411    Function     :       optimalExtractionFromCube()
00412    In           :       cube: input data cube
00413                         llx, lly: lower left edge points of the 2d Gaussian fitting box
00414                         halfbox_x, halfbox_y: half width of a box inside which
00415                                               a 2D-Gaussian fit is carried out
00416                         fwhm_factor: factor applied to the found fwhms of a 2D-Gaussian
00417                                      fit, defines the radius of the aperture inside which the 
00418                                      spectral extraction is carried out (default: 0.6).
00419                         backvariance: (readnoise^2 + dark current variance) needed to determine 
00420                                       the noise variance of the background. Must be given
00421                                       in counts/sec.
00422                         sky:         estimated sky variance in counts/sec
00423                         gain:        conversion factor electrons/count
00424                         exptime:     total exposure time
00425    Out          :       resulting spectrum stored in a 1D-image
00426    Job          :       does the optimal extraction of a standard star spectrum
00427                         according to the equation:
00428                         S = sum { (P^2 / V) * (I - B) / P } / sum{ P^2 / V } 
00429                         S: spectral flux at a particular wavelength
00430                         P: normalized PSF (determined by a 2D-Gaussian fit)
00431                         I: pixel value
00432                         B: background pixel value determined by the background parameter of the 2D-Gaussian fit
00433                         V: estimated variance of a pixel: V = [R^2 + D + sky + I,c/exptime]/gain
00434                            where R is the read noise, and D the dark current variance.
00435                            backvariance is R^2 + D in counts/sec. I,c is the source intensity in counts
00436                         Remember: sigma,e[e-] = gain[e/count] * sigma,c [counts] = sqrt(I,e) = sqrt(gain*I,c)
00437                         => V,c = sigma,c^2 = sigma,e^2/gain^2 
00438                         => sigma,c = sqrt(I,c/gain) => V,c = I,c/gain
00439  ---------------------------------------------------------------------------*/
00440 
00441 OneImage * optimalExtractionFromCube( OneCube * cube, 
00442                                       int       llx,
00443                                       int       lly,
00444                                       int       halfbox_x, 
00445                                       int       halfbox_y,
00446                                       float     fwhm_factor,
00447                                       float     backvariance,
00448                                       float     sky,
00449                                       float     gain,
00450                                       float     exptime,
00451                                       const char* name,
00452                                       cpl_table** spectrum,
00453                       int       qc_info,
00454                                       int*      check2)
00455 {
00456     int col, row, z ;
00457     OneImage * averagedIm ; 
00458     OneImage * retIm ;
00459     double fit_par[7] ;
00460     double derv_par[7] ;
00461     int mpar[7] ;
00462     double gfit_par[7] ;
00463     double gderv_par[7] ;
00464     int gmpar[7] ;
00465     int fitInd ;
00466     int i ;
00467     double sum ;
00468     double weight[cube->lx][cube->ly] ;
00469     double psf[cube->lx][cube->ly] ;
00470 
00471     double variance ;
00472     double xdat[2] ;
00473     float weighted_sum ;
00474     float counts_tot ;
00475     float counts_bkg ;
00476     float bkg_tot ;
00477 
00478 
00479     int first_col, last_col ;
00480     int first_row, last_row ;
00481     float norm ;
00482      float sum_psf=0;
00483     float sum_wgt=0;
00484       float cenpix = 0;
00485      float cenLambda = 0;
00486      float dispersion = 0;
00487      float lambda=0;
00488      float lambda_start=0;
00489 
00490 
00491     const char* _id="optimalExtractionFromCube";
00492 
00493     if ( NullCube == cube )
00494     {
00495         cpl_msg_error("optimalExtractionFromCube:"," no cube given!\n") ;
00496         return NullImage ;
00497     }
00498     if ( llx < 0 || llx + 2*halfbox_x >= cube->lx || lly < 0 || lly+2*halfbox_y >= cube->ly )
00499     {
00500         cpl_msg_error("optimalExtractionFromCube:"," lower left edge points wrong position!\n") ;
00501         return NullImage ;
00502     }
00503     if ( halfbox_x <= 0 || halfbox_y <= 0 || 2*halfbox_x > cube->lx || 2*halfbox_y > cube->ly )
00504     {
00505         cpl_msg_error("optimalExtractionFromCube:"," wrong halfbox width given!\n") ;
00506         return NullImage ;
00507     }
00508     if ( fwhm_factor <= 0. )
00509     {
00510         cpl_msg_error("optimalExtractionFromCube:"," wrong fwhm_factor given!\n") ;
00511         return NullImage ;
00512     }
00513     if ( backvariance < 0. )
00514     {
00515         cpl_msg_error("optimalExtractionFromCube:"," wrong backvariance given!\n") ;
00516         return NullImage ;
00517     }
00518     if ( exptime <= 0. || exptime == FLAG )
00519     {
00520         cpl_msg_error("optimalExtractionFromCube:"," wrong exposure time given!\n") ;
00521         return NullImage ;
00522     }
00523 
00524     /* allocate memory for spectrum */
00525     if ( NULL == (retIm = new_image(1, cube->np)) )
00526     {
00527         cpl_msg_error("optimalExtractionFromCube:"," memory allocation failed!\n") ;
00528         return NullImage ;
00529     }
00530     /* collapse the cube to be able to compute the weights for optimal extraction */
00531     if ( NullImage == (averagedIm = averageCubeToImage(cube)) )
00532     {
00533         cpl_msg_error("optimalExtractionFromCube:"," averageCubeToImage failed!\n") ;
00534         destroy_image(retIm) ;
00535         return NullImage ;
00536     }
00537 
00538     /* call the 2D-Gaussian fit routine */
00539     for ( i = 0 ; i < 7 ; i++ )
00540     {
00541         mpar[i] = 1 ;
00542     }
00543     /* save_image_to_fits(averagedIm,"out_av_im.fits",-32); */
00544     if ( -1 == (fitInd = fit2dGaussian(averagedIm, fit_par, derv_par, mpar, llx, lly, halfbox_x, halfbox_y, check2 )) )
00545     {
00546         cpl_msg_error("optimalExtractionFromCube:"," fit2dGaussian failed!\n") ;
00547         destroy_image(retIm) ;
00548         destroy_image(averagedIm) ;
00549         return NullImage ;
00550     }
00551 
00552     /* determine the PSF by using the found 2D-Gaussian */
00553     sum = 0 ;
00554     for ( row = 0 ; row < cube->ly ; row++ )
00555     {
00556         for ( col = 0 ; col < cube->lx ; col++ )
00557         {
00558             xdat[0] = (double) col ;
00559             xdat[1] = (double) row ;
00560             psf[col][row] = gaussianEllipse(xdat , fit_par) - fit_par[3] ;
00561             sum += psf[col][row] ;
00562         }
00563     }
00564     /* Scale the PSF and determine the pixel variances and the normalization factor */
00565     norm = 0. ;
00566     variance = 0. ;
00567     sum_psf=0;
00568     for ( row = 0 ; row < cube->ly ; row++ )
00569     {
00570         for ( col = 0 ; col < cube->lx ; col++ )
00571         {
00572             psf[col][row] /= sum ;
00573         sum_psf +=  psf[col][row];
00574             if ( !qfits_isnan(averagedIm->data[col+row*cube->lx]) )
00575             {
00576           /*
00577                 variance = (backvariance + sky + averagedIm -> data[col+row*cube->lx] / exptime) / gain ;
00578           */
00579                 variance = averagedIm -> data[col+row*cube->lx] / gain ;
00580 
00581             }
00582             else 
00583             {
00584                 weight[col][row] = 0. ;
00585             }
00586             if (variance == 0.)
00587             {
00588                 weight[col][row] = 0. ;
00589             }
00590             else
00591             {
00592          
00593                 weight[col][row] = psf[col][row]/variance ;
00594            
00595                 norm += weight[col][row] * weight[col][row] * variance ;
00596 
00597             }
00598         
00599         }
00600     }
00601 
00602     sum_wgt=0;
00603     for ( row = 0 ; row < cube->ly ; row++ )
00604     {
00605         for ( col = 0 ; col < cube->lx ; col++ )
00606         {
00607       weight[col][row] /= norm;
00608           sum_wgt += weight[col][row]*psf[col][row];
00609     }
00610     }
00611     cpl_msg_info(_id,"sum_psf=%f sum_wgt=%f norm=%f",sum_psf,sum_wgt,norm);
00612     destroy_image(averagedIm) ;
00613     if ( norm == 0. )
00614     {
00615         cpl_msg_error("optimalExtractionFromCube:"," normalization sum is zero\n") ;
00616         destroy_image(retIm) ;
00617         return NullImage ;
00618     }
00619 
00620     /* limit the extraction region to the Gaussian, center +- fwhmx/y * cos(theta)  */
00621     first_col = (int) (fit_par[0] - fwhm_factor*fit_par[4]*cos(fit_par[6])) ;
00622     first_col = (first_col > 2 ) ? first_col : 2;
00623     last_col =  (int) (fit_par[0] + fwhm_factor*fit_par[4]*cos(fit_par[6])) ;
00624     last_col = (last_col < 63 ) ? last_col : 63;
00625 
00626     first_row = (int) (fit_par[1] - fwhm_factor*fit_par[5]*cos(fit_par[6])) ;
00627     first_row = (first_row > 2 ) ? first_row : 2;
00628     last_row =  (int) (fit_par[1] + fwhm_factor*fit_par[5]*cos(fit_par[6])) ;
00629     last_row = (last_row < 63 ) ? last_row : 63;
00630  
00631     cpl_msg_info(_id,"fit_par: %f %f %f %f %f %f %f", fit_par[0],fit_par[1],fit_par[2],fit_par[3],
00632                                  fit_par[4],fit_par[5],fit_par[6]);
00633 
00634     if ( first_col < 0 || first_row < 0 || last_col >= cube->lx || last_row >= cube->ly )
00635     {
00636         cpl_msg_error("optimalExtractionFromCube:"," star badly centered in FOV!\n") ;
00637         destroy_image(retIm) ;
00638         return NullImage ;
00639     }
00640 
00641 
00642     cpl_table_new_column(*spectrum,"wavelength", CPL_TYPE_FLOAT);
00643     /* cpl_table_new_column(*spectrum,"intensity" , CPL_TYPE_FLOAT); */
00644     cpl_table_new_column(*spectrum,"counts_tot" , CPL_TYPE_FLOAT);
00645     cpl_table_new_column(*spectrum,"counts_bkg" , CPL_TYPE_FLOAT);
00646     cpl_table_new_column(*spectrum,"bkg_tot" , CPL_TYPE_FLOAT);
00647 
00648     if(qc_info==1) {
00649        cpl_table_new_column(*spectrum,"AMP" , CPL_TYPE_FLOAT);
00650        cpl_table_new_column(*spectrum,"XC" , CPL_TYPE_FLOAT);
00651        cpl_table_new_column(*spectrum,"YC" , CPL_TYPE_FLOAT);
00652        cpl_table_new_column(*spectrum,"BKG" , CPL_TYPE_FLOAT);
00653        cpl_table_new_column(*spectrum,"FWHMX" , CPL_TYPE_FLOAT);
00654        cpl_table_new_column(*spectrum,"FWHMY" , CPL_TYPE_FLOAT);
00655        cpl_table_new_column(*spectrum,"ANGLE" , CPL_TYPE_FLOAT);
00656     }
00657 
00658     cenpix = spiffi_get_cenpix((char*)name);
00659     cenLambda = spiffi_get_cenLambda((char*)name);
00660     dispersion = spiffi_get_dispersion((char*)name);
00661     lambda_start=cenLambda-cenpix*dispersion;
00662 
00663     cpl_msg_info(_id,"frow %d lrow %d fcol %d lcol %d",first_row, last_row, first_col, last_col);
00664     /* go through the planes */
00665     for ( z = 0 ; z < cube->np ; z++ )
00666     {
00667         weighted_sum = 0. ;
00668         counts_tot=0.;
00669         counts_bkg=0.;
00670 
00671         bkg_tot=0.;
00672 
00673         if(qc_info==1) {
00674            fit2dGaussian(cube->plane[z],gfit_par,gderv_par,gmpar,llx,lly,halfbox_x,halfbox_y,check2);
00675     }
00676 
00677         for ( row = first_row ; row <= last_row ; row++ )
00678         {
00679             for ( col = first_col ; col < last_col ; col++ )
00680             {
00681                 if ( !qfits_isnan(cube->plane[z]->data[col+row*cube->lx]) )
00682                 {
00683           
00684                     weighted_sum += weight[col][row] * (cube->plane[z]->data[col+row*cube->lx] - fit_par[3]);
00685             
00686                     counts_bkg += (cube->plane[z]->data[col+row*cube->lx] - fit_par[3]);
00687                     counts_tot += (cube->plane[z]->data[col+row*cube->lx]);
00688                     bkg_tot += fit_par[3];
00689             
00690                 }
00691             }
00692         }    
00693 
00694         if (weighted_sum == 0.)
00695         {
00696             weighted_sum = ZERO ;
00697             counts_tot = ZERO;
00698             counts_bkg = ZERO;
00699             bkg_tot = ZERO;
00700 
00701         }
00702         else
00703         {
00704       /*
00705             weighted_sum /= norm ;
00706       */
00707       
00708       
00709         }
00710 
00711         retIm->data[z] = weighted_sum ;
00712         lambda=lambda_start+z*dispersion;
00713         cpl_table_set_float(*spectrum,"wavelength" ,z,lambda);
00714         /* cpl_table_set_float(*spectrum,"intensity" ,z,weighted_sum); */
00715         cpl_table_set_float(*spectrum,"counts_tot" ,z,counts_tot);
00716         cpl_table_set_float(*spectrum,"counts_bkg" ,z,counts_bkg);
00717         cpl_table_set_float(*spectrum,"bkg_tot" ,z,bkg_tot);
00718 
00719         if(qc_info==1) {
00720            cpl_table_set_float(*spectrum,"AMP" ,z,gfit_par[0]);
00721            cpl_table_set_float(*spectrum,"XC" ,z,gfit_par[1]);
00722            cpl_table_set_float(*spectrum,"YC" ,z,gfit_par[2]);
00723            cpl_table_set_float(*spectrum,"BKG" ,z,gfit_par[3]);
00724            cpl_table_set_float(*spectrum,"FWHMX" ,z,gfit_par[4]);
00725            cpl_table_set_float(*spectrum,"FWHMY" ,z,gfit_par[5]);
00726            cpl_table_set_float(*spectrum,"ANGLE" ,z,gfit_par[6]);
00727     }
00728 
00729     }
00730 
00731     return retIm ;
00732 }
00733 
00734 /*----------------------------------------------------------------------------
00735    Function     :       extractSkyFromCube()
00736    In           :       cube: reduced cube from sky spider observation
00737                         loReject, hiReject: fraction (percentage) of the extreme high and low
00738                                             sky spectrum values that are rejected before 
00739                                             averaging the found sky spectra.
00740                         position: end pixel positions of the straight line in the image
00741                                   dividing the sky from the object pixels.
00742                         tolerance: pixel tolerance which are not considered and subtracted from
00743                                    the diagonal line to be sure to get a clean sky, default: 2  
00744                         posindicator: indicates in which edge of the field of view the sky is projected.  
00745                                       output of spiffi_get_spiderposindex() in fitshead.c
00746    Out          :       resulting averaged sky spectrum
00747    Job          :       extracts a sky spectrum from a reduced sky spider observation, that
00748                         means from a data cube. Therefore, the position of the sky within the
00749                         field of view must be first read from the fits header. A pixel tolerance is
00750                         subtracted.
00751                         The found sky spectra are averaged by rejecting the extreme high and low values.
00752  ---------------------------------------------------------------------------*/
00753 
00754 Vector * extractSkyFromCube( OneCube * cube,
00755                              float     loReject,
00756                              float     hiReject,
00757                              int     * position,
00758                              int       tolerance,
00759                              int       posindicator )
00760 {
00761     Vector * spectrum ;
00762     int x, y, z ;
00763     int n ;
00764     int n_sky ;
00765     int x_low , x_high ;
00766     int y_low , y_high ;
00767     int hi_x, lo_x ;
00768     float * to_average ;
00769     float   cleanMean ;
00770 
00771     if ( NullCube == cube )
00772     {
00773         cpl_msg_error("extractSkyFromCube:"," no cube given!\n") ;
00774         return NullVector ;
00775     }
00776     if ( loReject < 0. || hiReject < 0. || loReject + hiReject >= 90. )
00777     {
00778         cpl_msg_error("extractSkyFromCube:"," wrong or unrealistic loReject and hiReject values!\n") ;
00779         return NullVector ;
00780     }
00781     if ( position == NULL)
00782     {
00783         cpl_msg_error("extractSkyFromCube:"," no position array given!\n") ;
00784         return NullVector ;
00785     }
00786     if ( position[0] < 0 || position[1] < 0 || position[0] > cube->lx  || position[1] > cube->ly )
00787     {
00788         cpl_msg_error("extractSkyFromCube:"," wrong position of sky spider!\n") ;
00789         return NullVector ;
00790     }
00791     if ( tolerance < 0 || tolerance >= cube->lx )
00792     {
00793         cpl_msg_error("extractSkyFromCube:"," wrong tolerance given!\n") ;
00794         return NullVector ;
00795     }
00796     if ( posindicator == 0 )
00797     {
00798         cpl_msg_error("extractSkyFromCube:"," no edge indicator given!\n") ;
00799         return NullVector ;
00800     }
00801 
00802     /* determine the edge of the image where the sky spectra are placed */
00803     switch(posindicator)
00804     {
00805         /* lower right edge */
00806         case 1:
00807             x_low  = position[0] + tolerance ;
00808             x_high = cube->lx ;
00809             y_low  = 0 ;
00810             y_high = position[1] - tolerance ;
00811             break ;
00812         /* upper right edge */
00813         case 2:
00814             x_low  = position[0] + tolerance ;
00815             x_high = cube->lx ;
00816             y_low  = position[1] + tolerance ;
00817             y_high = cube->ly ;
00818             break ;
00819         /* upper left edge */
00820         case 3:
00821             x_low  = 0 ;
00822             x_high = position[0] - tolerance ;
00823             y_low  = position [1] + tolerance ;
00824             y_high = cube->ly ;
00825             break ;
00826         default:
00827             cpl_msg_error("extractSkyFromCube:"," wrong position indicator index!\n") ;
00828             return NullVector ;
00829             break ;
00830     }
00831     if ( x_low >= cube -> lx || x_high < 1 || y_low >= cube->ly || y_high < 1 )
00832     {
00833         cpl_msg_error("extractSkyFromCube:"," tolerance too high!\n") ;
00834         return NullVector ;
00835     }
00836     if ( x_high - x_low != y_high - y_low )
00837     {
00838         cpl_msg_error("extractSkyFromCube:"," sky edge is not a diagonal line!\n") ;
00839         return NullVector ;
00840     }
00841 
00842     /* determine the number of sky pixels in one image plane, take only the full sky pixels
00843        which are not cut by the diagonal line */
00844     n_sky = (x_high - x_low) * (x_high - x_low - 1) / 2 ;
00845     if ( n_sky <= 0 )
00846     {
00847         cpl_msg_error("extractSkyFromCube:"," no sky spectrum in found in cube!\n") ;
00848         return NullVector ;
00849     }
00850     if ( n_sky == 1 )
00851     {
00852         cpl_msg_warning("extractSkyFromCube:"," only one sky spectrum is taken, no averaging!\n") ;
00853     }
00854 
00855     /* allocate memory for the output spectrum */
00856     if ( NullVector == (spectrum = newVector(cube->np)) )
00857     {
00858         cpl_msg_error("extractSkyFromCube:"," could not allocate memory!\n") ;
00859         return NullVector ;
00860     }
00861 
00862     /* go through the image planes */
00863     for ( z = 0 ; z < cube->np ; z++ )
00864     {
00865         /* allocate memory for the sky pixels in one image plane */
00866         if ( NULL == ( to_average = (float*) cpl_calloc(n_sky, sizeof (float)) ) )
00867         {
00868             cpl_msg_error("extractSkyFromCube:"," could not allocate memory!\n") ;
00869             destroyVector(spectrum) ;
00870             return NullVector ;
00871         }
00872         n = 0 ;
00873         switch(posindicator)
00874         {
00875             /* lower right edge */
00876             case 1:
00877                 lo_x = x_low ;
00878                 for ( y = y_low ; y < y_high - 1 ; y++ )
00879                 {
00880                     lo_x++ ;
00881                     for ( x = lo_x ; x < x_high ; x++ )
00882                     {
00883                         to_average[n] = cube->plane[z]->data[x+y*cube->lx] ;
00884                         n++ ;
00885                     }
00886                 }
00887                 break ;
00888             /* lower left edge */
00889             case 2:
00890                 hi_x = x_high ;
00891                 for ( y = y_low ; y < y_high - 1 ; y++ )
00892                 {
00893                     hi_x-- ;
00894                     for ( x = x_low ; x < hi_x ; x++ )
00895                     {
00896                         to_average[n] = cube->plane[z]->data[x+y*cube->lx] ;
00897                         n++ ;
00898                     }
00899                 }
00900                 break ;
00901             /* upper right edge */
00902             case 3:
00903                 lo_x = x_high ;
00904                 for ( y = y_low+1 ; y < y_high ; y++ )
00905                 {
00906                     lo_x-- ;
00907                     for ( x = lo_x ; x < x_high ; x++ )
00908                     {
00909                         to_average[n] = cube->plane[z]->data[x+y*cube->lx] ;
00910                         n++ ;
00911                     }
00912                 }
00913                 break ;
00914             /* upper left edge */
00915             case 4:
00916                 hi_x = x_low ;
00917                 for ( y = y_low+1 ; y < y_high ; y++ )
00918                 {
00919                     hi_x++ ;
00920                     for ( x = x_low ; x < hi_x ; x++ )
00921                     {
00922                         to_average[n] = cube->plane[z]->data[x+y*cube->lx] ;
00923                         n++ ;
00924                     }
00925                 }
00926                 break ;
00927             default:
00928                 cpl_msg_error("extractSkyFromCube:"," wrong position indicator index!\n") ;
00929                 return NullVector ;
00930                 break ;
00931         }
00932         if ( n != n_sky )
00933         {
00934             cpl_msg_warning("extractSkyFromCube:","number of stored sky image pixels does not equal number of computed sky pixels!\n") ;
00935         }
00936 
00937         /* now take a clean mean of the sky "image" */
00938         cleanMean = clean_mean (to_average, n, loReject, hiReject) ;
00939         if (cleanMean == FLT_MAX)
00940         {
00941             cpl_msg_error("extractSkyFromCube:"," could not take a clean mean!\n") ;
00942             destroyVector(spectrum) ;
00943             cpl_free(to_average) ;
00944             return NullVector ;
00945         }
00946         spectrum->data[z] = cleanMean ;
00947         cpl_free (to_average) ;
00948     }
00949      
00950     return spectrum ;
00951 }
00952 
00953 /*---------------------------------------------------------------------------
00954    Function     :       sumRectangleOfCubeSpectra()
00955    In           :       cube: 1 allocated cube, 
00956                         llx, lly, urx, ury: lower left and upper right
00957                                             position of rectangle in x-y plane ,
00958                                             image coordinates 0...
00959    Out          :       result spectrum vector
00960    Job          :       summing routine for a reduced data to get a better spectral S/N
00961                         only for a rectangular aperture.
00962  ---------------------------------------------------------------------------*/
00963 
00964 Vector * sumRectangleOfCubeSpectra( OneCube * cube,
00965                                      int llx,
00966                                      int lly,
00967                                      int urx,
00968                                      int ury )
00969 {
00970     Vector          * sum ;
00971     pixelvalue      *rectangle ;
00972     int             i, j, k, l, m ;
00973     int             recsize ;
00974 
00975     if ( cube == NULL || cube -> np < 1 )
00976     {
00977         cpl_msg_error ("sumRectangleOfCubeSpectra:"," no cube to take the mean of his spectra\n") ;
00978         return NullVector ;
00979     }
00980 
00981     if ((llx<0) || (llx>=cube->lx) ||
00982         (urx<0) || (urx>=cube->lx) ||
00983         (lly<0) || (lly>=cube->ly) ||
00984         (ury<0) || (ury>=cube->ly) ||
00985         (llx>=urx) || (lly>=ury))
00986     {
00987         cpl_msg_error("sumRectangleOfCubeSpectra:"," invalid rectangle coordinates:") ;
00988         cpl_msg_error("sumRectangleOfCubeSpectra:","lower left is [%d %d] upper right is [%d %d]", llx, lly, urx, ury) ;
00989         return NullVector ;
00990     }
00991 
00992     recsize = (urx - llx + 1) * (ury - lly + 1) ;
00993 
00994     /* allocate a new vector to store the average spectral values */
00995     if (NULL == (sum = newVector (cube -> np)) )
00996     {
00997         cpl_msg_error ("sumRectangleOfCubeSpectra:"," cannot allocate a new vector \n") ;
00998         return NullVector ;
00999     }
01000 
01001     /*------------------------------------------------------------------------
01002      *  loop through the cube planes, through the x axis and the y-axis of the
01003      *  plane rectangle and store pixel values in a buffer.
01004      */
01005     for ( i = 0 ; i < cube->np ; i++ )
01006     {
01007         m = 0 ;
01008         rectangle = (pixelvalue *) cpl_calloc (recsize, sizeof (pixelvalue*));
01009 
01010         for ( j = lly ; j <= ury ; j++ )
01011         {
01012             for ( k = llx ; k <= urx ; k++ )
01013             {
01014                 rectangle[m] = cube -> plane[i] -> data[k + j * cube -> lx] ;
01015                 m ++ ;
01016             }
01017         }
01018         for ( l = 0 ; l < recsize ; l++ )
01019         {
01020             if ( qfits_isnan(rectangle[l]) )
01021             {
01022                 continue ;
01023             }
01024             sum -> data[i] += rectangle[l] ;
01025         }
01026         cpl_free ( rectangle ) ;
01027     }
01028     return sum ;
01029 }
01030 
01031 /*---------------------------------------------------------------------------
01032    Function     :       sumCircleOfCubeSpectra()
01033    In           :       cube: 1 allocated cube, 
01034                         centerx, centery: center pixel of circular aperture 
01035                                           in image coordinates
01036                         radius: integer radius of circular aperture
01037    Out          :       result spectrum vector
01038    Job          :       summing routine for a reduced data to get a better spectral S/N
01039                         only for a circular aperture.
01040  ---------------------------------------------------------------------------*/
01041 
01042 Vector * sumCircleOfCubeSpectra( OneCube * cube,
01043                                   int       centerx,
01044                                   int       centery,
01045                                   int       radius )
01046 {
01047     Vector          * sum ;
01048     pixelvalue      * circle ;
01049     int             i, j, k, l, m, n ;
01050     int             circsize ;
01051 
01052     if ( cube == NULL || cube -> np < 1 )
01053     {
01054         cpl_msg_error ("sumCircleOfCubeSpectra:"," no cube to take the mean of his spectra\n") ;
01055         return NullVector ;
01056     }
01057 
01058     if ((centerx+radius>=cube->lx) ||
01059         (centery+radius>=cube->ly) ||
01060         (centerx-radius<0) ||
01061         (centery-radius<0))
01062     {
01063         cpl_msg_error("sumCircleOfCubeSpectra:"," invalid circular coordinates") ;
01064         return NullVector ;
01065     }
01066 
01067     n = 0 ;
01068     for ( j = centery - radius ; j <= centery + radius ; j++ )
01069     {
01070         for ( k = centerx - radius ; k <= centerx + radius ; k++ )
01071         {
01072             if ( (k-centerx)*(k-centerx)+(j-centery)*(j-centery) <= radius*radius )
01073             {
01074                 n ++ ;
01075             }
01076         }
01077     }
01078     if (n == 0)
01079     {
01080         cpl_msg_error ("sumCircleOfCubeSpectra:"," no data points found!\n") ;
01081         return NullVector ;
01082     }
01083     circsize = n ;
01084 
01085     /* allocate a new vector to store the average spectral values */
01086     if (NULL == (sum = newVector (cube -> np)) )
01087     {
01088         cpl_msg_error ("sumCircleOfCubeSpectra:","  cannot allocate a new vector \n") ;
01089         return NullVector ;
01090     }
01091 
01092     /*------------------------------------------------------------------------
01093      *  loop through the cube planes, through the x axis and the y-axis of the
01094      *  plane circle and store pixel values in a buffer.
01095      */
01096     for ( i = 0 ; i < cube->np ; i++ )
01097     {
01098         m = 0 ;
01099         circle = (pixelvalue *) cpl_calloc (circsize, sizeof (pixelvalue*));
01100 
01101         for ( j = centery - radius ; j <= centery + radius ; j++ )
01102         {
01103             for ( k = centerx - radius ; k <= centerx + radius ; k++ )
01104             {
01105                 if ( (k-centerx)*(k-centerx)+(j-centery)*(j-centery) <= radius*radius )
01106                 {
01107                     circle[m] = cube -> plane[i] -> data[k + j * cube -> lx] ;
01108                     m ++ ;
01109                 }
01110             }
01111         }
01112 
01113         for ( l = 0 ; l < circsize ; l++ )
01114         {
01115             if ( qfits_isnan(circle[l]) )
01116             {
01117                 continue ;
01118             }
01119             sum -> data[i] += circle[l] ;
01120         }
01121         cpl_free (circle) ;
01122     }
01123     return sum ;
01124 }
01125 
01126 /*---------------------------------------------------------------------------
01127    Function     :       meanRectangleOfCubeSpectra()
01128    In           :       cube: 1 allocated cube, 
01129                         llx, lly, urx, ury: lower left and upper right
01130                                             position of rectangle in x-y plane ,
01131                                             image coordinates 0...
01132    Out          :       result spectrum vector
01133    Job          :       averaging routine for a reduced data to get a better spectral S/N
01134                         only for a rectangular aperture.
01135  ---------------------------------------------------------------------------*/
01136 
01137 Vector * meanRectangleOfCubeSpectra( OneCube * cube,
01138                                      int llx,
01139                                      int lly,
01140                                      int urx,
01141                                      int ury )
01142 {
01143     Vector          * mean ;
01144     pixelvalue      *rectangle ;
01145     int             i, j, k, l, m ;
01146     int             recsize, nv ;
01147 
01148     if ( cube == NULL || cube -> np < 1 )
01149     {
01150         cpl_msg_error ("meanRectangleOfCubeSpectra:","  no cube to take the mean of his spectra\n") ;
01151         return NullVector ;
01152     }
01153 
01154     if ((llx<0) || (llx>=cube->lx) ||
01155         (urx<0) || (urx>=cube->lx) ||
01156         (lly<0) || (lly>=cube->ly) ||
01157         (ury<0) || (ury>=cube->ly) ||
01158         (llx>=urx) || (lly>=ury))
01159     {
01160         cpl_msg_error("meanRectangleOfCubeSpectra:","  invalid rectangle coordinates:") ;
01161         cpl_msg_error("meanRectangleOfCubeSpectra:","lower left is [%d %d] upper right is [%d %d]",
01162                  llx, lly, urx, ury) ;
01163         return NullVector ;
01164     }
01165 
01166     recsize = (urx - llx + 1) * (ury - lly + 1) ;
01167 
01168     /* allocate a new vector to store the average spectral values */
01169     if (NULL == (mean = newVector (cube -> np)) )
01170     {
01171         cpl_msg_error ("meanRectangleOfCubeSpectra:"," cannot allocate a new vector \n") ;
01172         return NullVector ;
01173     }
01174 
01175     /*------------------------------------------------------------------------
01176      *  loop through the cube planes, through the x axis and the y-axis of the
01177      *  plane rectangle and store pixel values in a buffer.
01178      */
01179     for ( i = 0 ; i < cube->np ; i++ )
01180     {
01181         m = 0 ;
01182         rectangle = (pixelvalue *) cpl_calloc (recsize, sizeof (pixelvalue*));
01183 
01184         for ( j = lly ; j <= ury ; j++ )
01185         {
01186             for ( k = llx ; k <= urx ; k++ )
01187             {
01188                 rectangle[m] = cube -> plane[i] -> data[k + j * cube -> lx] ;
01189                 m ++ ;
01190             }
01191         }
01192         nv = 0 ;
01193         for ( l = 0 ; l < recsize ; l++ )
01194         {
01195             if ( qfits_isnan(rectangle[l]) )
01196             {
01197                 continue ;
01198             }
01199             mean -> data[i] += rectangle[l] ;
01200             nv ++;
01201         }
01202         if ( nv == 0 )
01203         {
01204             mean -> data[i] = ZERO ;
01205         }
01206         else
01207         {
01208             mean -> data[i] /= nv ;
01209         }
01210         cpl_free ( rectangle ) ;
01211     }
01212     return mean ;
01213 }
01214 
01215 /*---------------------------------------------------------------------------
01216    Function     :       meanCircleOfCubeSpectra()
01217    In           :       cube: 1 allocated cube, 
01218                         centerx, centery: center pixel of circular aperture 
01219                                           in image coordinates
01220                         radius: integer radius of circular aperture
01221    Out          :       result spectrum vector
01222    Job          :       averaging routine for a reduced data to get a better spectral S/N
01223                         only for a circular aperture.
01224  ---------------------------------------------------------------------------*/
01225 
01226 Vector * meanCircleOfCubeSpectra( OneCube * cube,
01227                                   int       centerx,
01228                                   int       centery,
01229                                   int       radius )
01230 {
01231     Vector          * mean ;
01232     pixelvalue      * circle ;
01233     int             i, j, k, l, m, n ;
01234     int             circsize, nv ;
01235 
01236     if ( cube == NULL || cube -> np < 1 )
01237     {
01238         cpl_msg_error ("meanCircleOfCubeSpectra:"," no cube to take the mean of his spectra\n") ;
01239         return NullVector ;
01240     }
01241 
01242     if ((centerx+radius>=cube->lx) ||
01243         (centery+radius>=cube->ly) ||
01244         (centerx-radius<0) ||
01245         (centery-radius<0))
01246     {
01247         cpl_msg_error("meanCircleOfCubeSpectra:"," invalid circular coordinates") ;
01248         return NullVector ;
01249     }
01250 
01251     n = 0 ;
01252     for ( j = centery - radius ; j <= centery + radius ; j++ )
01253     {
01254         for ( k = centerx - radius ; k <= centerx + radius ; k++ )
01255         {
01256             if ( (k-centerx)*(k-centerx)+(j-centery)*(j-centery) <= radius*radius )
01257             {
01258                 n ++ ;
01259             }
01260         }
01261     }
01262     if (n == 0)
01263     {
01264         cpl_msg_error ("meanCircleOfCubeSpectra:"," no data points found!\n") ;
01265         return NullVector ;
01266     }
01267     circsize = n ;
01268 
01269     /* allocate a new vector to store the average spectral values */
01270     if (NULL == (mean = newVector (cube -> np)) )
01271     {
01272         cpl_msg_error ("meanCircleOfCubeSpectra:"," cannot allocate a new vector \n") ;
01273         return NullVector ;
01274     }
01275 
01276     /*------------------------------------------------------------------------
01277      *  loop through the cube planes, through the x axis and the y-axis of the
01278      *  plane circle and store pixel values in a buffer.
01279      */
01280     for ( i = 0 ; i < cube->np ; i++ )
01281     {
01282         m = 0 ;
01283         circle = (pixelvalue *) cpl_calloc (circsize, sizeof (pixelvalue*));
01284 
01285         for ( j = centery - radius ; j <= centery + radius ; j++ )
01286         {
01287             for ( k = centerx - radius ; k <= centerx + radius ; k++ )
01288             {
01289                 if ( (k-centerx)*(k-centerx)+(j-centery)*(j-centery) <= radius*radius )
01290                 {
01291                     circle[m] = cube -> plane[i] -> data[k + j * cube -> lx] ;
01292                     m ++ ;
01293                 }
01294             }
01295         }
01296 
01297         nv = 0 ;
01298         for ( l = 0 ; l < circsize ; l++ )
01299         {
01300             if ( qfits_isnan(circle[l]) )
01301             {
01302                 continue ;
01303             }
01304             mean -> data[i] += circle[l] ;
01305             nv ++;
01306         }
01307         if ( nv == 0 )
01308         {
01309             mean -> data[i] = ZERO ;
01310         }
01311         else
01312         {
01313             mean -> data[i] /= nv ;
01314         }
01315 
01316         cpl_free (circle) ;
01317     }
01318     return mean ;
01319 }
01320 
01321 /*---------------------------------------------------------------------------
01322    Function     :       blackbodySpectrum()
01323    In           :       templateSpec: spectrum of a standard star
01324                         (1-d image with fits header)
01325                         temp: blackbody temperature in Kelvin (standard Star temp), 
01326    Out          :       resulting spectrum vector
01327    Job          :       computes a blackbody spectral intensity distribution
01328                         (W/(m^2 lambda ster)) 
01329  ---------------------------------------------------------------------------*/
01330 
01331 Vector * blackbodySpectrum( char * templateSpec, double temp )
01332 {
01333     Vector * retSpec ;
01334     fits_header * hdr ;
01335     int n ;
01336     double cenpix ;
01337     int npix ;
01338     double cenLambda ;
01339     double firstLambda ;
01340     double disp ;
01341     double lambda ;
01342     double intens ;
01343     double denom ;
01344     double norm ;
01345    
01346     if ( NULL == templateSpec )
01347     {
01348         cpl_msg_error ("blackbodySpectrum:"," now input image given!\n") ;
01349         return NULL ;
01350     }
01351     if ( temp < 0. )
01352     {
01353         cpl_msg_error ("blackbodySpectrum:"," wrong temperature given!\n") ;
01354         return NULL ;
01355     }
01356 
01357     /* get the fits header information needed */
01358     if ( NULL == (hdr = fits_read_header( templateSpec )) )
01359     {
01360         cpl_msg_error ("blackbodySpectrum:"," cannot read fits header and allocate data structure\n") ;
01361         return NULL ;
01362     }
01363 
01364     if ( -1. == (cenpix = fits_header_getdouble(hdr, "CRPIX2", -1.)) )
01365     {
01366         cpl_msg_error ("blackbodySpectrum:"," cannot get CRPIX2\n") ;
01367         fits_header_destroy(hdr) ;
01368         return NULL ;
01369     }
01370     if ( -1. == (cenLambda = fits_header_getdouble(hdr, "CRVAL2", -1.)) )
01371     {
01372         cpl_msg_error ("blackbodySpectrum:"," cannot get CRVAL2\n") ;
01373         fits_header_destroy(hdr) ;
01374         return NULL ;
01375     }
01376     if ( -1. == (disp = fits_header_getdouble(hdr, "CDELT2", -1.)) )
01377     {
01378         cpl_msg_error ("blackbodySpectrum:"," cannot get CDELT2\n") ;
01379         fits_header_destroy(hdr) ;
01380         return NULL ;
01381     }
01382     if ( -1. == (npix = fits_header_getint(hdr, "NAXIS2", -1.)) )
01383     {
01384         cpl_msg_error ("blackbodySpectrum:"," cannot get NAXIS2\n") ;
01385         fits_header_destroy(hdr) ;
01386         return NULL ;
01387     }
01388     fits_header_destroy(hdr) ;
01389 
01390     if (NULL == (retSpec = newVector (npix)))
01391     {
01392         cpl_msg_error ("blackbodySpectrum:"," could not allocate memory!\n") ;
01393         return NULL ;
01394     }
01395     
01396     /* shift from fits to image coordinates */
01397     cenpix-- ;
01398   
01399     firstLambda = cenLambda - cenpix * disp ;
01400     for ( n = 0 ; n < npix ; n++ )
01401     {
01402         lambda = firstLambda + disp * (double)n ;
01403         /* convert from microns to m */
01404 
01405         lambda /= 1.0e6 ;
01406         denom = 1./(exp(PLANCK*SPEED_OF_LIGHT/(lambda*BOLTZMANN*temp)) - 1.) ;
01407         intens = 2.*PI_NUMB*PLANCK*SPEED_OF_LIGHT*SPEED_OF_LIGHT / pow(lambda, 5) * denom ;
01408         retSpec->data[n] = intens ;   
01409     }
01410     norm = retSpec->data[npix/2] ;
01411     for ( n = 0 ; n < npix ; n++ )
01412     {
01413         retSpec->data[n] /= norm ;
01414     }
01415     
01416     return retSpec ;
01417 }
01418 
01419 
01420 /*---------------------------------------------------------------------------
01421    Function     :       medianRectangleOfCubeSpectra()
01422    In           :       cube: 1 allocated cube, 
01423                         llx, lly, urx, ury: lower left and upper right
01424                                             position of rectangle in x-y plane ,
01425                                             image coordinates 0...
01426    Out          :       result spectrum vector
01427    Job          :       median routine for a reduced data to get a better spectral S/N
01428                         only for a rectangular aperture.
01429  ---------------------------------------------------------------------------*/
01430 
01431 Vector * medianRectangleOfCubeSpectra( OneCube * cube,
01432                                              int llx,
01433                                              int lly,
01434                                              int urx,
01435                                              int ury )
01436 {
01437     Vector          * med ;
01438     pixelvalue      *rectangle ;
01439     int             i, j, k, m ;
01440     int             recsize ;
01441 
01442     if ( cube == NULL || cube -> np < 1 )
01443     {
01444         cpl_msg_error ("medianRectangleOfCubeSpectra:"," no cube to take the mean of his spectra\n") ;
01445         return NullVector ;
01446     }
01447 
01448     if ((llx<0) || (llx>=cube->lx) ||
01449         (urx<0) || (urx>=cube->lx) ||
01450         (lly<0) || (lly>=cube->ly) ||
01451         (ury<0) || (ury>=cube->ly) ||
01452         (llx>=urx) || (lly>=ury))
01453     {
01454         cpl_msg_error("medianRectangleOfCubeSpectra:"," invalid rectangle coordinates:") ;
01455         cpl_msg_error("medianRectangleOfCubeSpectra:","lower left is [%d %d] upper right is [%d %d]", llx, lly, urx, ury) ;
01456         return NullVector ;
01457     }
01458 
01459     recsize = (urx - llx + 1) * (ury - lly + 1) ;
01460 
01461     /* allocate a new vector to store the average spectral values */
01462     if (NULL == (med = newVector (cube -> np)) )
01463     {
01464         cpl_msg_error ("medianRectangleOfCubeSpectra:"," cannot allocate a new vector \n") ;
01465         return NullVector ;
01466     }
01467 
01468     /*------------------------------------------------------------------------
01469      *  loop through the cube planes, through the x axis and the y-axis of the
01470      *  plane rectangle and store pixel values in a buffer.
01471      */
01472     for ( i = 0 ; i < cube->np ; i++ )
01473     {
01474         m = 0 ;
01475         rectangle = (pixelvalue *) cpl_calloc (recsize, sizeof (pixelvalue*));
01476 
01477         for ( j = lly ; j <= ury ; j++ )
01478         {
01479             for ( k = llx ; k <= urx ; k++ )
01480             {
01481                 if ( qfits_isnan(cube->plane[i]->data[k+j*cube->lx]) )
01482                 {
01483                     continue ;
01484                 }
01485                 else
01486                 {
01487                     rectangle[m] = cube -> plane[i] -> data[k + j * cube -> lx] ;
01488                     m ++ ;
01489                 }
01490             }
01491         }
01492         if ( m == 0 )
01493         {
01494             med->data[i] = 0. ;
01495         }
01496         else
01497         {
01498             med->data[i] = median(rectangle, m) ;
01499         }
01500         cpl_free ( rectangle ) ;
01501     }
01502     return med ;
01503 }
01504 
01505 /*---------------------------------------------------------------------------
01506    Function     :       medianCircleOfCubeSpectra()
01507    In           :       cube: 1 allocated cube, 
01508                         centerx, centery: center pixel of circular aperture 
01509                                           in image coordinates
01510                         radius: integer radius of circular aperture
01511    Out          :       result spectrum vector
01512    Job          :       median routine for a reduced data to get a better spectral S/N
01513                         only for a circular aperture.
01514  ---------------------------------------------------------------------------*/
01515 
01516 Vector * medianCircleOfCubeSpectra( OneCube * cube,
01517                                   int       centerx,
01518                                   int       centery,
01519                                   int       radius )
01520 {
01521     Vector          * med ;
01522     pixelvalue      * circle ;
01523     int             i, j, k, l, m, n ;
01524     int             circsize, nv ;
01525 
01526     if ( cube == NULL || cube -> np < 1 )
01527     {
01528         cpl_msg_error ("medianCircleOfCubeSpectra:"," no cube to take the mean of his spectra\n") ;
01529         return NullVector ;
01530     }
01531 
01532     if ((centerx+radius>=cube->lx) ||
01533         (centery+radius>=cube->ly) ||
01534         (centerx-radius<0) ||
01535         (centery-radius<0))
01536     {
01537         cpl_msg_error("medianCircleOfCubeSpectra:"," invalid circular coordinates") ;
01538         return NullVector ;
01539     }
01540 
01541     n = 0 ;
01542     for ( j = centery - radius ; j <= centery + radius ; j++ )
01543     {
01544         for ( k = centerx - radius ; k <= centerx + radius ; k++ )
01545         {
01546             if ( (k-centerx)*(k-centerx)+(j-centery)*(j-centery) <= radius*radius )
01547             {
01548                 n ++ ;
01549             }
01550         }
01551     }
01552     if (n == 0)
01553     {
01554         cpl_msg_error ("medianCircleOfCubeSpectra:"," no data points found!\n") ;
01555         return NullVector ;
01556     }
01557     circsize = n ;
01558 
01559     /* allocate a new vector to store the average spectral values */
01560     if (NULL == (med = newVector (cube -> np)) )
01561     {
01562         cpl_msg_error ("medianCircleOfCubeSpectra:"," cannot allocate a new vector \n") ;
01563         return NullVector ;
01564     }
01565 
01566     /*------------------------------------------------------------------------
01567      *  loop through the cube planes, through the x axis and the y-axis of the
01568      *  plane circle and store pixel values in a buffer.
01569      */
01570     for ( i = 0 ; i < cube->np ; i++ )
01571     {
01572         m = 0 ;
01573         circle = (pixelvalue *) cpl_calloc (circsize, sizeof (pixelvalue*));
01574 
01575         for ( j = centery - radius ; j <= centery + radius ; j++ )
01576         {
01577             for ( k = centerx - radius ; k <= centerx + radius ; k++ )
01578             {
01579                 if ( (k-centerx)*(k-centerx)+(j-centery)*(j-centery) <= radius*radius )
01580                 {
01581                     circle[m] = cube -> plane[i] -> data[k + j * cube -> lx] ;
01582                     m ++ ;
01583                 }
01584             }
01585         }
01586 
01587         nv = 0 ;
01588         for ( l = 0 ; l < circsize ; l++ )
01589         {
01590             if ( qfits_isnan(circle[l]) )
01591             {
01592                 continue ;
01593             }
01594             med -> data[i] += circle[l] ;
01595             nv ++;
01596         }
01597         if ( nv == 0 )
01598         {
01599             med->data[i] = 0. ;
01600         }
01601         else
01602         {
01603             med->data[i] = median(circle, nv) ; 
01604         }
01605         cpl_free (circle) ;
01606     }
01607     return med ;
01608 }
01609 
01610 /*---------------------------------------------------------------------------
01611    Function     :       cleanmeanRectangleOfCubeSpectra()
01612    In           :       cube: 1 allocated cube, 
01613                         llx, lly, urx, ury: lower left and upper right
01614                                             position of rectangle in x-y plane ,
01615                                             image coordinates 0...
01616    Out          :       result spectrum vector
01617    Job          :       clean averaging routine for a reduced data to get a better spectral S/N
01618                         only for a rectangular aperture.
01619  ---------------------------------------------------------------------------*/
01620 
01621 Vector * cleanmeanRectangleOfCubeSpectra( OneCube * cube,
01622                                           int llx,
01623                                           int lly,
01624                                           int urx,
01625                                           int ury,
01626                                           float lo_reject,
01627                                           float hi_reject )
01628 {
01629     Vector          * clean ;
01630     pixelvalue      *rectangle ;
01631     int             i, j, k, m ;
01632     int             recsize ;
01633 
01634     if ( cube == NULL || cube -> np < 1 )
01635     {
01636         cpl_msg_error ("cleanmeanRectangleOfCubeSpectra:"," no cube to take the mean of his spectra\n") ;
01637         return NullVector ;
01638     }
01639 
01640     if ((llx<0) || (llx>=cube->lx) ||
01641         (urx<0) || (urx>=cube->lx) ||
01642         (lly<0) || (lly>=cube->ly) ||
01643         (ury<0) || (ury>=cube->ly) ||
01644         (llx>=urx) || (lly>=ury))
01645     {
01646         cpl_msg_error("cleanmeanRectangleOfCubeSpectra:"," invalid rectangle coordinates:") ;
01647         cpl_msg_error("cleanmeanRectangleOfCubeSpectra:",
01648                        "lower left is [%d %d] upper right is [%d %d]",
01649                  llx, lly, urx, ury) ;
01650         return NullVector ;
01651     }
01652 
01653     recsize = (urx - llx + 1) * (ury - lly + 1) ;
01654 
01655     /* allocate a new vector to store the average spectral values */
01656     if (NULL == (clean = newVector (cube -> np)) )
01657     {
01658         cpl_msg_error ("cleanmeanRectangleOfCubeSpectra:"," cannot allocate a new vector") ;
01659         return NullVector ;
01660     }
01661 
01662     /*------------------------------------------------------------------------
01663      *  loop through the cube planes, through the x axis and the y-axis of the
01664      *  plane rectangle and store pixel values in a buffer.
01665      */
01666     for ( i = 0 ; i < cube->np ; i++ )
01667     {
01668         m = 0 ;
01669         rectangle = (pixelvalue *) cpl_calloc (recsize, sizeof (pixelvalue*));
01670 
01671         for ( j = lly ; j <= ury ; j++ )
01672         {
01673             for ( k = llx ; k <= urx ; k++ )
01674             {
01675                 if ( qfits_isnan(cube->plane[i]->data[k+j*cube->lx]) )
01676                 {
01677                     continue ;
01678                 }
01679                 else
01680                 {
01681                     rectangle[m] = cube -> plane[i] -> data[k + j * cube -> lx] ;
01682                     m ++ ;
01683                 }
01684             }
01685         }
01686         if ( m == 0 )
01687         {
01688             clean->data[i] = 0. ;
01689         }
01690         else
01691         {
01692             clean->data[i] = clean_mean(rectangle, m, lo_reject, hi_reject) ;
01693         }
01694         cpl_free ( rectangle ) ;
01695     }
01696     return clean ;
01697 }
01698 
01699 /*---------------------------------------------------------------------------
01700    Function     :       cleanmeanCircleOfCubeSpectra()
01701    In           :       cube: 1 allocated cube, 
01702                         centerx, centery: center pixel of circular aperture 
01703                                           in image coordinates
01704                         radius: integer radius of circular aperture
01705    Out          :       result spectrum vector
01706    Job          :       clean averaging routine for a reduced data to get a better spectral S/N
01707                         only for a circular aperture.
01708  ---------------------------------------------------------------------------*/
01709 
01710 Vector * cleanmeanCircleOfCubeSpectra( OneCube * cube,
01711                                   int       centerx,
01712                                   int       centery,
01713                                   int       radius,
01714                                   float     lo_reject,
01715                                   float     hi_reject )
01716 {
01717     Vector          * clean ;
01718     pixelvalue      * circle ;
01719     int             i, j, k, l, m, n ;
01720     int             circsize, nv ;
01721 
01722     if ( cube == NULL || cube -> np < 1 )
01723     {
01724         cpl_msg_error ("cleanmeanCircleOfCubeSpectra:"," no cube to take the mean of his spectra\n") ;
01725         return NullVector ;
01726     }
01727 
01728     if ((centerx+radius>=cube->lx) ||
01729         (centery+radius>=cube->ly) ||
01730         (centerx-radius<0) ||
01731         (centery-radius<0))
01732     {
01733         cpl_msg_error("cleanmeanCircleOfCubeSpectra:"," invalid circular coordinates") ;
01734         return NullVector ;
01735     }
01736 
01737     n = 0 ;
01738     for ( j = centery - radius ; j <= centery + radius ; j++ )
01739     {
01740         for ( k = centerx - radius ; k <= centerx + radius ; k++ )
01741         {
01742             if ( (k-centerx)*(k-centerx)+(j-centery)*(j-centery) <= radius*radius )
01743             {
01744                 n ++ ;
01745             }
01746         }
01747     }
01748     if (n == 0)
01749     {
01750         cpl_msg_error ("cleanmeanCircleOfCubeSpectra:"," no data points found!\n") ;
01751         return NullVector ;
01752     }
01753     circsize = n ;
01754 
01755     /* allocate a new vector to store the average spectral values */
01756     if (NULL == (clean = newVector (cube -> np)) )
01757     {
01758         cpl_msg_error ("cleanmeanCircleOfCubeSpectra:"," cannot allocate a new vector \n") ;
01759         return NullVector ;
01760     }
01761 
01762     /*------------------------------------------------------------------------
01763      *  loop through the cube planes, through the x axis and the y-axis of the
01764      *  plane circle and store pixel values in a buffer.
01765      */
01766     for ( i = 0 ; i < cube->np ; i++ )
01767     {
01768         m = 0 ;
01769         circle = (pixelvalue *) cpl_calloc (circsize, sizeof (pixelvalue*));
01770 
01771         for ( j = centery - radius ; j <= centery + radius ; j++ )
01772         {
01773             for ( k = centerx - radius ; k <= centerx + radius ; k++ )
01774             {
01775                 if ( (k-centerx)*(k-centerx)+(j-centery)*(j-centery) <= radius*radius )
01776                 {
01777                     circle[m] = cube -> plane[i] -> data[k + j * cube -> lx] ;
01778                     m ++ ;
01779                 }
01780             }
01781         }
01782 
01783         nv = 0 ;
01784         for ( l = 0 ; l < circsize ; l++ )
01785         {
01786             if ( qfits_isnan(circle[l]) )
01787             {
01788                 continue ;
01789             }
01790             clean -> data[i] += circle[l] ;
01791             nv ++;
01792         }
01793         if ( nv == 0 )
01794         {
01795             clean->data[i] = 0. ;
01796         }
01797         else
01798         {
01799             clean->data[i] = clean_mean(circle, nv, lo_reject, hi_reject) ; 
01800         }
01801         cpl_free (circle) ;
01802     }
01803     return clean ;
01804 }
01805 
01806 /*---------------------------------------------------------------------------
01807    Function     :       shiftArray()
01808    In           :       input: input array, 
01809                         n_elements: number of elements in input array
01810                         shift: sub-pixel shift value (must be < 1.)
01811                         ker: interpolation kernel
01812    Out          :       resulting float array 
01813    Job          :       shifts an array by a sub-pixel shift value using a tanh
01814                         interpolation kernel
01815  ---------------------------------------------------------------------------*/
01816 
01817 float * shiftArray ( float * input, int n_elements, float shift, double * ker ) 
01818 {
01819     float  *         shifted ;
01820     int              samples = KERNEL_SAMPLES ;
01821     int           i ;
01822     float            fx ;
01823     float            rx ;
01824     int              px ;
01825     int              tabx ;
01826     float            value ;
01827     /*size_t           pos ;*/
01828     register float * pix ;
01829     int              mid;
01830     float            norm ;
01831 
01832     /* error handling: test entries */
01833     if (input==NULL) 
01834     {
01835         cpl_msg_error("shiftArray:"," no input array given!\n") ;
01836         return NULL ;
01837     }
01838     if (n_elements<=0) 
01839     {
01840         cpl_msg_error("shiftArray:"," wrong number of elements in input array given!\n") ;
01841         return NULL ;
01842     }
01843 
01844     shifted    = (float*) cpl_calloc(n_elements, sizeof(float)) ;
01845 
01846     /* Shifting by a zero offset returns a copy of the input image */
01847     if ((fabs(shift)<1e-2))
01848     {
01849         for (i = 0 ; i <  n_elements ; i++ )
01850         {
01851             shifted[i] = input[i] ;
01852         }
01853         return shifted ;
01854     }
01855     
01856     mid = (int)samples/(int)2 ;
01857 
01858     for (i=1 ; i<  n_elements-2 ; i++) 
01859     {
01860         fx = (float)i+shift ;
01861         px = nint(fx) ;
01862         rx = fx - (float)px ;
01863         pix = input ;
01864 
01865         if ((px>=1) && (px<(n_elements-2))) 
01866         {
01867             tabx = (int)(fabs((float)mid * rx)) ;
01868             /* exclude blank (ZERO) pixels from interpolation */
01869             if (qfits_isnan(pix[i]))
01870             {
01871                 value = ZERO ;
01872             }
01873             else
01874             {
01875                 if (qfits_isnan(pix[i-1]))
01876                 {
01877                     pix[i-1] = 0. ;
01878                 }
01879                 if (qfits_isnan(pix[i+1]))
01880                 {
01881                     pix[i+1] = 0. ;
01882                 }
01883                 if (qfits_isnan(pix[i+2]))
01884                 {
01885                     pix[i+2] = 0. ;
01886                 }
01887 
01888                 /*
01889                  * Sum up over 4 closest pixel values,
01890                  * weighted by interpolation kernel values
01891                  */
01892                 value =     pix[i-1] * ker[mid+tabx] +
01893                             pix[i] *   ker[tabx] +
01894                             pix[i+1] * ker[mid-tabx] +
01895                             pix[i+2] * ker[samples-tabx-1] ;
01896                 /*
01897                  * Also sum up interpolation kernel coefficients
01898                  * for further normalization
01899                  */
01900                 norm =      ker[mid+tabx] +
01901                             ker[tabx] +
01902                             ker[mid-tabx] +
01903                             ker[samples-tabx-1] ;
01904                 if (fabs(norm) > 1e-4) 
01905                 {
01906                     value /= norm ;
01907                 }
01908             }
01909         }   
01910         else 
01911         {
01912             value = 0.0 ;
01913         }
01914         if ( qfits_isnan(value) )
01915         {
01916             shifted[i] = ZERO ;
01917         }
01918         else
01919         {
01920             shifted[i] = value ;
01921         }
01922     }  
01923     return shifted ;
01924 }
01925 
01926 
01927 /*--------------------------------------------------------------------------*/
01928 
01929 /*----------------------------------------------------------------------------
01930    Function     :       divImageBySpectrum()
01931    In           :       image: resampled image
01932                         spectrum: resampled spectrum in image format
01933    Out          :       resulting image
01934    Job          :       divides a resampled image with a resampled spectrum
01935                         in the same spectral range
01936                         that means all image columns are multiplied with the same spectrum
01937  ---------------------------------------------------------------------------*/
01938 OneImage * divImageBySpectrum( OneImage * image, OneImage * spectrum )
01939 {
01940     int col, row ;
01941     OneImage * retImage ;
01942     if ( image == NullImage )
01943     {
01944         cpl_msg_error("divImageBySpectrum:","no image given!") ;
01945         return NullImage ;
01946     }
01947     if ( spectrum == NullImage )
01948     {
01949         cpl_msg_error("divImageBySpectrum:","no spectrum image given!") ;
01950         return NullImage ;
01951     }
01952     if ( spectrum -> ly != image -> ly )
01953     {
01954         cpl_msg_error("divImageBySpectrum:","images are not compatible in pixel length!") ;
01955         return NullImage ;
01956     }
01957     if ( NullImage == (retImage = copy_image(image)) )
01958     {
01959         cpl_msg_error("divImageBySpectrum:","could not copy original image!") ;
01960         return NullImage ;
01961     }
01962     for ( col = 0 ; col < image -> lx ; col++ )
01963     {
01964         for ( row = 0 ; row < image -> ly ; row++ )
01965         {
01966             if ( !qfits_isnan(image->data[col+row*image->lx]) &&
01967                  !qfits_isnan(spectrum->data[col+row*image->lx]))
01968             {
01969                 retImage->data[col+row*image->lx] =
01970                     image->data[col+row*image->lx] / spectrum->data[row] ;
01971             }
01972          }
01973     }
01974     return retImage ;
01975 }
01976 
01977 /*--------------------------------------------------------------------------*/

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