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 /*--------------------------------------------------------------------------*/
1.2.13.1 written by Dimitri van Heesch,
© 1997-2001