boltzmann.c

00001 /*******************************************************************************
00002 * E.S.O. - VLT project
00003 *
00004 *
00005 *
00006 * who       when      what
00007 * --------  --------  ----------------------------------------------
00008 * schreib  27/02/01  created
00009 */
00010 
00011 /************************************************************************
00012 *   NAME
00013 *        boltzmann.c -
00014 *        routines to determine the absolute positions of the slitlets out
00015 *        of an emission line frame
00016 *
00017 *   SYNOPSIS
00018 *   #include "absolute.h"
00019 *
00020 *   1) float boltz ( float * xdat, float * parlist )
00021 *   2) void boltz_deriv( float * xdat, float * parlist, float * dervs )
00022 *   3) static int inv_mat (void)
00023 *   4) static void get_mat ( float * xdat,
00024 *                            int   * xdim,
00025 *                            float * ydat,
00026 *                            float * wdat,
00027 *                            int   * ndat,
00028 *                            float * fpar,
00029 *                            float * epar,
00030 *                            int   * npar )
00031 *   5) static int get_vec ( float * xdat,
00032 *                           int   * xdim,
00033 *                           float * ydat,
00034 *                           float * wdat,
00035 *                           int   * ndat,
00036 *                           float * fpar,
00037 *                           float * epar,
00038 *                           int   * npar )
00039 *   6) int lsqfit ( float * xdat,
00040 *                   int   * xdim,
00041 *                   float * ydat,
00042 *                   float * wdat,
00043 *                   int   * ndat,
00044 *                   float * fpar,
00045 *                   float * epar,
00046 *                   int   * mpar,
00047 *                   int   * npar,
00048 *                   float * tol ,
00049 *                   int   * its ,
00050 *                   float * lab  )
00051 *   7) int fitSlitsBoltz( OneImage   * lineImage, 
00052 *                         FitParams ** par,
00053 *                         float     ** slit_pos,
00054 *                         int          box_length,
00055 *                         float        y_box,
00056 *                         float        diff_tol )
00057 *   8) int fitSlitsBoltzSingleLine ( OneImage   * lineImage, 
00058 *                                    float     ** slit_pos,
00059 *                                    int          box_length,
00060 *                                    float        y_box,
00061 *                                int          low_pos,
00062 *                                    int          high_pos )
00063 *   9) int fitSlitsBoltzWithEstimate ( OneImage   * lineImage, 
00064 *                                      float     ** slit_pos,
00065 *                                      int          box_length,
00066 *                                      float        y_box,
00067 *                                      float        diff_tol,
00068 *                                      int          low_pos,
00069 *                          int          high_pos )
00070 *
00071 *   DESCRIPTION
00072 *   1) calculates the value of a Boltzmann function with parameters 
00073 *      parlist at the position xdat 
00074 *   2) calculates the partial derivatives for a Boltzmann function with
00075 *      parameters parlist at position xdat 
00076 *   3) calculates the inverse of matrix2. The algorithm used 
00077 *      is the Gauss-Jordan algorithm described in Stoer,
00078 *      Numerische Mathematik, 1. Teil.
00079 *   4) builds the matrix 
00080 *   5) calculates the correction vector. The matrix has been
00081 *      built by get_mat(), we only have to rescale it for the 
00082 *      current value of labda. The matrix is rescaled so that
00083 *      the diagonal gets the value 1 + labda.
00084 *      Next we calculate the inverse of the matrix and then
00085 *      the correction vector.
00086 *   6) this is a routine for making a least-squares fit of a
00087 *      function to a set of data points. The method used is
00088 *      described in: Marquardt, J.Soc.Ind.Appl.Math. 11. 431 (1963).
00089 *      This method is a mixture of the steepest descent method 
00090 *      and the Taylor method.
00091 *   7) fits the beginning and end position of the slitlets
00092 *      by using non-linear least square fitting of Boltzmann function
00093 *      fits a Boltzmann function to the slitlet edges exposed and indicated
00094 *      by the brightest emission lines. To achieve this, the fit
00095 *      parameters are used to find the brightest emission line
00096 *      and to get its position for each column.
00097 *      The least squares fit is done by using a box smaller than
00098 *      the size of two slitlets
00099 *   8) fits the beginning and end position of the slitlets
00100 *      by using non-linear least square fitting of a Boltzmann function
00101 *      fits a Boltzmann function to the slitlet edges exposed and indicated
00102 *      by the brightest emission lines. The slitlet is searched within
00103 *      user given positions.
00104 *      The least squares fit is done by using a box smaller than
00105 *      the size of two slitlets 
00106 *   9) fits the beginning and end position of the slitlets
00107 *      by using non-linear least square fitting of a Boltzmann function
00108 *      fits a Boltzmann function to the slitlet edges exposed and indicated
00109 *      by the brightest emission lines. The slitlet is searched within
00110 *      user given positions.
00111 *      The least squares fit is done by using a box smaller than
00112 *      the size of two slitlets 
00113 *
00114 *   FILES
00115 *
00116 *   ENVIRONMENT
00117 *
00118 *   RETURN VALUES
00119 *
00120 *   CAUTIONS
00121 *
00122 *   EXAMPLES
00123 *
00124 *   SEE ALSO
00125 *
00126 *   BUGS
00127 *
00128 *------------------------------------------------------------------------
00129 */
00130 
00131 #include "vltPort.h"
00132 
00133 /*
00134  * System Headers
00135  */
00136 
00137 /*
00138  * Local Headers
00139  */
00140 
00141 #include "absolute.h"
00142 #include "recipes.h"
00143 
00144 /*----------------------------------------------------------------------------
00145  *                              Defines
00146  *--------------------------------------------------------------------------*/
00147 
00148 #define XDIMA         1         /* dimension of the x values */
00149 #define TOLA          0.001     /* fitting tolerance */
00150 #define LABA          0.1       /* labda parameter */
00151 #define ITSA          200       /* maximum number of iterations */
00152 #define LABFACA       10.0      /* labda step factor */
00153 #define LABMAXA       1.0e+10   /* maximum value for labda */
00154 #define LABMINA       1.0e-10   /* minimum value for labda */
00155 #define NPAR          4         /* number of fit parameters */
00156 
00157 /*----------------------------------------------------------------------------
00158  *                                 Local variables
00159  *--------------------------------------------------------------------------*/
00160 
00161 static double chi1 ;                    /* old reduced chi-squared */
00162 static double chi2 ;                    /* new reduced chi-squared */
00163 static double labda ;                   /* mixing parameter */
00164 static double vec[NPAR] ;               /* correction vector */
00165 static double matrix1[NPAR][NPAR] ;     /* original matrix */
00166 static double matrix2[NPAR][NPAR] ;     /* inverse of matrix1 */
00167 static int    nfree ;                   /* number of free parameters */
00168 static int    parptr[NPAR] ;            /* parameter pointer */
00169 
00170 /*----------------------------------------------------------------------------
00171  *                  Functions private to this module
00172  *--------------------------------------------------------------------------*/
00173 
00174 static int inv_mat (void) ;
00175 
00176 static void get_mat ( float * xdat,
00177                       int   * xdim,
00178                       float * ydat,
00179                       float * wdat,
00180                       int   * ndat,
00181                       float * fpar,
00182                       float * epar/*,
00183                       int   * npar*/ ) ;
00184 
00185 static int get_vec ( float * xdat,
00186                      int   * xdim,
00187                      float * ydat,
00188                      float * wdat,
00189                      int   * ndat,
00190                      float * fpar,
00191                      float * epar,
00192                      int   * npar ) ;
00193 
00194 /*----------------------------------------------------------------------------
00195  *                          Function codes
00196  *--------------------------------------------------------------------------*/
00197 
00198 /*----------------------------------------------------------------------------
00199    Function     :       boltz()
00200    In           :       position array xdat, parameter list parlist
00201                         The parameters are:
00202                         parlist(0): background1
00203                         parlist(1): background2
00204                         parlist(2): central position
00205                         parlist(3): width
00206    Out          :       function value of a Boltzmann function
00207                         that is 
00208                         y = (parlist(0) - parlist(1)) / (1+exp((x-parlist(2))/parlist(3))) + parlist(1)
00209    Job          :       calculates the value of a Boltzmann function with parameters 
00210                         parlist at the position xdat 
00211  ---------------------------------------------------------------------------*/
00212 
00213 float boltz ( float * xdat, float * parlist )
00214 {
00215     float return_value ;
00216 
00217     /* now build the boltzman function out of the parameters */
00218     return_value = 
00219     (parlist[0] - parlist[1]) / (1 + exp(( xdat[0] - parlist[2] ) / parlist[3])) + parlist[1] ;
00220     
00221     return return_value ;
00222 }
00223        
00224 /*----------------------------------------------------------------------------
00225    Function     :       boltz_deriv()
00226    In           :       position array xdat, parameter list parlist 
00227                         The parameters are:
00228                         parlist(0): background1
00229                         parlist(1): background2
00230                         parlist(2): central position
00231                         parlist(3): width
00232                         derivative value of a Boltzmann function at position xdat: dervs
00233                         dervs[0]: partial derivative by background1
00234                         dervs[1]: partial derivative by background2 
00235                         dervs[2]: partial derivative by central position
00236                         dervs[3]: partial derivative by the width
00237    Out          :       nothing, void
00238    Job          :       calculates the partial derivatives for a Boltzmann function with
00239                         parameters parlist at position xdat 
00240 ----------------------------------------------------------------------------*/
00241 
00242 void boltz_deriv( float * xdat, float * parlist, float * dervs )
00243 {
00244     float subst ;
00245  
00246     subst = (xdat[0] - parlist[2]) / parlist[3] ;
00247 
00248     dervs[0] = 1. / ( 1. + exp(subst) ) ;
00249 
00250     dervs[1] = -1. / ( 1. + exp(subst) ) + 1. ;
00251 
00252     dervs[2] = ( (parlist[0] - parlist[1]) / parlist[3] * exp(subst) ) /
00253                ( (1. + exp(subst)) * (1. + exp(subst)) ) ;
00254 
00255     dervs[3] = ( (parlist[0] - parlist[1]) * (xdat[0] - parlist[2]) /
00256                (parlist[3]*parlist[3]) * exp(subst) ) /
00257                ( (1. + exp(subst)) * (1. + exp(subst)) ) ;
00258 }
00259 
00260 /*----------------------------------------------------------------------------
00261    Function     :       inv_mat()
00262    In           :       void
00263    Out          :       integer (0 if it worked, -6 if determinant is zero) 
00264    Job          :       calculates the inverse of matrix2. The algorithm used 
00265                         is the Gauss-Jordan algorithm described in Stoer,
00266                         Numerische Mathematik, 1. Teil.
00267  ---------------------------------------------------------------------------*/
00268 
00269 static int inv_mat (void)
00270 {
00271     double even ;
00272     double hv[NPAR] ;
00273     double mjk ;
00274     double rowmax ;
00275     int evin ;
00276     int i, j, k, row ;
00277     int per[NPAR] ;
00278    
00279     /* set permutation array */
00280     for ( i = 0 ; i < nfree ; i++ )
00281     {
00282         per[i] = i ;
00283     }
00284     
00285     for ( j = 0 ; j < nfree ; j++ ) /* in j-th column */
00286     {
00287         /* determine largest element of a row */                                   
00288         rowmax = fabs ( matrix2[j][j] ) ;
00289         row = j ;                         
00290 
00291         for ( i = j + 1 ; i < nfree ; i++ )
00292         {
00293             if ( fabs ( matrix2[i][j] ) > rowmax )
00294             {
00295                 rowmax = fabs( matrix2[i][j] ) ;
00296                 row = i ;
00297             }
00298         }
00299 
00300         /* determinant is zero! */
00301         if ( matrix2[row][j] == 0.0 )
00302         {
00303             return -6 ;
00304         }
00305         
00306         /* if the largest element is not on the diagonal, then permutate rows */
00307         if ( row > j )
00308         {
00309             for ( k = 0 ; k < nfree ; k++ )
00310             {
00311                 even = matrix2[j][k] ;
00312                 matrix2[j][k] = matrix2[row][k] ;
00313                 matrix2[row][k] = even ;
00314             }
00315             /* keep track of permutation */
00316             evin = per[j] ;
00317             per[j] = per[row] ;
00318             per[row] = evin ;
00319         }
00320         
00321         /* modify column */
00322         even = 1.0 / matrix2[j][j] ;
00323         for ( i = 0 ; i < nfree ; i++ )
00324         {
00325             matrix2[i][j] *= even ;
00326         }
00327         matrix2[j][j] = even ;
00328         
00329         for ( k = 0 ; k < j ; k++ )
00330         {
00331             mjk = matrix2[j][k] ;
00332             for ( i = 0 ; i < j ; i++ )
00333             {
00334                 matrix2[i][k] -= matrix2[i][j] * mjk ;
00335             }
00336             for ( i = j + 1 ; i < nfree ; i++ )
00337             {
00338                 matrix2[i][k] -= matrix2[i][j] * mjk ;
00339             }
00340             matrix2[j][k] = -even * mjk ;
00341         }
00342     
00343         for ( k = j + 1 ; k < nfree ; k++ )
00344         {
00345             mjk = matrix2[j][k] ;
00346             for ( i = 0 ; i < j ; i++ )
00347             {
00348                 matrix2[i][k]  -= matrix2[i][j] * mjk ;
00349             }
00350             for ( i = j + 1 ; i < nfree ; i++ )
00351             {
00352                 matrix2[i][k]  -= matrix2[i][j] * mjk ;
00353             }
00354             matrix2[j][k] = -even * mjk ;
00355         }
00356     }
00357     
00358     /* finally, repermute the columns */
00359     for ( i = 0 ; i < nfree ; i++ )
00360     {
00361         for ( k = 0 ; k < nfree ; k++ )
00362         {
00363             hv[per[k]] = matrix2[i][k] ;
00364         }
00365         for ( k = 0 ; k < nfree ; k++ )
00366         {
00367             matrix2[i][k] = hv[k] ;
00368         }
00369     }
00370         
00371     /* all is well */
00372     return 0 ;
00373 }
00374     
00375 /*----------------------------------------------------------------------------
00376    Function     :       get_mat()
00377    In           :       xdat: position
00378                         xdim: factor of the indices of the position array
00379                         ydat: real data
00380                         wdat: weights
00381                         ndat: number of data points
00382                         fpar: function parameters
00383                         epar: partial derivatives of the function
00384                         npar: number of function parameters
00385    Out          :       void
00386    Job          :       builds the matrix 
00387  ---------------------------------------------------------------------------*/
00388 
00389 static void get_mat ( float * xdat,
00390                       int   * xdim,
00391                       float * ydat,
00392                       float * wdat,
00393                       int   * ndat,
00394                       float * fpar,
00395                       float * epar/*,
00396                       int   * npar*/ )
00397 {
00398     double wd ;
00399     double wn ;
00400     double yd ;
00401     int i, j, n ;
00402 
00403     for ( j = 0 ; j < nfree ; j++ )
00404     {
00405         vec[j] = 0.0 ; /* zero vector */
00406         for ( i = 0 ; i<= j ; i++ )   /* zero matrix only on and below diagonal */
00407         {
00408             matrix1[j][i] = 0.0 ;
00409         }
00410     }
00411     chi2 = 0.0 ;  /* reset reduced chi-squared */
00412     
00413     /* loop through data points */
00414     for ( n = 0 ; n < (*ndat) ; n++ )
00415     {
00416         wn = wdat[n] ;
00417         if ( wn > 0.0 )  /* legal weight ? */
00418         {
00419             yd = ydat[n] - boltz( &xdat[(*xdim) * n], fpar ) ;
00420             boltz_deriv( &xdat[(*xdim) * n], fpar, epar ) ;
00421             chi2 += yd * yd * wn ; /* add to chi-squared */
00422             for ( j = 0 ; j < nfree ; j++ )
00423             {
00424                 wd = epar[parptr[j]] * wn ;  /* weighted derivative */
00425                 vec[j] += yd * wd ;       /* fill vector */
00426                 for ( i = 0 ; i <= j ; i++ ) /* fill matrix */
00427                 {
00428                     matrix1[j][i] += epar[parptr[i]] * wd ;
00429                 }
00430             }
00431         }
00432     }                   
00433 }  
00434                 
00435             
00436 /*----------------------------------------------------------------------------
00437    Function     :       get_vec()
00438    In           :       xdat: position
00439                         xdim: factor of the indices of the position array
00440                         ydat: real data
00441                         wdat: weights
00442                         ndat: number of data points
00443                         fpar: function parameters
00444                         epar: partial derivatives of the function
00445                         npar: number of function parameters
00446    Out          :       integer (0 if it had worked, -5 or -7 
00447                                  if diagonal element is wrong, or -6,
00448                                  if determinant is zero )
00449    Job          :       calculates the correction vector. The matrix has been
00450                         built by get_mat(), we only have to rescale it for the 
00451                         current value of labda. The matrix is rescaled so that
00452                         the diagonal gets the value 1 + labda.
00453                         Next we calculate the inverse of the matrix and then
00454                         the correction vector.
00455  ---------------------------------------------------------------------------*/
00456             
00457 static int get_vec ( float * xdat,
00458                      int   * xdim,
00459                      float * ydat,
00460                      float * wdat,
00461                      int   * ndat,
00462                      float * fpar,
00463                      float * epar,
00464                      int   * npar )
00465 {
00466     double dj ;
00467     double dy ;
00468     double mii ;
00469     double mji ;
00470     double mjj ;
00471     double wn ;
00472     int i, j, n, r ;
00473 
00474     /* loop to modify and scale the matrix */
00475     for ( j = 0 ; j < nfree ; j++ )
00476     {
00477         mjj = matrix1[j][j] ;
00478         if ( mjj <= 0.0 )             /* diagonal element wrong */
00479         {
00480             return -5 ;
00481         }
00482         mjj = sqrt( mjj ) ;
00483         for ( i = 0 ; i < j ; i++ )
00484         {
00485             mji = matrix1[j][i] / mjj / sqrt( matrix1[i][i] ) ;
00486             matrix2[i][j] = matrix2[j][i] = mji ;
00487         }
00488         matrix2[j][j] = 1.0 + labda ;  /* scaled value on diagonal */
00489     }    
00490     
00491     if ( (r = inv_mat()) )                /* invert matrix inlace */
00492     {
00493         return r ;
00494     }
00495     
00496     for ( i = 0 ; i < (*npar) ; i ++ )
00497     {
00498         epar[i] = fpar[i] ;
00499     }
00500     
00501     /* loop to calculate correction vector */
00502     for ( j = 0 ; j < nfree ; j++ )
00503     {
00504         dj = 0.0 ;
00505         mjj = matrix1[j][j] ;
00506         if ( mjj <= 0.0)               /* not allowed */
00507         {
00508             return -7 ;
00509         }
00510         mjj = sqrt ( mjj ) ;
00511         for ( i = 0 ; i < nfree ; i++ )
00512         {
00513             mii = matrix1[i][i] ;
00514             if ( mii <= 0.0 )
00515             {
00516                 return -7 ;
00517             }
00518             mii = sqrt( mii ) ;
00519             dj += vec[i] * matrix2[j][i] / mjj / mii ;
00520         }
00521         epar[parptr[j]] += dj ;       /* new parameters */
00522     }    
00523     chi1 = 0.0 ;                      /* reset reduced chi-squared */
00524  
00525     /* loop through the data points */
00526     for ( n = 0 ; n < (*ndat) ; n++ )
00527     {
00528         wn = wdat[n] ;               /* get weight */
00529         if ( wn > 0.0 )              /* legal weight */
00530         {
00531             dy = ydat[n] - boltz( &xdat[(*xdim) * n], epar) ;
00532             chi1 += wdat[n] * dy * dy ;
00533         }
00534     }
00535     return 0 ;
00536 }   
00537     
00538         
00539 
00540 /*----------------------------------------------------------------------------
00541    Function     :       lsqfit()
00542    In           :       xdat: position, coordinates of data points.
00543                               xdat is 2 dimensional: XDAT ( XDIM, NDAT )
00544                         xdim: dimension of fit
00545                         ydat: data points
00546                         wdat: weights for data points
00547                         ndat: number of data points
00548                         fpar: on input contains initial estimates of the 
00549                               parameters for non-linear fits, on output the
00550                               fitted parameters.
00551                         epar: contains estimates of the errors in fitted parameters
00552                         mpar: logical mask telling which parameters are free (non-zero)
00553                               and which parameters are fixed (0)
00554                         npar: number of function parameters ( free + fixed )
00555                         tol:  relative tolerance. lsqfit stops when successive iterations
00556                               fail to produce a decrement in reduced chi-squared less
00557                               than tol. If tol is less than the minimum tolerance 
00558                               possible, tol will be set to this value. This means
00559                               that maximum accuracy can be obtained by setting
00560                               tol = 0.0.
00561                         its:  maximum number of iterations
00562                         lab:  mixing parameter, lab determines the initial weight
00563                               of steepest descent method relative to the Taylor method
00564                               lab should be a small value (i.e. 0.01). lab can only
00565                               be zero when the partial derivatives are independent
00566                               of the parameters. In fact in this case lab should be
00567                               exactly equal to zero.
00568    Out          :       returns number of iterations needed to achieve convergence
00569                         according to tol. When this number is negative, the fitting
00570                         was not continued because a fatal error occurred:
00571                         -1 too many free parameters, maximum is 32
00572                         -2 no free parameters
00573                         -3 not enough degrees of freedom
00574                         -4 maximum number of iterations too small to obtain
00575                            a solution which satisfies tol.
00576                         -5 diagonal of matrix contains elements which are zero
00577                         -6 determinant of the coefficient matrix is zero
00578                         -7 square root of a negative number 
00579    Job          :       this is a routine for making a least-squares fit of a
00580                         function to a set of data points. The method used is
00581                         described in: Marquardt, J.Soc.Ind.Appl.Math. 11. 431 (1963).
00582                         This method is a mixture of the steepest descent method 
00583                         and the Taylor method.
00584  ---------------------------------------------------------------------------*/
00585 
00586 int lsqfit ( float * xdat,
00587              int   * xdim,
00588              float * ydat,
00589              float * wdat,
00590              int   * ndat,
00591              float * fpar,
00592              float * epar,
00593              int   * mpar,
00594              int   * npar,
00595              float * tol ,
00596              int   * its ,
00597              float * lab  )
00598 {
00599     int i, n, r ;
00600     int itc ;                      /* fate of fit */
00601     int found ;                    /* fit converged: 1, not yet converged: 0 */
00602     int  nuse ;                    /* number of useable data points */
00603     double tolerance ;             /* accuracy */
00604 
00605     itc   = 0 ;                    /* fate of fit */
00606     found = 0 ;                    /* reset */
00607     nfree = 0 ;                    /* number of free parameters */
00608     nuse  = 0 ;                    /* number of legal data points */
00609 
00610     if ( *tol < (FLT_EPSILON * 10.0 ) )
00611     {
00612         tolerance = FLT_EPSILON * 10.0 ;  /* default tolerance */
00613     }
00614     else
00615     {
00616         tolerance = *tol ;                /* tolerance */
00617     }
00618     
00619     labda = fabs( *lab ) * LABFACA ;       /* start value for mixing parameter */
00620     for ( i = 0 ; i < (*npar) ; i++ )
00621     {
00622         if ( mpar[i] )
00623         {
00624             if ( nfree > NPAR )         /* too many free parameters */
00625             {
00626                 return -1 ;
00627             }
00628             parptr[nfree++] = i ;         /* a free parameter */
00629         }
00630     }
00631     
00632     if (nfree == 0)                       /* no free parameters */     
00633     {
00634         return -2 ;
00635     }
00636     
00637     for ( n = 0 ; n < (*ndat) ; n++ )
00638     {
00639         if ( wdat[n] > 0.0 )              /* legal weight */
00640         {
00641             nuse ++ ;
00642         }
00643     }
00644     
00645     if ( nfree >= nuse )
00646     {
00647         return -3 ;                       /* no degrees of freedom */
00648     }
00649     if ( labda == 0.0 )                   /* linear fit */
00650     {
00651         for ( i = 0 ; i < nfree ; fpar[parptr[i++]] = 0.0 ) ;  /* initialize fpar array */
00652         get_mat ( xdat, xdim, ydat, wdat, ndat, fpar, epar/*, npar*/ ) ;
00653         r =  get_vec ( xdat, xdim, ydat, wdat, ndat, fpar, epar, npar ) ;
00654         if ( r )                         /* error */
00655         {
00656             return r ;
00657         }
00658         for ( i = 0 ; i < (*npar) ; i++ )
00659         {
00660             fpar[i] = epar[i] ;           /* save new parameters */
00661             epar[i] = 0.0 ;               /* and set errors to zero */
00662         }
00663         chi1 = sqrt( chi1 / (double) (nuse - nfree) ) ;
00664         for ( i = 0 ; i < nfree ; i++ )
00665         {
00666             if ( (matrix1[i][i] <= 0.0 ) || (matrix2[i][i] <= 0.0) )
00667             {
00668                 return -7 ;
00669             }
00670             epar[parptr[i]] = chi1 * sqrt( matrix2[i][i] ) / sqrt( matrix1[i][i] ) ;
00671         }
00672     }
00673     else                                  /* non-linear fit */
00674     {
00675         /*----------------------------------------------------------------
00676          * the non-linear fit uses the steepest descent method in combination
00677          * with the Taylor method. The mixing of these methods is controlled
00678          * by labda. In the outer loop ( called the iteration loop ) we build
00679          * the matrix and calculate the correction vector. In the inner loop
00680          * (called the interpolation loop) we check whether we have obtained a
00681          * better solution than the previous one. If so, we leave the inner loop
00682          * else we increase labda ( give more weight to the steepest descent method)
00683          * calculate the correction vector and check again. After the inner loop
00684          * we do a final check on the goodness of the fit and if this satisfies
00685          * the tolerance we calculate the errors of the fitted parameters.
00686          */
00687         while ( !found )                  /* iteration loop */
00688         {      
00689             if ( itc++ == (*its) )        /* increase iteration counter */
00690             {
00691                 return -4 ;               
00692             }
00693             get_mat( xdat, xdim, ydat, wdat, ndat, fpar, epar/*, npar */) ;
00694             
00695             /*-------------------------------------------------------------
00696              * here we decrease labda since we may assume that each iteration
00697              * brings us closer to the answer.
00698              */
00699             if ( labda > LABMINA )
00700             {
00701                 labda = labda / LABFACA ;         /* decrease labda */
00702             }
00703             r = get_vec ( xdat, xdim, ydat, wdat, ndat, fpar, epar, npar ) ;
00704 
00705             if ( r )                      /* error */
00706             {
00707                 return r ;
00708             }
00709 
00710             while ( chi1 >= chi2 )        /* interpolation loop */
00711             {
00712                 /*-----------------------------------------------------------
00713                  * The next statement is based on experience, not on the mathematics
00714                  * of the problem. It is assumed that we have reached convergence 
00715                  * when the pure steepest descent method does not produce a better
00716                  * solution.
00717                  */
00718                 if ( labda > LABMAXA )    /* assume solution found */
00719                 {
00720                     break ;
00721                 }
00722                 labda = labda * LABFACA ;         /* increase mixing parameter */
00723                 r = get_vec ( xdat, xdim, ydat, wdat, ndat, fpar, epar, npar ) ;
00724 
00725                 if ( r )                  /* error */
00726                 {
00727                     return r ;
00728                 }
00729             }
00730 
00731             if ( labda <= LABMAXA )        /* save old parameters */
00732             {
00733                 for ( i = 0 ; i < *npar ; i++ )
00734                 {
00735                     fpar[i] = epar[i] ;
00736                 }
00737             }
00738             if ( (fabs( chi2 - chi1 ) <= (tolerance * chi1)) || (labda > LABMAXA) )
00739             {
00740                 /*---------------------------------------------------------------------
00741                  * we have a satisfying solution, so now we need to calculate the correct
00742                  * errors of the fitted parameters. This we do by using the pure Taylor
00743                  * method because we are very close to the real solution.
00744                  */
00745                 labda = LABMINA ;              /* for Taylor solution */
00746                 get_mat ( xdat, xdim, ydat, wdat, ndat, fpar, epar/*, npar */) ;
00747                 r = get_vec ( xdat, xdim, ydat, wdat, ndat, fpar, epar, npar ) ;
00748 
00749                 if ( r )                    /* error */
00750                 {
00751                     return r ;
00752                 }
00753                 for ( i = 0 ; i < (*npar) ; i++ )
00754                 {
00755                     epar[i] = 0.0 ;          /* set error to zero */
00756                 }
00757                 chi2 = sqrt ( chi2 / (double) (nuse - nfree) ) ;
00758 
00759                 for ( i = 0 ; i < nfree ; i++ )
00760                 {
00761                     if ( (matrix1[i][i] <= 0.0) || (matrix2[i][i] <= 0.0) )
00762                     {
00763                         return -7 ;
00764                     }
00765                     epar[parptr[i]] = chi2 * sqrt( matrix2[i][i] ) / sqrt( matrix1[i][i] ) ;
00766                 }
00767                 found = 1 ;                  /* we found a solution */
00768             }
00769         }
00770     }
00771     return itc ;                             /* return number of iterations */
00772 }
00773 
00774 /*----------------------------------------------------------------------------
00775    Function     :       fitSlitsBoltz()
00776    In           :       lineImage:  emission line frame
00777                         par:        fit parameter data structure of fitted lines
00778                         slit_pos:   allocated dummy array for the slitlet positions [32][2]
00779                         box_length: pixel length of the row box within the fit is done
00780                         y_box:      small box in spectral direction within the slitlet
00781                                     may lie.
00782                         diff_tol:   maximum tolerable difference of the resulting fit position
00783                                     with respect to the expected position. If difference is 
00784                                     greater the expected position is taken.
00785    Out          :       slit_pos:  beginning and end position of the slitlets to
00786                                    sub-pixel accuracy
00787                          0  if it worked,
00788                         -1  if there was no line image given,
00789                         -2  if there were no line fit parameters given,
00790                         -3  if there was no dummy array for the slit positions
00791                             allocated
00792                         -4  if the given box length is impossible
00793                         -5  if the given y box length is impossible
00794                         -6  if the given difference tolerance is too small
00795                         -7  if there were no emission lines found in the first image columns
00796                         -8  if not all slitlets could be found
00797    Job          :       fits the beginning and end position of the slitlets
00798                         by using non-linear least square fitting of a Boltzmann function
00799                         fits a Boltzmann function to the slitlet edges exposed and indicated
00800                         by the brightest emission lines. To achieve this, the fit
00801                         parameters are used to find the brightest emission line
00802                         and to get its position for each column.
00803                         The least squares fit is done by using a box smaller than
00804                         the size of two slitlets 
00805  ---------------------------------------------------------------------------*/
00806 
00807 int fitSlitsBoltz ( OneImage   * lineImage, 
00808                     FitParams ** par,
00809                     float     ** slit_pos,
00810                     int          box_length,
00811                     float        y_box,
00812                     float        diff_tol )
00813 {
00814     float position[lineImage -> lx] ;
00815     int   * edge, * edgeclean ;
00816     int   * dummyedge ;
00817     int   * pos_row, * pos_rowclean ;
00818     Vector * box_buffer ;
00819     Vector * half_buffer ;
00820     Vector * in_buffer ;
00821     float max_intensity ;
00822     float row_pos ;
00823     int   row, col ;
00824     int   i, j, k, m, n, ed ;
00825     int   found, init1 ;
00826     int   line ; 
00827     int   nel, n_right, left_right ;
00828     int   n_buf, edge_ind, shift ;
00829     int   column ;
00830     int   slit_length ;
00831     int   agreed ;
00832     int   bad_line ;
00833     int   margin ;
00834     int   iters, xdim, ndat ;
00835     int   numpar, its ;
00836     int   * mpar ;
00837     float * xdat, * wdat ;
00838     float tol, lab ;
00839     float fitpar[NPAR] ;
00840     float dervpar[NPAR] ;
00841     float minval, maxval ;
00842     float min ;
00843     float pos, last_pos ;
00844     int old_col=0;
00845     int old_pos=0;
00846 
00847 
00848     slit_length = SLITLENGTH ;
00849     if ( NullImage == lineImage )
00850     {
00851         cpl_msg_error( "fitSlitsBoltz:"," no line image given!\n" ) ;
00852         return -1 ;
00853     }
00854 
00855     if ( NULL == par )
00856     {
00857         cpl_msg_error( "fitSlitsBoltz:"," no line fit parameters given!\n" ) ;
00858         return -2 ;
00859     }
00860 
00861     if ( NULL == slit_pos )
00862     {
00863         cpl_msg_error( "fitSlitsBoltz:"," no position array allocated!\n" ) ;
00864         return -3 ;
00865     }
00866 
00867     if ( box_length <  4 ||
00868          box_length > 2*slit_length )
00869     {
00870         cpl_msg_error( "fitSlitsBoltz:"," wrong fitting box length given!\n" ) ;
00871         return -4 ;
00872     }
00873 
00874     if ( y_box <= 0.  || y_box > 6. )
00875     {
00876         cpl_msg_error( "fitSlitsBoltz:"," wrong y box length given!\n" ) ;
00877         return -5 ;
00878     }
00879 
00880     if ( diff_tol < 1. )
00881     {
00882         cpl_msg_error( "fitSlitsBoltz:"," diff_tol too small!\n" ) ;
00883         return -6 ;
00884     }
00885  
00886     /* allocate memory for the edges and the row positon of the slitlets */
00887     edge         = (int*) cpl_calloc( 3*slit_length, sizeof(int) ) ;
00888     dummyedge    = (int*) cpl_calloc( 3*slit_length, sizeof(int) ) ;
00889     edgeclean    = (int*) cpl_calloc( slit_length-1, sizeof(int) ) ;
00890     pos_row      = (int*) cpl_calloc( 3*slit_length, sizeof(int) ) ;
00891     pos_rowclean = (int*) cpl_calloc( slit_length, sizeof(int) ) ;
00892 
00893 
00894     /* --------------------------------------------------------------------------
00895      * go through the first image columns and the fit parameters and find the line 
00896      * with the highest intensity 
00897      */
00898     agreed = -1 ;
00899     bad_line = -1 ;
00900     while( agreed == -1 )
00901     {
00902         found = -1 ;
00903         max_intensity = -FLT_MAX ;
00904         for ( col = 0 ; col < slit_length ; col++ )
00905         {
00906             for ( i = 0 ; i < par[0]->n_params ; i++ )
00907             {
00908                 if ( par[i]->column == col && par[i]->line != bad_line )
00909                 {
00910                     if ( par[i]->fit_par[0] > max_intensity )
00911                     {
00912                         if ( par[i]->fit_par[1] >= 1. && par[i]->fit_par[2] > 0. )
00913                         {
00914                             max_intensity = par[i]->fit_par[0] ;
00915                             found = i ;
00916                         }
00917                     }
00918                 }
00919             }  
00920         }
00921 
00922         /* --------------------------------------------------------------------
00923          * check if the found line is usable and if the neighbouring line 
00924          * have intensity on near rows in neighbouring slitlets 
00925          */
00926         line    = par[found]->line ;
00927         column  = par[found]->column ;
00928         row_pos = par[found]->fit_par[2] ;
00929         if ( found >= 0 && max_intensity > 0. )
00930         {
00931             for ( i = 0 ; i < par[0]->n_params ; i++ )
00932             {
00933                 if ( par[i]->line == line-1 && 
00934                      par[i]->column == column + slit_length )
00935                 {
00936                     if ( par[i]->fit_par[2] <= (row_pos + y_box) &&
00937                          par[i]->fit_par[2] >= (row_pos - y_box) )
00938                     {
00939                         bad_line = line ;
00940                     } 
00941                 }
00942             }
00943             if ( bad_line != line )
00944             {
00945                 agreed = 1 ;
00946                 break ;
00947             }
00948         }
00949         else 
00950         {
00951             cpl_msg_error("fitSlitsBoltz:"," no emission line found in the first image columns\n") ;
00952             cpl_free( edge ) ;
00953             cpl_free( pos_row ) ;
00954             cpl_free( edgeclean ) ;
00955             cpl_free( dummyedge ) ;
00956             cpl_free( pos_rowclean ) ;
00957             return -7 ;
00958         }    
00959     }
00960 
00961 
00962 
00963     if ( agreed == -1 )
00964     {
00965         cpl_msg_error("fitSlitsBoltz:"," no emission line found in the first image columns\n") ;
00966         cpl_free( edge ) ;
00967         cpl_free( pos_row ) ;
00968         cpl_free( edgeclean ) ;
00969         cpl_free( dummyedge ) ;
00970         cpl_free( pos_rowclean ) ;
00971         return -7 ;
00972     }    
00973 
00974  
00975     /* now find and store the raw edge positions of the found slitlet */ 
00976     n  = 0 ;
00977     ed = 0 ;
00978     /* was for ( col = 0 ; col < lineImage -> lx - slit_length/2 ; col++ ) */
00979     for ( col = slit_length/2 ; col < lineImage -> lx - slit_length/2 ; col++ )
00980     {
00981         for ( i = 0 ; i < par[0]->n_params ; i++ )
00982       {   
00983         /*
00984               printf("p1=%f p2=%f p3=%f\n",
00985               par[i]->fit_par[0],par[i]->fit_par[1],par[i]->fit_par[2]);
00986         */
00987             if ( par[i]->column == col && par[i]->line == line )
00988             {
00989                 if ( par[i]->fit_par[0] > 0.  && 
00990                      par[i]->fit_par[1] >= 1. && 
00991                      par[i]->fit_par[2] > 0. )
00992                 {
00993                     position[n] = par[i]->fit_par[2] ;
00994                     old_pos=position[n];
00995                     if ( n > 0 && 
00996                          fabs(position[n] - position[n-1]) > y_box && 
00997                          (col-old_col) > (slit_length-SLIT_POS_ERR) )
00998                     {
00999 
01000               old_col=col;
01001               edge[ed] = col ; 
01002               pos_row[ed] = nint( position[n-1] ) ;
01003                       /* printf("edge[%d]=%d , pos_row=%d\n",ed,edge[ed],pos_row[ed]); */ 
01004               ed++ ;
01005               if ( col >= lineImage -> lx - slit_length - SLIT_POS_ERR ) {
01006               pos_row[ed] =  nint( position[n] ) ;
01007               } 
01008             } else if ( ((col-old_col) > (slit_length+SLIT_POS_ERR)) &&                                 (col>120) ) {
01009               old_col=col;
01010               edge[ed] = col ; 
01011               pos_row[ed] = nint( position[n-1] ) ;
01012                       cpl_msg_warning("fitSlitsBoltz:","added1 slitlet edge[%d]=%d, pos_row=%d",
01013                      ed,edge[ed],pos_row[ed]);
01014               ed++ ;
01015               if ( col >= lineImage -> lx - slit_length - SLIT_POS_ERR ) {
01016               pos_row[ed] =  nint( position[n] ) ;
01017               } 
01018             }
01019                     n++ ;
01020                 }
01021             } else if ( ((col-old_col) > (slit_length+SLIT_POS_ERR)) && 
01022                         (col>120) ) {
01023           /*
01024               printf("check col=%d col-old_col=%d check=%d\n",
01025                      col,(col-old_col),(slit_length+SLIT_POS_ERR));
01026           */
01027           position[n] = old_pos ;
01028              
01029           old_col+=slit_length;
01030           edge[ed] = old_col; ; 
01031           pos_row[ed] = nint( position[n-1] ) ;
01032 
01033 
01034           cpl_msg_warning("fitSlitsBoltz:","added2 slitlet edge[%d]=%d, pos_row=%d",
01035                   ed,edge[ed],pos_row[ed]);
01036           ed++ ;
01037           if ( old_col >= lineImage -> lx - slit_length - SLIT_POS_ERR ) {
01038         pos_row[ed] =  old_pos ;
01039           }
01040           n++;
01041         }
01042         }
01043     }
01044 
01045 
01046     if ( ed < (N_SLITLETS - 1) )
01047     {
01048         cpl_msg_error("fitSlitsBoltz:"," not enough slitlets, found: %d\n", ed) ;
01049         cpl_free( edge ) ;
01050         cpl_free( pos_row ) ;
01051         cpl_free( edgeclean ) ;
01052         cpl_free( dummyedge ) ;
01053         cpl_free( pos_rowclean ) ;
01054         return -8 ;
01055     } 
01056 
01057     /* now find the clean edge and row positions of the slitlets */
01058     /* printf("ed=%d\n",ed); */
01059     for ( i = 1 ; i <= ed ; i ++ )
01060     {
01061         if ( i == ed )
01062         {
01063             if ( (edge[i-1] - edge[i-2]) < slit_length - SLIT_LEN_ERR ||
01064                  (edge[i-1] - edge[i-2]) > slit_length + SLIT_LEN_ERR )
01065             {
01066         /* printf("e(i-1)=%d e(i-2)=%d i=%d\n",edge[i-1], edge[i-2],i); */
01067                 dummyedge[i-1]   = -1 ;
01068             }
01069         }
01070         if (dummyedge[i-1] != -1)
01071         {
01072             dummyedge[i-1] = edge[i-1] ;
01073         }
01074         else
01075         {
01076             continue ;
01077         }
01078         if ( i < ed )
01079         {
01080             if ( (edge[i] - edge[i-1]) < slit_length - SLIT_LEN_ERR ||
01081                  (edge[i] - edge[i-1]) > slit_length + SLIT_LEN_ERR )
01082             {
01083         /* printf("e(i)=%d e(i-1)=%d i=%d\n",edge[i], edge[i-1],i); */
01084                 dummyedge[i]   = -1 ;
01085             }
01086         }
01087         if ( i+1 < ed && dummyedge[i] != -1 )
01088         {
01089             if ( (edge[i+1] - edge[i]) < slit_length - SLIT_LEN_ERR ||
01090                  (edge[i+1] - edge[i]) > slit_length + SLIT_LEN_ERR )
01091             {
01092         /* printf("e(i+1)=%d e(i)=%d i=%d\n",edge[i+1], edge[i],i); */
01093                 dummyedge[i+1] = -1 ; 
01094             }
01095         }
01096     }
01097 
01098     k = 0 ;
01099     for ( i = 0 ; i < ed ; i++ )
01100     {
01101         if ( dummyedge[i] != -1 && dummyedge[i] != 0 )
01102         {
01103             edgeclean[k] = dummyedge[i] ;
01104             pos_rowclean[k] = pos_row[i] ;
01105             k++ ;
01106             if( edgeclean[k-1] > (lineImage -> lx - slit_length -2*SLIT_LEN_ERR ) )
01107             {
01108                 pos_rowclean[k] = pos_row[ed] ;
01109             } 
01110         }
01111     }
01112     /*
01113     for ( i = 0 ; i < k ; i++ )
01114     {
01115       cpl_msg_warning("fitSlitsBoltz:","%d %d", edgeclean[i], pos_rowclean[i]);
01116     }
01117     */
01118     if ( k != N_SLITLETS - 1 )
01119     {
01120         cpl_msg_error("fitSlitsBoltz:"," wrong number of clean slitlets found: %d\n", k+1) ;
01121         cpl_free( edge ) ;
01122         cpl_free( pos_row ) ;
01123         cpl_free( edgeclean ) ;
01124         cpl_free( dummyedge ) ;
01125         cpl_free( pos_rowclean ) ;
01126         return -7 ;
01127     } 
01128 
01129     /* determine the margins of the fitting box outside the slitlets */
01130     margin = box_length / 2 ;
01131 
01132     /* --------------------------------------------------------------------------
01133      * now go through the slitlets, search along each column within a box with 
01134      * half width y_box the maximum value and store these found values in a buffer
01135      */
01136      if(
01137          ( (pos_rowclean[0]-nint(y_box)) < 0 ) ||
01138          ( (pos_rowclean[0]+nint(y_box)) >lineImage->ly )
01139        ) {
01140 
01141              cpl_msg_error("fitSlitsBoltz:",
01142               "pos_rowclean[0] <0 something wrong!\n") ;
01143              cpl_free( edge ) ;
01144              cpl_free( pos_row ) ;
01145              cpl_free( edgeclean ) ;
01146              cpl_free( dummyedge ) ;
01147              cpl_free( pos_rowclean ) ;
01148              return -7 ;
01149 
01150     }
01151 
01152     for ( j = 0 ; j <= k ; j++ )
01153     {
01154         m = 0 ;
01155         if ( j == 0 )
01156         {
01157             box_buffer = newVector( edgeclean[0] + margin ) ;
01158             for( col = 0 ; col < edgeclean[0] + margin ; col++ )
01159             {
01160                 maxval = -FLT_MAX ;
01161                 for ( row = pos_rowclean[0] - nint(y_box) ; row <= pos_rowclean[0] + nint(y_box) ; row++ )
01162                 {
01163                     if ( maxval < lineImage->data[col + lineImage->lx*row] )
01164                     {
01165                         maxval = lineImage -> data[col + lineImage->lx*row] ;
01166                     }
01167                 }
01168                 box_buffer->data[m] = maxval ;
01169                 m++ ;
01170             }
01171         }
01172         else if ( j < k )
01173         {
01174             box_buffer = newVector( edgeclean[j] - edgeclean[j-1] + 2*margin ) ;
01175             for ( col = edgeclean[j - 1] - margin ; col < edgeclean[j] + margin ; col++ )
01176             {
01177                 maxval = -FLT_MAX ;
01178                 for ( row = pos_rowclean[j] - nint(y_box) ; row <= pos_rowclean[j] + nint(y_box) ; row++ )
01179                 {
01180                     if ( maxval < lineImage->data[col + lineImage->lx*row] )
01181                     {
01182                         maxval = lineImage -> data[col + lineImage->lx*row] ;
01183                     }
01184                 }
01185                 box_buffer->data[m] = maxval ;
01186                 m++ ;
01187             }
01188         }
01189         else 
01190         {
01191             box_buffer = newVector( lineImage -> lx - edgeclean[k-1] + margin ) ;
01192             for ( col = edgeclean[k - 1] - margin ; col < lineImage -> lx ; col++ )
01193             {
01194                 maxval = -FLT_MAX ;
01195                 for ( row = pos_rowclean[k-2] - nint(y_box) ; row <= pos_rowclean[k-2] + nint(y_box) ; row++ )
01196                 {
01197                     if ( maxval < lineImage->data[col + lineImage->lx*row] )
01198                     {
01199                         maxval = lineImage -> data[col + lineImage->lx*row] ;
01200                     }
01201                 }
01202                 if(maxval>0) box_buffer->data[m] = maxval ;
01203         else box_buffer->data[m] = 0;
01204         m++ ;
01205             }
01206         }
01207 
01208         /* determine the minimum value in the box to get background1 value for the edge slitlets */
01209         min = FLT_MAX ;
01210         for ( i = 0 ; i < box_buffer->n_elements ; i++ )
01211         {
01212             if ( box_buffer -> data[i] < min )
01213             {
01214                 min = box_buffer -> data[i] ;
01215             }
01216         }
01217 
01218         for ( left_right = 0 ; left_right <= 1 ; left_right++ )
01219         { 
01220             nel = 0 ;
01221             if ( left_right == 0 )
01222             {
01223                 nel = box_buffer -> n_elements / 2 ;
01224             }
01225             else
01226             {
01227                 if ( box_buffer -> n_elements % 2 == 0 )
01228                 {
01229                     nel = box_buffer -> n_elements / 2 ;
01230                 }
01231                 else
01232                 {
01233                     nel = box_buffer -> n_elements / 2 + 1 ;
01234                 }
01235             }
01236 
01237             /* now split the buffer in the midth in a left and right part for fitting */
01238             half_buffer = newVector( nel ) ;
01239             if ( left_right == 0 )
01240             {
01241                 for ( i = 0 ; i < nel ; i++ )
01242                 {
01243                     half_buffer -> data[i] = box_buffer -> data[i] ;
01244                 }
01245             }
01246             else
01247             {
01248                 n_right = 0 ;
01249                 for ( i = box_buffer -> n_elements - 1 ; i >= box_buffer -> n_elements - nel ; i-- )
01250                 {
01251                     half_buffer -> data[n_right] = box_buffer -> data[i] ;
01252                     n_right++ ;
01253                 }
01254             }
01255 
01256             xdat = (float *) cpl_calloc( nel, sizeof (float) ) ;
01257             wdat = (float *) cpl_calloc( nel, sizeof (float) ) ;
01258             mpar = (int *)   cpl_calloc( NPAR, sizeof (int) ) ;
01259 
01260             /* set initial values for the fitting routine */
01261             minval =  FLT_MAX ;
01262             maxval = -FLT_MAX ;
01263             for ( i = 0 ; i < nel ; i++ )
01264             {
01265                 xdat[i] = i ;
01266                 wdat[i] = 1.0 ;
01267                 if ( half_buffer -> data[i] < minval )
01268                 {
01269                     minval = half_buffer -> data[i] ;
01270                 }
01271                 if ( half_buffer -> data[i] > maxval )
01272                 {
01273                     maxval = half_buffer -> data[i] ;
01274                 }
01275             }
01276 
01277             fitpar[0] = minval ;
01278             fitpar[1] = maxval ; 
01279 
01280             /* search for both positions of the half intensity of the hat within the buffer */
01281             init1 = -1 ; 
01282             for ( i = 0 ; i < nel ; i++ )
01283             {
01284                 if ( half_buffer -> data[i] >= ( maxval + minval ) / 2. )
01285                 {
01286                     init1 = i ;
01287                     break ;
01288                 }
01289             }
01290 
01291             /*------------------------------------------------------------------
01292              * if we have too few left background values (at the image edges)
01293              * the left margin of the buffer to fit is filled with the minimal
01294              * values in order to get a good fit
01295              */
01296               
01297             edge_ind = 0 ;
01298             if ( init1 < 3 )
01299             {
01300                 n_buf = half_buffer->n_elements + margin ;
01301                 in_buffer = newVector( n_buf ) ;     
01302                 for ( i = 0 ; i < margin ; i++ )
01303                 {
01304                     in_buffer -> data[i] = min ;
01305                 }
01306                 shift = 0 ;
01307                 for ( i = margin ; i < n_buf ; i++ )
01308                 {
01309                     in_buffer -> data[i] = half_buffer -> data[shift] ;
01310                     shift++ ;
01311                 }
01312                 destroyVector ( half_buffer ) ;
01313                 half_buffer = newVector ( n_buf ) ;
01314                 for ( i = 0 ; i < n_buf ; i++ )
01315                 {
01316                     half_buffer -> data[i] = in_buffer -> data[i] ;
01317                 }
01318                 edge_ind = 1 ;
01319                 init1 += margin ;
01320                 destroyVector ( in_buffer ) ;
01321             }
01322 
01323             /* determine the initial positions from the found values */
01324             if ( init1 != -1 )
01325             {
01326                 fitpar[2] = (float)init1 ;
01327             }
01328             fitpar[3] = 1. ;
01329 
01330             for ( i = 0 ; i < NPAR ; i++ )
01331             {
01332                 mpar[i] = 1 ;
01333                 dervpar[i] = 0. ;
01334             }
01335       
01336             xdim     = XDIMA ;
01337             ndat     = nel ;
01338             numpar   = NPAR ;
01339             tol      = TOLA ;
01340             lab      = LABA ;
01341             its      = ITSA ;
01342          
01343             /* finally, do the least squares fit over the buffer data */
01344             if ( 0 > ( iters = lsqfit( xdat, &xdim, half_buffer -> data, wdat, &ndat, fitpar,
01345                                        dervpar, mpar, &numpar, &tol, &its, &lab )) )
01346             { 
01347                 /* if the fit doesn't succeed the initial values are taken */
01348                 cpl_msg_warning ("lsqfit:"," least squares fit failed, error no.: %d in slitlet: %d\n", iters, j) ;
01349                 fitpar[2] = (float)init1 ;
01350             }
01351 
01352             pos = fitpar[2] ;
01353             if ( edge_ind == 1 )
01354             {
01355                 pos -= (float)margin ;
01356             }
01357 
01358             /*------------------------------------------------------------------------------------------ 
01359              * now discern the left and the right edge fit of the slitlets and associate the fit results
01360              * with the absolute positions in the image 
01361              * consider the difference of the fitted slit position to the expected position
01362              * and decide wether the fit is taken or the expected value is taken.
01363              */
01364             if ( left_right == 0 )
01365             {
01366                 /* take care of the column position of the fit boxes to get the absolute positions */
01367                 if ( j == 0 )
01368                 {
01369                     if ( fabs(pos - ((float)edgeclean[0] - 1. - (float)slit_length)) < diff_tol )
01370                     {
01371                         slit_pos[0][0] = pos ;
01372                     }
01373                     else
01374                     {
01375                         cpl_msg_warning("lsqfit:","something wrong with fitted left position of slitlet 0\n") ;
01376                         if ( (float) edgeclean[0] - 1. - (float)slit_length < 0. )
01377                         {
01378                             slit_pos[0][0] = 0. ;
01379                         }
01380                         else
01381                         {
01382                             slit_pos[0][0] = (float)edgeclean[0] - 1. - (float)slit_length ;
01383                         }
01384                     }
01385                 }
01386                 else if ( j < k )
01387                 {
01388                     if ( fabs( pos - (float)margin ) < diff_tol )
01389                     {
01390                         slit_pos[j][0] = pos + (float)edgeclean[j-1] - (float)margin ;
01391                     }
01392                     else
01393                     {
01394                         cpl_msg_warning("lsqfit:","something wrong with fitted left position of slitlet %d\n", j) ;
01395                         slit_pos[j][0] = (float)edgeclean[j-1] - 1. ;
01396                     }
01397                 }
01398                 else
01399                 {
01400                     if ( fabs( pos - (float)margin ) < diff_tol )
01401                     {
01402                         slit_pos[k][0] = pos + (float)edgeclean[k-1] - (float)margin ;
01403                     }
01404                     else
01405                     {
01406                         cpl_msg_warning("lsqfit:","something wrong with fitted left position of slitlet %d\n", j) ;
01407                         slit_pos[k][0] = (float)edgeclean[k-1] - 1. ;
01408                     }
01409                 }
01410             }
01411             else
01412             {
01413                 /* take care of the column position of the fit boxes to get the absolute positions */
01414                 if ( j == 0 )
01415                 {
01416                     if ( fabs( (float)box_buffer->n_elements - pos - (float)edgeclean[0] ) < diff_tol )
01417                     {
01418                         slit_pos[0][1] = (float)(box_buffer->n_elements - 1) - pos ;
01419                     }
01420                     else
01421                     {
01422                         cpl_msg_warning("lsqfit:","something wrong with fitted right position of slitlet 0\n") ;
01423                         slit_pos[0][1] = (float)edgeclean[0] - 1. ;
01424                     }
01425                 }
01426                 else if ( j < k )
01427                 {
01428                     if ( fabs( (float)box_buffer->n_elements - pos
01429                              + (float)edgeclean[j-1] - (float)margin - (float)edgeclean[j] ) < diff_tol )
01430                     {
01431                         slit_pos[j][1] = (float)(box_buffer->n_elements - 1) - pos
01432                                        + (float)edgeclean[j-1] - (float)margin ;
01433                     }
01434                     else
01435                     {
01436                         cpl_msg_warning("lsqfit:","something wrong with fitted right position of slitlet %d\n", j) ;
01437                         slit_pos[j][1] = (float)edgeclean[j] - 1. ;
01438                     }
01439                 }
01440                 else
01441                 {
01442                     if ( edgeclean[k-1] + slit_length > lineImage -> lx )
01443                     {
01444                         last_pos = (float)(lineImage -> lx - 1) ;
01445                     }
01446                     else
01447                     {
01448                         last_pos = (float)(edgeclean[k-1] - 1 + slit_length) ;
01449                     }
01450                     if ( fabs( (float)(box_buffer->n_elements - 1) - pos
01451                              + (float)edgeclean[k-1] - (float)margin - last_pos ) < diff_tol )
01452                     {
01453                         slit_pos[k][1] = (float)(box_buffer->n_elements - 1) - pos
01454                                        + (float)edgeclean[k-1] - (float)margin ;
01455                     }
01456                     else
01457                     {
01458                         cpl_msg_warning("lsqfit:","something wrong with fitted right position of slitlet %d\n", j) ;
01459                         slit_pos[k][1] = last_pos ;
01460                     }
01461                 }
01462             }
01463 
01464             destroyVector ( half_buffer ) ;
01465             cpl_free( xdat ) ;
01466             cpl_free( wdat ) ;
01467             cpl_free( mpar ) ;
01468         }
01469         destroyVector ( box_buffer ) ;
01470     }
01471      
01472 
01473     cpl_free( edge ) ;
01474     cpl_free( pos_row ) ;
01475     cpl_free( edgeclean ) ;
01476     cpl_free( dummyedge ) ;
01477     cpl_free( pos_rowclean ) ;
01478     return 0 ;
01479 }
01480 
01481 /*----------------------------------------------------------------------------
01482    Function     :       fitSlitsBoltzSingleLine()
01483    In           :       lineImage:  emission line frame
01484                         slit_pos:   allocated dummy array for the slitlet positions [min32][2]
01485                         box_length: pixel length of the row box within the fit is done
01486                         y_box:      small box in spectral direction within the slitlet
01487                                     may lie.
01488                         low_pos, high_pos: pixel positions in spectral direction between which the line
01489                        should be located.
01490    Out          :       slit_pos:  beginning and end position of the slitlets to
01491                                    sub-pixel accuracy
01492                          0  if it worked,
01493                         -1  if it failed,
01494    Job          :       fits the beginning and end position of the slitlets
01495                         by using non-linear least square fitting of a Boltzmann function
01496                         fits a Boltzmann function to the slitlet edges exposed and indicated
01497                         by the brightest emission lines. The slitlet is searched within
01498             user given positions.
01499                         The least squares fit is done by using a box smaller than
01500                         the size of two slitlets
01501  ---------------------------------------------------------------------------*/
01502 
01503 int fitSlitsBoltzSingleLine ( OneImage   * lineImage, 
01504                               float     ** slit_pos,
01505                               int          box_length,
01506                               float        y_box,
01507                   int          low_pos,
01508                   int          high_pos )
01509 {
01510     int     position[lineImage->lx] ;
01511     int   * edge, * edgeclean ;
01512     int   * dummyedge ;
01513     int   * pos_row, * pos_rowclean ;
01514     Vector * box_buffer ;
01515     Vector * half_buffer ;
01516     Vector * in_buffer ;
01517     int   found_row ;
01518     int   row, col ;
01519     int   i, j, k, m, ed ;
01520     int   init1 ;
01521     int   nel, n_right, left_right ;
01522     int   n_buf, edge_ind, shift ;
01523     int   slit_length ;
01524     int   margin ;
01525     int   iters, xdim, ndat ;
01526     int   numpar, its ;
01527     int   * mpar ;
01528     float * xdat, * wdat ;
01529     float tol, lab ;
01530     float fitpar[NPAR] ;
01531     float dervpar[NPAR] ;
01532     float minval, maxval ;
01533     float min ;
01534     float pos, last_pos ;
01535 
01536     slit_length = SLITLENGTH ;
01537 
01538     if ( NullImage == lineImage )
01539     {
01540         cpl_msg_error( "fitSlitsBoltzSingleLine:"," no line image given!\n" ) ;
01541         return -1 ;
01542     }
01543 
01544     if ( NULL == slit_pos )
01545     {
01546         cpl_msg_error( "fitSlitsBoltzSingleLine:"," no position array allocated!\n" ) ;
01547         return -1 ;
01548     }
01549 
01550     if ( box_length <  4 ||
01551          box_length >  2*slit_length )
01552     {
01553         cpl_msg_error( "fitSlitsBoltzSingleLine:"," wrong fitting box length given!\n" ) ;
01554         return -1 ;
01555     }
01556 
01557     if ( y_box <= 0. || y_box > 6. )
01558     {
01559         cpl_msg_error( "fitSlitsBoltzSingleLine:"," wrong y box length given!\n" ) ;
01560         return -1 ;
01561     }
01562 
01563     if ( low_pos >= high_pos || low_pos < 0 || high_pos <= 0 || high_pos >= lineImage->lx )
01564     {
01565         cpl_msg_error( "fitSlitsBoltzSingleLine:"," wrong user given search positions!\n" ) ;
01566         return -1 ;
01567     }
01568 
01569     /* allocate memory for the edges and the row position of the slitlets */
01570     edge         = (int*) cpl_calloc( lineImage->lx/2, sizeof(int) ) ;
01571     dummyedge    = (int*) cpl_calloc( lineImage->lx/2, sizeof(int) ) ;
01572     edgeclean    = (int*) cpl_calloc( lineImage->lx/2, sizeof(int) ) ;
01573     pos_row      = (int*) cpl_calloc( lineImage->lx/2, sizeof(int) ) ;
01574     pos_rowclean = (int*) cpl_calloc( lineImage->lx/2, sizeof(int) ) ;
01575 
01576     /* now search for the maximum between the given positions for each col */
01577     for ( col = 0 ; col < lineImage -> lx ; col++ )
01578     {
01579         found_row = -1 ;
01580         maxval = -FLT_MAX ;
01581     for ( row = low_pos ; row <= high_pos ; row++ )
01582     {
01583         if ( maxval < lineImage->data[col+row*lineImage->lx] )
01584             {
01585         maxval = lineImage->data[col+row*lineImage->lx] ;
01586                 found_row = row ;
01587         }
01588     }
01589     if ( maxval > -FLT_MAX && found_row > low_pos )
01590     {
01591             position[col] = found_row ;
01592         }
01593     else
01594     {
01595         position[col] = 0 ;
01596         }
01597     }
01598 
01599     /* now find and store the raw edge positions of the found slitlet */ 
01600     ed = 0 ;
01601     for ( col = 0 ; col < (lineImage->lx) - 1 ; col++ )
01602     {
01603         if ( position[col] > 0 && position[col+1] > 0 &&
01604          abs(position[col+1] - position[col]) > 10 ) 
01605         {
01606             edge[ed] = col ; 
01607             pos_row[ed] = position[col] ;
01608             ed++ ;
01609         }
01610 
01611     }
01612     if (ed <= 1)
01613     {
01614         cpl_msg_error( "fitSlitsBoltzSingleLine:"," no slitlets found!\n" ) ;
01615         cpl_free( edge ) ;
01616         cpl_free( pos_row ) ;
01617         cpl_free( edgeclean ) ;
01618         cpl_free( dummyedge ) ;
01619         cpl_free( pos_rowclean ) ;
01620         return -1 ;
01621     }
01622 
01623     /* now find the clean edge and row positions of the slitlets */
01624     for ( i = 1 ; i <= ed ; i ++ )
01625     {
01626         if ( i == ed )
01627         {
01628             if ( (edge[i-1] - edge[i-2]) < slit_length - SLIT_LEN_ERR ||
01629                  (edge[i-1] - edge[i-2]) > slit_length + SLIT_LEN_ERR )
01630             {
01631                 dummyedge[i-1]   = -1 ;
01632             }
01633         }
01634         if (dummyedge[i-1] != -1)
01635         {
01636             dummyedge[i-1] = edge[i-1] ;
01637         }
01638         else
01639         {
01640             continue ;
01641         }
01642         if ( i < ed )
01643         {
01644             if ( (edge[i] - edge[i-1]) < slit_length - SLIT_LEN_ERR ||
01645                  (edge[i] - edge[i-1]) > slit_length + SLIT_LEN_ERR )
01646             {
01647                 dummyedge[i]   = -1 ;
01648             }
01649         }
01650         if ( i+1 < ed && dummyedge[i] != -1 )
01651         {
01652             if ( (edge[i+1] - edge[i]) < slit_length - SLIT_LEN_ERR ||
01653                  (edge[i+1] - edge[i]) > slit_length + SLIT_LEN_ERR )
01654             {
01655                 dummyedge[i+1] = -1 ; 
01656             }
01657         }
01658     }
01659     
01660     k = 0 ;
01661     for ( i = 0 ; i < ed ; i++ )
01662     {
01663         if ( dummyedge[i] != -1 && dummyedge[i] != 0 )
01664         {
01665             edgeclean[k] = dummyedge[i] ;
01666             pos_rowclean[k] = pos_row[i] ;
01667             k++ ;
01668             if( edgeclean[k-1] > (lineImage -> lx - slit_length - 2*SLIT_LEN_ERR ) )
01669             {
01670                 pos_rowclean[k] = pos_row[ed] ;
01671             }
01672         }
01673     }
01674 
01675     /* determine the margins of the fitting box outside the slitlets */
01676     margin = box_length / 2 ;
01677 
01678     /* --------------------------------------------------------------------------
01679      * now go through the slitlets, search along each column within a box with 
01680      * half width y_box the maximum value and store these found values in a buffer
01681      */
01682     for ( j = 0 ; j <= k ; j++ )
01683     {
01684         m = 0 ;
01685         if ( j == 0 )
01686         {
01687             box_buffer = newVector( edgeclean[0] + margin ) ;
01688             for( col = 0 ; col < edgeclean[0] + margin ; col++ )
01689             {
01690                 maxval = -FLT_MAX ;
01691                 for ( row = pos_rowclean[0] - nint(y_box) ; row <= pos_rowclean[0] + nint(y_box) ; row++ )
01692                 {
01693                     if ( maxval < lineImage->data[col + lineImage->lx*row] )
01694                     {
01695                         maxval = lineImage -> data[col + lineImage->lx*row] ;
01696                     }
01697                 }
01698                 box_buffer->data[m] = maxval ;
01699                 m++ ;
01700             }
01701         }
01702         else if ( j < k )
01703         {
01704             box_buffer = newVector( edgeclean[j] - edgeclean[j-1] + 2*margin ) ;
01705             for ( col = edgeclean[j - 1] - margin ; col < edgeclean[j] + margin ; col++ )
01706             {
01707                 maxval = -FLT_MAX ;
01708                 for ( row = pos_rowclean[j] - nint(y_box) ; row <= pos_rowclean[j] + nint(y_box) ; row++ )
01709                 {
01710                     if ( maxval < lineImage->data[col + lineImage->lx*row] )
01711                     {
01712                         maxval = lineImage -> data[col + lineImage->lx*row] ;
01713                     }
01714                 }
01715                 box_buffer->data[m] = maxval ;
01716                 m++ ;
01717             }
01718         }
01719         else 
01720         {
01721             box_buffer = newVector( lineImage -> lx - edgeclean[k-1] + margin ) ;
01722             for ( col = edgeclean[k - 1] - margin ; col < lineImage -> lx ; col++ )
01723             {
01724         if ( col < 0 )
01725         {
01726             col = 0 ;
01727                 }
01728 
01729                 maxval = -FLT_MAX ;
01730                 for ( row = pos_rowclean[k] - nint(y_box) ; row <= pos_rowclean[k] + nint(y_box) ; row++ )
01731                 {
01732             if ( row < 0 )
01733                     {
01734             continue ;
01735                     }
01736                     if ( maxval < lineImage->data[col + row * lineImage->lx] )
01737                     {
01738                         maxval = lineImage -> data[col + row * lineImage->lx] ;
01739                     }
01740                 }
01741                 box_buffer->data[m] = maxval ;
01742                 m++ ;
01743             }
01744         }
01745 
01746         /* determine the minimum value in the box to get background1 value for the edge slitlets */
01747         min = FLT_MAX ;
01748         for ( i = 0 ; i < box_buffer->n_elements ; i++ )
01749         {
01750             if ( box_buffer -> data[i] < min )
01751             {
01752                 min = box_buffer -> data[i] ;
01753             }
01754         }
01755 
01756         for ( left_right = 0 ; left_right <= 1 ; left_right++ )
01757         { 
01758             nel = 0 ;
01759             if ( left_right == 0 )
01760             {
01761                 nel = box_buffer -> n_elements / 2 ;
01762             }
01763             else
01764             {
01765                 if ( box_buffer -> n_elements % 2 == 0 )
01766                 {
01767                     nel = box_buffer -> n_elements / 2 ;
01768                 }
01769                 else
01770                 {
01771                     nel = box_buffer -> n_elements / 2 + 1 ;
01772                 }
01773             }
01774 
01775             /* now split the buffer in the midth in a left and right part for fitting */
01776             half_buffer = newVector( nel ) ;
01777             if ( left_right == 0 )
01778             {
01779                 for ( i = 0 ; i < nel ; i++ )
01780                 {
01781                     half_buffer -> data[i] = box_buffer -> data[i] ;
01782                 }
01783             }
01784             else
01785             {
01786                 n_right = 0 ;
01787                 for ( i = box_buffer -> n_elements - 1 ; i >= box_buffer -> n_elements - nel ; i-- )
01788                 {
01789                     half_buffer -> data[n_right] = box_buffer -> data[i] ;
01790                     n_right++ ;
01791                 }
01792             }
01793 
01794             xdat = (float *) cpl_calloc( nel, sizeof (float) ) ;
01795             wdat = (float *) cpl_calloc( nel, sizeof (float) ) ;
01796             mpar = (int *)   cpl_calloc( NPAR, sizeof (int) ) ;
01797 
01798             /* set initial values for the fitting routine */
01799             minval =  FLT_MAX ;
01800             maxval = -FLT_MAX ;
01801             for ( i = 0 ; i < nel ; i++ )
01802             {
01803                 xdat[i] = i ;
01804                 wdat[i] = 1.0 ;
01805                 if ( half_buffer -> data[i] < minval )
01806                 {
01807                     minval = half_buffer -> data[i] ;
01808                 }
01809                 if ( half_buffer -> data[i] > maxval )
01810                 {
01811                     maxval = half_buffer -> data[i] ;
01812                 }
01813             }
01814             fitpar[0] = minval ;
01815             fitpar[1] = maxval ; 
01816 
01817             /* search for both positions of the half intensity of the hat within the buffer */
01818             init1 = -1 ; 
01819             for ( i = 0 ; i < nel ; i++ )
01820             {
01821                 if ( half_buffer -> data[i] >= ( maxval + minval ) / 2. )
01822                 {
01823                     init1 = i ;
01824                     break ;
01825                 }
01826             }
01827 
01828             /*------------------------------------------------------------------
01829              * if we have too few left background values (at the image edges)
01830              * the left margin of the buffer to fit is filled with the minimal
01831              * values in order to get a good fit
01832              */
01833               
01834             edge_ind = 0 ;
01835             if ( init1 < 3 )
01836             {
01837                 n_buf = half_buffer->n_elements + margin ;
01838                 in_buffer = newVector( n_buf ) ;     
01839                 for ( i = 0 ; i < margin ; i++ )
01840                 {
01841                     in_buffer -> data[i] = min ;
01842                 }
01843                 shift = 0 ;
01844                 for ( i = margin ; i < n_buf ; i++ )
01845                 {
01846                     in_buffer -> data[i] = half_buffer -> data[shift] ;
01847                     shift++ ;
01848                 }
01849                 destroyVector ( half_buffer ) ;
01850                 half_buffer = newVector ( n_buf ) ;
01851                 for ( i = 0 ; i < n_buf ; i++ )
01852                 {
01853                     half_buffer -> data[i] = in_buffer -> data[i] ;
01854                 }
01855                 edge_ind = 1 ;
01856                 init1 += margin ;
01857                 destroyVector ( in_buffer ) ;
01858             }
01859 
01860             /* determine the initial positions from the found values */
01861             if ( init1 != -1 )
01862             {
01863                 fitpar[2] = (float)init1 ;
01864             }
01865             fitpar[3] = 1. ;
01866 
01867             for ( i = 0 ; i < NPAR ; i++ )
01868             {
01869                 mpar[i] = 1 ;
01870                 dervpar[i] = 0. ;
01871             }
01872       
01873             xdim     = XDIMA ;
01874             ndat     = nel ;
01875             numpar   = NPAR ;
01876             tol      = TOLA ;
01877             lab      = LABA ;
01878             its      = ITSA ;
01879          
01880             /* finally, do the least squares fit over the buffer data */
01881             if ( 0 > ( iters = lsqfit( xdat, &xdim, half_buffer -> data, wdat, &ndat, fitpar,
01882                                        dervpar, mpar, &numpar, &tol, &its, &lab )) )
01883             { 
01884                 cpl_msg_warning ("lsqfit:"," least squares fit failed, error no.: %d in slitlet: %d\n", iters, j) ;
01885                 fitpar[2] = 0. ;
01886             }
01887         if ( fitpar[3] <=0. )
01888             { 
01889                 cpl_msg_warning ("fitSlitsBoltzSingleLine:"," fit failed due to negative width of boltzmann function in slitlet: %d\n", j) ;
01890                 fitpar[2] = 0. ;
01891             }
01892 
01893             pos = fitpar[2] ;
01894             if ( edge_ind == 1 )
01895             {
01896                 pos -= (float)margin ;
01897             }
01898 
01899             /*------------------------------------------------------------------------------------------ 
01900              * now discern the left and the right edge fit of the slitlets and associate the fit results
01901              * with the absolute positions in the image 
01902              * consider the difference of the fitted slit position to the expected position
01903              * and decide wether the fit is taken or the expected value is taken.
01904              */
01905             if ( left_right == 0 )
01906             {
01907                 /* take care of the column position of the fit boxes to get the absolute positions */
01908                 if ( j == 0 )
01909                 {
01910                     slit_pos[0][0] = pos ;
01911             if ( slit_pos[0][0] - (int) slit_pos[0][0] == 0.)
01912             {
01913                         slit_pos[0][0] = 0. ;
01914                     }
01915                 }
01916                 else if ( j < k )
01917                 {
01918                     slit_pos[j][0] = pos + (float)edgeclean[j-1] - (float)margin ;
01919             if ( slit_pos[j][0] - (int) slit_pos[j][0] == 0.)
01920             {
01921                         slit_pos[j][0] = 0. ;
01922                     }
01923                 }
01924                 else
01925                 {
01926                     slit_pos[k][0] = pos + (float)edgeclean[k-1] - (float)margin ;
01927             if ( slit_pos[k][0] - (int) slit_pos[k][0] == 0.)
01928             {
01929                         slit_pos[k][0] = 0. ;
01930                     }
01931                 }
01932             }
01933             else
01934             {
01935                 /* take care of the column position of the fit boxes to get the absolute positions */
01936                 if ( j == 0 )
01937                 {
01938                     slit_pos[0][1] = (float)(box_buffer->n_elements - 1) - pos ;
01939             if ( slit_pos[0][1] - (int) slit_pos[0][1] == 0.)
01940             {
01941                         slit_pos[0][1] = 0. ;
01942                     }
01943                 }
01944                 else if ( j < k )
01945                 {
01946                     slit_pos[j][1] = (float)(box_buffer->n_elements - 1) - pos
01947                                      + (float)edgeclean[j-1] - (float)margin ;
01948             if ( slit_pos[j][1] - (int) slit_pos[j][1] == 0.)
01949             {
01950                         slit_pos[j][1] = 0. ;
01951                     }
01952                 }
01953                 else
01954                 {
01955                     last_pos = (float)(edgeclean[k-1] - 1 + slit_length) ;
01956                     slit_pos[k][1] = (float)(box_buffer->n_elements - 1) - pos
01957                                      + (float)edgeclean[k-1] - (float)margin ;
01958             if ( slit_pos[k][1] - (int) slit_pos[k][1] == 0.)
01959             {
01960                         slit_pos[k][1] = 0. ;
01961                     }
01962                 }
01963             }
01964 
01965             destroyVector ( half_buffer ) ;
01966             cpl_free( xdat ) ;
01967             cpl_free( wdat ) ;
01968             cpl_free( mpar ) ;
01969         }
01970         destroyVector ( box_buffer ) ;
01971     }
01972         
01973     cpl_free( edge ) ;
01974     cpl_free( pos_row ) ;
01975     cpl_free( edgeclean ) ;
01976     cpl_free( dummyedge ) ;
01977     cpl_free( pos_rowclean ) ;
01978     return 0 ;
01979 }
01980 
01981 /*----------------------------------------------------------------------------
01982    Function     :       fitSlitsBoltzWithEstimate()
01983    In           :       lineImage:  emission line frame
01984                         slit_pos:   estimation array for the slitlet positions [min32][2]
01985                         box_length: pixel length of the row box within the fit is done
01986                         y_box:      small box in spectral direction within the slitlet
01987                                     may lie.
01988                         low_pos, high_pos: pixel positions in spectral direction between which the line
01989                        should be located.
01990    Out          :       slit_pos:  beginning and end position of the slitlets to
01991                                    sub-pixel accuracy
01992                          0  if it worked,
01993                         -1  if it failed,
01994    Job          :       fits the beginning and end position of the slitlets
01995                         by using non-linear least square fitting of a Boltzmann function
01996                         fits a Boltzmann function to the slitlet edges exposed and indicated
01997                         by the brightest emission lines. The slitlet is searched within
01998             user given positions.
01999                         The least squares fit is done by using a box smaller than
02000                         the size of two slitlets 
02001  ---------------------------------------------------------------------------*/
02002 
02003 int fitSlitsBoltzWithEstimate ( OneImage   * lineImage, 
02004                                 float     ** slit_pos,
02005                                 int          box_length,
02006                                 float        y_box,
02007                                 float        diff_tol,
02008                     int          low_pos,
02009                     int          high_pos )
02010 {
02011     int     position[lineImage->lx] ;
02012     Vector * box_buffer ;
02013     Vector * in_buffer ;
02014     int   found_row ;
02015     int   row, col ;
02016     int   col_first, col_last ;
02017     int   row_first, row_last ;
02018     int   i, j, m, n ;
02019     int   init1 ;
02020     int   left_right ;
02021     int   n_buf, shift ;
02022     int   slit_length ;
02023     int   iters, xdim, ndat ;
02024     int   numpar, its ;
02025     int   * mpar ;
02026     float * xdat, * wdat ;
02027     float tol, lab ;
02028     float fitpar[NPAR] ;
02029     float dervpar[NPAR] ;
02030     float minval, maxval ;
02031     float pos ;
02032     float new_pos ;
02033     int   slitposition[SLITLENGTH] ;
02034     pixelvalue rowpos[SLITLENGTH] ;
02035 
02036     slit_length = SLITLENGTH ; /* this setting is too much 64 */
02037     slit_length = N_SLITLETS ; /* this setting is better: 32 */
02038 
02039     if ( NullImage == lineImage )
02040     {
02041         cpl_msg_error( "fitSlitsBoltzWithEstimate:"," no line image given!\n" ) ;
02042         return -1 ;
02043     }
02044 
02045     if ( NULL == slit_pos )
02046     {
02047         cpl_msg_error( "fitSlitsBoltzWithEstimate:"," no position array allocated!\n" ) ;
02048         return -1 ;
02049     }
02050 
02051     if ( box_length < 4 ||
02052          box_length > 2*slit_length )
02053     {
02054         cpl_msg_error( "fitSlitsBoltzWithEstimate:"," wrong fitting box length given!\n" ) ;
02055         return -1 ;
02056     }
02057 
02058     if ( y_box <= 0. || y_box > 6. )
02059     {
02060         cpl_msg_error( "fitSlitsBoltzWithEstimate:"," wrong y box length given!\n" ) ;
02061         return -1 ;
02062     }
02063     if ( diff_tol <= 0. )
02064     {
02065         cpl_msg_error( "fitSlitsBoltzWithEstimate:"," wrong diff_tol given!\n" ) ;
02066         return -1 ;
02067     }
02068 
02069     if ( low_pos >= high_pos || low_pos < 0 || high_pos <= 0 || high_pos > lineImage->ly )
02070     {
02071         cpl_msg_error( "fitSlitsBoltzWithEstimate:"," wrong user given search positions!\n" ) ;
02072         return -1 ;
02073     }
02074 
02075     /* now search for the maximum between the given positions for each col */
02076     for ( col = 0 ; col < lineImage -> lx ; col++ )
02077     {
02078         found_row = -1 ;
02079         maxval = -FLT_MAX ;
02080     for ( row = low_pos ; row <= high_pos ; row++ )
02081     {
02082         if ( maxval < lineImage->data[col+row*lineImage->lx] )
02083             {
02084         maxval = lineImage->data[col+row*lineImage->lx] ;
02085                 found_row = row ;
02086         }
02087     }
02088     if ( maxval > -FLT_MAX && found_row > low_pos )
02089     {
02090             position[col] = found_row ;
02091         }
02092     else
02093     {
02094         position[col] = 0 ;
02095         }
02096     }
02097     
02098     /* --------------------------------------------------------------------------
02099      * now go through the slitlets, search along each column within a box with 
02100      * half width y_box the maximum value and store these found values in a buffer
02101      */
02102     for ( j = 0 ; j < slit_length ; j++ )
02103     {
02104         /* now go through the columns and determine the slitlet positions by
02105          * calculating the median of the found positions 
02106          */
02107         n = 0 ;
02108         
02109         for ( col = nint(slit_pos[j][0])+ 1 ; col < nint(slit_pos[j][1]) -1 ; col++ )
02110         {
02111             rowpos[n] = (pixelvalue)position[col] ;
02112             n++ ;
02113         }
02114 
02115         slitposition[j] = (int)median(rowpos, n) ;
02116         for ( left_right = 0 ; left_right <= 1 ; left_right++ )
02117         {
02118             init1 = 0 ;
02119             col_first = nint( slit_pos[j][left_right] ) - box_length/2 ;
02120             col_last  = nint( slit_pos[j][left_right] ) + box_length/2 ;
02121             if ( col_first < 0 )
02122             {
02123                 col_first = 0 ;
02124                 init1 = 1 ;
02125             }
02126             if ( col_last > lineImage->lx )
02127             {
02128                 col_last = lineImage->lx ;
02129                 init1 = 1 ;
02130             }
02131             if ( col_last - col_first <= 0 )
02132             {
02133                 cpl_msg_error( "fitSlitsBoltzWithEstimate:"," first and last column wrong!\n" ) ;
02134                 return -1 ;
02135             }
02136             box_buffer = newVector( col_last - col_first ) ;
02137             m = 0 ;
02138 
02139 
02140             if ( left_right == 0 )
02141             {
02142                 for( col = col_first ; col < col_last ; col++ )
02143                 {
02144                     row_first = slitposition[j] - nint(y_box) ;
02145                     row_last  = slitposition[j] + nint(y_box) ;
02146                     if ( row_first < 0 )
02147                     {
02148                         row_first = 0 ;
02149                     }
02150                     if ( row_last >= lineImage->ly  )
02151                     {
02152                         row_last = lineImage->ly - 1 ;
02153                     }
02154                     maxval = -FLT_MAX ;
02155                     for ( row = row_first ; row <= row_last ; row++ )
02156                     {
02157                         if ( maxval < lineImage->data[col + lineImage->lx*row] )
02158                         {
02159                             maxval = lineImage -> data[col + lineImage->lx*row] ;
02160                         }
02161                     }
02162                     box_buffer->data[m] = maxval ;
02163                     m++ ;
02164                 }
02165 
02166             }
02167             else 
02168             {
02169 
02170                 for( col = col_last-1 ; col >= col_first ; col-- )
02171                 {
02172                     row_first = slitposition[j] - nint(y_box) ;
02173                     row_last  = slitposition[j] + nint(y_box) ;
02174                     if ( row_first < 0 )
02175                     {
02176                         row_first = 0 ;
02177                     }
02178                     if ( row_last >= lineImage->ly  )
02179                     {
02180                         row_last = lineImage->ly - 1 ;
02181                     }
02182                     maxval = -FLT_MAX ;
02183                     for ( row = row_first ; row <= row_last ; row++ )
02184                     {
02185                         if ( maxval < lineImage->data[col + lineImage->lx*row] )
02186                         {
02187                             maxval = lineImage -> data[col + lineImage->lx*row] ;
02188                         }
02189                     }
02190                     box_buffer->data[m] = maxval ;
02191                     m++ ;
02192                 }
02193 
02194             }
02195 
02196             xdat = (float *) cpl_calloc( box_buffer->n_elements, sizeof (float) ) ;
02197             wdat = (float *) cpl_calloc( box_buffer->n_elements, sizeof (float) ) ;
02198             mpar = (int *)   cpl_calloc( NPAR, sizeof (int) ) ;
02199 
02200             /* set initial values for the fitting routine */
02201             minval =  FLT_MAX ;
02202             maxval = -FLT_MAX ;
02203 
02204             for ( i = 0 ; i < box_buffer->n_elements ; i++ )
02205             {
02206                 xdat[i] = i ;
02207                 wdat[i] = 1.0 ;
02208                 if ( box_buffer -> data[i] < minval )
02209                 {
02210                     minval = box_buffer -> data[i] ;
02211                 }
02212                 if ( box_buffer -> data[i] > maxval )
02213                 {
02214                     maxval = box_buffer -> data[i] ;
02215                 }
02216             }
02217             fitpar[0] = minval ;
02218             fitpar[1] = maxval ; 
02219 
02220             /*------------------------------------------------------------------
02221              * if we have too few left background values (at the image edges)
02222              * the left margin of the buffer to fit is filled with the minimal
02223              * values in order to get a good fit
02224              */
02225             
02226   
02227             if ( init1 == 1 )
02228             {
02229                 n_buf = box_buffer->n_elements + box_length/2 ;
02230                 in_buffer = newVector( n_buf ) ;     
02231                 for ( i = 0 ; i < box_length/2 ; i++ )
02232                 {
02233                     in_buffer -> data[i] = minval ;
02234                 }
02235                 shift = 0 ;
02236                 for ( i = box_length/2 ; i < n_buf ; i++ )
02237                 {
02238                     in_buffer -> data[i] = box_buffer -> data[shift] ;
02239                     shift++ ;
02240                 }
02241                 destroyVector ( box_buffer ) ;
02242                 box_buffer = newVector ( n_buf ) ;
02243                 for ( i = 0 ; i < n_buf ; i++ )
02244                 {
02245                     box_buffer -> data[i] = in_buffer -> data[i] ;
02246                 }
02247                 destroyVector ( in_buffer ) ;
02248             }
02249 
02250             /* determine the initial positions from the found values */
02251             fitpar[2] = (float)box_buffer->n_elements/2.  ;
02252             fitpar[3] = 1. ;
02253 
02254             for ( i = 0 ; i < NPAR ; i++ )
02255             {
02256                 mpar[i] = 1 ;
02257                 dervpar[i] = 0. ;
02258             }
02259       
02260             xdim     = XDIMA ;
02261             ndat     = box_buffer->n_elements ;
02262             numpar   = NPAR ;
02263             tol      = TOLA ;
02264             lab      = LABA ;
02265             its      = ITSA ;
02266          
02267             /* finally, do the least squares fit over the buffer data */
02268             if ( 0 > ( iters = lsqfit( xdat, &xdim, box_buffer -> data, wdat, &ndat, fitpar,
02269                                        dervpar, mpar, &numpar, &tol, &its, &lab )) )
02270             { 
02271                 cpl_msg_warning ("lsqfit:"," least squares fit failed, error no.: %d in slitlet: %d\n", iters, j) ;
02272                 destroyVector(box_buffer) ;
02273                 cpl_free( xdat ) ;
02274                 cpl_free( wdat ) ;
02275                 cpl_free( mpar ) ;
02276                 continue ;
02277             }
02278 
02279         if ( fitpar[3] <=0. )
02280             { 
02281                 cpl_msg_warning ("fitSlitsBoltzWithEstimate:"," fit failed due to negative width of boltzmann function in slitlet: %d\n", j) ;
02282                 destroyVector(box_buffer) ;
02283                 cpl_free( xdat ) ;
02284                 cpl_free( wdat ) ;
02285                 cpl_free( mpar ) ;
02286                 continue ;
02287             }
02288             pos = fitpar[2] ;
02289             if ( init1 == 1 )
02290             {
02291                 pos -= (float)box_length/2. ;
02292             }
02293 
02294             /*-------------------------------------------------------------
02295              * now compute the real slit positions using the guess positions
02296              * if the fit did not work the guess positions are taken
02297              * the same is done if the deviations are too big.
02298              */
02299             if ( pos != 0. )  
02300             {
02301                 if ( left_right == 0 )
02302                 {
02303                     new_pos = (float)col_first + pos ;
02304                 }
02305                 else
02306                 {
02307                     new_pos = (float)col_last-1 - pos ;
02308                 }
02309                 if ( fabs(new_pos - slit_pos[j][left_right]) < diff_tol )
02310                 {
02311                     slit_pos[j][left_right] = new_pos ;
02312                 }
02313                 else
02314                 {
02315                     cpl_msg_warning ("fitSlitsBoltzWithEstimate:"," deviation bigger than tolerance, take the estimated slitlet position in slitlet: %d\n", j) ;
02316                 }
02317             }
02318 
02319             cpl_free( xdat ) ;
02320             cpl_free( wdat ) ;
02321             cpl_free( mpar ) ;
02322             destroyVector ( box_buffer ) ;
02323 
02324         }
02325     }
02326         
02327     return 0 ;
02328 }
02329                               
02330 /*--------------------------------------------------------------------------*/

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