absolute.c

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

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