wavecal.c

00001 /*******************************************************************************
00002 * E.S.O. - VLT project
00003 *
00004 * 
00005 *
00006 * who       when      what
00007 * --------  --------  ----------------------------------------------
00008 * schreib  22/01/02  created
00009 */
00010 
00011 /************************************************************************
00012 *   NAME
00013 *     wavecal.c -
00014 *     routines needed for wavelength calibration with smoothing only
00015 *     within the slitlets
00016 *
00017 *   SYNOPSIS
00018 *  1) Bcoeffs * newBcoeffs( int n_slitlets,
00019 *                           int n_acoeffs,
00020 *                           int n_bcoeffs )
00021 *  2) void destroyBcoeffs ( Bcoeffs * bco )
00022 *  3) int   coeffsCrossSlitFit ( int      n_columns,
00023 *                                float ** acoefs,
00024 *                                float ** dacoefs,
00025 *                                Bcoeffs* bco,
00026 *                                float    sigma_factor,
00027 *                                float    dispersion,
00028 *                                float    pixel_dist,
00029 *                                float  * chisq )
00030 *  4) OneImage * waveMapSlit ( float ** acoefs,
00031 *                              int      n_acoefs,
00032 *                              int      n_rows,
00033 *                              int      n_columns )
00034 *  5) OneImage * waveCal( OneImage   * image, 
00035 *                         FitParams ** par ,
00036 *                         float     ** abuf,
00037 *                         int          n_slitlets,
00038 *                         int       ** row_clean,
00039 *                         float     ** wavelength_clean,
00040 *                         int        * n_found_lines,
00041 *                         float        dispersion,
00042 *                         int          halfWidth,
00043 *                         float        minAmplitude,
00044 *                         float        max_residual,
00045 *                         float        fwhm,
00046 *                         int          n_a_fitcoefs,
00047 *                         int          n_b_fitcoefs,
00048 *                         float        sigmaFactor,
00049 *                         float        pixel_dist )
00050 *  6) int checkForFakeLines ( FitParams ** par,
00051 *                             float        dispersion,
00052 *                             float     ** wavelength_clean,
00053 *                             int       ** row_clean,
00054 *                             int        * n_found_lines,
00055 *                             int          n_columns,
00056 *                             float        pixel_tolerance )
00057 *  7) OneImage * createShiftedSlitWavemap ( OneImage * lineIm,
00058 *                                           float    ** coeffs,
00059 *                                           int      n_fitcoeffs,
00060 *                                           float  * wavelength,
00061 *                                           float  * intensity,
00062 *                                           int      n_lines,
00063 *                                           int      magFactor )
00064 *   
00065 *   DESCRIPTION
00066 *  1) allocates memory for a new array of 
00067 *     Bcoeffs data structures
00068 *  2) frees memory of an array of Bcoeffs data structures
00069 *  3) Fits each single polynomial coefficient acoefs resulting from polyfit 
00070 *     across the columns of each slitlet and use the result of this fit to
00071 *     smooth the acoefs.
00072 *  4) builds a new wavelength calibration map as fits image
00073 *     by using the fit coeficients.
00074 *  5) this routine takes an image from a calibration
00075 *     emission lamp and delivers the smoothed fit coefficients of  
00076 *     a polynomial fit along the columns of the line positions as output. 
00077 *     This routine expects Nyquist sampled spectra 
00078 *     (either an interleaved image or an image convolved with an 
00079 *      appropriate function in spectral direction)
00080 *  6) this routine searches for successfully fitted fake lines like
00081 *     bad pixels by comparing the found line positons with 
00082 *     estimated template positions. This routine should be
00083 *     inserted in the wavelength calibration routine just after
00084 *     the fitLines() routine.
00085 *  7) This routine cross-correlates a shifted emission line frames 
00086 *     and determines the shift to the old one which is given by
00087 *     its polynomial coefficients.
00088 *     Then the a0 coefficients is recalculated and afterwards
00089 *     a new wavelength calibration map is generated using the 
00090 *     already calculated smoothed polynomial coefficients.
00091 *
00092 *   FILES
00093 *
00094 *   ENVIRONMENT
00095 *
00096 *   RETURN VALUES 
00097 *
00098 *   CAUTIONS 
00099 *
00100 *   EXAMPLES
00101 *
00102 *   SEE ALSO
00103 *
00104 *   BUGS   
00105 *
00106 *------------------------------------------------------------------------
00107 */
00108 
00109 #include "vltPort.h"
00110 #include <math.h>
00111 
00112 /* 
00113  * System Headers
00114  */
00115 /* 
00116  * Local Headers
00117  */
00118 
00119 #include "wavecal.h"
00120 #include "wave_calibration.h"
00121 #include "solve_poly_root.h"
00122 
00123 /*
00124  * Private functions prototype
00125  */
00126 Bcoeffs * newBcoeffs( int n_slitlets,
00127                       int n_acoeffs,
00128                       int n_bcoeffs ) ;
00129 
00130 void destroyBcoeffs ( Bcoeffs * bco ) ;
00131 
00132 int   coeffsCrossSlitFit ( int      n_columns,
00133                            float ** acoefs,
00134                            float ** dacoefs,
00135                            Bcoeffs* bco,
00136                            float    sigma_factor,
00137                            float    dispersion,
00138                            float    pixel_dist,
00139                            float  * chisq ) ;
00140 
00141 
00142 /*
00143  * function definitions
00144  */
00145 
00146 /*----------------------------------------------------------------------------
00147    Function     :       newBcoeffs()
00148    In           :       n_slitlets: number of slitlets
00149                         n_acoeffs: number of acoeffs
00150                         n_bcoeffs: number of bcoeffs
00151    Out          :       allocated array of Bcoeffs data structures
00152    Job          :       allocates memory for a new array of 
00153                         Bcoeffs data structures
00154  ---------------------------------------------------------------------------*/
00155 
00156 Bcoeffs * newBcoeffs( int n_slitlets,
00157                       int n_acoeffs,
00158                       int n_bcoeffs )
00159 {
00160     int i, n ;
00161     Bcoeffs * returnbco ;
00162 
00163     if ( NULL == (returnbco = (Bcoeffs*) cpl_calloc(n_slitlets, sizeof(Bcoeffs))) )
00164     {
00165         cpl_msg_error ("newBcoeffs:","could not allocate memory") ;
00166         return NULL ;
00167     }
00168     returnbco -> n_slitlets = n_slitlets ;
00169     returnbco -> n_acoeffs  = n_acoeffs ;
00170     returnbco -> n_bcoeffs  = n_bcoeffs ;
00171     for ( i = 0 ; i < n_slitlets ; i++ )
00172     {
00173         returnbco[i].slitlet = i ;
00174         if ( NULL == (returnbco[i].b = (float**)cpl_calloc(n_acoeffs, sizeof(float*)) ) ) 
00175         {
00176             cpl_msg_error ("newBcoeffs:","could not allocate memory") ;
00177             return NULL ;
00178         }
00179         for ( n = 0 ; n < n_acoeffs ; n++ )
00180         {
00181             if ( NULL == (returnbco[i].b[n] = (float*)cpl_calloc(n_bcoeffs, sizeof(float))) )
00182             {
00183                 cpl_msg_error ("newBcoeffs:","could not allocate memory") ;
00184                 return NULL ;
00185             }
00186         }  
00187     }
00188     return returnbco ;
00189 }
00190 
00191 /*----------------------------------------------------------------------------
00192    Function     :       destroyBcoeffs()
00193    In           :       bcoeffs to destroy
00194    Out          :       nothing
00195    Job          :       frees memory of an array of Bcoeffs data structures
00196  ---------------------------------------------------------------------------*/
00197 
00198 void destroyBcoeffs ( Bcoeffs * bco )
00199 {
00200     int i, n ;
00201   
00202     for ( i = 0 ; i < bco->n_slitlets ; i++ )
00203     {
00204         for ( n = 0 ; n < bco->n_acoeffs ; n++ )
00205         {
00206             cpl_free (bco[i].b[n]) ;
00207         }
00208         cpl_free(bco[i].b) ;
00209     }
00210     
00211     cpl_free (bco) ;  
00212 }
00213 
00214 /*----------------------------------------------------------------------------
00215    Function     :       coeffsCrossSlitFit()
00216    In           :       n_columns:    number of image columns
00217                         acoefs:       coeffs fitted in polyfit
00218                                       note: this is a matrix of the polynomial
00219                                       coefficients, each index for each column:
00220                                       acoefs[index][column]
00221                         dacoefs:      fit errors of the corresponding acoefs
00222                         bco:          the fitted Bcoeffs data structure for each slitlet
00223                         sigma_factor: factor of sigma beyond which the
00224                                       column coefficients are discarded for the fit
00225                         dispersion:   estimate of the dispersion
00226                         pixel_dist:   estimate of minimal pixel distance in spectral direction between slitle
00227 ts
00228    Out          :       0 if all went o.k
00229                         -1 if an error occurred
00230                         the found fit coefficients bco data structure to calculate
00231                         the smoothed acoefs
00232                         chisq:        list of the chi squared of each fit for each slitlet
00233    Job          :       Fits each single polynomial coefficient acoefs resulting from polyfit
00234                         across the columns of each slitlet and use the result of this fit to
00235                         smooth the acoefs.
00236  ---------------------------------------------------------------------------*/
00237 
00238 int   coeffsCrossSlitFit ( int      n_columns,
00239                            float ** acoefs,
00240                            float ** dacoefs,
00241                            Bcoeffs* bco,
00242                            float    sigma_factor,
00243                            float    dispersion,
00244                            float    pixel_dist,
00245                            float  * chisq )
00246 {
00247     float col_index, sub_col_index[n_columns] ;
00248     float sub_acoefs[n_columns], sub_dacoefs[n_columns] ;
00249     float wcoefs[bco->n_bcoeffs] ;
00250     float ** ucoefs, **vcoefs, **covar ;
00251     float * acoefsclean ;
00252     double sum, sumq, mean ;
00253     double sigma ;
00254     double cliphi, cliplo ;
00255     float offset ;
00256     float threshold ;
00257     int edge[bco->n_slitlets] ;
00258     int ed1, ed2 ;
00259     int i, n, num, ndata ;
00260     int nc, ns ;
00261     int index ;
00262     int last_i=PIXEL;
00263     
00264     if ( n_columns < 1 )
00265     {
00266         cpl_msg_error("coeffsCrossSlitFit:","wrong number of image columns given") ;
00267         return -1 ;
00268     }
00269     if ( acoefs == NULL || dacoefs == NULL )
00270     {
00271         cpl_msg_error("coeffsCrossSlitFit:","acoeffs or errors of coefficients are not given") ;
00272         return -1 ;
00273     }
00274     if ( bco == NULL )
00275     {
00276         cpl_msg_error("coeffsCrossSlitFit:","bcoeffs are not allocated") ;
00277         return -1 ;
00278     }
00279     if ( sigma_factor <= 0. )
00280     {
00281         cpl_msg_error("coeffsCrossSlitFit:","impossible sigma_factor given!") ;
00282         return -1 ;
00283     }
00284     if ( dispersion == 0. )
00285     {
00286         cpl_msg_error("coeffsCrossSlitFit:","impossible dispersion given!") ;
00287         return -1 ;
00288     }
00289 
00290     /*---------------------------------------------------------------------------------
00291      * search for the slitlet edges by comparing the a0 coefficients along the columns
00292      * if a bigger deviation occurrs it is assumed that there is an edge.
00293      */
00294     n = 0 ;
00295     threshold = pixel_dist * fabs(dispersion) ;
00296 
00297     for ( i = PIXEL ; i < n_columns - PIXEL ; )
00298     {
00299         if ( !qfits_isnan(acoefs[0][i+1]) && acoefs[0][i+1] != 0. && acoefs[0][i] != 0.
00300                   && dacoefs[0][i+1] != 0.)
00301         {
00302             if ( qfits_isnan(acoefs[0][i]) || acoefs[0][i] == 0. )
00303             {
00304                 if (fabs(acoefs[0][i+1] - acoefs[0][i-1]) >= threshold )
00305                 {
00306            if( (i-last_i) < 60 && (i > 80) ) {
00307              cpl_msg_info("wavecal","skip1 i=%d diff=%d\n",i,i-last_i);
00308              goto skip;
00309                    } else {  
00310                      /*
00311                  printf("diff1=%f i=%d threshold-%f size=%d\n", 
00312                      fabs(acoefs[0][i+1] - acoefs[0][i-1]),i,threshold,i-last_i);
00313              */
00314                      edge[n] = i+1 ;
00315                      n++ ;
00316                      last_i = i;
00317                      i += PIXEL ;
00318            }
00319                 }
00320             }
00321             else
00322             {
00323                 if (fabs(acoefs[0][i+1] - acoefs[0][i]) >= threshold )
00324                 {
00325            if( (i-last_i) < 60 && (i> 80)) {
00326              cpl_msg_info("wavecal","skip2 i=%d diff=%d\n",i,i-last_i);
00327              goto skip;
00328                    } else {  
00329                       /*
00330               printf("diff2=%f i=%d threshold-%f size=%d\n",
00331                       fabs(acoefs[0][i+1] - acoefs[0][i]),i,threshold,i-last_i);
00332               */
00333  
00334                   
00335                       edge[n] = i+1 ;
00336                       n++ ;
00337                       last_i = i;
00338                       i += PIXEL ;
00339            }
00340                 }
00341             }
00342         /* sometimes a slitlet may be lost due to divergences in acoeffs[0]
00343                we try to recover it */
00344             if( ( (i-last_i) > 63 ) && 
00345                 ( qfits_isnan(fabs(acoefs[0][i+1] - acoefs[0][i])) ||
00346                   qfits_isnan(fabs(acoefs[0][i+1] - acoefs[0][i-1])) ) ) 
00347           {
00348                    edge[n] = i+1 ;
00349                     n++ ;
00350                     last_i = i;
00351                  cpl_msg_warning("wavecal:","recovered slitlet edge i=%d",i);
00352                     i += PIXEL ;
00353 
00354           }
00355     }
00356       skip:
00357         i++ ;
00358     }
00359     /*
00360     printf("X min %d max %d last %d\n", PIXEL, n_columns - PIXEL, i);
00361     printf("n=%d check=%d\n",n,bco->n_slitlets - 1);
00362     */
00363     if ( n != bco->n_slitlets - 1 )
00364     {
00365         cpl_msg_error("coeffsCrossSlitFit:","could not find the right number of slitlets, found: %d\n",n+1) ;
00366         return -1 ;
00367     }
00368  
00369     /* go through the coefficents indices */
00370     for ( index = 0 ; index < bco->n_acoeffs ; index++ )
00371     {
00372         /* go through the single slitlets */
00373         for ( ns = 0 ; ns < bco->n_slitlets ; ns++ )
00374         {
00375             /* determine the slitlet edges */
00376             if ( ns == 0 )
00377             {
00378                 ed1 = 0 ;
00379                 ed2 = edge[0] ;
00380             }
00381             else if ( ns == bco->n_slitlets - 1 )
00382             {
00383                 ed1 = edge[bco->n_slitlets - 2] ;
00384                 ed2 = n_columns ;
00385             }
00386             else
00387             {
00388                 ed1 = edge[ns-1] ;
00389                 ed2 = edge[ns] ;
00390             }
00391 
00392             nc = 0 ;
00393             for ( i = ed1 ; i < ed2 ; i++ )
00394             {
00395                 if ( qfits_isnan(acoefs[index][i]) || acoefs[index][i] == 0. || dacoefs[index][i] == 0. )
00396                 {
00397                     continue ;
00398                 }
00399                 else
00400                 {
00401                     nc++ ;
00402                 }
00403             }
00404             if ( NULL == (acoefsclean = (float*) cpl_calloc(nc , sizeof(float))) )
00405             {
00406                 cpl_msg_error("coeffsCrossSlitFit:","could not allocate memory for acoefsclean!") ;
00407                 return -1 ;
00408             }
00409             nc = 0 ;
00410             for ( i = ed1 ; i < ed2 ; i++ )
00411             {
00412                 if ( qfits_isnan(acoefs[index][i]) || acoefs[index][i] == 0. || dacoefs[index][i] == 0. )
00413                 {
00414                     continue ;
00415                 }
00416                 else
00417                 {
00418                     acoefsclean[nc] = acoefs[index][i] ;
00419                     nc++ ;
00420                 }
00421             }
00422 
00423             /* ----------------------------------------------------------
00424              * determine the clean mean and sigma value of the coefficients,
00425              * that means reject 10 % of the extreme low and high values
00426              */
00427             pixel_qsort(acoefsclean, nc) ;
00428 
00429             sum   = 0. ;
00430             sumq  = 0. ;
00431             mean  = 0. ;
00432             sigma = 0. ;
00433             n     = 0 ;
00434             for ( i = (int)((float)nc*LOW_REJECT) ; i < (int)((float)nc*HIGH_REJECT) ; i++ )
00435             {
00436                 sum  += (double)acoefsclean[i] ;
00437                 sumq += ((double)acoefsclean[i] * (double)acoefsclean[i]) ;
00438                 n ++ ;
00439             }
00440             mean          = sum/(double)n ;
00441             sigma         = sqrt( sumq/(double)n - (mean * mean) ) ;
00442             cliphi        = mean + sigma * (double)sigma_factor ;
00443             cliplo        = mean - sigma * (double)sigma_factor ;
00444             /* fit only the reasonnable values */
00445             num = 0 ;
00446             col_index = 0 ;
00447             /*
00448             printf("ed1=%d ed2=%d\n",ed1,ed2);
00449         */
00450             for ( i = ed1 ; i < ed2 ; i++ )
00451             {
00452                 /* take only the reasonnable coefficients */
00453                  /*
00454         printf("acoeffs=%f dacoefs=%f cliphi=%f cliplo=%f\n",
00455                    acoefs[index][i],dacoefs[index][i],cliphi,cliplo);
00456          */
00457                 if ( !qfits_isnan(acoefs[index][i])     && 
00458                      (acoefs[index][i] <= cliphi) && 
00459                      (acoefs[index][i] >= cliplo) &&
00460                      (dacoefs[index][i] != 0. )   && 
00461                      (acoefs[index][i] != 0.)     )
00462                 {
00463                     sub_acoefs[num]    = acoefs[index][i] ;
00464                     sub_dacoefs[num]   = dacoefs[index][i] ;
00465                     sub_col_index[num] = col_index ;
00466                     num ++ ;
00467                 }
00468                 col_index++ ;
00469             }
00470             ndata = num ;
00471             offset = (float)(col_index-1) / 2. ;
00472             /* printf("ndata=%d bco->n_bcoeffs=%d\n",ndata,bco->n_bcoeffs); */
00473 
00474             if ( ndata < bco->n_bcoeffs )
00475             {
00476                 cpl_msg_error("coefsCrossSlitFit:","not enough data found in slitlet %d to determine the fit coefficients.", ns) ;
00477                 cpl_free(acoefsclean) ;
00478                 return -1 ;
00479             }
00480 
00481             /* allocate coefficient matrices */
00482             ucoefs = matrix(1, ndata, 1, bco->n_bcoeffs) ;
00483             vcoefs = matrix(1, ndata, 1, bco->n_bcoeffs) ;
00484             covar  = matrix(1, bco->n_bcoeffs, 1, bco->n_bcoeffs) ;
00485 
00486             /* scale the x-values for the fit */
00487             for ( i = 0 ; i < ndata ; i++ )
00488             {
00489                 sub_col_index[i] = (sub_col_index[i] - offset) / offset ;
00490             }
00491 
00492             /* finally, do the singular value decomposition fit */
00493             svd_fitting ( sub_col_index-1, sub_acoefs-1, sub_dacoefs-1, ndata, bco[ns].b[index]-1,
00494                      bco->n_bcoeffs, ucoefs, vcoefs, wcoefs-1, covar, &chisq[ns], fpol ) ;
00495 
00496             /* scale the found coefficients */
00497             for ( i = 0 ; i < bco->n_bcoeffs ; i ++ )
00498             {
00499                 bco[ns].b[index][i] /= pow( offset, i ) ;
00500             }
00501 
00502             /* free memory */
00503             cpl_free (acoefsclean) ;
00504             free_matrix( ucoefs, 1/*, ndata*/, 1/*, bco->n_bcoeffs */) ;
00505             free_matrix( vcoefs, 1/*, ndata*/, 1/*, bco->n_bcoeffs */) ;
00506             free_matrix( covar, 1/*, bco->n_bcoeffs*/, 1/*, bco->n_bcoeffs */) ;
00507 
00508             /* now calculate the smoothed acoefs for each column */
00509             col_index = 0 ;
00510             for ( i = ed1 ; i < ed2  ; i++ )
00511             {
00512                 acoefs[index][i] = 0. ;
00513                 for ( n = 0 ; n < bco->n_bcoeffs ; n++ )
00514                 {
00515                     acoefs[index][i] += bco[ns].b[index][n] * pow(col_index - offset, n) ;
00516                 }
00517                 col_index++ ;
00518             }
00519 
00520         }
00521     }
00522     return 0 ;
00523 }
00524 
00525 
00526 /*----------------------------------------------------------------------------
00527    Function     :       waveMapSlit()
00528    In           :       acoefs: fit coefficient matrix: output of coeffsCrossSlitFit()
00529                         n_acoefs: number of acoefs, polynomial order + 1
00530                         n_rows: number of final image rows
00531                         n_columns: number of final image columns
00532    Out          :       wavemap calibration map image
00533    Job          :       builds a new wavelength calibration map as fits image
00534                         by using the fit coeficients.
00535  ---------------------------------------------------------------------------*/
00536 
00537 OneImage * waveMapSlit ( float ** acoefs,
00538                          int      n_acoefs,
00539                          int      n_rows,
00540                          int      n_columns )
00541 {
00542     OneImage * newIm ;
00543     float lambda ;
00544     float offset ;
00545     int col, row ;
00546     int i ;
00547     float row_index ;
00548 
00549     if ( NULL == acoefs )
00550     {
00551         cpl_msg_error ("waveMapSlit:"," no coefficient matrix given!") ;
00552         return NullImage ;
00553     }
00554 
00555     /* allocate new image */
00556     if ( NullImage == (newIm = new_image(n_columns , n_rows)) )
00557     {
00558         cpl_msg_error ("waveMapSlit:","could not allocate new image!") ;
00559         return NullImage ;
00560     }
00561 
00562     /* make the parabola symmetric to the image */
00563     offset = (float)(n_rows - 1) / 2. ; 
00564 
00565     /* go through the rows */
00566     for ( col = 0 ; col < n_columns ; col++ )
00567     {
00568         /* go through the columns */
00569         for ( row = 0 ; row < n_rows ; row++ )
00570         {
00571             lambda = 0. ;
00572             row_index = (float)row - offset ;
00573             for ( i = 0 ; i < n_acoefs ; i++ )
00574             {
00575                 lambda += acoefs[i][col] * pow(row_index, i) ;
00576             }
00577             newIm -> data[col+row*n_columns] = lambda ;
00578         }
00579     }
00580     return newIm ;
00581 }
00582 
00583 
00584 /*----------------------------------------------------------------------------
00585    Function     :       waveCal()
00586    In           :       image:        merged image from a calibration emission lamp,
00587                         par:          fit parameters data structure storage
00588                         abuf:         buffer array for fit coefficients abuf[index][column]
00589                         row_clean: resulting list of the row indices but without the
00590                                    lines that are too close to each other for the fit
00591                                    output of findLines()
00592                         wavelength_clean: corrected wavelength list corresponding to
00593                                           the row_clean array
00594                                           output of findLines()
00595                         n_found_lines: output of findLines(): total number of found emission lines
00596                         dispersion:   dispersion of spectrum: micron per pixel
00597                         halfWidth:    half width of the box where the line must sit
00598                         minAmplitude: minimum amplitude of the Gaussian to do the fit
00599                         max_residual: maximum residual value, beyond that value
00600                                       the polynomial lambda-position fit is rejected.
00601                         fwhm:         first guess for the full width of half maximum
00602                                       of the gaussian line fit
00603                         n_a_fitcoefs: number of fit coefficients for the single
00604                                       column fits: lambda-position
00605                         n_b_fitcoefs: number of fit coefficients for the fits of
00606                                       the single a coefficients across the columns
00607                         sigmaFactor:  factor of the standard deviation of the determined
00608                                       polynomial coefficients of the columns beyond
00609                                       which these coefficients are not used to carry out
00610                                       the polynomial fit across the columns.
00611                         pixel_dist:   estimate of minimal pixel distance in spectral direction between slitle
00612 ts
00613                         pixel_tolerance: maximum tolerated difference between estimated
00614                                          and fitted line positions.
00615    Out          :       wavelength map image
00616                         abuf: array of smoothed coefficients of the polynomial fit along the columns:
00617                               abuf[index][column].
00618                         par: array of the resulting FitParams data structure
00619    Job          :       this routine takes an image from a calibration
00620                         emission lamp and delivers the smoothed fit coefficients of
00621                         a polynomial fit along the columns of the line positions as output.
00622                         This routine expects Nyquist sampled spectra
00623                         (either an interleaved image or an image convolved with an
00624                          appropriate function in spectral direction)
00625  ---------------------------------------------------------------------------*/
00626 OneImage * waveCal( OneImage   * image,
00627                     FitParams ** par ,
00628                     float     ** abuf,
00629                     int          n_slitlets,
00630                     int       ** row_clean,
00631                     float     ** wavelength_clean,
00632                     int        * n_found_lines,
00633                     float        dispersion,
00634                     int          halfWidth,
00635                     float        minAmplitude,
00636                     float        max_residual,
00637                     float        fwhm,
00638                     int          n_a_fitcoefs,
00639                     int          n_b_fitcoefs,
00640                     float        sigmaFactor,
00641                     float        pixel_dist,
00642                     float        pixel_tolerance )
00643 
00644 {
00645     int          i, j, k ;
00646     int          n_fit ;
00647     int          n_reject ;
00648     float     *  acoefs ;
00649     float     *  dacoefs ;
00650     float     ** dabuf ;
00651     float        chisq_poly ;
00652     float     *  chisq_cross ;
00653     int          zeroind ;
00654     int          crossInd ;
00655     Bcoeffs   *  bco ;
00656     OneImage  *  wavemap ;
00657 
00658 
00659     if (  NullImage == image )
00660     {
00661         cpl_msg_error("waveCal:","no image given") ;
00662         return NullImage ;
00663     }
00664     if ( par == NULL )
00665     {
00666         cpl_msg_error("waveCal:","no fit parameter data structure given") ;
00667         return NullImage ;
00668     }
00669     if ( abuf == NULL )
00670     {
00671         cpl_msg_error("waveCal:","no buffer for fit coefficients given") ;
00672         return NullImage ;
00673     }
00674     if ( n_slitlets <= 0 )
00675     {
00676         cpl_msg_error("waveCal:","impossible number of slitlets given") ;
00677         return NullImage ;
00678     }
00679     if ( row_clean == NULL )
00680     {
00681         cpl_msg_error("waveCal:","no row_clean array given") ;
00682         return NullImage ;
00683     }
00684     if ( wavelength_clean == NULL )
00685     {
00686         cpl_msg_error("waveCal:","no wavelength_clean array given") ;
00687         return NullImage ;
00688     }
00689 
00690     if ( dispersion == 0. )
00691     {
00692         cpl_msg_error("waveCal:","impossible dispersion given") ;
00693         return NullImage ;
00694     }
00695 
00696     if ( halfWidth <= 0 || halfWidth > image->ly/2 )
00697     {
00698         cpl_msg_error("waveCal:","impossible half width of the fitting box given") ;
00699         return NullImage ;
00700     }
00701     if ( minAmplitude < 1. )
00702     {
00703         cpl_msg_error("waveCal:","impossible minimal amplitude") ;
00704         return NullImage ;
00705     }
00706 
00707     if ( max_residual <= 0. || max_residual > 1. )
00708     {
00709         cpl_msg_error("waveCal:","impossible max_residual given") ;
00710         return NullImage ;
00711     }
00712     if ( fwhm <= 0. || fwhm > 10. )
00713     {
00714         cpl_msg_error("waveCal:","impossible fwhm given") ;
00715         return NullImage ;
00716     }
00717 
00718     if ( n_a_fitcoefs <= 0 || n_a_fitcoefs > 9 )
00719     {
00720         cpl_msg_error("waveCal:","unrealistic n_a_fitcoefs given") ;
00721         return NullImage ;
00722     }
00723 
00724     if ( n_b_fitcoefs <= 0 || n_b_fitcoefs > 9 )
00725     {
00726         cpl_msg_error("waveCal:","unrealistic n_b_fitcoefs given") ;
00727         return NullImage ;
00728     }
00729     if ( sigmaFactor <= 0. )
00730     {
00731         cpl_msg_error("waveCal:","impossible sigmaFactor given") ;
00732         return NullImage ;
00733     }
00734 
00735     /* initialize the variables */
00736     n_reject = 0 ;
00737     n_fit = 0 ;
00738 
00739     /* fit each found line by using a gaussian function and determine the exact position */
00740     if ( 0 > (n_fit = fitLines( image , par, fwhm, n_found_lines, row_clean, wavelength_clean,
00741                                 halfWidth, minAmplitude )) )
00742     {
00743         cpl_msg_error("waveCal:","cannot fit the lines, error code of fitLines: %d\n", n_fit) ;
00744         return NullImage ;
00745     }
00746 
00747     /* first check for faked lines like bad pixels */
00748     if ( -1 == checkForFakeLines (par, dispersion, wavelength_clean, row_clean, n_found_lines,
00749                                   image->lx, pixel_tolerance) )
00750     {
00751         cpl_msg_error("waveCal:","cannot fit the lines, error code of fitLines: %d\n", n_fit) ;
00752         return NullImage ;
00753     }
00754 
00755 
00756     /* allocate memory */
00757     if ( NULL == (acoefs = (float*) cpl_calloc (n_a_fitcoefs, sizeof(float))) ||
00758          NULL == (dacoefs = (float*) cpl_calloc (n_a_fitcoefs, sizeof(float))) ||
00759          NULL == (dabuf = (float**) cpl_calloc (n_a_fitcoefs, sizeof(float*))) ||
00760          NULL == (chisq_cross = (float*) cpl_calloc(n_slitlets, sizeof(float))) )
00761     {
00762         cpl_msg_error("waveCal:","cannot allocate memory\n") ;
00763         return NullImage ;
00764     }
00765     for ( i = 0 ; i < n_a_fitcoefs ; i++ )
00766     {
00767         if (  NULL == (dabuf[i] = (float*) cpl_calloc(image -> lx, sizeof(float))) )
00768         {
00769             cpl_msg_error("wavelengthCalibration:","cannot allocate memory") ;
00770             cpl_free ( acoefs ) ;
00771             cpl_free ( dacoefs ) ;
00772             cpl_free ( chisq_cross ) ;
00773             cpl_free(dabuf) ;
00774             return NullImage ;
00775         }
00776     }
00777 
00778     /* fit wavelengths to the corresponding found positions for each column */
00779     k = 0 ;
00780     for ( i = 0 ; i < image -> lx ; i++ )
00781     {
00782         zeroind = 0 ;
00783         if ( FLT_MAX == (chisq_poly = polyfit( par, i, n_found_lines[i], image -> ly, dispersion,
00784                                                max_residual, acoefs, dacoefs, &n_reject, n_a_fitcoefs)) )
00785         {
00786       /*
00787             cpl_msg_warning ("waveCal:","error in polyfit in column: %d\n", i) ;
00788       */
00789             for ( j = 0 ; j < n_a_fitcoefs ; j++ )
00790             {
00791                 acoefs[j] = ZERO ;
00792                 dacoefs[j] = ZERO ;
00793             }
00794         }
00795 
00796         for ( j = 0 ; j < n_a_fitcoefs ; j++ )
00797         {
00798 
00799             if ( acoefs[0] <= 0. || acoefs[1] ==0. ||
00800                  dacoefs[j] == 0. || qfits_isnan(acoefs[j]) )
00801             {
00802                 zeroind = 1 ;
00803             }
00804         }
00805         for ( j = 0 ; j < n_a_fitcoefs ; j++ )
00806         {
00807             if ( zeroind == 0 )
00808             {
00809                 abuf[j][i]  = acoefs[j] ;
00810                 dabuf[j][i] = dacoefs[j] ;
00811             }
00812             else
00813             {
00814                 abuf[j][i]  = ZERO ;
00815                 dabuf[j][i] = ZERO ;
00816             }
00817         }
00818     }
00819 
00820     /* allocate memory for the fitting coefficients */
00821     if ( NULL == ( bco = newBcoeffs( n_slitlets, n_a_fitcoefs, n_b_fitcoefs)) )
00822     {
00823         cpl_msg_error ("waveCal:","cannot allocate memory for the bcoeffs") ;
00824         for ( i = 0 ; i < n_a_fitcoefs ; i++ )
00825         {
00826             cpl_free (dabuf[i]) ;
00827         }
00828         cpl_free (dabuf) ;
00829         cpl_free ( acoefs ) ;
00830         cpl_free ( dacoefs ) ;
00831         cpl_free ( chisq_cross ) ;
00832         return NullImage ;
00833     }
00834 
00835     /* fit each acoefs across the slitlets to smooth the result */
00836     if ( -1 == ( crossInd = coeffsCrossSlitFit( image -> lx, 
00837                                                 abuf, 
00838                                                 dabuf,
00839                                                 bco, 
00840                                                 sigmaFactor, 
00841                                                 dispersion, 
00842                                                 pixel_dist, 
00843                                                 chisq_cross )) )
00844     {
00845         cpl_msg_error ("waveCal:","cannot carry out the fitting of coefficients across the columns") ;
00846         for ( i = 0 ; i < n_a_fitcoefs ; i++ )
00847         {
00848             cpl_free (dabuf[i]) ;
00849         }
00850 
00851         cpl_free (dabuf) ;
00852         cpl_free ( acoefs ) ;
00853         cpl_free ( dacoefs ) ;
00854         destroyBcoeffs(bco) ;
00855         cpl_free ( chisq_cross ) ;
00856         return NullImage ;
00857     }
00858 
00859     if ( NullImage == (wavemap = waveMapSlit (abuf, n_a_fitcoefs, image->ly, image->lx)) )
00860     {
00861         cpl_msg_error ("waveCal:","cannot carry out wavemap creation") ;
00862         for ( i = 0 ; i < n_a_fitcoefs ; i++ )
00863         {
00864             cpl_free (dabuf[i]) ;
00865         }
00866 
00867         cpl_free (dabuf) ;
00868         cpl_free ( acoefs ) ;
00869         cpl_free ( dacoefs ) ;
00870         destroyBcoeffs(bco) ;
00871         cpl_free ( chisq_cross ) ;
00872         return NullImage ;
00873     }
00874 
00875     /* free all allocated memory */
00876     for ( i = 0 ; i < n_a_fitcoefs ; i++ )
00877     {
00878         cpl_free (dabuf[i]) ;
00879     }
00880     cpl_free (dabuf) ;
00881     cpl_free ( acoefs ) ;
00882     cpl_free ( dacoefs ) ;
00883     destroyBcoeffs(bco) ;
00884     cpl_free ( chisq_cross ) ;
00885 
00886     return wavemap ;
00887 }
00888 
00889 
00890 
00891 /*----------------------------------------------------------------------------
00892    Function     :       checkForFakeLines()
00893    In           :       par: array of the fit parameter data structure FitParams
00894             dispersion: estimated dispersion
00895             wavelength_clean: corrected wavelength list
00896             row_clean: corrected row list corresponding to the wavelength list
00897             n_found_lines: array of numbers of found lines of each column
00898                         n_columns: total number of image columns
00899             pixel_tolerance: maximum tolerated difference between estimated 
00900                      and fitted line positions.
00901    Out          :       par: corrected FitParams 
00902             0 in case of success
00903             -1 in case of error
00904    Job          :       this routine searches for successfully fitted fake lines like
00905             bad pixels by comparing the found line positons with 
00906             estimated template positions. This routine should be
00907             inserted in the wavelength calibration routine just after
00908             the fitLines() routine.
00909  ---------------------------------------------------------------------------*/
00910 
00911 int checkForFakeLines ( FitParams ** par,
00912                         float        dispersion,
00913             float     ** wavelength_clean,
00914             int       ** row_clean,
00915             int        * n_found_lines,
00916             int          n_columns,
00917             float        pixel_tolerance )
00918 {
00919     int i,  k ;
00920     int col ;
00921     int found ;
00922     float row ;
00923     float * beginWave ;
00924     float firstWave ;
00925 
00926     if ( par == NULL )
00927     {
00928         cpl_msg_error("checkForFakeLines:","no fit parameter data structure given") ;
00929         return -1 ;
00930     }
00931     if ( dispersion == 0. )
00932     {
00933         cpl_msg_error("checkForFakeLines:","dispersion zero given!") ;
00934         return -1 ;
00935     }
00936     if ( wavelength_clean == NULL )
00937     {
00938         cpl_msg_error("checkForFakeLines:","no wavelength array given!") ;
00939         return -1 ;
00940     }
00941     if ( row_clean == NULL )
00942     {
00943         cpl_msg_error("checkForFakeLines:","no row array given!") ;
00944         return -1 ;
00945     }
00946     if ( n_found_lines == NULL )
00947     {
00948         cpl_msg_error("checkForFakeLines:","no number of lines given!") ;
00949         return -1 ;
00950     }
00951     if ( n_columns < 200 )
00952     {
00953         cpl_msg_error("checkForFakeLines:","wrong number of columns given!") ;
00954         return -1 ;
00955     }
00956 
00957     /* first determine the estimated beginning wavelength of the first row */
00958     for ( col = 0 ; col < n_columns ; col++ )
00959     {
00960         if ( n_found_lines[col] == 0 )
00961         {
00962         continue ;
00963         } 
00964         if ( NULL == (beginWave = (float*) cpl_calloc( n_found_lines[col], sizeof(float) ) ) )
00965         {
00966             cpl_msg_error("checkForFakeLines:","could not allocate memory!") ;
00967             return -1 ;
00968     }
00969     for ( k = 0 ; k < n_found_lines[col] ; k++ ) 
00970     {
00971         beginWave[k] = wavelength_clean[col][k] - (float)row_clean[col][k] * dispersion ;
00972         }
00973     /* determine the clean mean of the estimated beginning wavelengths of one column */
00974     if ( FLT_MAX == (firstWave = clean_mean (beginWave, n_found_lines[col], 10., 10.) ) )
00975     {
00976             cpl_msg_error("checkForFakeLines:","clean mean did not work!") ;
00977             return -1 ;
00978     }
00979 
00980         cpl_free (beginWave) ;
00981         /* go through the lines in that column and select the correct FitParam structure */
00982         for ( k = 0 ; k < n_found_lines[col] ; k++ ) 
00983         {    
00984         /* compute the estimated line position */
00985         row = ( wavelength_clean[col][k] - firstWave ) / dispersion ;
00986 
00987             /* go through all fit parameters and find the corresponding one */
00988         found = -1 ;
00989             for ( i = 0 ; i < (par[0] -> n_params) ; i ++ )
00990             {   
00991             /* find the given column and go through the lines in that column */
00992         if ( (par[i] -> column == col) && (par[i] -> line == k) && 
00993              (par[i] -> wavelength == wavelength_clean[col][k]) )
00994         {
00995                 found = i ;
00996             break ;
00997                 }
00998             }
00999         if ( found != -1 )
01000         {
01001             /* set fit params to zero where the fitted row position and the estimated 
01002             row positions are outside the tolerance */
01003                 if ( fabs(row - par[found]->fit_par[2]) > pixel_tolerance ) 
01004             {
01005                 cpl_msg_warning("checkForFakeLines:","found bad line in col: %d line: %d in row: %f difference: %f", col, k, par[found]->fit_par[2],row - par[found]->fit_par[2])  ;
01006                 par[found]->fit_par[2] = 0. ;
01007                 }
01008             }
01009         else
01010         {
01011         cpl_msg_warning("checkForFakeLines:","fit parameter of col %d and line no %d not found!\n", col, k ) ;
01012         }
01013         }
01014     }
01015 
01016     return 0 ;
01017 }
01018 
01019 /*----------------------------------------------------------------------------
01020    Function     :       createShiftedSlitWavemap()
01021    In           :       lineIm:       new shifted emission line frame
01022                         coeffs:       calculated polynomial coefficients,
01023                                       output of coeffsCrossSlitFit()
01024                         n_fitcoeffs:  number of polynomial coefficients, order + 1
01025                         wavelength:   wavelength list from line list
01026                         intensity:    corresponding line intensity from line list
01027                         n_lines:      number of lines in the list
01028    Out          :       wavelength map  
01029    Job          :       This routine cross-correlates a shifted emission line frames 
01030                         and determines the shift to the old one which is given by
01031                         its polynomial coefficients.
01032                         Then the a0 coefficients is recalculated and afterwards
01033                         a new wavelength calibration map is generated using the 
01034                         already calculated smoothed polynomial coefficients.
01035  ---------------------------------------------------------------------------*/
01036 
01037 OneImage * createShiftedSlitWavemap ( OneImage * lineIm,
01038                               float    ** coeffs,
01039                           int      n_fitcoeffs,
01040                                       float  * wavelength,
01041                                       float  * intensity,
01042                                       int      n_lines,
01043                                       int      magFactor )
01044 {
01045     OneImage * wavemap ;
01046     float emline[lineIm -> ly] ;
01047     float spec[lineIm -> ly] ;
01048     double  * result ;
01049     float * filter_spec ;
01050     float centreval ;
01051     float centrepix ;
01052     float cenpos, cenpix ;
01053     float pixvalue ;
01054     float wavelag ;
01055     float angst ;
01056     float wave[n_lines] ;
01057     float a_initial ;
01058     double a[n_fitcoeffs] ;
01059     float par[MAXPAR] ;
01060     float derv_par[MAXPAR] ;
01061     int        numpar, its ;
01062     int        * mpar ;
01063     float      tol, lab ;
01064     float      * xdat, * wdat ;
01065     Vector     * peak;
01066     int   iters, xdim, ndat ;
01067     int   row , col ;
01068     int   i, j, k/*, l, m*/ ;
01069     int   sign, found, line, width ;
01070     int var, maxlag, cmin, cmax ;
01071     double z[2*(n_fitcoeffs - 1)] ;
01072     gsl_poly_complex_workspace * w ;
01073     double xcorr_max ;
01074     int delta ;
01075 
01076 
01077     if ( lineIm == NullImage )
01078     {
01079         cpl_msg_error ("createShiftedSlitWavemap:","no input image given!") ;
01080         return NullImage ;
01081     }
01082     if ( coeffs == NULL )
01083     {
01084         cpl_msg_error ("createShiftedSlitWavemap:","no coefficient matrix given!") ;
01085         return NullImage ;
01086     }
01087     if ( n_fitcoeffs < 2 )
01088     {
01089         cpl_msg_error ("createShiftedSlitWavemap:","wrong number of polynomial coefficients given!") ;
01090         return NullImage ;
01091     }
01092     if ( wavelength == NULL || intensity == NULL )
01093     {
01094         cpl_msg_error ("createShiftedSlitWavemap:","no input image given!") ;
01095         return NullImage ;
01096     }
01097     if ( n_lines < 1 )
01098     {
01099         cpl_msg_error ("createShiftedSlitWavemap:","no input image given!") ;
01100         return NullImage ;
01101     }
01102 
01103     /* find out if Angstroem or microns are used */
01104     if ( wavelength[0] > 10000. )
01105     {
01106         /* Angstroem */
01107         angst = 10000. ;
01108     }
01109     else if ( wavelength[0] > 1000. && wavelength[0] < 10000. )
01110     {
01111         /* nanometers */
01112         angst = 1000. ;
01113     }
01114     else
01115     {
01116         /* microns */
01117         angst = 1. ;
01118     }
01119 
01120     /* allocate memory */
01121     if ( NullImage == (wavemap = new_image ( lineIm->lx, lineIm->ly)) )
01122     {
01123         cpl_msg_error ("createShiftedSlitWavemap:","could not allocate memory!") ;
01124         return NullImage ;
01125     }
01126     var = (magFactor-1)*(magFactor-1) ;
01127     
01128     /* first store each spectrum in a buffer */
01129     for ( col = 0 ; col < lineIm->lx ; col++ )
01130     {
01131         /* initialize the emline array for each column */
01132         for ( i = 0 ; i < lineIm -> ly ; i++ )
01133         {
01134             emline[i] = 0. ;
01135         }
01136         /* determine the coefficients by using the given bcoefs */
01137         for ( i = 0 ; i < n_fitcoeffs ; i++ )
01138         {
01139             /* initialize coefficients and solution */
01140             if (i < n_fitcoeffs-1)
01141             {
01142                 z[2*i] = 0. ;
01143                 z[2*i+1] = 0. ;
01144             }
01145             a[i] = coeffs[i][col] ;
01146         }
01147          
01148         a_initial = coeffs[0][col] ;
01149         /* go through the lines and generate an artificial spectrum */
01150         for ( line = 0 ; line < n_lines ; line++ )
01151         {
01152             /* go from Angstroem to micron */
01153             wave[line] = wavelength[line]/angst ;
01154 
01155             /* ---------------------------------------------------------------------
01156              * solve the polynomial for the exact offset of the line that means
01157              * find the root of the polynomial of order n_fitcoefs - 1
01158              */
01159             a[0] = a_initial - wave[line] ;
01160 
01161             if (NULL == (w = gsl_poly_complex_workspace_alloc(n_fitcoeffs))) 
01162             {
01163                 cpl_msg_error("createShiftedSlitWavemap:","could not allocate complex workspace!") ;
01164                 destroy_image(wavemap) ;
01165                 return NullImage ;
01166             }
01167             if (-1 == gsl_poly_complex_solve(a, n_fitcoeffs, w, z)) 
01168             {
01169                 cpl_msg_error("createShiftedSlitWavemap:","gsl_poly_complex_solve did not work!") ;
01170                 destroy_image(wavemap) ;
01171                 return NullImage ;
01172             }
01173             gsl_poly_complex_workspace_free(w) ;           
01174 
01175             j = 0 ;
01176             found = -1 ;
01177             for ( i = 0 ; i < n_fitcoeffs - 1 ; i++ )
01178             {
01179                 /* test for appropriate solution */
01180                 if( (z[2*i] > (-1.)*(float) lineIm->ly/2. &&
01181                      z[2*i] < (float)lineIm -> ly/2.) && z[2*i+1] == 0. )
01182                 {
01183                     found = 2*i ;
01184                     j ++ ;
01185                 }
01186                 else
01187                 {
01188                     continue ;
01189                 }
01190             }
01191             if ( j == 0 )
01192             {
01193                 cpl_msg_warning("createShiftedSlitWavemap:","no offset solution found for line %d in column %d", line, col) ;
01194                 continue ;
01195             }
01196             else if ( j == 1 )
01197             {
01198                 cenpos = z[found] + (float) lineIm->ly /2. ;
01199             }
01200             else
01201             {
01202                 cpl_msg_warning("createShiftedSlitWavemap:","two or more offset solutions found for line %d in column %d", line, col) ;
01203                 continue ;
01204             }
01205 
01206             /*---------------------------------------------------------------------------------
01207              * magnify image by the given factor add an additional offset
01208              */
01209             cenpix = cenpos ;
01210 
01211             /* determine max and min pixel limits over which line should be convolved */
01212             cmin = (nint(cenpix) - (var-1)) > 0 ? nint(cenpix) - (var-1) : 0 ;
01213             cmax = (nint(cenpix) + (var-1)) < lineIm -> ly ?
01214                     nint(cenpix) + (var-1) : lineIm -> ly ;
01215 
01216             /* convolve neon lines with Gaussian function */
01217             for ( j = cmin ; j < cmax ; j++ )
01218             {
01219                 emline[j] += intensity[line] * exp((double)(-0.5*(j-cenpix)*(j-cenpix))/(double)var) ;
01220             }
01221         }
01222 
01223         /*-----------------------------------------------------------------------------------------
01224          * for each column, map the image data points onto an magFactor times bigger element grid for
01225 FFT
01226          * in the cross correlation, first initialize the two helping arrays for each new column.
01227          */
01228         for ( k = 0 ; k < lineIm -> ly ; k++ )
01229         {
01230             spec[k] = 0. ;
01231         }
01232 
01233         /* now take the image data points of the column and put them into the spec array */
01234         for ( row = 0 ; row < lineIm -> ly ; row++ ) /* go through the column */
01235         {
01236             /* set bad pixels or negative values to zero */
01237             if (!qfits_isnan(lineIm -> data[col + row*lineIm -> lx]) &&
01238                 (lineIm -> data[col + row*lineIm -> lx] > 0.))
01239             {
01240                 spec[row] = lineIm -> data[col + row*lineIm -> lx] ;
01241             }
01242             else
01243             {
01244                 spec[row] = 0. ;
01245             }
01246         }
01247         /* convolve the spectrum by Gaussian */
01248         filter_spec = function1d_filter_lowpass(spec, lineIm->ly, LOW_PASS_GAUSSIAN, magFactor) ; 
01249 
01250         /* now call the cross correlation routine */
01251         result = xcorrel( filter_spec, lineIm->ly, emline, lineIm->ly,
01252                           lineIm->ly/2, &delta, &maxlag, &xcorr_max) ;
01253 
01254         if ( xcorr_max <= 0. )
01255         {
01256             cpl_msg_warning("createShiftedSlitWavemap:","no positive cross correlation sum , col %d set to ZERO \n", col) ;
01257             for ( row = 0 ; row < lineIm -> ly ; row++ )
01258             {
01259                 wavemap -> data[col + row*lineIm -> lx] = ZERO ;
01260             }
01261             function1d_del(filter_spec) ;
01262             cpl_free(result) ;
01263             continue ;
01264         }
01265 
01266        /* in this section, we fit the correlation function with a gauss, and find its peak, th
01267            us getting subpixel-accuracy */
01268 
01269         i = maxlag; j = i+1;
01270         while (result[j] < result[i])
01271         {
01272             i++; j++;
01273         }
01274         i = maxlag; k = i-1;
01275         while (result[k] < result[i])
01276         {
01277             i--; k--;
01278         }
01279         width = j-k+1;
01280         /* allocate memory for the spectral vector */
01281         if ( NULL == (peak = newVector (width)) )
01282         {
01283             cpl_msg_error ("determineShiftByCorrelation:","cannot allocate new Vector \n") ;
01284             function1d_del(filter_spec) ;
01285             cpl_free(result) ;
01286             return NullImage ;
01287         }
01288 
01289 
01290         /* allocate memory */
01291         xdat = (float *) cpl_calloc( peak -> n_elements, sizeof (float) ) ;
01292         wdat = (float *) cpl_calloc( peak -> n_elements, sizeof (float) ) ;
01293         mpar = (int *)   cpl_calloc( MAXPAR, sizeof (int) ) ;
01294 
01295         /* determine the values of the spectral vector given as input */
01296         /* go through the chosen column */
01297 
01298         for ( i = 0 ; i < width ; i++ )
01299         {
01300             peak -> data[i] = result[k+i]/xcorr_max * 100. ;
01301             xdat[i] = i;
01302             wdat[i] = 1.0;
01303         }
01304 
01305         /* set initial values for the fitting routine */
01306         xdim     = XDIM;
01307         ndat     = peak -> n_elements ;
01308         numpar   = MAXPAR ;
01309         tol      = TOL ;
01310         lab      = LAB ;
01311         its      = ITS ;
01312         par[1] = width/2.0 ;
01313         par[2] = (float) (maxlag - k) ;
01314         par[3] = (peak -> data[0] + peak -> data[peak->n_elements - 1]) / 2.0 ;
01315         par[0]  = result[maxlag]/xcorr_max * 100. - (par[3]) ;
01316 
01317         for ( i = 0 ; i < MAXPAR ; i++ )
01318         {
01319             derv_par[i] = 0.0 ;
01320             mpar[i] = 1 ;
01321         }
01322 
01323         /* finally, do the least square fit using a gaussian */
01324         if ( 0 > ( iters = lsqfit_c( xdat, &xdim, peak -> data, wdat, &ndat, par, derv_par, mpar,
01325                                      &numpar, &tol, &its, &lab )) )
01326         {
01327             cpl_msg_warning ("linefit:","lsqfit_c: least squares fit failed in col: %d, error no.: %d", col, iters) ;
01328             destroyVector ( peak ) ;
01329             cpl_free ( xdat ) ;
01330             cpl_free ( wdat ) ;
01331             cpl_free ( mpar ) ;
01332             function1d_del(filter_spec) ;
01333             cpl_free(result) ;
01334             continue ;
01335         }
01336 
01337         destroyVector ( peak ) ;
01338         cpl_free (xdat) ;
01339         cpl_free (wdat) ;
01340         cpl_free (mpar) ;
01341         function1d_del(filter_spec) ;
01342         cpl_free(result) ;
01343 
01344         wavelag =((float)lineIm->ly/2 - (float)k - par[2]) ;
01345 
01346         if ( fabs(wavelag) > (float)lineIm->ly/20. )
01347         {
01348             cpl_msg_warning("linefit","wavelag very big , col %d set to ZERO \n", col) ;
01349             for ( row = 0 ; row < lineIm -> ly ; row++ )
01350             {
01351                 wavemap -> data[col + row*lineIm -> lx] = ZERO ;
01352             }
01353             continue ;
01354         }
01355 
01356         /*-----------------------------------------------------------------------------
01357          * determine new zero order coefficient centreval, of which the formula is
01358          * determined by setting equal a polynomial shifted by wavelag with the same
01359          * higher order coefficients and set the new zero order coefficient to
01360          * get both sides of the equation approximately equal.
01361          */
01362         centreval = a_initial ;
01363         for ( i = 1 ; i < n_fitcoeffs ; i++ )
01364         {
01365             if ( i%2 == 0 ) 
01366             {
01367                 sign = -1 ;
01368             }
01369             else
01370             {
01371                 sign = 1 ;
01372             }
01373             centreval += (float)sign * coeffs[i][col]*pow(wavelag, i) ;
01374         }
01375 
01376         /* prepare to write out wavelength as pixel values */
01377         for ( row = 0 ; row < wavemap -> ly ; row++ )
01378         {
01379             centrepix = (float)row - ((float)wavemap->ly - 1.)/2. ;
01380             pixvalue = 0. ;
01381             for ( i = 1 ; i < n_fitcoeffs ; i++ )
01382             {
01383                 pixvalue += coeffs[i][col]*pow(centrepix, i) ;
01384             }
01385             wavemap -> data[col+row*wavemap->lx] =
01386                 centreval + pixvalue ;
01387         }
01388     }
01389     return wavemap ;
01390 }
01391 
01392 /*----------------------------------------------------------------------------
01393    Function     :       createShiftedSlitWavemap2()
01394    In           :       lineIm:       new shifted emission line frame
01395                         coeffs:       calculated polynomial coefficients,
01396                                       output of coeffsCrossSlitFit()
01397                         n_fitcoeffs:  number of polynomial coefficients, order + 1
01398                         wavelength:   wavelength list from line list
01399                         intensity:    corresponding line intensity from line list
01400                         n_lines:      number of lines in the list
01401                         magFactor:    magnification factor for help arrays
01402                         dispersion:   estimate of the dispersion
01403                         pixel_dist:   estimate of minimal pixel distance in spectral direction between sli
01404 tlets
01405    Out          :       wavelength map
01406    Job          :       This routine cross-correlates a shifted emission line frames
01407                         and determines the shift to the old one which is given by
01408                         its polynomial coefficients.
01409                         Then the a0 coefficients is recalculated and afterwards
01410                         a new wavelength calibration map is generated using the
01411                         already calculated smoothed polynomial coefficients.
01412  ---------------------------------------------------------------------------*/
01413 
01414 OneImage * createShiftedSlitWavemap2 ( OneImage * lineIm,
01415                                       float    ** coeffs,
01416                                       int      n_fitcoeffs,
01417                                       float  * wavelength,
01418                                       float  * intensity,
01419                                       int      n_lines,
01420                                       int      magFactor,
01421                                       float    dispersion,
01422                                       float    pixel_dist )
01423 {
01424     OneImage * wavemap ;
01425     float emline[lineIm -> ly] ;
01426     float spec[lineIm -> ly] ;
01427     double * result ;
01428     float * filter_spec ;
01429     float centreval ;
01430     float centrepix ;
01431     float cenpos, cenpix ;
01432     float pixvalue ;
01433     float wavelag ;
01434     /*float maxres ;*/
01435     float angst ;
01436     float wave[n_lines] ;
01437     /*float temp ;*/
01438     float a_initial ;
01439     double a[n_fitcoeffs] ;
01440     float par[MAXPAR] ;
01441     float derv_par[MAXPAR] ;
01442     int        numpar, its ;
01443     int        * mpar ;
01444     float      tol, lab ;
01445     float      * xdat, * wdat ;
01446     Vector     * peak;
01447     int   iters, xdim, ndat ;
01448     int   row , col ;
01449     int   i, j, k/*, l, m*/, n ;
01450     int   sign, found, line, width ;
01451     int var, maxlag, cmin, cmax ;
01452     float offset2 ;
01453     float threshold ;
01454     int edge[N_SLITLETS] ;
01455     int ed1, ed2 ;
01456     float a0[lineIm->lx] ;
01457     float a0_clean[lineIm->lx] ;
01458     float bcoef[N_SLITLETS][n_fitcoeffs] ;
01459     int ns, nc ;
01460     float * acoefsclean ;
01461     float col_index, sub_col_index[lineIm->lx] ;
01462     float sub_acoefs[lineIm->lx], sub_dacoefs[lineIm->lx] ;
01463     float wcoefs[n_fitcoeffs] ;
01464     float ** ucoefs, **vcoefs, **covar ;
01465     double sum, sumq, mean ;
01466     double sigma ;
01467     double cliphi, cliplo ;
01468     float chisq ;
01469     int num, ndata ;
01470     double z[2*(n_fitcoeffs - 1)] ;
01471     gsl_poly_complex_workspace * w ;
01472     double xcorr_max ;
01473     int delta ;
01474 
01475     if ( lineIm == NullImage )
01476     {
01477         cpl_msg_error ("createShiftedSlitWavemap2:"," no input image given!\n") ;
01478         return NullImage ;
01479     }
01480     if ( coeffs == NULL )
01481     {
01482         cpl_msg_error ("createShiftedSlitWavemap2:"," no coefficient matrix given!\n") ;
01483         return NullImage ;
01484     }
01485     if ( n_fitcoeffs < 2 )
01486     {
01487         cpl_msg_error ("createShiftedSlitWavemap2:"," wrong number of polynomial coefficients given!\n") ;
01488         return NullImage ;
01489     }
01490     if ( wavelength == NULL || intensity == NULL )
01491     {
01492         cpl_msg_error ("createShiftedSlitWavemap2:"," no input image given!\n") ;
01493         return NullImage ;
01494     }
01495     if ( n_lines < 1 || magFactor < 1 )
01496     {
01497         cpl_msg_error ("createShiftedSlitWavemap2:"," no input image given!\n") ;
01498         return NullImage ;
01499     }
01500     var    = (magFactor - 1)*(magFactor - 1) ;
01501     /* find out if Angstroem or microns are used */
01502     if ( wavelength[0] > 10000. )
01503     {
01504         /* Angstroem */
01505         angst = 10000. ;
01506     }
01507     else if ( wavelength[0] > 1000. && wavelength[0] < 10000. )
01508     {
01509         /* nanometers */
01510         angst = 1000. ;
01511     }
01512     else
01513     {
01514         /* microns */
01515         angst = 1. ;
01516     }
01517 
01518 
01519     /* find the slitlet edges from the a0 coefficient */
01520     n = 0 ;
01521     threshold = pixel_dist * fabs(dispersion) ;
01522     for ( i = PIXEL ; i < lineIm->lx - PIXEL ; )
01523     {
01524         if (fabs(coeffs[0][i+1] - coeffs[0][i]) >= threshold )
01525         {
01526             edge[n] = i+1 ;
01527             n++ ;
01528             i += PIXEL ;
01529         }
01530         i++ ;
01531     }
01532 
01533 
01534     /* allocate memory */
01535     if ( NullImage == (wavemap = new_image ( lineIm->lx, lineIm->ly)) )
01536     {
01537         cpl_msg_error ("createShiftedSlitWavemap2:"," could not allocate memory!\n") ;
01538         return NullImage ;
01539     }
01540 
01541     /* first store each spectrum in a buffer */
01542     for ( col = 0 ; col < lineIm->lx ; col++ )
01543     {
01544         /* initialize the emline array for each column */
01545         for ( i = 0 ; i < lineIm -> ly ; i++ )
01546         {
01547             emline[i] = 0. ;
01548         }
01549         /* determine the coefficients by using the given bcoefs */
01550         for ( i = 0 ; i < n_fitcoeffs ; i++ )
01551         {
01552             /* initialize coefficients and solution */
01553             if (i < n_fitcoeffs-1)
01554             {
01555                 z[2*i] = 0. ;
01556                 z[2*i+1] = 0. ;
01557             }
01558             a[i] = coeffs[i][col] ;
01559         }
01560          
01561         a_initial = coeffs[0][col] ;
01562         /* go through the lines and generate an artificial spectrum */
01563         for ( line = 0 ; line < n_lines ; line++ )
01564         {
01565             /* go from Angstroem to micron */
01566             wave[line] = wavelength[line]/angst ;
01567 
01568             /* ---------------------------------------------------------------------
01569              * solve the polynomial for the exact offset of the line that means
01570              * find the root of the polynomial of order n_fitcoefs - 1
01571              */
01572             a[0] = a_initial - wave[line] ;
01573 
01574             if (NULL == (w = gsl_poly_complex_workspace_alloc(n_fitcoeffs))) 
01575             {
01576                 cpl_msg_error("createShiftedSlitWavemap2:"," could not allocate complex workspace!") ;
01577                 destroy_image(wavemap) ;
01578                 return NullImage ;
01579             }
01580             if (-1 == gsl_poly_complex_solve(a, n_fitcoeffs, w, z)) 
01581             {
01582                 cpl_msg_error("createShiftedSlitWavemap2:"," gsl_poly_complex_solve did not work!") ;
01583                 destroy_image(wavemap) ;
01584                 return NullImage ;
01585             }
01586             gsl_poly_complex_workspace_free(w) ;           
01587 
01588             j = 0 ;
01589             found = -1 ;
01590             for ( i = 0 ; i < n_fitcoeffs - 1 ; i++ )
01591             {
01592                 /* test for appropriate solution */
01593                 if( (z[2*i] > (-1.)*(float) lineIm->ly/2. &&
01594                      z[2*i] < (float)lineIm -> ly/2.) && z[2*i+1] == 0. )
01595                 {
01596                     found = 2*i ;
01597                     j ++ ;
01598                 }
01599                 else
01600                 {
01601                     continue ;
01602                 }
01603             }
01604             if ( j == 0 )
01605             {
01606                 cpl_msg_warning("createShiftedSlitWavemap2:"," no offset solution found for line %d in column %d\n", line, col) ;
01607                 continue ;
01608             }
01609             else if ( j == 1 )
01610             {
01611                 cenpos = z[found] + (float) lineIm->ly/2. ;
01612             }
01613             else
01614             {
01615                 cpl_msg_warning("createShiftedSlitWavemap2:"," two or more offset solutions found for line %d in column %d\n", line, col) ;
01616                 continue ;
01617             }
01618 
01619             /*---------------------------------------------------------------------------------
01620              * magnify image by the given factor add an additional offset
01621              */
01622             cenpix = cenpos ;
01623 
01624             /* determine max and min pixel limits over which line should be convolved */
01625             cmin = (nint(cenpix) - (var-1)) > 0 ? nint(cenpix) - (var-1) : 0 ;
01626             cmax = (nint(cenpix) + (var-1)) < lineIm -> ly ?
01627                     nint(cenpix) + (var-1) : lineIm -> ly ;
01628 
01629             /* convolve neon lines with Gaussian function */
01630             for ( j = cmin ; j < cmax ; j++ )
01631             {
01632                 emline[j] += intensity[line] * exp((double)(-0.5*(j-cenpix)*(j-cenpix))/(double)var) ;
01633             }
01634         }
01635 
01636         /*-----------------------------------------------------------------------------------------
01637          * for each column, map the image data points onto an magFactor times bigger element grid for
01638 FFT
01639          * in the cross correlation, first initialize the two helping arrays for each new column.
01640          */
01641         for ( k = 0 ; k < lineIm -> ly ; k++ )
01642         {
01643             spec[k] = 0. ;
01644         }
01645 
01646         /* now take the image data points of the column and put them into the spec array */
01647         for ( row = 0 ; row < lineIm -> ly ; row++ ) /* go through the column */
01648         {
01649             /* set bad pixels or negative values to zero */
01650             if (!qfits_isnan(lineIm -> data[col + row*lineIm -> lx]) &&
01651                 (lineIm -> data[col + row*lineIm -> lx] > 0.))
01652             {
01653                 spec[row] = lineIm -> data[col + row*lineIm -> lx] ;
01654             }
01655             else
01656             {
01657                 spec[row] = 0. ;
01658             }
01659         }
01660         /* convolve the spectrum by Gaussian */
01661         filter_spec = function1d_filter_lowpass(spec, lineIm->ly, LOW_PASS_GAUSSIAN, magFactor) ; 
01662 
01663         /* now call the cross correlation routine */
01664         result = xcorrel( filter_spec, lineIm->ly, emline, lineIm->ly,
01665                           lineIm->ly/2, &delta, &maxlag, &xcorr_max) ;
01666 
01667         if ( xcorr_max <= 0. )
01668         {
01669             cpl_msg_warning("cross-correl","no positive cross correlation sum , col %d set to ZERO \n", col) ;
01670             for ( row = 0 ; row < lineIm -> ly ; row++ )
01671             {
01672                 wavemap -> data[col + row*lineIm -> lx] = ZERO ;
01673             }
01674             function1d_del(filter_spec) ;
01675             cpl_free(result) ;
01676             continue ;
01677         }
01678 
01679        /* in this section, we fit the correlation function with a gauss, and find its peak, th
01680            us getting subpixel-accuracy */
01681 
01682         i = maxlag; j = i+1;
01683         while (result[j] < result[i])
01684         {
01685             i++; j++;
01686         }
01687         i = maxlag; k = i-1;
01688         while (result[k] < result[i])
01689         {
01690             i--; k--;
01691         }
01692         width = j-k+1;
01693         /* allocate memory for the spectral vector */
01694         if ( NULL == (peak = newVector (width)) )
01695         {
01696             cpl_msg_error ("createShiftedWavemap2:"," cannot allocate new Vector \n") ;
01697             function1d_del(filter_spec) ;
01698             cpl_free(result) ;
01699             return NullImage ;
01700         }
01701 
01702 
01703         /* allocate memory */
01704         xdat = (float *) cpl_calloc( peak -> n_elements, sizeof (float) ) ;
01705         wdat = (float *) cpl_calloc( peak -> n_elements, sizeof (float) ) ;
01706         mpar = (int *)   cpl_calloc( MAXPAR, sizeof (int) ) ;
01707 
01708         /* determine the values of the spectral vector given as input */
01709         /* go through the chosen column */
01710 
01711         for ( i = 0 ; i < width ; i++ )
01712         {
01713             peak -> data[i] = result[k+i]/xcorr_max * 100. ;
01714             xdat[i] = i;
01715             wdat[i] = 1.0;
01716         }
01717 
01718         /* set initial values for the fitting routine */
01719         xdim     = XDIM;
01720         ndat     = peak -> n_elements ;
01721         numpar   = MAXPAR ;
01722         tol      = TOL ;
01723         lab      = LAB ;
01724         its      = ITS ;
01725         par[1] = width/2.0 ;
01726         par[2] = (float) (maxlag - k) ;
01727         par[3] = (peak -> data[0] + peak -> data[peak->n_elements - 1]) / 2.0 ;
01728         par[0]  = result[maxlag]/xcorr_max * 100. - (par[3]) ;
01729 
01730         for ( i = 0 ; i < MAXPAR ; i++ )
01731         {
01732             derv_par[i] = 0.0 ;
01733             mpar[i] = 1 ;
01734         }
01735 
01736         /* finally, do the least square fit using a gaussian */
01737         if ( 0 > ( iters = lsqfit_c( xdat, &xdim, peak -> data, wdat, &ndat, par, derv_par, mpar,
01738                                      &numpar, &tol, &its, &lab )) )
01739         {
01740             cpl_msg_warning ("linefit:"," lsqfit_c: least squares fit failed in col: %d, error no.: %d\n", col, iters) ;
01741             destroyVector ( peak ) ;
01742             cpl_free ( xdat ) ;
01743             cpl_free ( wdat ) ;
01744             cpl_free ( mpar ) ;
01745             function1d_del(filter_spec) ;
01746             cpl_free(result) ;
01747             continue ;
01748         }
01749 
01750         destroyVector ( peak ) ;
01751         cpl_free (xdat) ;
01752         cpl_free (wdat) ;
01753         cpl_free (mpar) ;
01754         function1d_del(filter_spec) ;
01755         cpl_free(result) ;
01756 
01757         wavelag =((float)lineIm->ly/2 - (float)k - par[2]) ;
01758 
01759         if ( fabs(wavelag) > (float)lineIm->ly/20. )
01760         {
01761             cpl_msg_warning("lfit","wavelag very big , col %d set to ZERO \n", col) ;
01762             for ( row = 0 ; row < lineIm -> ly ; row++ )
01763             {
01764                 wavemap -> data[col + row*lineIm -> lx] = ZERO ;
01765             }
01766             continue ;
01767         }
01768 
01769         /*-----------------------------------------------------------------------------
01770          * determine new zero order coefficient centreval, of which the formula is
01771          * determined by setting equal a polynomial shifted by wavelag with the same
01772          * higher order coefficients and set the new zero order coefficient to
01773          * get both sides of the equation approximately equal.
01774          */
01775         centreval = a_initial ;
01776         for ( i = 1 ; i < n_fitcoeffs ; i++ )
01777         {
01778             if ( i%2 == 0 )
01779             {
01780                 sign = -1 ;
01781             }
01782             else
01783             {
01784                 sign = 1 ;
01785             }
01786             centreval += (float)sign * coeffs[i][col]*pow(wavelag, i) ;
01787         }
01788         a0[col] = centreval ;
01789     }
01790 
01791     /* go through the single slitlets */
01792     for ( ns = 0 ; ns < N_SLITLETS ; ns++ )
01793     {
01794         /* determine the slitlet edges */
01795         if ( ns == 0 )
01796         {
01797             ed1 = 0 ;
01798             ed2 = edge[0] ;
01799         }
01800         else if ( ns == N_SLITLETS - 1 )
01801         {
01802             ed1 = edge[N_SLITLETS - 2] ;
01803             ed2 = lineIm->lx ;
01804         }
01805         else
01806         {
01807             ed1 = edge[ns-1] ;
01808             ed2 = edge[ns] ;
01809         }
01810 
01811         nc = 0 ;
01812         for ( i = ed1 ; i < ed2 ; i++ )
01813         {
01814             if ( qfits_isnan(a0[i]) || a0[i] == 0. )
01815             {
01816                 continue ;
01817             }
01818             else
01819             {
01820                 nc++ ;
01821             }
01822         }
01823         if ( NULL == (acoefsclean = (float*) cpl_calloc(nc , sizeof(float))) )
01824         {
01825             cpl_msg_error("createShiftedWavemap2:","could not allocate memory for acoefsclean!\n") ;
01826             return NullImage ;
01827         }
01828         nc = 0 ;
01829         for ( i = ed1 ; i < ed2 ; i++ )
01830         {
01831             if ( qfits_isnan(a0[i]) || a0[i] == 0. )
01832             {
01833                 continue ;
01834             }
01835             else
01836             {
01837                 acoefsclean[nc] = a0[i] ;
01838                 nc++ ;
01839             }
01840         }
01841 
01842         /* ----------------------------------------------------------
01843          * determine the clean mean and sigma value of the coefficients,
01844          * that means reject 10 % of the extreme low and high values
01845          */
01846         pixel_qsort(acoefsclean, nc) ;
01847         sum   = 0. ;
01848         sumq  = 0. ;
01849         mean  = 0. ;
01850         sigma = 0. ;
01851         n     = 0 ;
01852         for ( i = (int)((float)nc*LOW_REJECT) ; i < (int)((float)nc*HIGH_REJECT) ; i++ )
01853         {
01854             sum  += (double)acoefsclean[i] ;
01855             sumq += ((double)acoefsclean[i] * (double)acoefsclean[i]) ;
01856             n ++ ;
01857         }
01858         mean          = sum/(double)n ;
01859         sigma         = sqrt( sumq/(double)n - (mean * mean) ) ;
01860         cliphi        = mean + sigma * (double)3. ;
01861         cliplo        = mean - sigma * (double)3. ;
01862         /* fit only the reasonnable values */
01863         num = 0 ;
01864         col_index = 0 ;
01865         for ( i = ed1 ; i < ed2 ; i++ )
01866         {
01867             /* take only the reasonnable coefficients */
01868             if ( !qfits_isnan(a0[i]) && (a0[i] <= cliphi) && (a0[i] >= cliplo) &&
01869                  (a0[i] != 0.) )
01870             {
01871                 sub_acoefs[num]    = a0[i] ;
01872                 sub_dacoefs[num]   = 0.0000005 ;
01873                 sub_col_index[num] = col_index ;
01874                 num ++ ;
01875             }
01876             col_index++ ;
01877         }
01878         ndata = num ;
01879         offset2 = (float)(col_index-1) / 2. ;
01880 
01881         if ( ndata < n_fitcoeffs )
01882         {
01883             cpl_msg_error("createShiftedWavemap2:"," not enough data found in slitlet %d\
01884                     to determine the fit coefficients.\n", ns) ;
01885             cpl_free(acoefsclean) ;
01886             return NullImage ;
01887         }
01888 
01889         /* allocate coefficient matrices, see numerical recipe function matrix */
01890         ucoefs = matrix(1, ndata, 1, n_fitcoeffs) ;
01891         vcoefs = matrix(1, ndata, 1, n_fitcoeffs) ;
01892         covar  = matrix(1, n_fitcoeffs, 1, n_fitcoeffs) ;
01893 
01894         /* scale the x-values for the fit */
01895         for ( i = 0 ; i < ndata ; i++ )
01896         {
01897             sub_col_index[i] = (sub_col_index[i] - offset2) / offset2 ;
01898         }
01899 
01900         /* finally, do the singular value decomposition fit */
01901         svd_fitting ( sub_col_index-1, sub_acoefs-1, sub_dacoefs-1, ndata, bcoef[ns]-1,
01902                  n_fitcoeffs, ucoefs, vcoefs, wcoefs-1, covar, &chisq, fpol ) ;
01903 
01904         /* scale the found coefficients */
01905         for ( i = 0 ; i < n_fitcoeffs ; i ++ )
01906         {
01907             bcoef[ns][i] /= pow( offset2, i ) ;
01908         }
01909 
01910         /* free memory */
01911         cpl_free (acoefsclean) ;
01912         free_matrix( ucoefs, 1/*, ndata*/, 1/*, n_fitcoeffs */) ;
01913         free_matrix( vcoefs, 1/*, ndata*/, 1/*, n_fitcoeffs */) ;
01914         free_matrix( covar, 1/*, n_fitcoeffs*/, 1/*, n_fitcoeffs*/) ;
01915 
01916         /* now calculate the smoothed acoefs for each column */
01917         col_index = 0 ;
01918         for ( i = ed1 ; i < ed2 ; i++ )
01919         {
01920             a0_clean[i] = 0. ;
01921             for ( n = 0 ; n < n_fitcoeffs ; n++ )
01922             {
01923                 a0_clean[i] += bcoef[ns][n] * pow((float)col_index - offset2, n) ;
01924             }
01925             col_index++ ;
01926         }
01927 
01928     }
01929 
01930     for ( col = 0 ; col < lineIm->lx ; col++ )
01931     {
01932         /* prepare to write out wavelength as pixel values */
01933         for ( row = 0 ; row < wavemap -> ly ; row++ )
01934         {
01935             centrepix = (float)row - ((float)wavemap->ly - 1.)/2. ;
01936             pixvalue = 0. ;
01937             for ( i = 1 ; i < n_fitcoeffs ; i++ )
01938             {
01939                 pixvalue += coeffs[i][col]*pow(centrepix, i) ;
01940             }
01941             wavemap -> data[col+row*wavemap->lx] =
01942                 a0_clean[col] + pixvalue ;
01943         }
01944     }
01945     return wavemap ;
01946 }
01947 
01948 /*----------------------------------------------------------------------------
01949    Function     :       createShiftedSlitWavemap3()
01950    In           :       lineIm:       new shifted emission line frame
01951                         coeffs:       calculated polynomial coefficients,
01952                                       output of coeffsCrossSlitFit()
01953                         n_fitcoeffs:  number of polynomial coefficients, order + 1
01954                         wavelength:   wavelength list from line list
01955                         intensity:    corresponding line intensity from line list
01956                         n_lines:      number of lines in the list
01957                         magFactor:    magnification factor for help arrays
01958    Out          :       wavelength map  
01959    Job          :       This routine cross-correlates a shifted emission line frames 
01960                         and determines the shift to the old one which is given by
01961                         its polynomial coefficients.
01962                         Then the a0 coefficients is recalculated by using a clean mean
01963                         of the determined offset (shift) over the whole frame
01964                         and afterwards a new wavelength calibration map is generated using the 
01965                         already calculated smoothed polynomial coefficients.
01966  ---------------------------------------------------------------------------*/
01967 
01968 OneImage * createShiftedSlitWavemap3 ( OneImage * lineIm,
01969                               float    ** coeffs,
01970                           int      n_fitcoeffs,
01971                                       float  * wavelength,
01972                                       float  * intensity,
01973                                       int      n_lines,
01974                                       int      magFactor )
01975 {
01976 
01977     OneImage * wavemap ;
01978     float emline[lineIm -> ly] ;
01979     float spec[lineIm -> ly] ;
01980     double * result ;
01981     float * filter_spec ;
01982     float centreval ;
01983     float centrepix ;
01984     float cenpos, cenpix ;
01985     float pixvalue ;
01986     float wavelag[lineIm->lx] ;
01987     float wavelag_mean ;
01988     /*float maxres ;*/
01989     float angst ;
01990     float wave[n_lines] ;
01991     /*float temp ;*/
01992     float a_initial ;
01993     /*float solution[n_fitcoeffs-1], im_solution[n_fitcoeffs-1] ;*/
01994     double a[n_fitcoeffs] ;
01995     float par[MAXPAR] ;
01996     float derv_par[MAXPAR] ;
01997     int        numpar, its ;
01998     int        * mpar ;
01999     float      tol, lab ;
02000     float      * xdat, * wdat ;
02001     Vector     * peak;
02002     int   iters, xdim, ndat ;
02003     int   row , col ;
02004     int   i, j, k/*, l, m, n*/ ;
02005     int   sign, found, line, width ;
02006     int var, maxlag, cmin, cmax ;
02007     double z[2*(n_fitcoeffs - 1)] ;
02008     gsl_poly_complex_workspace * w ;
02009     double xcorr_max ;
02010     int delta ;
02011     
02012     if ( lineIm == NullImage )
02013     {
02014         cpl_msg_error ("createShiftedSlitWavemap3:"," no input image given!\n") ;
02015         return NullImage ;
02016     }
02017     if ( coeffs == NULL )
02018     {
02019         cpl_msg_error ("createShiftedSlitWavemap3:"," no coefficient matrix given!\n") ;
02020         return NullImage ;
02021     }
02022     if ( n_fitcoeffs < 2 )
02023     {
02024         cpl_msg_error ("createShiftedSlitWavemap3:"," wrong number of polynomial coefficients given!\n") ;
02025         return NullImage ;
02026     }
02027 
02028     if ( wavelength == NULL || intensity == NULL )
02029     {
02030         cpl_msg_error ("createShiftedSlitWavemap3:"," no wavelength list given!\n") ;
02031         return NullImage ;
02032     }
02033     if ( n_lines < 1 || magFactor < 1 )
02034     {
02035         cpl_msg_error ("createShiftedSlitWavemap3:"," wrong n_lines or magFactor given!\n") ;
02036         return NullImage ;
02037     }
02038     
02039     var    = (magFactor - 1)*(magFactor - 1) ;
02040     /* find out if Angstroem or microns are used */
02041     if ( wavelength[0] > 10000. )
02042     {
02043         /* Angstroem */
02044         angst = 10000. ;
02045     }
02046     else if ( wavelength[0] > 1000. && wavelength[0] < 10000. )
02047     {
02048         /* nanometers */
02049         angst = 1000. ;
02050     }
02051     else
02052     {
02053         /* microns */
02054         angst = 1. ;
02055     }
02056 
02057 
02058 
02059     /* allocate memory */
02060     if ( NullImage == (wavemap = new_image ( lineIm->lx, lineIm->ly)) )
02061     {
02062         cpl_msg_error ("createShiftedSlitWavemap3:"," could not allocate memory!\n") ;
02063         return NullImage ;
02064     }
02065 
02066     
02067     /* first store each spectrum in a buffer */
02068     for ( col = 0 ; col < lineIm->lx ; col++ )
02069     {
02070         /* initialize the emline array for each column */
02071         for ( i = 0 ; i < lineIm -> ly ; i++ )
02072         {
02073             emline[i] = 0. ;
02074         }
02075         /* determine the coefficients by using the given bcoefs */
02076         for ( i = 0 ; i < n_fitcoeffs ; i++ )
02077         {
02078             /* initialize coefficients and solution */
02079             if (i < n_fitcoeffs-1)
02080             {
02081                 z[2*i] = 0. ;
02082                 z[2*i+1] = 0. ;
02083             }
02084             a[i] = coeffs[i][col] ;
02085         }
02086          
02087         a_initial = coeffs[0][col] ;
02088         /* go through the lines and generate an artificial spectrum */
02089         for ( line = 0 ; line < n_lines ; line++ )
02090         {
02091             /* go from Angstroem to micron */
02092             wave[line] = wavelength[line]/angst ;
02093 
02094             /* ---------------------------------------------------------------------
02095              * solve the polynomial for the exact offset of the line that means
02096              * find the root of the polynomial of order n_fitcoefs - 1
02097              */
02098             a[0] = a_initial - wave[line] ;
02099 
02100             if (NULL == (w = gsl_poly_complex_workspace_alloc(n_fitcoeffs))) 
02101             {
02102                 cpl_msg_error("createShiftedSlitWavemap3:"," could not allocate complex workspace!") ;
02103                 destroy_image(wavemap) ;
02104                 return NullImage ;
02105             }
02106             if (-1 == gsl_poly_complex_solve(a, n_fitcoeffs, w, z)) 
02107             {
02108                 cpl_msg_error("createShiftedSlitWavemap3:"," gsl_poly_complex_solve did not work!") ;
02109                 destroy_image(wavemap) ;
02110                 return NullImage ;
02111             }
02112             gsl_poly_complex_workspace_free(w) ;           
02113 
02114             j = 0 ;
02115             found = -1 ;
02116             for ( i = 0 ; i < n_fitcoeffs - 1 ; i++ )
02117             {
02118                 /* test for appropriate solution */
02119                 if( (z[2*i] > (-1.)*(float) lineIm->ly/2. &&
02120                      z[2*i] < (float)lineIm -> ly/2.) && z[2*i+1] == 0. )
02121                 {
02122                     found = 2*i ;
02123                     j ++ ;
02124                 }
02125                 else
02126                 {
02127                     continue ;
02128                 }
02129             }
02130             if ( j == 0 )
02131             {
02132                 cpl_msg_warning("createShiftedSlitWavemap3:"," no offset solution found for line %d in column %d\n", line, col) ;
02133                 continue ;
02134             }
02135             else if ( j == 1 )
02136             {
02137                 cenpos = z[found] + (float) lineIm->ly /2. ;
02138             }
02139             else
02140             {
02141                 cpl_msg_warning("createShiftedSlitWavemap3:"," two or more offset solutions found for line %d in column %d\n", line, col) ;
02142                 continue ;
02143             }
02144 
02145             /*---------------------------------------------------------------------------------
02146              * magnify image by the given factor add an additional offset
02147              */
02148             cenpix = cenpos ;
02149 
02150             /* determine max and min pixel limits over which line should be convolved */
02151             cmin = (nint(cenpix) - (var-1)) > 0 ? nint(cenpix) - (var-1) : 0 ;
02152             cmax = (nint(cenpix) + (var-1)) < lineIm -> ly ?
02153                     nint(cenpix) + (var-1) : lineIm -> ly ;
02154 
02155             /* convolve neon lines with Gaussian function */
02156             for ( j = cmin ; j < cmax ; j++ )
02157             {
02158                 emline[j] += intensity[line] * exp((double)(-0.5*(j-cenpix)*(j-cenpix))/(double)var) ;
02159             }
02160         }
02161 
02162         /*-----------------------------------------------------------------------------------------
02163          * for each column, map the image data points onto an magFactor times bigger element grid for
02164 FFT
02165          * in the cross correlation, first initialize the two helping arrays for each new column.
02166          */
02167         for ( k = 0 ; k < lineIm -> ly ; k++ )
02168         {
02169             spec[k] = 0. ;
02170         }
02171 
02172         /* now take the image data points of the column and put them into the spec array */
02173         for ( row = 0 ; row < lineIm -> ly ; row++ ) /* go through the column */
02174         {
02175             /* set bad pixels or negative values to zero */
02176             if (!qfits_isnan(lineIm -> data[col + row*lineIm -> lx]) &&
02177                 (lineIm -> data[col + row*lineIm -> lx] > 0.))
02178             {
02179                 spec[row] = lineIm -> data[col + row*lineIm -> lx] ;
02180             }
02181             else
02182             {
02183                 spec[row] = 0. ;
02184             }
02185         }
02186         /* convolve the spectrum by Gaussian */
02187         filter_spec = function1d_filter_lowpass(spec, lineIm->ly, LOW_PASS_GAUSSIAN, magFactor) ; 
02188 
02189         /* now call the cross correlation routine */
02190         result = xcorrel( filter_spec, lineIm->ly, emline, lineIm->ly,
02191                           lineIm->ly/2, &delta, &maxlag, &xcorr_max) ;
02192 
02193         if ( xcorr_max <= 0. )
02194         {
02195             cpl_msg_warning("cross-corr","no positive cross correlation sum , col %d set to ZERO \n", col) ;
02196             for ( row = 0 ; row < lineIm -> ly ; row++ )
02197             {
02198                 wavemap -> data[col + row*lineIm -> lx] = ZERO ;
02199             }
02200             function1d_del(filter_spec) ;
02201             cpl_free(result) ;
02202             continue ;
02203         }
02204 
02205        /* in this section, we fit the correlation function with a gauss, and find its peak, th
02206            us getting subpixel-accuracy */
02207 
02208         i = maxlag; j = i+1;
02209         while (result[j] < result[i])
02210         {
02211             i++; j++;
02212         }
02213         i = maxlag; k = i-1;
02214         while (result[k] < result[i])
02215         {
02216             i--; k--;
02217         }
02218         width = j-k+1;
02219         /* allocate memory for the spectral vector */
02220         if ( NULL == (peak = newVector (width)) )
02221         {
02222             cpl_msg_error ("createShiftedWavemap3:"," cannot allocate new Vector \n") ;
02223             function1d_del(filter_spec) ;
02224             cpl_free(result) ;
02225             return NullImage ;
02226         }
02227 
02228 
02229         /* allocate memory */
02230         xdat = (float *) cpl_calloc( peak -> n_elements, sizeof (float) ) ;
02231         wdat = (float *) cpl_calloc( peak -> n_elements, sizeof (float) ) ;
02232         mpar = (int *)   cpl_calloc( MAXPAR, sizeof (int) ) ;
02233 
02234         /* determine the values of the spectral vector given as input */
02235         /* go through the chosen column */
02236 
02237         for ( i = 0 ; i < width ; i++ )
02238         {
02239             peak -> data[i] = result[k+i]/xcorr_max * 100. ;
02240             xdat[i] = i;
02241             wdat[i] = 1.0;
02242         }
02243 
02244         /* set initial values for the fitting routine */
02245         xdim     = XDIM;
02246         ndat     = peak -> n_elements ;
02247         numpar   = MAXPAR ;
02248         tol      = TOL ;
02249         lab      = LAB ;
02250         its      = ITS ;
02251         par[1] = width/2.0 ;
02252         par[2] = (float) (maxlag - k) ;
02253         par[3] = (peak -> data[0] + peak -> data[peak->n_elements - 1]) / 2.0 ;
02254         par[0]  = result[maxlag]/xcorr_max * 100. - (par[3]) ;
02255 
02256         for ( i = 0 ; i < MAXPAR ; i++ )
02257         {
02258             derv_par[i] = 0.0 ;
02259             mpar[i] = 1 ;
02260         }
02261 
02262         /* finally, do the least square fit using a gaussian */
02263         if ( 0 > ( iters = lsqfit_c( xdat, &xdim, peak -> data, wdat, &ndat, par, derv_par, mpar,
02264                                      &numpar, &tol, &its, &lab )) )
02265         {
02266             cpl_msg_warning ("linefit:"," lsqfit_c: least squares fit failed in col: %d, error no.: %d\n", col, iters) ;
02267             destroyVector ( peak ) ;
02268             cpl_free ( xdat ) ;
02269             cpl_free ( wdat ) ;
02270             cpl_free ( mpar ) ;
02271             function1d_del(filter_spec) ;
02272             cpl_free(result) ;
02273             continue ;
02274         }
02275 
02276         destroyVector ( peak ) ;
02277         cpl_free (xdat) ;
02278         cpl_free (wdat) ;
02279         cpl_free (mpar) ;
02280         function1d_del(filter_spec) ;
02281         cpl_free(result) ;
02282 
02283         wavelag[col] =((float)lineIm->ly/2 - (float)k - par[2]) ;
02284 
02285     }
02286     
02287     if (FLT_MAX == (wavelag_mean = clean_mean(wavelag, lineIm->lx, 10., 10.)) )
02288     {
02289         cpl_msg_error("createShiftedWavemap3:"," could not determine a mean offset\n") ;
02290         return NullImage ;
02291     }
02292 
02293     if ( fabs(wavelag_mean) > (float)lineIm->ly/20. )
02294     {
02295         cpl_msg_error("createShiftedWavemap3:"," wavelag too big \n") ;
02296         return NullImage ;
02297     }
02298 
02299     
02300 
02301     for ( col = 0 ; col < lineIm->lx ; col++ )
02302     {
02303         /*-----------------------------------------------------------------------------
02304          * determine new zero order coefficient centreval, of which the formula is
02305          * determined by setting equal a polynomial shifted by wavelag with the same
02306          * higher order coefficients and set the new zero order coefficient to
02307          * get both sides of the equation approximately equal.
02308          */
02309         a_initial = coeffs[0][col] ;
02310     centreval = a_initial ;
02311     for ( i = 1 ; i < n_fitcoeffs ; i++ )
02312         {
02313             if ( i%2 == 0 )
02314             {
02315                 sign = -1 ;
02316             }
02317             else
02318             {
02319                 sign = 1 ;
02320             }
02321             centreval += (float)sign * coeffs[i][col]*pow(wavelag_mean, i) ;
02322         }
02323 
02324     
02325         /* prepare to write out wavelength as pixel values */
02326         for ( row = 0 ; row < wavemap -> ly ; row++ )
02327         {
02328             centrepix = (float)row - ((float)wavemap->ly - 1.)/2. ;
02329             pixvalue = 0. ;
02330             for ( i = 1 ; i < n_fitcoeffs ; i++ )
02331             {
02332                 pixvalue += coeffs[i][col]*pow(centrepix, i) ;
02333             }
02334             wavemap -> data[col+row*wavemap->lx] =
02335                 centreval + pixvalue ;
02336         }
02337     }
02338     return wavemap ;
02339 }
02340 
02341 /*----------------------------------------------------------------------------
02342    Function     :       checkLinePositions()
02343    In           :       lineIm:       new shifted emission line frame
02344                         coeffs:       calculated polynomial coefficients,
02345                                       output of coeffsCrossSlitFit()
02346                         n_fitcoeffs:  number of polynomial coefficients, order + 1
02347                         par:          fit parameters
02348    Out          :       clean averaged position error (shift) of the brightest lines in the emission
02349                         line image with respect to the expected positions from the wavelength
02350                         calibration (error in the wavelength regime (microns)).
02351    Job          :       This routine determines the clean averaged error (shift) of the brightest
02352                         lines in the emission line image with respect to the expected positions
02353                         from the wavelength calibration. The error is given in the wavelength
02354                         regime (microns).
02355                         It should give the user an impression of the quality 
02356                         of the wavelength calibration.  
02357  ---------------------------------------------------------------------------*/
02358 
02359 float checkLinePositions ( OneImage     * lineIm,
02360                            float       ** coeffs,
02361                            int            n_fitcoeffs,
02362                            FitParams   ** par )
02363 {
02364     float wave_shift ;
02365     float amp[100] ;
02366     float sort_amp[100] ;
02367     float offset ;
02368     float shift ;
02369     float shift_col[lineIm->lx] ;
02370     float position, lambda=0, wave ;
02371     int i, j, k, l, m, n ;
02372     int col, firstj ;
02373     int foundit[par[0]->n_params] ;
02374     int n_lines ;
02375     int lin, found=0 ;
02376 
02377     if ( lineIm == NullImage )
02378     {
02379         cpl_msg_error ("checkLinePositions:"," no input image given!\n") ;
02380         return FLAG ;
02381     }
02382     if ( coeffs == NULL )
02383     {
02384         cpl_msg_error ("checkLinePositions:"," no coefficient matrix given!\n") ;
02385         return FLAG ;
02386     }
02387     if ( par == NULL )
02388     {
02389         cpl_msg_error ("checkLinePositions:"," no fit parameters given!\n") ;
02390         return FLAG ;
02391     }
02392     if ( n_fitcoeffs < 2 )
02393     {
02394         cpl_msg_error ("checkLinePositions:"," wrong number of polynomial coefficients given!\n") ;
02395         return FLAG ;
02396     }
02397 
02398     offset = (float) (lineIm->ly -1.) / 2. ;
02399     n_lines = par[0]->n_params/lineIm->lx ;
02400     
02401     /*search for the brightest 5 lines in each column and compute the wavelength difference*/
02402     for ( col = 0 ; col < lineIm->lx ; col++ )
02403     {
02404         n = 0 ;
02405         for ( i = 0 ; i < par[0]->n_params ; i++ )
02406         {
02407             if (par[i]->column == col && par[i]->fit_par[2] != 0. && 
02408                 par[i]->fit_par[1] > 1. && par[i]->fit_par[1] < 7. )
02409             {
02410                 foundit[n] = i ;
02411                 amp[n] = par[i]->fit_par[0] ;
02412                 sort_amp[n] = amp[n] ;
02413                 n++ ;
02414             }
02415         }
02416         pixel_qsort(sort_amp, n) ;
02417 
02418         if ( n > 5 )
02419         {
02420             firstj = n - 5 ;
02421         }
02422         else
02423         {
02424             firstj = 0 ;
02425         }
02426         l = 0 ; 
02427         shift = 0 ;
02428         for ( j = firstj ; j < n ; j++ )
02429         {
02430             for ( m = 0 ; m < n ; m++ )
02431             {
02432                 if ( sort_amp[j] == amp[m] )
02433                 {
02434                     position = par[foundit[m]]->fit_par[2] ; 
02435                     lambda   = par[foundit[m]]->wavelength ;
02436                     wave = 0 ;
02437                     for ( k = 0 ; k < n_fitcoeffs ; k++ ) 
02438                     {
02439                         wave += coeffs[k][col]*pow(position-offset, k) ;
02440                     }
02441                     shift += lambda - wave ;
02442                     l++ ;
02443                 }
02444             }
02445         }
02446         if ( l == 0 ) continue ;
02447         shift_col[col] = shift/(float)l ; 
02448     }
02449     wave_shift = clean_mean(shift_col, lineIm->lx, 10., 10.) ;
02450     cpl_msg_info("checkLinePositions:", "overall positioning error in microns: %g", wave_shift) ;
02451 
02452 
02453     /* determine positioning error for each found line */
02454     for ( lin = 0 ; lin < n_lines ; lin++ )
02455     {
02456         for ( col = 0 ; col < lineIm->lx ; col++ )
02457         {
02458             shift_col[col] = 0. ;
02459             found = -1 ;
02460             for ( i = 0 ; i < par[0]->n_params ; i++ )
02461             {
02462                 if (par[i]->column == col && par[i]->fit_par[2] != 0. && 
02463                     par[i]->fit_par[1] > 1. && par[i]->fit_par[1] < 7. && 
02464                     par[i]->line == lin )
02465                 {
02466                     found = i ;
02467                 }
02468             }
02469             if (found == -1) break ;
02470 
02471             position = par[found]->fit_par[2] ; 
02472             lambda  = par[found]->wavelength ;
02473             wave = 0 ;
02474             for ( k = 0 ; k < n_fitcoeffs ; k++ ) 
02475             {
02476                 wave += coeffs[k][col]*pow(position-offset, k) ;
02477             }
02478             shift_col[col] = lambda - wave ;
02479         }
02480         if (found != -1 )
02481         {
02482             cpl_msg_info("checkLinePositions:", "shift in microns: %g at wavelength: %f\n", clean_mean(shift_col, lineIm->lx, 10., 10.), lambda) ;
02483         }
02484     }
02485 
02486     return wave_shift ;
02487 }
02488 
02489 
02490 /*----------------------------------------------------------------------------
02491    Function     :       checkCorrelatedLinePositions()
02492    In           :       lineIm:       new shifted emission line frame
02493                         coeffs:       calculated polynomial coefficients,
02494                                       output of coeffsCrossSlitFit()
02495                         n_fitcoeffs:  number of polynomial coefficients, order + 1
02496                         wavelength:   wavelength list from line list
02497                         intensity:    corresponding line intensity from line list
02498                         n_lines:      number of lines in list
02499                         fwhm:         guess value for full width of half maximum of gaussian
02500                         width:        half width of the box where the line must sit
02501                         min_amplitude: minimum line amplitude with respect to the background
02502                                        to do the fit
02503                         par:          fit parameters
02504    Out          :       clean averaged position error (shift) of the brightest lines in the emission
02505                         line image with respect to the expected positions from the wavelength
02506                         calibration (error in the wavelength regime (microns)).
02507    Job          :       This routine determines the clean averaged error (shift) of the brightest
02508                         lines in the emission line image with respect to the expected positions
02509                         from the wavelength calibration. The error is given in the wavelength
02510                         regime (microns).
02511                         It should give the user an impression of the quality 
02512                         of the wavelength calibration.  
02513  ---------------------------------------------------------------------------*/
02514 
02515 float checkCorrelatedLinePositions ( OneImage     * lineIm,
02516                                      float       ** coeffs,
02517                                      int            n_fitcoeffs,
02518                                      float        * wavelength,
02519                                      float        * intensity,
02520                                      int            n_lines,
02521                                      float          fwhm, 
02522                                      float          width,
02523                                      float          min_amplitude,
02524                                      float          dispersion,
02525                                      FitParams   ** par )
02526 {
02527     float wave_shift=0 ;
02528     float offset ;
02529     float shift ;
02530     float shift_col[lineIm->lx] ;
02531     float position, lambda, wave ;
02532     int i, j, k/*, l, m*/, n, c, z ;
02533     int col/*, firstj*/ ;
02534     int foundit[par[0]->n_params] ;
02535     int /*lin,*/ found ;
02536     int line ;
02537     int result ;
02538     double a[n_fitcoeffs] ;
02539     float cenpos ;
02540     float angst ;
02541     float wave_cor[n_lines] ;
02542     float a_initial ;
02543     double zroot[2*(n_fitcoeffs - 1)] ;
02544     gsl_poly_complex_workspace * w ;
02545     Vector * vline;
02546     int    * mpar;
02547     float  * xdat, * wdat;
02548 
02549     if ( lineIm == NullImage )
02550     {
02551         cpl_msg_error ("checkCorrelatedLinePositions:"," no input image given!\n") ;
02552         return FLAG ;
02553     }
02554     if ( coeffs == NULL )
02555     {
02556         cpl_msg_error ("checkCorrelatedLinePositions:"," no coefficient matrix given!\n") ;
02557         return FLAG ;
02558     }
02559     if ( par == NULL )
02560     {
02561         cpl_msg_error ("checkCorrelatedLinePositions:"," no fit parameters given!\n") ;
02562         return FLAG ;
02563     }
02564     if ( n_fitcoeffs < 2 )
02565     {
02566         cpl_msg_error ("checkCorrelatedLinePositions:"," wrong number of polynomial coefficients given!\n") ;
02567         return FLAG ;
02568     }
02569     if ( wavelength == NULL || intensity == NULL )
02570     {
02571         cpl_msg_error ("checkCorrelatedLinePositions:"," no line list given!\n") ;
02572         return FLAG ;
02573     }
02574     if ( fwhm <= 0 )
02575     {
02576         cpl_msg_error ("checkCorrelatedLinePositions:"," wrong guess fwhm given!\n") ;
02577         return FLAG ;
02578     }
02579     if ( width <= 0 )
02580     {
02581         cpl_msg_error ("checkCorrelatedLinePositions:"," wrong half width given!\n") ;
02582         return FLAG ;
02583     } 
02584     if ( min_amplitude <= 0 )
02585     {
02586         cpl_msg_error ("checkCorrelatedLinePositions:"," wrong guess amplitude given!\n") ;
02587         return FLAG ;
02588     } 
02589 
02590     /* allocate memory for the spectral vector */
02591     if ( NULL == (vline = newVector (2*width + 1)) )
02592     {
02593         cpl_msg_error ("linefit:"," cannot allocate new Vector \n") ;
02594         return -14 ;
02595     }
02596     /* allocate memory */
02597     xdat = (float *) cpl_calloc( vline -> n_elements, sizeof (float) ) ;
02598     wdat = (float *) cpl_calloc( vline -> n_elements, sizeof (float) ) ;
02599     mpar = (int *)   cpl_calloc( MAXPAR, sizeof (int) ) ;
02600 
02601 
02602     /* determine the approximate line positions using the line list and the coefficients */
02603     /* find out if Angstroem or microns are used */
02604     if ( wavelength[0] > 10000. )
02605     {
02606         /* Angstroem */
02607         angst = 10000. ;
02608     }
02609     else if ( wavelength[0] > 1000. && wavelength[0] < 10000. )
02610     {
02611         /* nanometers */
02612         angst = 1000. ;
02613     }
02614     else
02615     {
02616         /* microns */
02617         angst = 1. ;
02618     }
02619     offset = ((float) lineIm->ly -1.) / 2. ;
02620 
02621     k = 0 ;
02622     for ( col = 10 ; col < 25 ; col++ )
02623     {
02624         /* determine the coefficients by using the given bcoefs */
02625         for ( i = 0 ; i < n_fitcoeffs ; i++ )
02626         {
02627             /* initialize coefficients and solution */
02628             if (i < n_fitcoeffs-1)
02629             {
02630                 zroot[2*i] = 0. ;
02631                 zroot[2*i+1] = 0. ;
02632             }
02633             a[i] = coeffs[i][col] ;
02634         }
02635         a_initial = a[0] ;
02636 
02637         /* go through the lines */
02638         for ( line = 0 ; line < n_lines ; line++ )
02639         {
02640             /* go from Angstroem to micron */
02641             wave_cor[line] = wavelength[line]/angst ;
02642             if (line > 0 && line < n_lines-1)
02643             {
02644                 if (fabs((wave_cor[line] - wave_cor[line-1]) / dispersion ) < 2*width ||
02645                     fabs((wave_cor[line] - wave_cor[line+1]) / dispersion ) < 2*width )
02646                 {
02647                     continue ;
02648                 }
02649             }
02650 
02651             a[0] = a_initial - wave_cor[line] ;
02652 
02653             if (NULL == (w = gsl_poly_complex_workspace_alloc(n_fitcoeffs)))
02654             {
02655                 cpl_msg_error("checkCorrelatedLinePositions:"," could not allocate complex workspace!") ;
02656                 return FLAG ;
02657             }
02658             if (-1 == gsl_poly_complex_solve(a, n_fitcoeffs, w, zroot))
02659             {
02660                 cpl_msg_error("checkCorrelatedLinePositions:"," gsl_poly_complex_solve did not work!") ;
02661                 return FLAG ;
02662             }
02663             gsl_poly_complex_workspace_free(w) ;
02664 
02665             j = 0 ;
02666             found = -1 ;
02667             for ( i = 0 ; i < n_fitcoeffs - 1 ; i++ )
02668             {
02669                 /* test for appropriate solution */
02670                 if( (zroot[2*i] > (-1.)*(float) lineIm->ly/2. &&
02671                      zroot[2*i] < (float)lineIm -> ly/2.) && zroot[2*i+1] == 0. )
02672                 {
02673                     found = 2*i ;
02674                     j ++ ;
02675                 }
02676                 else
02677                 {
02678                     continue ;
02679                 }
02680             }
02681 
02682             if ( j == 0 )
02683             {
02684                 cpl_msg_warning("checkCorrelatedLinePositions:"," no offset solution found for line %d in column %d\n", line, col) ;
02685                 continue ;
02686             }
02687             else if ( j == 1 )
02688             {
02689                 cenpos = zroot[found] + (float)lineIm->ly / 2. ; ;
02690             }
02691             else
02692             {
02693                 cpl_msg_warning("checkCorrelatedLinePositions:"," two or more offset solutions found for \
02694                        line %d in column %d\n", line, col) ;
02695                 continue ;
02696             }
02697 
02698             if ( cenpos <= 0 )
02699             {
02700                 continue ;
02701             }
02702 
02703             /* --------------------------------------------------------------
02704              * fit the single lines using linefit and store the parameters in
02705              * an array of the FitParams data structure allParams[].
02706              */
02707             if ( (result = linefit ( lineIm, par[k], fwhm, line, col,
02708                                      width, cenpos, min_amplitude,vline,mpar,xdat,wdat ) ) < 0 )
02709             {
02710                 cpl_msg_debug ("fitLines:"," linefit failed, error no.: %d, column: %d,\
02711                             row: %f, line: %d\n", result, col, cenpos, line) ;
02712                 continue ;
02713             }
02714             if ( (par[k] -> fit_par[0] <= 0.) || (par[k] -> fit_par[1] <= 0.)
02715                   || (par[k] -> fit_par[2] <= 0.) )
02716             {
02717                 cpl_msg_warning ("linefit","negative fit parameters in column: %d, line: %d\n", col, line) ;
02718                 continue ;
02719             }
02720             par[k] -> wavelength = wavelength[line] ;
02721             k++ ;
02722         }
02723 
02724     }
02725     
02726     /* free memory */
02727     destroyVector(vline);
02728     cpl_free(xdat);
02729     cpl_free(wdat);
02730     cpl_free(mpar);
02731 
02732     
02733     c = 0 ;
02734     for ( col = 10 ; col < 25 ; col++ )
02735     {
02736         n = 0 ;
02737         for ( i = 0 ; i < k ; i++ )
02738         {
02739             if (par[i]->column == col && par[i]->fit_par[2] != 0. && 
02740                 par[i]->fit_par[1] > 1. && par[i]->fit_par[1] < 7. )
02741             {
02742                 foundit[n] = i ;
02743                 n++ ;
02744             }
02745         }
02746         if ( n == 0 ) continue ;
02747 
02748         shift = 0 ;
02749         z = 0 ;
02750         for ( j = 0 ; j < n ; j++ )
02751         {
02752             position = par[foundit[j]]->fit_par[2] ; 
02753             lambda   = par[foundit[j]]->wavelength ;
02754             line     = par[foundit[j]]->line ;
02755             if (line > 0 && line < n_lines-1)
02756             {
02757                 if (fabs((wave_cor[line] - wave_cor[line-1]) / dispersion ) < 2*width ||
02758                     fabs((wave_cor[line] - wave_cor[line+1]) / dispersion ) < 2*width )
02759                 {
02760                     continue ;
02761                 }
02762             }
02763             wave = 0 ;
02764             for ( i = 0 ; i < n_fitcoeffs ; i++ ) 
02765             {
02766                 wave += coeffs[i][col]*pow(position-offset, i) ;
02767             }
02768             shift += lambda - wave ;
02769             z++ ;
02770         }
02771         shift_col[c] = shift/(float)z ; 
02772         c++ ;
02773     }
02774     if ( c > 0 )
02775     {
02776         wave_shift = clean_mean(shift_col, c, 10., 10.) ;
02777         cpl_msg_info("checkCorrelatedLinePositions:", 
02778              "overall positioning error in microns: %g", wave_shift) ;
02779     }
02780 
02781     /* determine positioning error for each found line */
02782     for ( line = 0 ; line < n_lines ; line++ )
02783     {
02784         if (line > 0 && line < n_lines-1)
02785         {
02786             if (fabs((wave_cor[line] - wave_cor[line-1]) / dispersion ) < 2*width ||
02787                 fabs((wave_cor[line] - wave_cor[line+1]) / dispersion ) < 2*width )
02788             {
02789                 continue ;
02790             }
02791         }
02792 
02793         c = 0 ;
02794         for ( col = 10 ; col < 25 ; col++ )
02795         {
02796             shift_col[c] = 0. ;
02797             found = -1 ;
02798             for ( i = 0 ; i < k ; i++ )
02799             {
02800                 if (par[i]->column == col && par[i]->fit_par[2] != 0. && 
02801                     par[i]->fit_par[1] > 1. && par[i]->fit_par[1] < 7. && 
02802                     par[i]->line == line )
02803                 {
02804                     found = i ;
02805                 }
02806             }
02807             if (found == -1) break ;
02808 
02809             position = par[found]->fit_par[2] ; 
02810             lambda  = par[found]->wavelength ;
02811             wave = 0 ;
02812             for ( i = 0 ; i < n_fitcoeffs ; i++ ) 
02813             {
02814                 wave += coeffs[i][col]*pow(position-offset, i) ;
02815             }
02816             shift_col[c] = lambda - wave ;
02817             c++ ;
02818         }
02819         if (found != -1 && c > 0 )
02820         {
02821             cpl_msg_info("checkCorrelatedLinePositions:", 
02822          "shift in microns: %g at wavelength: %f\n", clean_mean(shift_col, c, 20., 20.), lambda) ;
02823         }
02824     }
02825     return wave_shift ;
02826 }
02827 
02828 
02829  int   spred_coeffsCrossSlitFit ( int      n_columns,
02830                             float ** acoefs,
02831                             float ** dacoefs,
02832                             Bcoeffs* bco,
02833                             float    sigma_factor,
02834                             float    dispersion,
02835                             float    pixel_dist,
02836                             float  * chisq,
02837                           float ** slit_pos) ;
02838 
02839  int   spred_coeffsCrossSlitFit ( int      n_columns,
02840                             float ** acoefs,
02841                             float ** dacoefs,
02842                             Bcoeffs* bco,
02843                             float    sigma_factor,
02844                             float    dispersion,
02845                             float    pixel_dist,
02846                             float  * chisq,
02847                           float ** slit_pos )
02848  {
02849      float col_index, sub_col_index[n_columns] ;
02850      float sub_acoefs[n_columns], sub_dacoefs[n_columns] ;
02851      float wcoefs[bco->n_bcoeffs] ;
02852      float ** ucoefs, **vcoefs, **covar ;
02853      float * acoefsclean ;
02854      double sum, sumq, mean ;
02855      double sigma ;
02856      double cliphi, cliplo ;
02857      float offset ;
02858      float threshold ;
02859      int edge[bco->n_slitlets] ;
02860      int ed1, ed2 ;
02861      int i, n, num, ndata ;
02862      int nc, ns ;
02863      int index ;
02864      int sl_index;
02865      int last_i=PIXEL;
02866  
02867      if ( n_columns < 1 )
02868      {
02869          cpl_msg_error("coeffsCrossSlitFit:"," wrong number of image columns given\n") ;
02870          return -1 ;
02871      }
02872      if ( acoefs == NULL || dacoefs == NULL )
02873      {
02874          cpl_msg_error("coeffsCrossSlitFit:"," acoeffs or errors of coefficients are not given\n") ;
02875          return -1 ;
02876      }
02877      if ( bco == NULL )
02878      {
02879          cpl_msg_error("coeffsCrossSlitFit:"," bcoeffs are not allocated\n") ;
02880          return -1 ;
02881      }
02882      if ( sigma_factor <= 0. )
02883      {
02884          cpl_msg_error("coeffsCrossSlitFit:"," impossible sigma_factor given!\n") ;
02885          return -1 ;
02886      }
02887      if ( dispersion == 0. )
02888      {
02889          cpl_msg_error("coeffsCrossSlitFit:"," impossible dispersion given!\n") ;
02890          return -1 ;
02891      }
02892 
02893      /*---------------------------------------------------------------------------------
02894       * search for the slitlet edges by comparing the a0 coefficients along the columns
02895       * if a bigger deviation occurrs it is assumed that there is an edge.
02896       */
02897      n = 0 ;
02898      threshold = pixel_dist * fabs(dispersion) ;
02899      slit_pos[0][0]=0 ;
02900      sl_index  = 0;
02901      for ( i = PIXEL ; i  < n_columns - PIXEL ; )
02902      {
02903          if ( !qfits_isnan(acoefs[0][i+1]) && acoefs[0][i+1] != 0. && acoefs[0][i] != 0.
02904                    && dacoefs[0][i+1] != 0.)
02905          {
02906              if ( qfits_isnan(acoefs[0][i]) || acoefs[0][i] == 0. )
02907              {
02908                  if (fabs(acoefs[0][i+1] - acoefs[0][i-1]) >= threshold )
02909              {
02910              /* printf("case a pos1 %d pos2 %d \n",i,i+1); */
02911                      edge[n] = i+1 ;
02912              slit_pos[sl_index][1] = i ;
02913              slit_pos[sl_index+1][0] = i + 1 ;
02914              sl_index++;
02915                      n++ ;
02916              last_i = i;
02917                      i += PIXEL ;
02918                  }
02919              }
02920              else
02921              {
02922                  if (fabs(acoefs[0][i+1] - acoefs[0][i]) >= threshold )
02923              {
02924              /* printf("case b pos1 %d pos2 %d \n",i,i+1); */
02925                      edge[n] = i+1 ;
02926              slit_pos[sl_index][1] = i ;
02927              slit_pos[sl_index+1][0] = i + 1 ;
02928              sl_index++;
02929                      n++ ;
02930              last_i = i;
02931                      i += PIXEL ;
02932                  }
02933              }
02934 
02935 
02936          /* sometimes a slitlet may be lost due to divergences in acoeffs[0]
02937                we try to recover it */
02938              if( ( (i-last_i) > 63 ) && 
02939                 ( qfits_isnan(fabs(acoefs[0][i+1] - acoefs[0][i])) ||
02940                   qfits_isnan(fabs(acoefs[0][i+1] - acoefs[0][i-1])) ) ) 
02941          {
02942                      edge[n] = i+1 ;
02943                      slit_pos[sl_index][1] = i ;
02944                      slit_pos[sl_index+1][0] = i + 1 ;
02945                      sl_index++;
02946                      n++ ;
02947 
02948                      last_i = i;
02949                      cpl_msg_warning("wavecal:","2 recovered slitlet edge i=%d",i);
02950                      i += PIXEL ;
02951 
02952          }
02953          }
02954          i++ ;
02955      }
02956      slit_pos[sl_index][1]  = 2047;
02957      /* printf("2 Found n slitlest: %d check %d\n", n,bco->n_slitlets - 1); */
02958      if ( n != bco->n_slitlets - 1 )
02959      {
02960          cpl_msg_error("coeffsCrossSlitFit:"," could not find the right number of slitlets, found: %d\n",n+1) ;
02961          return -1 ;
02962      }
02963 
02964      /* go through the coefficents indices */
02965      for ( index = 0 ; index < bco->n_acoeffs ; index++ )
02966      {
02967          /* go through the single slitlets */
02968          for ( ns = 0 ; ns < bco->n_slitlets ; ns++ )
02969          {
02970              /* determine the slitlet edges */
02971              if ( ns == 0 )
02972              {
02973                  ed1 = 0 ;
02974                  ed2 = edge[0] ;
02975              }
02976              else if ( ns == bco->n_slitlets - 1 )
02977              {
02978                  ed1 = edge[bco->n_slitlets - 2] ;
02979                  ed2 = n_columns ;
02980              }
02981              else
02982              {
02983                  ed1 = edge[ns-1] ;
02984                  ed2 = edge[ns] ;
02985              }
02986 
02987              nc = 0 ;
02988              for ( i = ed1 ; i < ed2 ; i++ )
02989              {
02990                  if ( qfits_isnan(acoefs[index][i]) || acoefs[index][i] == 0. || dacoefs[index][i] == 0. )
02991                  {
02992                      continue ;
02993                  }
02994                  else
02995                  {
02996                      nc++ ;
02997                  }
02998              }
02999              if ( NULL == (acoefsclean = (float*) cpl_calloc(nc , sizeof(float))) )
03000              {
03001                  cpl_msg_error("coeffsCrossSlitFit:"," could not allocate memory for acoefsclean!\n") ;
03002                  return -1 ;
03003              }
03004              nc = 0 ;
03005              for ( i = ed1 ; i  < ed2 ; i++ )
03006             {
03007                 if ( qfits_isnan(acoefs[index][i]) || acoefs[index][i] == 0. || dacoefs[index][i] == 0. )
03008                 {
03009                     continue ;
03010                 }
03011                 else
03012                 {
03013                     acoefsclean[nc] = acoefs[index][i] ;
03014                     nc++ ;
03015                 }
03016             }
03017 
03018             /* ----------------------------------------------------------
03019              * determine the clean mean and sigma value of the coefficients,
03020              * that means reject 10 % of the extreme low and high values
03021              */
03022             pixel_qsort(acoefsclean, nc) ;
03023 
03024             sum   = 0. ;
03025             sumq  = 0. ;
03026             mean  = 0. ;
03027             sigma = 0. ;
03028             n     = 0 ;
03029             for ( i = (int)((float)nc*LOW_REJECT) ; i < (int)((float)nc*HIGH_REJECT) ; i++ )
03030             {
03031                 sum  += (double)acoefsclean[i] ;
03032                 sumq += ((double)acoefsclean[i] * (double)acoefsclean[i]) ;
03033                 n ++ ;
03034             }
03035             mean          = sum/(double)n ;
03036             sigma         = sqrt( sumq/(double)n - (mean * mean) ) ;
03037             cliphi        = mean + sigma * (double)sigma_factor ;
03038             cliplo        = mean - sigma * (double)sigma_factor ;
03039             /* fit only the reasonnable values */
03040             num = 0 ;
03041             col_index = 0 ;
03042             for ( i = ed1 ; i < ed2 ; i++ )
03043             {
03044                 /* take only the reasonnable coefficients */
03045                 if ( !qfits_isnan(acoefs[index][i]) && (acoefs[index][i] <= cliphi) && (acoefs[index][i] >= cliplo) &&
03046                      (dacoefs[index][i] != 0. ) && (acoefs[index][i] != 0.) )
03047                 {
03048                     sub_acoefs[num]    = acoefs[index][i] ;
03049                     sub_dacoefs[num]   = dacoefs[index][i] ;
03050                     sub_col_index[num] = col_index ;
03051                     num ++ ;
03052                 }
03053                 col_index++ ;
03054             }
03055             ndata = num ;
03056             offset = (float)(col_index-1) / 2. ;
03057 
03058             if ( ndata < bco->n_bcoeffs )
03059             {
03060                 cpl_msg_error("coefsCrossSlitFit:"," not enough data found in slitlet %d to determine the fit coefficients.\n", ns) ;
03061                 cpl_free(acoefsclean) ;
03062                 return -1 ;
03063             }
03064 
03065             /* allocate coefficient matrices */
03066             ucoefs = matrix(1, ndata, 1, bco->n_bcoeffs) ;
03067             vcoefs = matrix(1, ndata, 1, bco->n_bcoeffs) ;
03068             covar  = matrix(1, bco->n_bcoeffs, 1, bco->n_bcoeffs) ;
03069 
03070             /* scale the x-values for the fit */
03071             for ( i = 0 ; i < ndata ; i++ )
03072             {
03073                 sub_col_index[i] = (sub_col_index[i] - offset) / offset ;
03074             }
03075 
03076             /* finally, do the singular value decomposition fit */
03077             svd_fitting ( sub_col_index-1, sub_acoefs-1, sub_dacoefs-1, ndata, bco[ns].b[index]-1,
03078                      bco->n_bcoeffs, ucoefs, vcoefs, wcoefs-1, covar, &chisq[ns], fpol ) ;
03079 
03080             /* scale the found coefficients */
03081             for ( i = 0 ; i < bco->n_bcoeffs ; i ++ )
03082             {
03083                 bco[ns].b[index][i] /= pow( offset, i ) ;
03084             }
03085 
03086             /* free memory */
03087             cpl_free (acoefsclean) ;
03088             free_matrix( ucoefs, 1/*, ndata*/, 1/*, bco->n_bcoeffs */) ;
03089             free_matrix( vcoefs, 1/*, ndata*/, 1/*, bco->n_bcoeffs */) ;
03090             free_matrix( covar, 1/*, bco->n_bcoeffs*/, 1/*, bco->n_bcoeffs */) ;
03091 
03092             /* now calculate the smoothed acoefs for each column */
03093             col_index = 0 ;
03094             for ( i = ed1 ; i < ed2  ; i++ )
03095             {
03096                 acoefs[index][i] = 0. ;
03097                 for ( n = 0 ; n < bco->n_bcoeffs ; n++ )
03098                 {
03099                     acoefs[index][i] += bco[ns].b[index][n] * pow(col_index - offset, n) ;
03100                 }
03101                 col_index++ ;
03102             }
03103 
03104         }
03105     }
03106     return 0 ;
03107 }
03108 
03109 
03110 
03111 
03112 OneImage * spred_waveCal( OneImage   * image,
03113                     FitParams ** par ,
03114                     float     ** abuf,
03115                     int          n_slitlets,
03116                     int       ** row_clean,
03117                     float     ** wavelength_clean,
03118                     int        * n_found_lines,
03119                     float        dispersion,
03120                     int          halfWidth,
03121                     float        minAmplitude,
03122                     float        max_residual,
03123                     float        fwhm,
03124                     int          n_a_fitcoefs,
03125                     int          n_b_fitcoefs,
03126                     float        sigmaFactor,
03127                     float        pixel_dist,
03128                     float        pixel_tolerance,
03129                   float ** slit_pos)
03130 
03131 
03132 {
03133     int          i, j, k ;
03134     int          n_fit ;
03135     int          n_reject ;
03136     float     *  acoefs ;
03137     float     *  dacoefs ;
03138     float     ** dabuf ;
03139     float        chisq_poly ;
03140     float     *  chisq_cross ;
03141     int          zeroind ;
03142     int          crossInd ;
03143     Bcoeffs   *  bco ;
03144     OneImage  *  wavemap ;
03145 
03146     if (  NullImage == image )
03147     {
03148         cpl_msg_error("waveCal:"," no image given\n") ;
03149         return NullImage ;
03150     }
03151     if ( par == NULL )
03152     {
03153         cpl_msg_error("waveCal:"," no fit parameter data structure given\n") ;
03154         return NullImage ;
03155     }
03156     if ( abuf == NULL )
03157     {
03158         cpl_msg_error("waveCal:"," no buffer for fit coefficients given\n") ;
03159         return NullImage ;
03160     }
03161     if ( n_slitlets <= 0 )
03162     {
03163         cpl_msg_error("waveCal:"," impossible number of slitlets given\n") ;
03164         return NullImage ;
03165     }
03166     if ( row_clean == NULL )
03167     {
03168         cpl_msg_error("waveCal:"," no row_clean array given\n") ;
03169         return NullImage ;
03170     }
03171     if ( wavelength_clean == NULL )
03172     {
03173         cpl_msg_error("waveCal:"," no wavelength_clean array given\n") ;
03174         return NullImage ;
03175     }
03176 
03177     if ( dispersion == 0. )
03178     {
03179         cpl_msg_error("waveCal:"," impossible dispersion given\n") ;
03180         return NullImage ;
03181     }
03182 
03183     if ( halfWidth <= 0 || halfWidth > image->ly/2 )
03184     {
03185         cpl_msg_error("waveCal:"," impossible half width of the fitting box given\n") ;
03186         return NullImage ;
03187     }
03188     if ( minAmplitude < 1. )
03189      {
03190          cpl_msg_error("waveCal:"," impossible minimal amplitude\n") ;
03191          return NullImage ;
03192      }
03193 
03194      if ( max_residual <= 0. || max_residual > 1. )
03195      {
03196          cpl_msg_error("waveCal:"," impossible max_residual given\n") ;
03197          return NullImage ;
03198      }
03199      if ( fwhm <= 0. || fwhm > 10. )
03200      {
03201          cpl_msg_error("waveCal:"," impossible fwhm given\n") ;
03202          return NullImage ;
03203      }
03204 
03205      if ( n_a_fitcoefs <= 0 || n_a_fitcoefs > 9 )
03206      {
03207          cpl_msg_error("waveCal:"," unrealistic n_a_fitcoefs given\n") ;
03208          return NullImage ;
03209      }
03210 
03211      if ( n_b_fitcoefs <= 0 || n_b_fitcoefs > 9 )
03212      {
03213          cpl_msg_error("waveCal:"," unrealistic n_b_fitcoefs given\n") ;
03214          return NullImage ;
03215      }
03216      if ( sigmaFactor <= 0. )
03217      {
03218          cpl_msg_error("waveCal:"," impossible sigmaFactor given\n") ;
03219          return NullImage ;
03220      }
03221 
03222      /* initialize the variables */
03223      n_reject = 0 ;
03224      n_fit = 0 ;
03225 
03226      /* fit each found line by using a gaussian function and determine the exact position */
03227      if ( 0 > (n_fit = fitLines( image , par, fwhm, n_found_lines, row_clean, wavelength_clean,
03228                                  halfWidth, minAmplitude )) )
03229      {
03230          cpl_msg_error("waveCal:"," cannot fit the lines, error code of fitLines: %d\n", n_fit) ;
03231          return NullImage ;
03232      }
03233 
03234      /* first check for faked lines like bad pixels */
03235      if ( -1 == checkForFakeLines (par, dispersion, wavelength_clean, row_clean, n_found_lines,
03236                                    image->lx, pixel_tolerance) )
03237      {
03238          cpl_msg_error("waveCal:"," cannot fit the lines, error code of fitLines: %d\n", n_fit) ;
03239          return NullImage ;
03240      }
03241 
03242 
03243      /* allocate memory */
03244      if ( NULL == (acoefs = (float*) cpl_calloc (n_a_fitcoefs, sizeof(float))) ||
03245           NULL == (dacoefs = (float*) cpl_calloc (n_a_fitcoefs, sizeof(float))) ||
03246           NULL == (dabuf = (float**) cpl_calloc (n_a_fitcoefs, sizeof(float*))) ||
03247           NULL == (chisq_cross = (float*) cpl_calloc(n_slitlets, sizeof(float))) )
03248      {
03249          cpl_msg_error("waveCal:"," cannot allocate memory\n") ;
03250          return NullImage ;
03251      }
03252      for ( i = 0 ; i < n_a_fitcoefs ; i++ )
03253      {
03254          if (  NULL == (dabuf[i] = (float*) cpl_calloc(image -> lx, sizeof(float))) )
03255          {
03256              cpl_msg_error("wavelengthCalibration:"," cannot allocate memory\n") ;
03257              cpl_free ( acoefs ) ;
03258              cpl_free ( dacoefs ) ;
03259              cpl_free ( chisq_cross ) ;
03260              cpl_free(dabuf) ;
03261              return NullImage ;
03262          }
03263      }
03264 
03265      /* fit wavelengths to the corresponding found positions for each column */
03266      k = 0 ;
03267      for ( i = 0 ; i < image -> lx ; i++ )
03268      {
03269          zeroind = 0 ;
03270          if ( FLT_MAX == (chisq_poly = polyfit( par, i, n_found_lines[i], image -> ly, dispersion,
03271                                                 max_residual, acoefs, dacoefs, &n_reject, n_a_fitcoefs)) )
03272          {
03273              cpl_msg_warning ("waveCal"," error in polyfit in column: %d\n", i) ;
03274              for ( j = 0 ; j < n_a_fitcoefs ; j++ )
03275              {
03276                  acoefs[j] = ZERO ;
03277                  dacoefs[j] = ZERO ;
03278              }
03279          }
03280 
03281          for ( j = 0 ; j < n_a_fitcoefs ; j++ )
03282          {
03283 
03284              if ( acoefs[0] <= 0. || acoefs[1] ==0. ||
03285                   dacoefs[j] == 0. || qfits_isnan(acoefs[j]) )
03286              {
03287                  zeroind = 1 ;
03288              }
03289          }
03290          for ( j = 0 ; j < n_a_fitcoefs ; j++ )
03291          {
03292              if ( zeroind == 0 )
03293              {
03294                  abuf[j][i]  = acoefs[j] ;
03295                  dabuf[j][i] = dacoefs[j] ;
03296              }
03297              else
03298              {
03299                  abuf[j][i]  = ZERO ;
03300                  dabuf[j][i] = ZERO ;
03301              }
03302          }
03303      }
03304 
03305      /* allocate memory for the fitting coefficients */
03306      if ( NULL == ( bco = newBcoeffs( n_slitlets, n_a_fitcoefs, n_b_fitcoefs)) )
03307      {
03308          cpl_msg_error ("waveCal:"," cannot allocate memory for the bcoeffs\n") ;
03309          for ( i = 0 ; i < n_a_fitcoefs ; i++ )
03310          {
03311              cpl_free (dabuf[i]) ;
03312          }
03313          cpl_free (dabuf) ;
03314          cpl_free ( acoefs ) ;
03315          cpl_free ( dacoefs ) ;
03316          cpl_free ( chisq_cross ) ;
03317          return NullImage ;
03318      }
03319 
03320      /* fit each acoefs across the slitlets to smooth the result */
03321      if ( -1 == ( crossInd = spred_coeffsCrossSlitFit( image -> lx, abuf, dabuf,
03322                                                  bco, sigmaFactor, dispersion, pixel_dist, chisq_cross,slit_pos )) )
03323      {
03324          cpl_msg_error ("waveCal:"," cannot carry out the fitting of coefficients across the columns\n") ;
03325          for ( i = 0 ; i < n_a_fitcoefs ; i++ )
03326          {
03327              cpl_free (dabuf[i]) ;
03328          }
03329 
03330          cpl_free (dabuf) ;
03331          cpl_free ( acoefs ) ;
03332          cpl_free ( dacoefs ) ;
03333          destroyBcoeffs(bco) ;
03334          cpl_free ( chisq_cross ) ;
03335          return NullImage ;
03336      }
03337 
03338      if ( NullImage == (wavemap = waveMapSlit (abuf, n_a_fitcoefs, image->ly, image->lx)) )
03339      {
03340          cpl_msg_error ("waveCal:"," cannot carry out wavemap creation\n") ;
03341          for ( i = 0 ; i < n_a_fitcoefs ; i++ )
03342          {
03343              cpl_free (dabuf[i]) ;
03344          }
03345 
03346          cpl_free (dabuf) ;
03347          cpl_free ( acoefs ) ;
03348          cpl_free ( dacoefs ) ;
03349          destroyBcoeffs(bco) ;
03350          cpl_free ( chisq_cross ) ;
03351          return NullImage ;
03352      }
03353 
03354      /* free all allocated memory */
03355      for ( i = 0 ; i < n_a_fitcoefs ; i++ )
03356      {
03357          cpl_free (dabuf[i]) ;
03358      }
03359      cpl_free (dabuf) ;
03360      cpl_free ( acoefs ) ;
03361      cpl_free ( dacoefs ) ;
03362      destroyBcoeffs(bco) ;
03363      cpl_free ( chisq_cross ) ;
03364 
03365      return wavemap ;
03366  }
03367 
03368 
03369 /*___oOo___*/

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