focus.c

00001 /*******************************************************************************
00002 * E.S.O. - VLT project
00003 *
00004 *
00005 *
00006 * who       when      what
00007 * --------  --------  ----------------------------------------------
00008 * schreib  16/01/02  created
00009 */
00010 
00011 /************************************************************************
00012 *   NAME
00013 *        focus.c -
00014 *        routines to determine the focus position of the detector
00015 *
00016 *   SYNOPSIS
00017 *   #include "focus.h"
00018 *
00019 *   1) double gaussianEllipse ( double * xdat, double * parlist )
00020 *   2) void gaussianEllipse_deriv( double * xdat, double * parlist, double * dervs )
00021 *   3) static int inv_mat (void)
00022 *   4) static void get_mat ( double * xdat,
00023 *                            int    * xdim,
00024 *                            double * ydat,
00025 *                            double * wdat,
00026 *                            int    * ndat,
00027 *                            double * fpar,
00028 *                            double * epar,
00029 *                            int    * npar )
00030 *   5) static int get_vec ( double * xdat,
00031 *                           int    * xdim,
00032 *                           double * ydat,
00033 *                           double * wdat,
00034 *                           int    * ndat,
00035 *                           double * fpar,
00036 *                           double * epar,
00037 *                           int    * npar )
00038 *   6) int lsqfit ( double * xdat,
00039 *                   int    * xdim,
00040 *                   double * ydat,
00041 *                   double * wdat,
00042 *                   int    * ndat,
00043 *                   double * fpar,
00044 *                   double * epar,
00045 *                   int    * mpar,
00046 *                   int    * npar,
00047 *                   double * tol ,
00048 *                   int    * its ,
00049 *                   double * lab  )
00050 *   7) int fit2dGaussian( OneImage   * lineImage, 
00051 *                         double     * fit_par, 
00052 *                         double     * derv_par   
00053 *                         int        * mpar,
00054 *                         int          lleftx,
00055 *                         int          llefty,
00056 *                         int          halfbox_x,
00057 *                         int          halfbox_y, int* check )
00058 *   8) OneImage * plotGaussian ( OneImage   * image, 
00059 *                                double     * parlist )
00060 *   9) static int gauss2Ellipse ( double     * parlist )
00061 *  10) float determineConversionFactor ( OneCube * cube, 
00062 *                                        float     mag,
00063 *                                        float     exptime,
00064 *                                        int       lleftx,
00065 *                                        int       llefty,
00066 *                                        int       halfbox_x,
00067 *                                        int       halfbox_y, 
00068 *                                        int* check )
00069 *
00070 *   DESCRIPTION
00071 *   1) Compute the value of a 2d Gaussian function at a given point.
00072 *      The ellptical 2D Gaussian is:
00073 *      F(x,y) = par(2) * EXP( -4.0*log(2.0)*[(xr/par(4))^2+(yr/par(5))^2]) + par(3),
00074 *      where: xr = xo * cos(par(6)) + yo * sin(par(6))
00075 *      yr = -xo * sin(par(6)) + yo * cos(par(6))
00076 *      and:   x0 = x - par(0)
00077 *             y0 = y - par(1)
00078 *   2) calculates the partial derivatives for a 2d Gaussian function with
00079 *      parameters parlist at position xdat 
00080 *   3) calculates the inverse of matrix2. The algorithm used 
00081 *      is the Gauss-Jordan algorithm described in Stoer,
00082 *      Numerische Mathematik, 1. Teil.
00083 *   4) builds the matrix 
00084 *   5) calculates the correction vector. The matrix has been
00085 *      built by get_mat(), we only have to rescale it for the 
00086 *      current value of labda. The matrix is rescaled so that
00087 *      the diagonal gets the value 1 + labda.
00088 *      Next we calculate the inverse of the matrix and then
00089 *      the correction vector.
00090 *   6) this is a routine for making a least-squares fit of a
00091 *      function to a set of data points. The method used is
00092 *      described in: Marquardt, J.Soc.Ind.Appl.Math. 11. 431 (1963).
00093 *      This method is a mixture of the steepest descent method 
00094 *      and the Taylor method.
00095 *   7) fits the image of a point source by using a 2-D Gaussian
00096 *      fit.
00097 *   8) plots an image of a given 2D-Gaussian 
00098 *   9) converts gauss parameters to ellipse parameters. 
00099 *  10) determines an intensity conversion factor for the instrument
00100 *      by fitting a 2D-Gaussian to an collapsed image of a standard star
00101 *      with known brightness (only for non-AO observations).
00102 *      Then the resulting Gaussian is integrated and the counts
00103 *      are divided by the exposure time (Fits header information) 
00104 *
00105 *   FILES
00106 *
00107 *   ENVIRONMENT
00108 *
00109 *   RETURN VALUES
00110 *
00111 *   CAUTIONS
00112 *
00113 *   EXAMPLES
00114 *
00115 *   SEE ALSO
00116 *
00117 *   BUGS
00118 *
00119 *------------------------------------------------------------------------
00120 */
00121 
00122 #include "vltPort.h"
00123 
00124 /*
00125  * System Headers
00126  */
00127 
00128 /*
00129  * Local Headers
00130  */
00131 
00132 #include "focus.h"
00133 #include "recipes.h"
00134 #include <float.h>
00135 
00136 /*----------------------------------------------------------------------------
00137  *                              Defines
00138  *--------------------------------------------------------------------------*/
00139 
00140 #define XDIMG          2         /* dimension of the x values */
00141 #define TOLG           0.001     /* fitting tolerance */
00142 #define LABG           0.1       /* labda parameter */
00143 #define ITSG           200       /* maximum number of iterations */
00144 #define LABFACG        10.0      /* labda step factor */
00145 #define LABMAXG        1.0e+10   /* maximum value for labda */
00146 #define LABMING        1.0e-10   /* minimum value for labda */
00147 #define NPAR           7         /* number of fit parameters */
00148 #define PI_NUMB        (3.1415926535897932384626433832795) /* pi */
00149 
00150 
00151 /*----------------------------------------------------------------------------
00152  *                                 Local variables
00153  *--------------------------------------------------------------------------*/
00154 
00155 static double chi1 ;                    /* old reduced chi-squared */
00156 static double chi2 ;                    /* new reduced chi-squared */
00157 static double labda ;                   /* mixing parameter */
00158 static double vec[NPAR] ;               /* correction vector */
00159 static double matrix1[NPAR][NPAR] ;     /* original matrix */
00160 static double matrix2[NPAR][NPAR] ;     /* inverse of matrix1 */
00161 static int    nfree ;                   /* number of free parameters */
00162 static int    parptr[NPAR] ;            /* parameter pointer */
00163 
00164 /*----------------------------------------------------------------------------
00165  *                  Functions private to this module
00166  *--------------------------------------------------------------------------*/
00167 
00168 static int inv_mat (void) ;
00169 
00170 static void get_mat ( double * xdat,
00171                       int    * xdim,
00172                       double * ydat,
00173                       double * wdat,
00174                       int    * ndat,
00175                       double * fpar,
00176                       double * epar/*,
00177                       int    * npar */) ;
00178 
00179 static int get_vec ( double * xdat,
00180                      int    * xdim,
00181                      double * ydat,
00182                      double * wdat,
00183                      int    * ndat,
00184                      double * fpar,
00185                      double * epar,
00186                      int    * npar ) ;
00187 
00188 static int gauss2Ellipse ( double  * parlist ) ;
00189 
00190 /*----------------------------------------------------------------------------
00191  *                          Function codes
00192  *--------------------------------------------------------------------------*/
00193 
00194 /*-------------------------------------------------------------------------*/
00215 /*--------------------------------------------------------------------------*/
00216 
00217 double gaussianEllipse(double * xdat, double * parlist)
00218 {
00219     double  result ;
00220     double x ;
00221     double y ;
00222     double fwhmx ;
00223     double fwhmy ;
00224     double costheta ;
00225     double sintheta ;
00226     double argX ;           /* arguments in the exponent */
00227     double argY ;
00228 
00229     /* some abbreviations */
00230     x =  xdat[0] -  parlist[0] ;
00231     y =  xdat[1] -  parlist[1] ;
00232 
00233     fwhmx = fabs(parlist[4]) ;
00234     fwhmy = fabs(parlist[5]) ;
00235     
00236     costheta = cos ( parlist[6] ) ;
00237     sintheta = sin ( parlist[6] ) ;
00238 
00239     argX = x * costheta + y * sintheta ;
00240     argY = -x * sintheta + y * costheta ;
00241     
00242     /* function */
00243     result =  parlist[2] * exp(-4.*log(2.0)*((argX/fwhmx)*(argX/fwhmx)+(argY/fwhmy)*(argY/fwhmy))) +
00244               parlist[3] ; 
00245 
00246     return result ;
00247 }
00248 
00249 /*----------------------------------------------------------------------------
00250                         gaussianEllipse_deriv()
00251    In           :       position array xdat, parameter list parlist 
00252                         The parameters are:
00253                         parlist(0): center of Gaussian in x direction
00254                         parlist(1): center of Gaussian in y direction
00255                         parlist(2): Amplitude of 2d Gaussian
00256                         parlist(3): Background level
00257                         parlist(4): FWHMx
00258                         parlist(5): FWHMy
00259                         parlist(6): theta
00260                         derivative value of a 2d Gaussian function at position xdat: dervs
00261                         dervs[0]: partial derivative by center x
00262                         dervs[1]: partial derivative by center y 
00263                         dervs[2]: partial derivative by the amplitude
00264                         dervs[3]: partial derivative by the background
00265                         dervs[4]: partial derivative by FWHMx 
00266                         dervs[5]: partial derivative by FWHMy
00267                         dervs[6]: partial derivative by theta
00268    Out          :       nothing, void
00269    Job          :       calculates the partial derivatives for a 2d Gaussian function with
00270                         parameters parlist at position xdat 
00271 ----------------------------------------------------------------------------*/
00272 
00273 void gaussianEllipse_deriv( double * xdat, double * parlist, double * dervs )
00274 {
00275     double x ;
00276     double y ;
00277     double fwhmx ;
00278     double fwhmy ;
00279     double argX ;
00280     double argY ;
00281     double expon ;
00282     double e8log2 ;
00283     double fwx2 ;
00284     double fwy2 ;
00285     double costheta ;
00286     double sintheta ;
00287 
00288     /* some abbreviations */
00289     x = xdat[0] - parlist[0] ;
00290     y = xdat[1] - parlist[1] ;
00291 
00292     fwhmx = fabs(parlist[4]) ;
00293     fwhmy = fabs(parlist[5]) ;
00294     fwx2 = fwhmx * fwhmx ;
00295     fwy2 = fwhmy * fwhmy ;
00296 
00297     costheta = cos ( parlist[6] ) ;
00298     sintheta = sin ( parlist[6] ) ;
00299 
00300     argX = x * costheta + y * sintheta ;
00301     argY = -x * sintheta + y * costheta ;
00302 
00303     expon = exp ( -4.0 * log(2.0) * ((argX/fwhmx)*(argX/fwhmx) + (argY/fwhmy)*(argY/fwhmy)) ) ;
00304     e8log2 = expon * 8.0 * log(2.0) ;
00305 
00306     /* determine the derivatives */
00307     /* partial derivative x-position */
00308     dervs[0] = -parlist[2]*e8log2 * (-argX*costheta/fwx2 + argY*sintheta/fwy2) ;
00309     /* partial derivative y-position */
00310     dervs[1] = -parlist[2]*e8log2 * (-argX*sintheta/fwx2 - argY*costheta/fwy2) ;
00311     /* partial derivative amplitude */
00312     dervs[2] = expon ;
00313     /* partial derivative background */
00314     dervs[3] = 1. ;
00315     /* partial derivative fwhmx */
00316     dervs[4] = parlist[2]*e8log2 * argX*argX/(fwx2*fwhmx) ;
00317     /* partial derivative fwhmy */
00318     dervs[5] = parlist[2]*e8log2 * argY*argY/(fwy2*fwhmy) ;
00319     /* partial derivative theta */
00320     dervs[6] = -parlist[2]*e8log2 * argY * argX * (1.0/fwx2 - 1.0/fwy2) ;
00321     
00322 }
00323 
00324 /*----------------------------------------------------------------------------
00325    Function     :       inv_mat()
00326    In           :       void
00327    Out          :       integer (0 if it worked, -6 if determinant is zero) 
00328    Job          :       calculates the inverse of matrix2. The algorithm used 
00329                         is the Gauss-Jordan algorithm described in Stoer,
00330                         Numerische Mathematik, 1. Teil.
00331  ---------------------------------------------------------------------------*/
00332 
00333 static int inv_mat (void)
00334 {
00335     double even ;
00336     double hv[NPAR] ;
00337     double mjk ;
00338     double rowmax ;
00339     int evin ;
00340     int i, j, k, row ;
00341     int per[NPAR] ;
00342    
00343     /* set permutation array */
00344     for ( i = 0 ; i < nfree ; i++ )
00345     {
00346         per[i] = i ;
00347     }
00348     
00349     for ( j = 0 ; j < nfree ; j++ ) /* in j-th column */
00350     {
00351         /* determine largest element of a row */                                   
00352         rowmax = fabs ( matrix2[j][j] ) ;
00353         row = j ;                         
00354 
00355         for ( i = j + 1 ; i < nfree ; i++ )
00356         {
00357             if ( fabs ( matrix2[i][j] ) > rowmax )
00358             {
00359                 rowmax = fabs( matrix2[i][j] ) ;
00360                 row = i ;
00361             }
00362         }
00363 
00364         /* determinant is zero! */
00365         if ( matrix2[row][j] == 0.0 )
00366         {
00367             return -6 ;
00368         }
00369         
00370         /* if the largest element is not on the diagonal, then permutate rows */
00371         if ( row > j )
00372         {
00373             for ( k = 0 ; k < nfree ; k++ )
00374             {
00375                 even = matrix2[j][k] ;
00376                 matrix2[j][k] = matrix2[row][k] ;
00377                 matrix2[row][k] = even ;
00378             }
00379             /* keep track of permutation */
00380             evin = per[j] ;
00381             per[j] = per[row] ;
00382             per[row] = evin ;
00383         }
00384         
00385         /* modify column */
00386         even = 1.0 / matrix2[j][j] ;
00387         for ( i = 0 ; i < nfree ; i++ )
00388         {
00389             matrix2[i][j] *= even ;
00390         }
00391         matrix2[j][j] = even ;
00392         
00393         for ( k = 0 ; k < j ; k++ )
00394         {
00395             mjk = matrix2[j][k] ;
00396             for ( i = 0 ; i < j ; i++ )
00397             {
00398                 matrix2[i][k] -= matrix2[i][j] * mjk ;
00399             }
00400             for ( i = j + 1 ; i < nfree ; i++ )
00401             {
00402                 matrix2[i][k] -= matrix2[i][j] * mjk ;
00403             }
00404             matrix2[j][k] = -even * mjk ;
00405         }
00406     
00407         for ( k = j + 1 ; k < nfree ; k++ )
00408         {
00409             mjk = matrix2[j][k] ;
00410             for ( i = 0 ; i < j ; i++ )
00411             {
00412                 matrix2[i][k]  -= matrix2[i][j] * mjk ;
00413             }
00414             for ( i = j + 1 ; i < nfree ; i++ )
00415             {
00416                 matrix2[i][k]  -= matrix2[i][j] * mjk ;
00417             }
00418             matrix2[j][k] = -even * mjk ;
00419         }
00420     }
00421     
00422     /* finally, repermute the columns */
00423     for ( i = 0 ; i < nfree ; i++ )
00424     {
00425         for ( k = 0 ; k < nfree ; k++ )
00426         {
00427             hv[per[k]] = matrix2[i][k] ;
00428         }
00429         for ( k = 0 ; k < nfree ; k++ )
00430         {
00431             matrix2[i][k] = hv[k] ;
00432         }
00433     }
00434         
00435     /* all is well */
00436     return 0 ;
00437 }
00438     
00439 /*----------------------------------------------------------------------------
00440    Function     :       get_mat()
00441    In           :       xdat: position
00442                         xdim: factor of the indices of the position array
00443                         ydat: real data
00444                         wdat: weights
00445                         ndat: number of data points
00446                         fpar: function parameters
00447                         epar: partial derivatives of the function
00448                         npar: number of function parameters
00449    Out          :       void
00450    Job          :       builds the matrix 
00451  ---------------------------------------------------------------------------*/
00452 
00453 static void get_mat ( double * xdat,
00454                       int    * xdim,
00455                       double * ydat,
00456                       double * wdat,
00457                       int    * ndat,
00458                       double * fpar,
00459                       double * epar/*,
00460                       int    * npar */)
00461 {
00462     double wd ;
00463     double wn ;
00464     double yd ;
00465     int i, j, n ;
00466 
00467     for ( j = 0 ; j < nfree ; j++ )
00468     {
00469         vec[j] = 0.0 ; /* zero vector */
00470         for ( i = 0 ; i<= j ; i++ )   /* zero matrix only on and below diagonal */
00471         {
00472             matrix1[j][i] = 0.0 ;
00473         }
00474     }
00475     chi2 = 0.0 ;  /* reset reduced chi-squared */
00476     
00477     /* loop through data points */
00478     for ( n = 0 ; n < (*ndat) ; n++ )
00479     {
00480         wn = wdat[n] ;
00481         if ( wn > 0.0 )  /* legal weight ? */
00482         {
00483             yd = ydat[n] - gaussianEllipse( &xdat[(*xdim) * n], fpar ) ;
00484             gaussianEllipse_deriv( &xdat[(*xdim) * n], fpar, epar ) ;
00485             chi2 += yd * yd * wn ; /* add to chi-squared */
00486             for ( j = 0 ; j < nfree ; j++ )
00487             {
00488                 wd = epar[parptr[j]] * wn ;  /* weighted derivative */
00489                 vec[j] += yd * wd ;       /* fill vector */
00490                 for ( i = 0 ; i <= j ; i++ ) /* fill matrix */
00491                 {
00492                     matrix1[j][i] += epar[parptr[i]] * wd ;
00493                 }
00494             }
00495         }
00496     }                   
00497 }  
00498                 
00499             
00500 /*----------------------------------------------------------------------------
00501    Function     :       get_vec()
00502    In           :       xdat: position
00503                         xdim: factor of the indices of the position array
00504                         ydat: real data
00505                         wdat: weights
00506                         ndat: number of data points
00507                         fpar: function parameters
00508                         epar: partial derivatives of the function
00509                         npar: number of function parameters
00510    Out          :       integer (0 if it had worked, -5 or -7 
00511                                  if diagonal element is wrong, or -6,
00512                                  if determinant is zero )
00513    Job          :       calculates the correction vector. The matrix has been
00514                         built by get_mat(), we only have to rescale it for the 
00515                         current value of labda. The matrix is rescaled so that
00516                         the diagonal gets the value 1 + labda.
00517                         Next we calculate the inverse of the matrix and then
00518                         the correction vector.
00519  ---------------------------------------------------------------------------*/
00520             
00521 static int get_vec ( double * xdat,
00522                      int    * xdim,
00523                      double * ydat,
00524                      double * wdat,
00525                      int    * ndat,
00526                      double * fpar,
00527                      double * epar,
00528                      int    * npar )
00529 {
00530     double dj ;
00531     double dy ;
00532     double mii ;
00533     double mji ;
00534     double mjj ;
00535     double wn ;
00536     int i, j, n, r ;
00537 
00538     /* loop to modify and scale the matrix */
00539     for ( j = 0 ; j < nfree ; j++ )
00540     {
00541         mjj = matrix1[j][j] ;
00542         if ( mjj <= 0.0 )             /* diagonal element wrong */
00543         {
00544             return -5 ;
00545         }
00546         mjj = sqrt( mjj ) ;
00547         for ( i = 0 ; i < j ; i++ )
00548         {
00549             mji = matrix1[j][i] / mjj / sqrt( matrix1[i][i] ) ;
00550             matrix2[i][j] = matrix2[j][i] = mji ;
00551         }
00552         matrix2[j][j] = 1.0 + labda ;  /* scaled value on diagonal */
00553     }    
00554     
00555     if ( (r = inv_mat()) )                /* invert matrix inlace */
00556     {
00557         return r ;
00558     }
00559     
00560     for ( i = 0 ; i < (*npar) ; i ++ )
00561     {
00562         epar[i] = fpar[i] ;
00563     }
00564     
00565     /* loop to calculate correction vector */
00566     for ( j = 0 ; j < nfree ; j++ )
00567     {
00568         dj = 0.0 ;
00569         mjj = matrix1[j][j] ;
00570         if ( mjj <= 0.0)               /* not allowed */
00571         {
00572             return -7 ;
00573         }
00574         mjj = sqrt ( mjj ) ;
00575         for ( i = 0 ; i < nfree ; i++ )
00576         {
00577             mii = matrix1[i][i] ;
00578             if ( mii <= 0.0 )
00579             {
00580                 return -7 ;
00581             }
00582             mii = sqrt( mii ) ;
00583             dj += vec[i] * matrix2[j][i] / mjj / mii ;
00584         }
00585         epar[parptr[j]] += dj ;       /* new parameters */
00586     }    
00587     chi1 = 0.0 ;                      /* reset reduced chi-squared */
00588  
00589     /* loop through the data points */
00590     for ( n = 0 ; n < (*ndat) ; n++ )
00591     {
00592         wn = wdat[n] ;               /* get weight */
00593         if ( wn > 0.0 )              /* legal weight */
00594         {
00595             dy = ydat[n] - gaussianEllipse( &xdat[(*xdim) * n], epar) ;
00596             chi1 += wdat[n] * dy * dy ;
00597         }
00598     }
00599     return 0 ;
00600 }   
00601     
00602         
00603 
00604 /*----------------------------------------------------------------------------
00605    Function     :       lsqfitd()
00606    In           :       xdat: position, coordinates of data points.
00607                               xdat is 2 dimensional: XDAT ( XDIMG, NDAT )
00608                         xdim: dimension of fit
00609                         ydat: data points
00610                         wdat: weights for data points
00611                         ndat: number of data points
00612                         fpar: on input contains initial estimates of the 
00613                               parameters for non-linear fits, on output the
00614                               fitted parameters.
00615                         epar: contains estimates of the errors in fitted parameters
00616                         mpar: logical mask telling which parameters are free (non-zero)
00617                               and which parameters are fixed (0)
00618                         npar: number of function parameters ( free + fixed )
00619                         tol:  relative tolerance. lsqfit stops when successive iterations
00620                               fail to produce a decrement in reduced chi-squared less
00621                               than tol. If tol is less than the minimum tolerance 
00622                               possible, tol will be set to this value. This means
00623                               that maximum accuracy can be obtained by setting
00624                               tol = 0.0.
00625                         its:  maximum number of iterations
00626                         lab:  mixing parameter, lab determines the initial weight
00627                               of steepest descent method relative to the Taylor method
00628                               lab should be a small value (i.e. 0.01). lab can only
00629                               be zero when the partial derivatives are independent
00630                               of the parameters. In fact in this case lab should be
00631                               exactly equal to zero.
00632    Out          :       returns number of iterations needed to achieve convergence
00633                         according to tol. When this number is negative, the fitting
00634                         was not continued because a fatal error occurred:
00635                         -1 too many free parameters, maximum is 32
00636                         -2 no free parameters
00637                         -3 not enough degrees of freedom
00638                         -4 maximum number of iterations too small to obtain
00639                            a solution which satisfies tol.
00640                         -5 diagonal of matrix contains elements which are zero
00641                         -6 determinant of the coefficient matrix is zero
00642                         -7 square root of a negative number 
00643    Job          :       this is a routine for making a least-squares fit of a
00644                         function to a set of data points. The method used is
00645                         described in: Marquardt, J.Soc.Ind.Appl.Math. 11. 431 (1963).
00646                         This method is a mixture of the steepest descent method 
00647                         and the Taylor method.
00648  ---------------------------------------------------------------------------*/
00649 
00650 int lsqfitd ( double * xdat,
00651               int    * xdim,
00652               double * ydat,
00653               double * wdat,
00654               int    * ndat,
00655               double * fpar,
00656               double * epar,
00657               int    * mpar,
00658               int    * npar,
00659               double * tol ,
00660               int    * its ,
00661               double * lab  )
00662 {
00663     int i, n, r ;
00664     int itc ;                      /* fate of fit */
00665     int found ;                    /* fit converged: 1, not yet converged: 0 */
00666     int  nuse ;                    /* number of useable data points */
00667     double tolerance ;             /* accuracy */
00668 
00669     itc   = 0 ;                    /* fate of fit */
00670     found = 0 ;                    /* reset */
00671     nfree = 0 ;                    /* number of free parameters */
00672     nuse  = 0 ;                    /* number of legal data points */
00673 
00674     if ( *tol < (DBL_EPSILON * 10.0 ) )
00675     {
00676         tolerance = DBL_EPSILON * 10.0 ;  /* default tolerance */
00677     }
00678     else
00679     {
00680         tolerance = *tol ;                /* tolerance */
00681     }
00682     
00683     labda = fabs( *lab ) * LABFACG ;       /* start value for mixing parameter */
00684     for ( i = 0 ; i < (*npar) ; i++ )
00685     {
00686         if ( mpar[i] )
00687         {
00688             if ( nfree > NPAR )         /* too many free parameters */
00689             {
00690                 return -1 ;
00691             }
00692             parptr[nfree++] = i ;         /* a free parameter */
00693         }
00694     }
00695     
00696     if (nfree == 0)                       /* no free parameters */     
00697     {
00698         return -2 ;
00699     }
00700     
00701     for ( n = 0 ; n < (*ndat) ; n++ )
00702     {
00703         if ( wdat[n] > 0.0 )              /* legal weight */
00704         {
00705             nuse ++ ;
00706         }
00707     }
00708     
00709     if ( nfree >= nuse )
00710     {
00711         return -3 ;                       /* no degrees of freedom */
00712     }
00713     if ( labda == 0.0 )                   /* linear fit */
00714     {
00715         for ( i = 0 ; i < nfree ; fpar[parptr[i++]] = 0.0 ) ;  /* initialize fpar array */
00716         get_mat ( xdat, xdim, ydat, wdat, ndat, fpar, epar/*, npar */) ;
00717         r =  get_vec ( xdat, xdim, ydat, wdat, ndat, fpar, epar, npar ) ;
00718         if ( r )                         /* error */
00719         {
00720             return r ;
00721         }
00722         for ( i = 0 ; i < (*npar) ; i++ )
00723         {
00724             fpar[i] = epar[i] ;           /* save new parameters */
00725             epar[i] = 0.0 ;               /* and set errors to zero */
00726         }
00727         chi1 = sqrt( chi1 / (double) (nuse - nfree) ) ;
00728         for ( i = 0 ; i < nfree ; i++ )
00729         {
00730             if ( (matrix1[i][i] <= 0.0 ) || (matrix2[i][i] <= 0.0) )
00731             {
00732                 return -7 ;
00733             }
00734             epar[parptr[i]] = chi1 * sqrt( matrix2[i][i] ) / sqrt( matrix1[i][i] ) ;
00735         }
00736     }
00737     else                                  /* non-linear fit */
00738     {
00739         /*----------------------------------------------------------------
00740          * the non-linear fit uses the steepest descent method in combination
00741          * with the Taylor method. The mixing of these methods is controlled
00742          * by labda. In the outer loop ( called the iteration loop ) we build
00743          * the matrix and calculate the correction vector. In the inner loop
00744          * (called the interpolation loop) we check whether we have obtained a
00745          * better solution than the previous one. If so, we leave the inner loop
00746          * else we increase labda ( give more weight to the steepest descent method)
00747          * calculate the correction vector and check again. After the inner loop
00748          * we do a final check on the goodness of the fit and if this satisfies
00749          * the tolerance we calculate the errors of the fitted parameters.
00750          */
00751         while ( !found )                  /* iteration loop */
00752         {      
00753             if ( itc++ == (*its) )        /* increase iteration counter */
00754             {
00755                 return -4 ;               
00756             }
00757             get_mat( xdat, xdim, ydat, wdat, ndat, fpar, epar/*, npar*/ ) ;
00758             
00759             /*-------------------------------------------------------------
00760              * here we decrease labda since we may assume that each iteration
00761              * brings us closer to the answer.
00762              */
00763             if ( labda > LABMING )
00764             {
00765                 labda = labda / LABFACG ;         /* decrease labda */
00766             }
00767             r = get_vec ( xdat, xdim, ydat, wdat, ndat, fpar, epar, npar ) ;
00768 
00769             if ( r )                      /* error */
00770             {
00771                 return r ;
00772             }
00773 
00774             while ( chi1 >= chi2 )        /* interpolation loop */
00775             {
00776                 /*-----------------------------------------------------------
00777                  * The next statement is based on experience, not on the mathematics
00778                  * of the problem. It is assumed that we have reached convergence 
00779                  * when the pure steepest descent method does not produce a better
00780                  * solution.
00781                  */
00782                 if ( labda > LABMAXG )    /* assume solution found */
00783                 {
00784                     break ;
00785                 }
00786                 labda = labda * LABFACG ;         /* increase mixing parameter */
00787                 r = get_vec ( xdat, xdim, ydat, wdat, ndat, fpar, epar, npar ) ;
00788 
00789                 if ( r )                  /* error */
00790                 {
00791                     return r ;
00792                 }
00793             }
00794 
00795             if ( labda <= LABMAXG )        /* save old parameters */
00796             {
00797                 for ( i = 0 ; i < *npar ; i++ )
00798                 {
00799                     fpar[i] = epar[i] ;
00800                 }
00801             }
00802             if ( (fabs( chi2 - chi1 ) <= (tolerance * chi1)) || (labda > LABMAXG) )
00803             {
00804                 /*---------------------------------------------------------------------
00805                  * we have a satisfying solution, so now we need to calculate the correct
00806                  * errors of the fitted parameters. This we do by using the pure Taylor
00807                  * method because we are very close to the real solution.
00808                  */
00809                 labda = LABMING ;              /* for Taylor solution */
00810                 get_mat ( xdat, xdim, ydat, wdat, ndat, fpar, epar/*, npar */) ;
00811                 r = get_vec ( xdat, xdim, ydat, wdat, ndat, fpar, epar, npar ) ;
00812 
00813                 if ( r )                    /* error */
00814                 {
00815                     return r ;
00816                 }
00817                 for ( i = 0 ; i < (*npar) ; i++ )
00818                 {
00819                     epar[i] = 0.0 ;          /* set error to zero */
00820                 }
00821                 chi2 = sqrt ( chi2 / (double) (nuse - nfree) ) ;
00822 
00823                 for ( i = 0 ; i < nfree ; i++ )
00824                 {
00825                     if ( (matrix1[i][i] <= 0.0) || (matrix2[i][i] <= 0.0) )
00826                     {
00827                         return -7 ;
00828                     }
00829                     epar[parptr[i]] = chi2 * sqrt( matrix2[i][i] ) / sqrt( matrix1[i][i] ) ;
00830                 }
00831                 found = 1 ;                  /* we found a solution */
00832             }
00833         }
00834     }
00835     return itc ;                             /* return number of iterations */
00836 }
00837 
00838 /*----------------------------------------------------------------------------
00839    Function     :       fit2dGaussian()
00840    In           :       image:      reconstructed image of a point source 
00841                         fit_par:    array of 7 free fit parameters of a 2D-Gauss
00842                         derv_par:   derivatives of the fit_par array
00843                         mpar:       mask to set free parameters, 1: free, 0: fixed. 
00844                         lleftx, llefty: lower left starting point of the box in which the fit is carried through
00845                                     to find the maximum of point source image
00846                         halfbox_x:  half box length in x-direction in pixels inside which the fit 
00847                                     is carried through
00848                         halfbox_y:  half box length in y-direction in pixels inside which the fit 
00849                                     is carried through
00850    Out          :       number of needed iterations
00851                         -1 if an error occurred.
00852    Job          :       fits the image of a point source by using a 2-D Gaussian
00853                         fit. 
00854                         Remark on the fit results for the fit parameters
00855                         (fwhmx and fwhmy and theta): theta will always be
00856                         between -PI/4 and +PI/4, exchange of the fwhmx and fwhmy
00857                         values corresponds to a shift of theta by PI/2.
00858                         Consequently, an expected theta > |PI/4| will result
00859                         in an exchange of the fwhm values and a shift of theta by
00860                         PI/2 to a value < |PI/4| what yields exactly the same image.
00861  ---------------------------------------------------------------------------*/
00862 
00863 int fit2dGaussian ( OneImage   * image, 
00864                     double     * fit_par, 
00865                     double     * derv_par,   
00866                     int        * mpar,
00867                     int          lleftx,
00868                     int          llefty,
00869                     int          halfbox_x,
00870                     int          halfbox_y, int* check )
00871 {
00872     int i, j, n ;
00873     int col, row ;
00874     int boxi, boxj ;
00875     int iters ;
00876     int ndata ;
00877     int xdim ;
00878     int npar ;
00879     int its ;
00880     double lab ;
00881     double tol ;
00882     double maxval ;
00883     double background ;
00884     double amplitude ;
00885     float * backarray ;
00886     double M, Mx, My ;
00887     double Mxx, Mxy, Myy ; 
00888     double X0, Y0 ;
00889     double xydat[4 *halfbox_x*halfbox_y][XDIMG] ;
00890     double zdat[4*halfbox_x*halfbox_y] ;
00891     double wdat[4*halfbox_x*halfbox_y] ;
00892     double xco, yco ;
00893     double value ;
00894     double denom ;
00895     double temp ;
00896     int llx, lly ;
00897     int foundrow ;
00898     int foundcol ;
00899     int k ;
00900 
00901     if ( NullImage == image )
00902     {
00903         cpl_msg_error("fit2dGaussian:"," no image given\n") ;
00904         return -1 ;
00905     }
00906     if ( NULL == fit_par )
00907     {
00908         cpl_msg_error("fit2dGaussian:"," no fit parameters given\n") ;
00909         return -1 ;
00910     }
00911     if ( NULL == derv_par )
00912     {
00913         cpl_msg_error("fit2dGaussian:"," no derivatives of fit parameters given\n") ;
00914         return -1 ;
00915     }
00916     if ( lleftx < 0 || lleftx + 2*halfbox_x >= image->lx ||
00917          llefty < 0 || llefty + 2*halfbox_y >= image->ly )
00918     {
00919         cpl_msg_error("fit2dGaussian:"," wrong lower left point of fitting box given! \n") ;
00920         return -1 ;
00921     }
00922     if ( halfbox_x <= 1 || halfbox_y <= 1 )
00923     {
00924         cpl_msg_error("fit2dGaussian:"," wrong box dimensions given\n") ;
00925         return -1 ;
00926     }
00927     /* allocate memory */
00928     if ( NULL == (backarray = (float*) cpl_calloc(4*halfbox_x+4*halfbox_y, sizeof(float) ) ) ) 
00929     {
00930         cpl_msg_error("fit2dGaussian:"," could not allocate memory\n") ;
00931         return -1 ;
00932     }
00933 
00934     /* -----------------------------------------------------------------------------
00935      * find the initial estimates for the free parameters
00936      */
00937 
00938     /* first search for the position of the maximum intensity */
00939     foundrow = 0 ;
00940     foundcol = 0 ;
00941     maxval   = -DBL_MAX ;
00942     for ( col = lleftx ; col < lleftx + 2*halfbox_x ; col++ )
00943     {
00944         for ( row = llefty ; row < llefty + 2*halfbox_y ; row++ )
00945         {
00946             if ( qfits_isnan(image->data[col+row*image->lx]) )
00947             {
00948                 continue ;
00949             }
00950             if ( maxval < image->data[col+row*image->lx] ) 
00951             {
00952                 maxval = image->data[col+row*image->lx] ;
00953                 foundrow = row ;
00954                 foundcol = col ;
00955             }
00956         }
00957     }
00958 
00959     if ( foundrow == 0 || foundcol == 0 || maxval <= 0. ||
00960          foundrow == image->lx-1 || foundcol == image->ly-1 )
00961     {
00962         cpl_msg_error("fit2dGaussian:"," no maximum found\n") ;
00963         cpl_free(backarray) ;
00964         return -1 ;
00965     }
00966 
00967     /* determine the lower left edge of the fitting box, center it 
00968        on the maximum value */
00969     llx = foundcol - halfbox_x ;
00970     lly = foundrow - halfbox_y ;
00971     if ((foundcol - halfbox_x) > 0) {
00972       llx = (foundcol - halfbox_x);
00973     } else {
00974       llx=1;
00975       check++;
00976     }
00977 
00978     if ((foundrow - halfbox_y) > 0) {
00979       lly = (foundrow - halfbox_y);
00980     } else {
00981       lly=1;
00982       check++;
00983     } 
00984 
00985     if ( ( llx + 2*halfbox_x) <  image->lx-1 ) {
00986       halfbox_x=halfbox_x;
00987     } else {
00988        halfbox_x=(int) (image->lx-2-llx)/2;
00989       check++;
00990     }
00991 
00992     if ( ( lly + 2*halfbox_y) <  image->ly-1 ) {
00993       halfbox_y= halfbox_y;
00994     } else {
00995       halfbox_y=(int) (image->ly-2-lly)/2;
00996       check++;
00997     }
00998 
00999     if ( llx <= 0 || lly < 0 || llx + 2*halfbox_x >= image->lx-1 ||
01000          lly + 2*halfbox_y >= image->ly )
01001     {
01002         cpl_msg_error("fit2dGaussian:"," box does not fit into image\n") ;
01003         cpl_free(backarray) ;
01004         return -1 ;
01005     }
01006     
01007     /* determine the zeroth and first order moments of the image  
01008        within the fitting box */ 
01009     M = Mx = My = 0. ;
01010     n = 0 ;
01011     boxi = boxj = 0 ;
01012     for ( j = lly ; j < lly + 2*halfbox_y ; j++ )
01013     {
01014         boxj = j - lly ;
01015         for ( i = llx ; i < llx + 2*halfbox_x ; i++ )
01016         {
01017             boxi = i - llx ;
01018             if ( !qfits_isnan(image->data[i+j*image->lx]) )
01019             {
01020                 M  += image->data[i+j*image->lx] ;  
01021                 Mx += (double)boxi * image->data[i+j*image->lx] ;
01022                 My += (double)boxj * image->data[i+j*image->lx] ;
01023                 /*-----------------------------------------------------------
01024                  * estimate the amplitude and the background 
01025                  * go through the margins of the fitting box and calculate the clean mean to
01026                  * determine the background 
01027                  */
01028                 if ( i == llx || i == llx + 2*halfbox_x -1 ||
01029                      j == lly || j == lly + 2*halfbox_y -1 )
01030                 {
01031                     backarray[n] = image->data[i+j*image->lx] ;
01032                     n++ ;
01033                 }
01034             }
01035         }
01036     }
01037     if ( M <= 0. )
01038     {
01039         cpl_msg_error("fit2dGaussian:"," only negative or zero values\n") ;
01040         cpl_free(backarray) ;
01041         return -1 ;
01042     }
01043     if ( n < 3 )
01044     {
01045         cpl_msg_error("fit2dGaussian:"," not enough data points to calculate background\n") ;
01046         cpl_free(backarray) ;
01047         return -1 ;
01048     }
01049     /* determine the background as median of the surrounding pixels */
01050     if ( FLT_MAX == (background = clean_mean(backarray, n, 10., 10.)) ) 
01051     {
01052         cpl_msg_error("fit2dGaussian:"," it was not possible to compute the clean mean of the background values\n") ;
01053         cpl_free(backarray) ;
01054         return -1 ;
01055     }
01056     cpl_free (backarray) ;
01057     /* now calculate the amplitude estimation */
01058     amplitude = maxval - background ;
01059     if ( amplitude < 1e-12 )
01060     {
01061         cpl_msg_error("fit2dGaussian:"," amplitude is too small\n") ;
01062         return -1 ;
01063     }
01064 
01065     /* determine the center of gravity = centroid */
01066     X0 = Mx / M ;
01067     Y0 = My / M ;
01068     /* if one of the values is outside the fitting box return with error */
01069     if ( X0 <= 0. || Y0 <= 0. || X0 >= 2.*(double)halfbox_x || 
01070          Y0 >= 2.*(double)halfbox_y )
01071     {
01072         cpl_msg_error("fit2dGaussian:"," center of gravity is outside the fitting box!\n") ;
01073         return -1 ;
01074     }
01075     
01076     /*------------------------------------------------------------------------------ 
01077      * put the data in the 2-d array xydat[][] (pixel position) and zdat[] (data values)
01078      * additionally, determine the second order momentum
01079      */
01080     n = 0 ;
01081     M = Mx = Mxx = My = Myy = Mxy = 0. ;
01082     boxi = boxj = 0 ;
01083     for ( j = lly ; j < lly + 2*halfbox_y ; j++ )
01084     {
01085         boxj = j - lly ;
01086         for ( i = llx ; i < llx + 2*halfbox_x ; i++ )
01087         {
01088             boxi = i - llx ;
01089             value = image->data[i+j*image->lx] ;
01090             if ( !qfits_isnan(value) )
01091             {
01092                 xydat[n][0] = (double) boxi ;
01093                 xydat[n][1] = (double) boxj ;
01094                 zdat[n]     = value ;
01095                 wdat[n]     = 1. ;
01096                 n++ ; 
01097                
01098                 /* now calculate the moments without background in the centroid coordinate system */
01099                 value -= background ;
01100                 xco = (double) boxi - X0 ;
01101                 yco = (double) boxj - Y0 ;
01102                 M   += value ;
01103                 Mx  += xco * value ;
01104                 My  += yco * value ;
01105                 Mxx += xco * xco * value ;
01106                 Myy += yco * yco * value ;
01107                 Mxy += xco * yco * value ;
01108             }
01109         }
01110     }
01111     if ( M <= 0. )
01112     {
01113         cpl_msg_error("fit2dGaussian:"," only negative or zero values\n") ;
01114         return -1 ;
01115     }
01116   
01117     /* ----------------------------------------------------------------
01118      * estimate the fwhm_x and fwhm_y and theta 
01119      */ 
01120   
01121     /* first scale the moments */
01122     Mx  /= M ;
01123     My  /= M ;
01124     Mxx /= M ;
01125     Myy /= M ;
01126     Mxy /= M ;
01127    
01128     denom = 2. * (Mxx*Myy - Mxy*Mxy) ;
01129     if ( denom == 0. )
01130     {
01131         cpl_msg_error("fit2dGaussian:"," denominator is zero!\n") ;
01132         return -1 ;
01133     }
01134     
01135     /* now associate the parameter list with the found estimates */
01136     fit_par[0] = X0 ;
01137     fit_par[1] = Y0 ;
01138     fit_par[2] = amplitude ;
01139     fit_par[3] = background ;
01140     fit_par[4] = Myy/denom ;
01141     fit_par[5] = Mxx/denom ;
01142     fit_par[6] = -Mxy/denom ;
01143 
01144     /* convert the moments to ellipse paramters */
01145     if ( 0 > gauss2Ellipse (fit_par) )
01146     {
01147         cpl_msg_error("fit2dGaussian:"," gauss2Ellipse does not run!\n") ;
01148         return -1 ;
01149     }
01150 
01151     /* total number of data points */
01152     ndata = 4 * halfbox_x * halfbox_y ;
01153     xdim = XDIMG ; /* dimension of xydat array */
01154     npar = NPAR ; /* number of parameters in the fit */
01155     its = ITSG ;
01156     lab = LABG ;
01157     tol = TOLG ;
01158     for ( i = 0 ; i < NPAR ; i++ )
01159     {
01160         derv_par[i] = 0. ;
01161     }
01162     
01163     if ( 0 > ( iters = lsqfitd ( &xydat[0][0],
01164                                  &xdim,
01165                                  zdat,
01166                                  wdat, 
01167                                  &ndata,
01168                                  fit_par,
01169                                  derv_par,
01170                                  mpar,
01171                                  &npar,
01172                                  &tol,
01173                                  &its,
01174                                  &lab )) ) 
01175     {
01176         cpl_msg_error("fit2dGaussian:"," least squares fit failed, error no: %d!\n", iters) ;
01177         return -1 ;
01178     }
01179 
01180     /* exclude impossible fit results */
01181     if ( fit_par[2] <= 0. || fit_par[4] < 0. || fit_par[5] < 0. )
01182     {
01183         cpl_msg_error("fit2dGaussian:"," sorry, some impossible negative fit results!\n") ;
01184         return -1 ;
01185     }
01186     fit_par[0] += llx ;
01187     fit_par[1] += lly ;
01188     if ( fit_par[0] < llx || fit_par[0] >= llx + 2*halfbox_x ||
01189          fit_par[1] < lly || fit_par[1] >= lly + 2*halfbox_y )
01190     {
01191         cpl_msg_error("fit2dGaussian:"," sorry, centroid after the fit outside the fitting box\n") ;
01192         return -1 ;
01193     }
01194 
01195     /* exchange fwhmx and fwhmy if |theta| is bigger than pi/4 and subtract pi/2 from theta */
01196     if ( fabs ( fit_par[6] ) > PI_NUMB / 4. )
01197     {
01198         /* first convert angle to smaller than 2 pi */
01199         if ( fabs (fit_par[6]) >= 2. * PI_NUMB )
01200         { 
01201             k = (int) (fit_par[6] / (2.*PI_NUMB)) ;
01202             if ( k > 0 ) 
01203             {
01204                 fit_par[6] -= k*2.*PI_NUMB ;
01205             }
01206             else
01207             {
01208                 fit_par[6] += k*2.*PI_NUMB ;
01209             }
01210         }
01211         /* first convert angle to smaller than pi/2 */
01212         if ( fabs (fit_par[6]) > PI_NUMB / 2. )
01213         {
01214             if ( fit_par[6] > 0. )
01215             {
01216                 fit_par[6] -= PI_NUMB ;
01217             }
01218             else
01219             {
01220                 fit_par[6] += PI_NUMB ;
01221             }
01222         }
01223 
01224         if ( fabs (fit_par[6]) > PI_NUMB / 4. )
01225         {
01226             temp       = fit_par[4] ;
01227             fit_par[4] = fit_par[5] ;
01228             fit_par[5] = temp ;
01229             if ( fit_par[6] < 0. )
01230             { 
01231                 fit_par[6] += PI_NUMB / 2. ;
01232             }
01233             else
01234             {
01235                 fit_par[6] -= PI_NUMB / 2. ;
01236             }  
01237         }
01238     }
01239     
01240     return iters ;
01241 }
01242 
01243 /*----------------------------------------------------------------------------
01244    Function     :       plotGaussian()
01245    In           :       image: image which should be fitted
01246                         parlist: parameters of 2D-Gaussian
01247    Out          :       image of the 2D-Gaussian
01248    Job          :       plots an image of a given 2D-Gaussian 
01249  ---------------------------------------------------------------------------*/
01250 
01251 OneImage * plotGaussian ( OneImage   * image, 
01252                           double     * parlist )
01253 {
01254     int col, row ;
01255     OneImage * retImage ;
01256     double xdat[2] ;
01257 
01258     if ( image == NullImage )
01259     {
01260         cpl_msg_error("plotGaussian","no input image given!\n") ;
01261         return NullImage ;
01262     }
01263     if ( parlist == NULL ) 
01264     {
01265         cpl_msg_error("plotGaussian","no Gaussian parameters given!\n") ;
01266         return NullImage ;
01267     }
01268  
01269     retImage = new_image (image->lx, image->ly) ;
01270 
01271     for ( row = 0 ; row < image->ly ; row++ )
01272     {
01273         for ( col = 0 ; col < image->lx ; col++ )
01274         {
01275             xdat[0] = (double) col ;
01276             xdat[1] = (double) row ;
01277             retImage->data[col+row*image->lx] = gaussianEllipse( xdat , parlist) ; 
01278         }
01279     }
01280 
01281     return retImage ;
01282 }
01283 
01284 /*----------------------------------------------------------------------------
01285    Function     :       gauss2Ellipse()
01286    In           :       parlist: parameters of 2D-Gaussian
01287    Out          :       parlist: transformed Gaussian parameters to ellipse parameters
01288    Job          :       converts gauss parameters to ellipse parameters. 
01289  ---------------------------------------------------------------------------*/
01290 
01291 static int gauss2Ellipse ( double     * parlist )
01292 {
01293     double a, b, c ;
01294     double ellipseconst ;
01295     double axisX, axisY, phi ;
01296     double p ;
01297     
01298     if ( parlist == NULL )
01299     {
01300         cpl_msg_error("gauss2Ellipse:"," no parameters given!\n") ;
01301         return -1 ;
01302     }
01303 
01304     a = parlist[4] ; /* fwhmx */
01305     b = parlist[5] ; /* fwhmy */
01306     c = parlist[6] ; /* theta */
01307 
01308     ellipseconst = 2. * log(2.) ;
01309 
01310     if ( a*b - c*c <= 0. )
01311     {
01312         cpl_msg_error("gauss2Ellipse:"," estimates of moments are unusable, they do not make an ellipse!\n") ;
01313         return -1 ;
01314     }
01315     
01316     if ( a == b )
01317     { 
01318         phi = 0. ;
01319     }
01320     else
01321     {
01322         phi = 0.5 * atan( 2. * c / (a-b) ) ;
01323     }
01324 
01325     p = sqrt ( (a-b) * (a-b) + 4. * c*c ) ;
01326 
01327     if ( a > b )
01328     {
01329         axisX = 2. * sqrt ( ellipseconst / (a+b+p) ) ;
01330         axisY = 2. * sqrt ( ellipseconst / (a+b-p) ) ;
01331     }
01332     else
01333     {
01334         axisX = 2. * sqrt ( ellipseconst / (a+b-p) ) ;
01335         axisY = 2. * sqrt ( ellipseconst / (a+b+p) ) ;
01336     }
01337 
01338     parlist[4] = axisX ;
01339     parlist[5] = axisY ;
01340     parlist[6] = phi ;
01341 
01342     return 0 ;
01343 }
01344 
01345 /*----------------------------------------------------------------------------
01346    Function     :       determineConversionFactor()
01347    In           :       cube: reduced data cube of a standard star
01348                         mag:  brightness of the standard star
01349                         exptime: exposure time read from the fits header
01350                         llx, lly:  lower left point of fitting box
01351                         halfbox_x, halfbox_y: half width of a box inside which
01352                                               a 2D-Gauss fit is carried out
01353                         check: if search box is outside image definition.
01354 
01355    Out          :       intensity conversion value: magnitude per counts/s 
01356                         -FLT_MAX if error occurred.
01357    Job          :       determines an intensity conversion factor for the instrument
01358                         by fitting a 2D-Gaussian to an collapsed image of a standard star
01359                         with known brightness (only for non-AO observations).
01360                         Then the resulting Gaussian is integrated and the counts
01361                         are divided by the exposure time (Fits header information) 
01362  ---------------------------------------------------------------------------*/
01363 
01364 float determineConversionFactor ( OneCube * cube, 
01365                                   float     mag,
01366                                   float     exptime,
01367                                   int       llx,
01368                                   int       lly,
01369                                   int       halfbox_x,
01370                                   int       halfbox_y, int* check )
01371 {
01372     int row, col, i ;
01373     int first_row, first_col ;
01374     int last_row, last_col ;
01375     OneImage * summedIm ;
01376     float factor ;
01377     int mpar[7] ;
01378     double fit_par[7] ;
01379     double derv_par[7] ;
01380     int fitInd ;
01381     double sum ;
01382     double xdat[2] ;
01383 
01384     if ( NullCube == cube )
01385     {
01386         cpl_msg_error("determineConversionFactor:"," no cube given!\n") ;
01387         return -FLT_MAX ;
01388     }
01389     if ( halfbox_x <= 0 || halfbox_y <= 0 || 2*halfbox_x > cube->lx || 2*halfbox_y > cube->ly)
01390     {
01391         cpl_msg_error("determineConversionFactor:"," wrong width of halfbox given!\n") ;
01392         return -FLT_MAX ;
01393     }
01394     if ( exptime <= 0. )
01395     {
01396         cpl_msg_error("determineConversionFactor:"," impossible exposure time given !\n") ;
01397         return -FLT_MAX ;
01398     }
01399 
01400     /* collapse the cube to be able to do 2D-Gaussian fitting */
01401     if ( NullImage == (summedIm = sumCubeToImage(cube)) )
01402     {
01403         cpl_msg_error("determineConversionFactor:"," averageCubeToImage failed!\n") ;
01404         return -FLT_MAX ;
01405     }
01406 
01407     /* call the 2D-Gaussian fit routine */
01408     for ( i = 0 ; i < 7 ; i++ )
01409     {
01410         mpar[i] = 1 ;
01411     }
01412     if ( -1 == (fitInd = fit2dGaussian(summedIm, fit_par, derv_par, mpar, llx, lly, halfbox_x, halfbox_y, check)) )
01413     {
01414         cpl_msg_error("determineConversionFactor:"," fit2dGaussian failed!\n") ;
01415         destroy_image( summedIm) ;
01416         return -FLT_MAX ;
01417     }
01418     destroy_image(summedIm) ;
01419 
01420     /* now integrate the found 2D Gaussian by first subtracting the background */
01421     if  ((fit_par[0] - halfbox_x) < 0) {
01422       first_col=0;
01423       check++;
01424     } else {
01425       first_col=(fit_par[0] - halfbox_x);
01426     }
01427 
01428     if ((fit_par[0] + halfbox_x) < cube->lx) {
01429       last_col  = (fit_par[0] + halfbox_x);
01430     } else {
01431       last_col = (cube->lx-1) ;
01432       check++;
01433     }
01434 
01435     if ((fit_par[1] - halfbox_y) < 0) {
01436       first_row=0;
01437       check++;
01438     } else {
01439       first_row=(fit_par[1] - halfbox_y) ;
01440     }
01441 
01442     if ((fit_par[1] + halfbox_y) < cube->ly) {
01443       last_row=(fit_par[1] + halfbox_y);
01444     } else {
01445       last_row= (cube->ly-1);
01446       check++;
01447     }
01448 
01449 
01450     if ( first_col < 0 || first_row < 0 || last_col >= cube->lx || last_row >= cube->ly )
01451     {
01452         cpl_msg_error("determineConversionFactor:"," star badly centered in FOV or fitting box too big!\n") ;
01453         return -FLT_MAX ;
01454     }
01455     sum = 0. ;
01456     for ( row = first_row ; row < last_row ; row++ )  
01457     {
01458         for( col = first_col ; col < last_col ; col++ )
01459         {
01460             xdat[0] = (double) col ;
01461             xdat[1] = (double) row ;
01462             sum += (gaussianEllipse( xdat, fit_par ) - fit_par[3]) ;
01463         }
01464     }
01465     if ( sum <= 0. )
01466     {
01467         cpl_msg_error("determineConversionFactor:"," zero or negative sum of counts!\n") ;
01468         return -FLT_MAX ;
01469     }
01470     factor = mag / (float)sum * exptime ;
01471     return factor ;
01472 }
01473                               
01474 /*--------------------------------------------------------------------------*/

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