00001 /*******************************************************************************
00002 * E.S.O. - VLT project
00003 *
00004 *
00005 *
00006 * who when what
00007 * -------- -------- ----------------------------------------------
00008 * schreib 05/06/00 created
00009 */
00010
00011 /************************************************************************
00012 * NAME
00013 * recipes.c -
00014 * some numerical recipes
00015 *
00016 * SYNOPSIS
00017 * #include "recipes.h"
00018 *
00019 * 1) void myfit(float x[], float y[], int ndata, float sig[], int mwt, float *a,
00020 * float *b, float *siga, float *sigb, float *chi2, float *q)
00021 *
00022 * 2) pixelvalue median(pixelvalue * array, int n)
00023 *
00024 * 3) float gaussian ( float * xdat, float * parlist, int * npar )
00025 *
00026 * 4) void gaussian_deriv( float * xdat, float * parlist, float * dervs, int * npar )
00027 *
00028 * 5) static int invmat (void)
00029 *
00030 * 6) static void getmat ( float * xdat,
00031 * int * xdim,
00032 * float * ydat,
00033 * float * wdat,
00034 * int * ndat,
00035 * float * fpar,
00036 * float * epar,
00037 * int * npar )
00038 *
00039 * 7) static int getvec ( float * xdat,
00040 * int * xdim,
00041 * float * ydat,
00042 * float * wdat,
00043 * int * ndat,
00044 * float * fpar,
00045 * float * epar,
00046 * int * npar )
00047 *
00048 * 8) int lsqfit_c ( float * xdat,
00049 * int * xdim,
00050 * float * ydat,
00051 * float * wdat,
00052 * int * ndat,
00053 * float * fpar,
00054 * float * epar,
00055 * int * mpar,
00056 * int * npar,
00057 * float * tol ,
00058 * int * its ,
00059 * float * lab )
00060 *
00061 * 9) int nint ( double x )
00062 *
00063 * 10) int correlation ( float * data1, float * data2, int ndata )
00064 *
00065 * 11) float clean_mean( float * array,
00066 * int n_elements,
00067 * float throwaway_low,
00068 * float throwaway_high )
00069 *
00070 * 12) void errorPrint(char * text)
00071 *
00072 * 13) void commentPrint(char * text)
00073 *
00074 * 14) void convertZEROsTo0ForImages(OneImage * im)
00075 *
00076 * 15) void convertZEROsTo0ForCubes(OneCube * cube)
00077 *
00078 * 16) void convert0ToZEROsForImages(OneImage * im)
00079 *
00080 * 17) void convert0ToZEROsForCubes(OneCube * cube)
00081 *
00082 * 18) void invert(OneImage * im)
00083 *
00084 * DESCRIPTION
00085 * 1) fits a straight line to a set of x, y data points by
00086 * minimizing chi-square. Numeric. Rec. page 665
00087 * 2) sorts an array into ascending numerical order by
00088 * using pixel_qsort(), and derives the median depending
00089 * on an odd or even number of array elements
00090 * 3) calculates the value of a gaussian with parameters
00091 * parlist at the position xdat
00092 * 4) calculates the partial derivatives for a gaussian with
00093 * parameters parlist at position xdat
00094 * 5) calculates the inverse of matrix2. The algorithm used
00095 * is the Gauss-Jordan algorithm described in Stoer,
00096 * Numerische Mathematik, 1. Teil.
00097 * 6) builds the matrix
00098 * 7) calculates the correction vector. The matrix has been
00099 * built by getmat(), we only have to rescale it for the
00100 * current value of labda. The matrix is rescaled so that
00101 * the diagonal gets the value 1 + labda.
00102 * Next we calculate the inverse of the matrix and then
00103 * the correction vector.
00104 * 8) this is a routine for making a least-squares fit of a
00105 * function to a set of data points. The method used is
00106 * described in: Marquardt, J.Soc.Ind.Appl.Math. 11. 431 (1963).
00107 * This method is a mixture of the steepest descent method
00108 * and the Taylor method.
00109 * 9) this routine determines the nearest integer to a
00110 * specified real value
00111 * 10) this routine computes the cross correlation between
00112 * two data arrays of the same size and returns the
00113 * integer shift between both arrays
00114 * 11) this routine computes the clean mean of a given data
00115 * array that means the array is first sorted and
00116 * a given percentage of the lowest and the highest values
00117 * is not considered for averaging
00118 * 12) this routine gives python the chance to print out
00119 * error messages
00120 * 13) this routine gives python the chance to print out
00121 * comment messages
00122 * 14) this function converts any ZEROs in an image to 0
00123 * 15) this function converts any ZEROs in a cube to 0
00124 * 16) this function converts any 0 in an image to ZERO
00125 * 17) this function converts any 0 in a cube to ZERO
00126 * 18) this function multiplies each image value by -1.
00127 *
00128 * FILES
00129 *
00130 * ENVIRONMENT
00131 *
00132 * RETURN VALUES
00133 *
00134 * CAUTIONS
00135 *
00136 * EXAMPLES
00137 *
00138 * SEE ALSO
00139 *
00140 * BUGS
00141 *
00142 *------------------------------------------------------------------------
00143 */
00144
00145 #include "vltPort.h"
00146
00147 /*
00148 * System Headers
00149 */
00150
00151 /*
00152 * Local Headers
00153 */
00154
00155 #include "recipes.h"
00156
00157 /*----------------------------------------------------------------------------
00158 * Local variables
00159 *--------------------------------------------------------------------------*/
00160
00161 static float sqrarg ;
00162
00163 static double chi1 ; /* old reduced chi-squared */
00164 static double chi2 ; /* new reduced chi-squared */
00165 static double labda ; /* mixing parameter */
00166 static double vec[MAXPAR] ; /* correction vector */
00167 static double matrix1[MAXPAR][MAXPAR] ; /* original matrix */
00168 static double matrix2[MAXPAR][MAXPAR] ; /* inverse of matrix1 */
00169 static int nfree ; /* number of free parameters */
00170 static int parptr[MAXPAR] ; /* parameter pointer */
00171
00172 /*----------------------------------------------------------------------------
00173 * Defines
00174 *--------------------------------------------------------------------------*/
00175
00176 #define SQR(a) (sqrarg = (a) , sqrarg*sqrarg)
00177
00178 /*----------------------------------------------------------------------------
00179 * Functions private to this module
00180 *--------------------------------------------------------------------------*/
00181
00182
00183 static int invmat (void) ;
00184
00185 static void getmat ( float * xdat,
00186 int * xdim,
00187 float * ydat,
00188 float * wdat,
00189 int * ndat,
00190 float * fpar,
00191 float * epar/*,
00192 int * npar */) ;
00193
00194 static int getvec ( float * xdat,
00195 int * xdim,
00196 float * ydat,
00197 float * wdat,
00198 int * ndat,
00199 float * fpar,
00200 float * epar,
00201 int * npar ) ;
00202
00203
00204
00205
00206 /*----------------------------------------------------------------------------
00207 * Function codes
00208 *--------------------------------------------------------------------------*/
00209
00210 /*----------------------------------------------------------------------------
00211 Function : myfit()
00212 In : set of data points, individual standard deviations,
00213 parameters of a straight line, chi-square, goodness-of-
00214 fit probability
00215 Out : void
00216 Job : fits a straight line to a set of x, y data points by
00217 minimizing chi-square. Numeric. Rec. page 665
00218 ---------------------------------------------------------------------------*/
00219 /*-------------------------------------------------------------------------*/
00231 /*--------------------------------------------------------------------------*/
00232
00233 OneCube *
00234 copy_cube_range(
00235 OneCube *src_cube,
00236 const int z_min,
00237 const int z_max
00238 )
00239 {
00240 OneCube * dest_cube ;
00241 int i,k ;
00242 const int z_siz=z_max-z_min;
00243
00244 if (src_cube==NULL) return NULL ;
00245 /* allocate data zone, fill up information structure */
00246 dest_cube = new_cube( src_cube->lx,
00247 src_cube->ly,
00248 z_siz) ;
00249 dest_cube ->orig_ptype = src_cube->orig_ptype ;
00250
00251 /* Then copy data zones from one cube to the other */
00252 if (src_cube->type==1){
00253 dest_cube->data = (pixelvalue *) calloc(z_siz * src_cube->lx * src_cube->ly, sizeof(pixelvalue));
00254 dest_cube->planes= (OneImage *) calloc(z_siz, sizeof(OneImage));
00255 dest_cube->type=1;
00256 k=0;
00257 for (i=z_min ; i<z_max ; i++)
00258 {
00259 fill_new_image(src_cube->lx, src_cube->ly, dest_cube->planes+k, dest_cube->data+(i*src_cube->lx*src_cube->ly));
00260 dest_cube->plane[k] = dest_cube->planes+k;
00261 k++;
00262 }
00263 memcpy(dest_cube->data, src_cube->data+(z_min * src_cube->lx * src_cube->ly)*sizeof(pixelvalue), (size_t) (z_siz * src_cube->lx * src_cube->ly) * sizeof(pixelvalue));
00264 }
00265 else {
00266 k=0;
00267 for (i=z_min ; i<z_max ; i++) {
00268 dest_cube->plane[k] = copy_image(src_cube->plane[i]) ;
00269 }
00270 }
00271 return dest_cube ;
00272 }
00273
00274
00275
00276
00277 /*==================================================================*/
00278 void myfit(float x[], float y[], int ndata, float sig[], int mwt, float *a,
00279 float *b, float *siga, float *sigb, float *chi2, float *q)
00280 {
00281 int i ;
00282 float wt, t, sxoss, sx=0., sy=0., st2=0., ss, sigdat ;
00283
00284 *b = 0. ; /*accumulate sums ...*/
00285 if ( mwt )
00286 {
00287 ss = 0. ;
00288 for ( i = 0 ; i < ndata ; i++ ) /*... with weights*/
00289 {
00290 wt = 1./SQR(sig[i]) ;
00291 ss += wt ;
00292 sx += x[i]*wt ;
00293 sy += y[i]*wt ;
00294 }
00295 }
00296 else
00297 {
00298 for ( i = 0 ; i < ndata ; i++ ) /*... or without weights*/
00299 {
00300 sx += x[i] ;
00301 sy += y[i] ;
00302 }
00303 ss = ndata ;
00304 }
00305 sxoss = sx/ss ;
00306
00307 if ( mwt )
00308 {
00309 for ( i = 0 ; i < ndata ; i ++ )
00310 {
00311 t = (x[i] - sxoss)/sig[i] ;
00312 st2 += t*t ;
00313 *b += t*y[i]/sig[i] ;
00314 }
00315 }
00316 else
00317 {
00318 for ( i = 0 ; i < ndata ; i++ )
00319 {
00320 t = x[i] - sxoss ;
00321 st2 += t*t ;
00322 *b += t*y[i] ;
00323 }
00324 }
00325
00326 *b /= st2 ;
00327 *a = (sy - sx*(*b))/ss ;
00328 *siga = sqrt ((1.0 + sx*sx/(ss*st2))/ss) ;
00329 *sigb = sqrt (1.0/st2) ;
00330 *chi2 = 0.0 ; /*calculate chi-square*/
00331 if ( mwt == 0 )
00332 {
00333 for ( i = 0 ; i < ndata ; i++ )
00334 {
00335 *chi2 += SQR (y[i] - (*a) - (*b)*x[i]) ;
00336 }
00337 *q = 1. ;
00338
00339 /*------------------------------------------------------------------
00340 * for unweighted data evaluate typical sig using chi2, and adjust
00341 * the standard deviation
00342 */
00343 sigdat = sqrt ((*chi2)/(ndata - 2)) ;
00344 *siga *= sigdat ;
00345 *sigb *= sigdat ;
00346 }
00347 else
00348 {
00349 for (i = 0 ; i < ndata ; i++)
00350 {
00351 *chi2 += SQR ((y[i] - (*a) - (*b) * x[i])/sig[i]) ;
00352 }
00353 *q = 1. ; /* delete rest of lines. q is not a good value */
00354 }
00355 }
00356
00357 /*----------------------------------------------------------------------------
00358 Function : median()
00359 In : array of pixelvalues starting at index 0,
00360 number of array elements
00361 Out : median of pixelvalues
00362 Job : sorts an array into ascending numerical order by
00363 using pixel_qsort(), and derives the median depending
00364 on an odd or even number of array elements
00365 ---------------------------------------------------------------------------*/
00366
00367 pixelvalue median(pixelvalue * array, int n)
00368 {
00369 pixelvalue med ;
00370
00371 if ( array == NULL || n <= 0 )
00372 {
00373 cpl_msg_warning ("median:","nothing in the pixelvalue array, ZERO returned\n") ;
00374 return ZERO ;
00375 }
00376
00377 if ( n == 1 )
00378 {
00379 return array[0] ;
00380 }
00381
00382 pixel_qsort((float*) array, n) ;
00383 if ( n % 2 == 1 )
00384 {
00385 med = array[n/2] ;
00386 }
00387 else
00388 {
00389 med = (array[n/2] + array[n/2 - 1])/2. ;
00390 }
00391 return med ;
00392 }
00393
00394 float f_median(float * array, int n)
00395 {
00396 pixelvalue p_array[100];
00397 int i;
00398
00399 for (i=0;i<n;i++)
00400 p_array[i]= (pixelvalue) array[i];
00401
00402 return (float) median(p_array, n);
00403 }
00404
00405 /*----------------------------------------------------------------------------
00406 Function : gaussian()
00407 In : position array xdat, parameter list parlist, number of
00408 parameters in the list npar
00409 The parameters are:
00410 parlist(0): Amplitude
00411 parlist(1): FWHM
00412 parlist(2): x0, center of gauss versus center of column
00413 parlist(3): Zero level
00414 Out : function value of a 1-d gaussian:
00415 F(x) = par(0) * EXP ( -4.0 * LN(2.0) *
00416 [((x - parlist(2))/parlist(1))^2]) + parlist(3)
00417 FWHM^2 = 8 * ln (2) * sigma^2
00418 Job : calculates the value of a gaussian with parameters
00419 parlist at the position xdat
00420 ---------------------------------------------------------------------------*/
00421
00422 float gaussian ( float * xdat, float * parlist/*, int * npar*/ )
00423 {
00424 double xd ; /* FWHM's of gauss function */
00425 double x ; /* position */
00426
00427 xd = fabs((double) parlist[1]) ;
00428 x = (double) xdat[0] - (double) parlist[2] ;
00429 return (float) (
00430 (double) parlist[0] * exp( -4.0 * log(2.0) * (x/xd) * (x/xd) )
00431 + (double) parlist[3] ) ;
00432 }
00433
00434
00435 /*----------------------------------------------------------------------------
00436 Function : gaussian_deriv()
00437 In : position array xdat, parameter list parlist, number of
00438 parameters in the list npar
00439 The parameters are:
00440 parlist(0): Amplitude
00441 parlist(1): Axis(FWHM)
00442 parlist(2): x0, center of gauss versus center of column
00443 parlist(3): Zero level
00444 derivative value of gaussian at position xdat: dervs
00445 dervs[0]: partial derivative by the amplitude
00446 dervs[1]: partial derivative by FWHM
00447 dervs[2]: partial derivative by the x position (parlist[2])
00448 dervs[3]: partial derivative by the zero level
00449 Out : nothing, void
00450 Job : calculates the partial derivatives for a gaussian with
00451 parameters parlist at position xdat
00452 ---------------------------------------------------------------------------*/
00453
00454 void gaussian_deriv( float * xdat, float * parlist, float * dervs/*, int * npar*/ )
00455 {
00456 double xd ; /* FWHM of gaussian */
00457 double x, expon ; /* position and exponent */
00458
00459 xd = fabs( (double) parlist[1] ) ;
00460
00461 /* offset from peak position */
00462 x = (double) xdat[0] - (double) parlist[2] ;
00463
00464 /* determine the derivatives: */
00465 expon = -4.0 * log(2.0) * (x/xd) * (x/xd) ;
00466 expon = exp( expon ) ;
00467
00468 /* partial derivative by the amplitude */
00469 dervs[0] = (float) expon ;
00470
00471 /* calculate a * exp(-arg) */
00472 expon = (double) parlist[0] * expon ;
00473
00474 /* partial derivative by FWHM */
00475 dervs[1] = (float) ( expon * 8.0 * log(2.0) * x*x / (xd*xd*xd) ) ;
00476
00477 /* partial derivative by the x position (parlist[2]) */
00478 dervs[2] = (float) (expon * 8.0 * log(2.0) * x/(xd*xd) ) ;
00479
00480 /* partial derivative by the zero level */
00481 dervs[3] = 1.0 ;
00482 }
00483
00484
00485 /*----------------------------------------------------------------------------
00486 Function : invmat()
00487 In : void
00488 Out : integer (0 if it worked, -6 if determinant is zero)
00489 Job : calculates the inverse of matrix2. The algorithm used
00490 is the Gauss-Jordan algorithm described in Stoer,
00491 Numerische Mathematik, 1. Teil.
00492 ---------------------------------------------------------------------------*/
00493
00494 static int invmat (void)
00495 {
00496 double even ;
00497 double hv[MAXPAR] ;
00498 double mjk ;
00499 double rowmax ;
00500 int evin ;
00501 int i, j, k, row ;
00502 int per[MAXPAR] ;
00503
00504 /* set permutation array */
00505 for ( i = 0 ; i < nfree ; i++ )
00506 {
00507 per[i] = i ;
00508 }
00509
00510 for ( j = 0 ; j < nfree ; j++ ) /* in j-th column */
00511 {
00512 /* determine largest element of a row */
00513 rowmax = fabs ( matrix2[j][j] ) ;
00514 row = j ;
00515
00516 for ( i = j + 1 ; i < nfree ; i++ )
00517 {
00518 if ( fabs ( matrix2[i][j] ) > rowmax )
00519 {
00520 rowmax = fabs( matrix2[i][j] ) ;
00521 row = i ;
00522 }
00523 }
00524
00525 /* determinant is zero! */
00526 if ( matrix2[row][j] == 0.0 )
00527 {
00528 return -6 ;
00529 }
00530
00531 /* if the largest element is not on the diagonal, then permutate rows */
00532 if ( row > j )
00533 {
00534 for ( k = 0 ; k < nfree ; k++ )
00535 {
00536 even = matrix2[j][k] ;
00537 matrix2[j][k] = matrix2[row][k] ;
00538 matrix2[row][k] = even ;
00539 }
00540 /* keep track of permutation */
00541 evin = per[j] ;
00542 per[j] = per[row] ;
00543 per[row] = evin ;
00544 }
00545
00546 /* modify column */
00547 even = 1.0 / matrix2[j][j] ;
00548 for ( i = 0 ; i < nfree ; i++ )
00549 {
00550 matrix2[i][j] *= even ;
00551 }
00552 matrix2[j][j] = even ;
00553
00554 for ( k = 0 ; k < j ; k++ )
00555 {
00556 mjk = matrix2[j][k] ;
00557 for ( i = 0 ; i < j ; i++ )
00558 {
00559 matrix2[i][k] -= matrix2[i][j] * mjk ;
00560 }
00561 for ( i = j + 1 ; i < nfree ; i++ )
00562 {
00563 matrix2[i][k] -= matrix2[i][j] * mjk ;
00564 }
00565 matrix2[j][k] = -even * mjk ;
00566 }
00567
00568 for ( k = j + 1 ; k < nfree ; k++ )
00569 {
00570 mjk = matrix2[j][k] ;
00571 for ( i = 0 ; i < j ; i++ )
00572 {
00573 matrix2[i][k] -= matrix2[i][j] * mjk ;
00574 }
00575 for ( i = j + 1 ; i < nfree ; i++ )
00576 {
00577 matrix2[i][k] -= matrix2[i][j] * mjk ;
00578 }
00579 matrix2[j][k] = -even * mjk ;
00580 }
00581 }
00582
00583 /* finally, repermute the columns */
00584 for ( i = 0 ; i < nfree ; i++ )
00585 {
00586 for ( k = 0 ; k < nfree ; k++ )
00587 {
00588 hv[per[k]] = matrix2[i][k] ;
00589 }
00590 for ( k = 0 ; k < nfree ; k++ )
00591 {
00592 matrix2[i][k] = hv[k] ;
00593 }
00594 }
00595
00596 /* all is well */
00597 return 0 ;
00598 }
00599
00600 /*----------------------------------------------------------------------------
00601 Function : getmat()
00602 In : xdat: position
00603 xdim: factor of the indices of the position array
00604 ydat: real data
00605 wdat: weights
00606 ndat: number of data points
00607 fpar: function parameters
00608 epar: partial derivatives of the function
00609 npar: number of function parameters
00610 Out : void
00611 Job : builds the matrix
00612 ---------------------------------------------------------------------------*/
00613
00614 static void getmat ( float * xdat,
00615 int * xdim,
00616 float * ydat,
00617 float * wdat,
00618 int * ndat,
00619 float * fpar,
00620 float * epar/*,
00621 int * npar */)
00622 {
00623 double wd ;
00624 double wn ;
00625 double yd ;
00626 int i, j, n ;
00627
00628 for ( j = 0 ; j < nfree ; j++ )
00629 {
00630 vec[j] = 0.0 ; /* zero vector */
00631 for ( i = 0 ; i<= j ; i++ ) /* zero matrix only on and below diagonal */
00632 {
00633 matrix1[j][i] = 0.0 ;
00634 }
00635 }
00636 chi2 = 0.0 ; /* reset reduced chi-squared */
00637
00638 /* loop through data points */
00639 for ( n = 0 ; n < (*ndat) ; n++ )
00640 {
00641 wn = wdat[n] ;
00642 if ( wn > 0.0 ) /* legal weight ? */
00643 {
00644 yd = ydat[n] - gaussian( &xdat[(*xdim) * n], fpar/*, npar*/ ) ;
00645 gaussian_deriv( &xdat[(*xdim) * n], fpar, epar/*, npar*/ ) ;
00646 chi2 += yd * yd * wn ; /* add to chi-squared */
00647 for ( j = 0 ; j < nfree ; j++ )
00648 {
00649 wd = epar[parptr[j]] * wn ; /* weighted derivative */
00650 vec[j] += yd * wd ; /* fill vector */
00651 for ( i = 0 ; i <= j ; i++ ) /* fill matrix */
00652 {
00653 matrix1[j][i] += epar[parptr[i]] * wd ;
00654 }
00655 }
00656 }
00657 }
00658 }
00659
00660
00661 /*----------------------------------------------------------------------------
00662 Function : getvec()
00663 In : xdat: position
00664 xdim: factor of the indices of the position array
00665 ydat: real data
00666 wdat: weights
00667 ndat: number of data points
00668 fpar: function parameters
00669 epar: partial derivatives of the function
00670 npar: number of function parameters
00671 Out : integer (0 if it had worked, -5 or -7
00672 if diagonal element is wrong, or -6,
00673 if determinant is zero )
00674 Job : calculates the correction vector. The matrix has been
00675 built by getmat(), we only have to rescale it for the
00676 current value of labda. The matrix is rescaled so that
00677 the diagonal gets the value 1 + labda.
00678 Next we calculate the inverse of the matrix and then
00679 the correction vector.
00680 ---------------------------------------------------------------------------*/
00681
00682 static int getvec ( float * xdat,
00683 int * xdim,
00684 float * ydat,
00685 float * wdat,
00686 int * ndat,
00687 float * fpar,
00688 float * epar,
00689 int * npar )
00690 {
00691 double dj ;
00692 double dy ;
00693 double mii ;
00694 double mji ;
00695 double mjj ;
00696 double wn ;
00697 int i, j, n, r ;
00698
00699 /* loop to modify and scale the matrix */
00700 for ( j = 0 ; j < nfree ; j++ )
00701 {
00702 mjj = matrix1[j][j] ;
00703 if ( mjj <= 0.0 ) /* diagonal element wrong */
00704 {
00705 return -5 ;
00706 }
00707 mjj = sqrt( mjj ) ;
00708 for ( i = 0 ; i < j ; i++ )
00709 {
00710 mji = matrix1[j][i] / mjj / sqrt( matrix1[i][i] ) ;
00711 matrix2[i][j] = matrix2[j][i] = mji ;
00712 }
00713 matrix2[j][j] = 1.0 + labda ; /* scaled value on diagonal */
00714 }
00715
00716 if ( (r = invmat()) ) /* invert matrix inlace */
00717 {
00718 return r ;
00719 }
00720
00721 for ( i = 0 ; i < (*npar) ; i ++ )
00722 {
00723 epar[i] = fpar[i] ;
00724 }
00725
00726 /* loop to calculate correction vector */
00727 for ( j = 0 ; j < nfree ; j++ )
00728 {
00729 dj = 0.0 ;
00730 mjj = matrix1[j][j] ;
00731 if ( mjj <= 0.0) /* not allowed */
00732 {
00733 return -7 ;
00734 }
00735 mjj = sqrt ( mjj ) ;
00736 for ( i = 0 ; i < nfree ; i++ )
00737 {
00738 mii = matrix1[i][i] ;
00739 if ( mii <= 0.0 )
00740 {
00741 return -7 ;
00742 }
00743 mii = sqrt( mii ) ;
00744 dj += vec[i] * matrix2[j][i] / mjj / mii ;
00745 }
00746 epar[parptr[j]] += dj ; /* new parameters */
00747 }
00748 chi1 = 0.0 ; /* reset reduced chi-squared */
00749
00750 /* loop through the data points */
00751 for ( n = 0 ; n < (*ndat) ; n++ )
00752 {
00753 wn = wdat[n] ; /* get weight */
00754 if ( wn > 0.0 ) /* legal weight */
00755 {
00756 dy = ydat[n] - gaussian( &xdat[(*xdim) * n], epar/*, npar*/ ) ;
00757 chi1 += wdat[n] * dy * dy ;
00758 }
00759 }
00760 return 0 ;
00761 }
00762
00763
00764
00765 /*----------------------------------------------------------------------------
00766 Function : lsqfit_c()
00767 In : xdat: position, coordinates of data points.
00768 xdat is 2 dimensional: XDAT ( XDIM, NDAT )
00769 xdim: dimension of fit
00770 ydat: data points
00771 wdat: weights for data points
00772 ndat: number of data points
00773 fpar: on input contains initial estimates of the
00774 parameters for non-linear fits, on output the
00775 fitted parameters.
00776 epar: contains estimates of the errors in fitted parameters
00777 mpar: logical mask telling which parameters are free (non-zero)
00778 and which parameters are fixed (0)
00779 npar: number of function parameters ( free + fixed )
00780 tol: relative tolerance. lsqfit stops when successive iterations
00781 fail to produce a decrement in reduced chi-squared less
00782 than tol. If tol is less than the minimum tolerance
00783 possible, tol will be set to this value. This means
00784 that maximum accuracy can be obtained by setting
00785 tol = 0.0.
00786 its: maximum number of iterations
00787 lab: mixing parameter, lab determines the initial weight
00788 of steepest descent method relative to the Taylor method
00789 lab should be a small value (i.e. 0.01). lab can only
00790 be zero when the partial derivatives are independent
00791 of the parameters. In fact in this case lab should be
00792 exactly equal to zero.
00793 Out : returns number of iterations needed to achieve convergence
00794 according to tol. When this number is negative, the fitting
00795 was not continued because a fatal error occurred:
00796 -1 too many free parameters, maximum is 32
00797 -2 no free parameters
00798 -3 not enough degrees of freedom
00799 -4 maximum number of iterations too small to obtain
00800 a solution which satisfies tol.
00801 -5 diagonal of matrix contains elements which are zero
00802 -6 determinant of the coefficient matrix is zero
00803 -7 square root of a negative number
00804 Job : this is a routine for making a least-squares fit of a
00805 function to a set of data points. The method used is
00806 described in: Marquardt, J.Soc.Ind.Appl.Math. 11. 431 (1963).
00807 This method is a mixture of the steepest descent method
00808 and the Taylor method.
00809 ---------------------------------------------------------------------------*/
00810
00811 int lsqfit_c ( float * xdat,
00812 int * xdim,
00813 float * ydat,
00814 float * wdat,
00815 int * ndat,
00816 float * fpar,
00817 float * epar,
00818 int * mpar,
00819 int * npar,
00820 float * tol ,
00821 int * its ,
00822 float * lab )
00823 {
00824 int i, n, r ;
00825 int itc ; /* fate of fit */
00826 int found ; /* fit converged: 1, not yet converged: 0 */
00827 int nuse ; /* number of useable data points */
00828 double tolerance ; /* accuracy */
00829
00830 itc = 0 ; /* fate of fit */
00831 found = 0 ; /* reset */
00832 nfree = 0 ; /* number of free parameters */
00833 nuse = 0 ; /* number of legal data points */
00834
00835 if ( *tol < (FLT_EPSILON * 10.0 ) )
00836 {
00837 tolerance = FLT_EPSILON * 10.0 ; /* default tolerance */
00838 }
00839 else
00840 {
00841 tolerance = *tol ; /* tolerance */
00842 }
00843
00844 labda = fabs( *lab ) * LABFAC ; /* start value for mixing parameter */
00845 for ( i = 0 ; i < (*npar) ; i++ )
00846 {
00847 if ( mpar[i] )
00848 {
00849 if ( nfree > MAXPAR ) /* too many free parameters */
00850 {
00851 return -1 ;
00852 }
00853 parptr[nfree++] = i ; /* a free parameter */
00854 }
00855 }
00856
00857 if (nfree == 0) /* no free parameters */
00858 {
00859 return -2 ;
00860 }
00861
00862 for ( n = 0 ; n < (*ndat) ; n++ )
00863 {
00864 if ( wdat[n] > 0.0 ) /* legal weight */
00865 {
00866 nuse ++ ;
00867 }
00868 }
00869
00870 if ( nfree >= nuse )
00871 {
00872 return -3 ; /* no degrees of freedom */
00873 }
00874 if ( labda == 0.0 ) /* linear fit */
00875 {
00876 for ( i = 0 ; i < nfree ; fpar[parptr[i++]] = 0.0 ) ; /* initialize fpar array */
00877 getmat ( xdat, xdim, ydat, wdat, ndat, fpar, epar/*, npar*/ ) ;
00878 r = getvec ( xdat, xdim, ydat, wdat, ndat, fpar, epar, npar ) ;
00879 if ( r ) /* error */
00880 {
00881 return r ;
00882 }
00883 for ( i = 0 ; i < (*npar) ; i++ )
00884 {
00885 fpar[i] = epar[i] ; /* save new parameters */
00886 epar[i] = 0.0 ; /* and set errors to zero */
00887 }
00888 chi1 = sqrt( chi1 / (double) (nuse - nfree) ) ;
00889 for ( i = 0 ; i < nfree ; i++ )
00890 {
00891 if ( (matrix1[i][i] <= 0.0 ) || (matrix2[i][i] <= 0.0) )
00892 {
00893 return -7 ;
00894 }
00895 epar[parptr[i]] = chi1 * sqrt( matrix2[i][i] ) / sqrt( matrix1[i][i] ) ;
00896 }
00897 }
00898 else /* non-linear fit */
00899 {
00900 /*----------------------------------------------------------------
00901 * the non-linear fit uses the steepest descent method in combination
00902 * with the Taylor method. The mixing of these methods is controlled
00903 * by labda. In the outer loop ( called the iteration loop ) we build
00904 * the matrix and calculate the correction vector. In the inner loop
00905 * (called the interpolation loop) we check whether we have obtained a
00906 * better solution than the previous one. If so, we leave the inner loop
00907 * else we increase labda ( give more weight to the steepest descent method)
00908 * calculate the correction vector and check again. After the inner loop
00909 * we do a final check on the goodness of the fit and if this satisfies
00910 * the tolerance we calculate the errors of the fitted parameters.
00911 */
00912 while ( !found ) /* iteration loop */
00913 {
00914 if ( itc++ == (*its) ) /* increase iteration counter */
00915 {
00916 return -4 ;
00917 }
00918 getmat( xdat, xdim, ydat, wdat, ndat, fpar, epar/*, npar */) ;
00919
00920 /*-------------------------------------------------------------
00921 * here we decrease labda since we may assume that each iteration
00922 * brings us closer to the answer.
00923 */
00924 if ( labda > LABMIN )
00925 {
00926 labda = labda / LABFAC ; /* decrease labda */
00927 }
00928 r = getvec ( xdat, xdim, ydat, wdat, ndat, fpar, epar, npar ) ;
00929 if ( r ) /* error */
00930 {
00931 return r ;
00932 }
00933
00934 while ( chi1 >= chi2 ) /* interpolation loop */
00935 {
00936 /*-----------------------------------------------------------
00937 * The next statement is based on experience, not on the mathematics
00938 * of the problem. It is assumed that we have reached convergence
00939 * when the pure steepest descent method does not produce a better
00940 * solution.
00941 */
00942 if ( labda > LABMAX ) /* assume solution found */
00943 {
00944 break ;
00945 }
00946 labda = labda * LABFAC ; /* increase mixing parameter */
00947 r = getvec ( xdat, xdim, ydat, wdat, ndat, fpar, epar, npar ) ;
00948 if ( r ) /* error */
00949 {
00950 return r ;
00951 }
00952 }
00953
00954 if ( labda <= LABMAX ) /* save old parameters */
00955 {
00956 for ( i = 0 ; i < *npar ; i++ )
00957 {
00958 fpar[i] = epar[i] ;
00959 }
00960 }
00961 if ( (fabs( chi2 - chi1 ) <= (tolerance * chi1)) || (labda > LABMAX) )
00962 {
00963 /*---------------------------------------------------------------------
00964 * we have a satisfying solution, so now we need to calculate the correct
00965 * errors of the fitted parameters. This we do by using the pure Taylor
00966 * method because we are very close to the real solution.
00967 */
00968 labda = 0.0 ; /* for Taylor solution */
00969 getmat ( xdat, xdim, ydat, wdat, ndat, fpar, epar/*, npar */) ;
00970 r = getvec ( xdat, xdim, ydat, wdat, ndat, fpar, epar, npar ) ;
00971
00972 if ( r ) /* error */
00973 {
00974 return r ;
00975 }
00976 for ( i = 0 ; i < (*npar) ; i++ )
00977 {
00978 epar[i] = 0.0 ; /* set error to zero */
00979 }
00980 chi2 = sqrt ( chi2 / (double) (nuse - nfree) ) ;
00981
00982 for ( i = 0 ; i < nfree ; i++ )
00983 {
00984 if ( (matrix1[i][i] <= 0.0) || (matrix2[i][i] <= 0.0) )
00985 {
00986 return -7 ;
00987 }
00988 epar[parptr[i]] = chi2 * sqrt( matrix2[i][i] ) / sqrt( matrix1[i][i] ) ;
00989 }
00990 found = 1 ; /* we found a solution */
00991 }
00992 }
00993 }
00994 return itc ; /* return number of iterations */
00995 }
00996
00997 /*----------------------------------------------------------------------------
00998 Function : nint()
00999 In : x : specified real value
01000 Out : nearest integer to specified real value
01001 Job : this routine determines the nearest integer to a
01002 specified real value
01003 ---------------------------------------------------------------------------*/
01004
01005 int nint ( double x )
01006 {
01007 int k ;
01008
01009 k = x ;
01010 if ( x >= 0. )
01011 {
01012 if ( (x - (double) k) <= 0.5 )
01013 {
01014 return k ;
01015 }
01016 else
01017 {
01018 return k + 1 ;
01019 }
01020 }
01021 else
01022 {
01023 if ( (x - (double) k) <= -0.5 )
01024 {
01025 return k - 1;
01026 }
01027 else
01028 {
01029 return k ;
01030 }
01031 }
01032 }
01033
01034 /*----------------------------------------------------------------------------
01035 Function : correlation()
01036 In : data1: first data array
01037 data2: second data array
01038 ndata: number of data elements in the arrays,
01039 both arrays must have the same number of
01040 elements.
01041 Out : integer shift of the second array with respect to
01042 the first array.
01043 2^32/2 if something went wrong.
01044 Job : this routine computes the cross correlation between
01045 two data arrays of the same size and returns the
01046 integer shift between both arrays
01047 ---------------------------------------------------------------------------*/
01048
01049 int correlation ( float * data1, float * data2, int ndata )
01050 {
01051 /*float help[3*ndata] ;
01052 float corsum[3*ndata] ;*/
01053 float help[ndata+300] ;
01054 float corsum[ndata+300] ;
01055 float maxval ;
01056 int i, j, k, position, shift ;
01057 int /*start,end,size,*/ndata3,limit;
01058
01059
01060 /*ndata3=3*ndata;*/
01061 ndata3=ndata+300;
01062
01063 if ( NULL == data1 || NULL == data2 || ndata <= 1 )
01064 {
01065 cpl_msg_error("correlation:"," wrong input for correlation\n") ;
01066 return INT32_MAX ;
01067 }
01068
01069 /* initialize the help arrays with zeros */
01070 for ( i = 0 ; i < ndata3 ; i++ )
01071 {
01072 help[i] = 0. ;
01073 }
01074
01075 /* shift the second data array by ndata in the help array */
01076 for ( i = 0 ; i < ndata ; i++ )
01077 {
01078 help[(300/2) + i] = data2[i] ;
01079 }
01080
01081 /* compute the cross correlation sum array */
01082 for ( j = 0 ; j < ndata3 ; j++ )
01083 {
01084 if ( ndata3-j <= ndata)
01085 limit = ndata3-j;
01086 else
01087 limit = ndata;
01088 corsum[j] = 0. ;
01089 for ( k = 0 ; k < limit ; k++ )
01090 {
01091 /*if ( k + j >= ndata3 )
01092 {
01093 break ;
01094 }*/
01095 corsum[j] += data1[k] * help[k + j] ;
01096 }
01097 }
01098
01099 /* search for the maximal corsum value and determine its position */
01100 maxval = -FLT_MAX ;
01101 position = -1 ;
01102 for ( i = 0 ; i < ndata3 ; i++ )
01103 {
01104 if ( maxval < corsum[i] )
01105 {
01106 maxval = corsum[i] ;
01107 position = i ;
01108 }
01109 }
01110
01111 /* determine shift of data2 relative to the data1 array */
01112 shift = position - 300/2 ;
01113 return shift ;
01114 }
01115
01116 /*----------------------------------------------------------------------------
01117 Function : clean_mean()
01118 In : array: data array to average
01119 n_elements: number of elements of the data array
01120 throwaway_low: percentage of low value elements to be thrown away before
01121 averaging
01122 throwaway_high: percentage of high value elements to be thrown away before
01123 averaging
01124 Out : the clean mean of a data array
01125 FLT_MAX in case of error
01126 Job : this routine computes the clean mean of a given data
01127 array that means the array is first sorted and
01128 a given percentage of the lowest and the highest values
01129 is not considered for averaging
01130 ---------------------------------------------------------------------------*/
01131 float clean_mean( float * array,
01132 int n_elements,
01133 float throwaway_low,
01134 float throwaway_high )
01135 {
01136 int i, n ;
01137 int lo_n, hi_n ;
01138 float sum ;
01139
01140 if ( array == NULL )
01141 {
01142 cpl_msg_error("clean_mean:"," no array given in clean_mean!\n") ;
01143 return FLT_MAX ;
01144 }
01145
01146 if ( n_elements <= 0 )
01147 {
01148 cpl_msg_error("clean_mean:","wrong number of elements given\n") ;
01149 return FLT_MAX ;
01150 }
01151
01152 if ( throwaway_low < 0. || throwaway_high < 0. ||
01153 throwaway_low + throwaway_high >= 100. )
01154 {
01155 cpl_msg_error("clean_mean:","wrong throw away percentage given!\n") ;
01156 return FLT_MAX ;
01157 }
01158
01159 lo_n = (int) (throwaway_low * (float)n_elements / 100.) ;
01160 hi_n = (int) (throwaway_high * (float)n_elements / 100.) ;
01161
01162 /* sort the array */
01163 pixel_qsort( array, n_elements ) ;
01164
01165 n = 0 ;
01166 sum = 0. ;
01167 for ( i = lo_n ; i < n_elements - hi_n ; i++ )
01168 {
01169 if ( !qfits_isnan(array[i]) )
01170 {
01171 sum += array[i] ;
01172 n++ ;
01173 }
01174 }
01175 if ( n == 0 )
01176 {
01177 return FLAG ;
01178 }
01179 else
01180 {
01181 return sum/(float)n ;
01182 }
01183 }
01184
01185 /*----------------------------------------------------------------------------
01186 Function : errorPrint()
01187 In : text: text for cpl_msg_error to print out
01188 Out : nothing
01189 Job : this routine gives python the chance to print out
01190 error messages
01191 ---------------------------------------------------------------------------*/
01192 void errorPrint(char * text)
01193 {
01194 if ( text == NULL )
01195 {
01196 cpl_msg_error ("errorPrint:","no string given\n") ;
01197 return ;
01198 }
01199 cpl_msg_error ( "errorPrint:","%s",text ) ;
01200 }
01201
01202 /*----------------------------------------------------------------------------
01203 Function : commentPrint()
01204 In : text: text for e_comment to print out
01205 level: comment indentation level
01206 Out : nothing
01207 Job : this routine gives python the chance to print out
01208 comment messages
01209 ---------------------------------------------------------------------------*/
01210 void commentPrint(char * text)
01211 {
01212 if ( text == NULL )
01213 {
01214 cpl_msg_error ("errorPrint:","no string given\n") ;
01215 return ;
01216 }
01217 cpl_msg_info ("commentPrint", text ) ;
01218 }
01219
01220 /*----------------------------------------------------------------------------
01221 Function : convertZEROsTo0ForImages()
01222 In : input image
01223 Out : nothing
01224 Job : this function converts any ZEROs in an image to 0
01225 ---------------------------------------------------------------------------*/
01226 void convertZEROsTo0ForImages(OneImage * im)
01227 {
01228 int i ;
01229 if ( im == NullImage )
01230 {
01231 cpl_msg_error("convertZEROsTo0ForImages:","no input image given!\n") ;
01232 return ;
01233 }
01234 for ( i = 0 ; i < (int) im -> nbpix ; i++ )
01235 {
01236 if( qfits_isnan(im -> data[i]) )
01237 {
01238 im -> data[i] = 0. ;
01239 }
01240 }
01241 return ;
01242 }
01243
01244 /*----------------------------------------------------------------------------
01245 Function : convertZEROsTo0ForCubes()
01246 In : input cube
01247 Out : nothing
01248 Job : this function converts any ZEROs in a cube to 0
01249 ---------------------------------------------------------------------------*/
01250 void convertZEROsTo0ForCubes(OneCube * cube)
01251 {
01252 int i ;
01253 if ( cube == NullCube )
01254 {
01255 cpl_msg_error("convertZEROsTo0ForCubes:","no input cube given!\n") ;
01256 return ;
01257 }
01258 for ( i = 0 ; i < cube -> np ; i++ )
01259 {
01260 convertZEROsTo0ForImages(cube->plane[i]) ;
01261 }
01262 return ;
01263 }
01264
01265
01266 /*----------------------------------------------------------------------------
01267 Function : convertZEROsTo0ForCubesRange()
01268 In : input cube
01269 Out : nothing
01270 Job : this function converts any ZEROs in a cube to 0
01271 ---------------------------------------------------------------------------*/
01272 void convertZEROsTo0ForCubesRange(OneCube * cube,const int z_min,const int z_max)
01273 {
01274 int i ;
01275 if ( cube == NullCube )
01276 {
01277 cpl_msg_error("convertZEROsTo0ForCubes:","no input cube given!\n") ;
01278 return ;
01279 }
01280 for ( i = z_min ; i < z_max ; i++ )
01281 {
01282 convertZEROsTo0ForImages(cube->plane[i]) ;
01283 }
01284 return ;
01285 }
01286 /*----------------------------------------------------------------------------
01287 Function : convert0ToZEROsForImages()
01288 In : input image
01289 Out : nothing
01290 Job : this function converts any 0 in an image to ZEROs
01291 ---------------------------------------------------------------------------*/
01292 void convert0ToZEROsForImages(OneImage * im)
01293 {
01294 int i ;
01295 if ( im == NullImage )
01296 {
01297 cpl_msg_error("convert0ToZEROsForImages:","no input image given!\n") ;
01298 return ;
01299 }
01300 for ( i = 0 ; i < (int) im -> nbpix ; i++ )
01301 {
01302 if( im -> data[i] == 0. )
01303 {
01304 im -> data[i] = ZERO ;
01305 }
01306 }
01307 return ;
01308 }
01309
01310 /*----------------------------------------------------------------------------
01311 Function : convert0ToZEROsForCubes()
01312 In : input cube
01313 Out : nothing
01314 Job : this function converts any 0 in a cube to ZERO
01315 ---------------------------------------------------------------------------*/
01316 void convert0ToZEROForCubes(OneCube * cube)
01317 {
01318 int i ;
01319 if ( cube == NullCube )
01320 {
01321 cpl_msg_error("convert0ToZEROsForCubes:","no input cube given!\n") ;
01322 return ;
01323 }
01324 for ( i = 0 ; i < cube -> np ; i++ )
01325 {
01326 convert0ToZEROsForImages(cube->plane[i]) ;
01327 }
01328 return ;
01329 }
01330
01331
01332 /*----------------------------------------------------------------------------
01333 Function : convert0ToZEROsForCubes()
01334 In : input cube
01335 z_min min z to be considered
01336 z_max max z to be considered
01337 Out : nothing
01338 Job : this function converts any 0 in a cube to ZERO
01339 ---------------------------------------------------------------------------*/
01340 void convert0ToZEROForCubesRange(OneCube * cube,const int z_min,const int z_max)
01341 {
01342 int i ;
01343 if ( cube == NullCube )
01344 {
01345 cpl_msg_error("convert0ToZEROsForCubes:","no input cube given!\n") ;
01346 return ;
01347 }
01348 for ( i = z_min ; i < z_max ; i++ )
01349 {
01350 convert0ToZEROsForImages(cube->plane[i]) ;
01351 }
01352 return ;
01353 }
01354 /*----------------------------------------------------------------------------
01355 Function : invertImage()
01356 In : input image
01357 Out : inverted input image (overwritten)
01358 Job : this function multiplies each image value by -1.
01359 ---------------------------------------------------------------------------*/
01360 void invert(OneImage * im)
01361 {
01362 int i ;
01363 for ( i = 0 ; i < (int) im->nbpix ; i++ )
01364 {
01365 im->data[i] = -im->data[i] ;
01366 }
01367 return ;
01368 }
01369
01370 /*----------------------------------------------------------------------------
01371 Function : xcorrel ()
01372 In : line_i: the reference signal
01373 width_i: number of samples in reference signal
01374 line_t: candidate signal to compare
01375 width_t: number of samples in candidate signal
01376 half_search: half-size of the search domain
01377 delta: output correlation offset.
01378 maxpos: position index of the maximum of the output array
01379 xcorr_max: maximum value of the output array
01380 Out : correlation function, must be freed.
01381 Job : simple cross-correlation of two 1d-signals without FFT
01382 ---------------------------------------------------------------------------*/
01383 #define STEP_MIN (-half_search)
01384 #define STEP_MAX (half_search)
01385
01386 double * xcorrel(
01387 float * line_i,
01388 int width_i,
01389 float * line_t,
01390 int width_t,
01391 int half_search,
01392 int * delta,
01393 int * maxpos,
01394 double * xcorr_max
01395
01396 )
01397 {
01398 double * xcorr ;
01399 double mean_i, mean_t ;
01400 double rms_i, rms_t ;
01401 double sum, sqsum ;
01402 double norm ;
01403 int nsteps ;
01404 int i ;
01405 int step ;
01406 int nval ;
01407 /*double r;*/
01408
01409
01410 /* Compute normalization factors */
01411 sum = sqsum = 0.00 ;
01412 for (i=0 ; i<width_i ; i++) {
01413 sum += (double)line_i[i] ;
01414 sqsum += (double)line_i[i] * (double)line_i[i];
01415 }
01416 mean_i = sum / (double)width_i ;
01417 sqsum /= (double)width_i ;
01418 rms_i = sqsum - mean_i*mean_i ;
01419
01420 sum = sqsum = 0.00 ;
01421 for (i=0 ; i<width_t ; i++) {
01422 sum += (double)line_t[i] ;
01423 sqsum += (double)line_t[i] * (double)line_t[i];
01424 }
01425 mean_t = sum / (double)width_t ;
01426 sqsum /= (double)width_t ;
01427 rms_t = sqsum - mean_t*mean_t ;
01428
01429 norm = 1.00 / sqrt(rms_i * rms_t);
01430
01431 nsteps = (STEP_MAX - STEP_MIN) ;
01432 xcorr = cpl_malloc(nsteps * sizeof(double));
01433 for (step=STEP_MIN ; step<STEP_MAX ; step++) {
01434 xcorr[step-STEP_MIN] = 0.00 ;
01435 nval = 0 ;
01436 for (i=0 ; i<width_t ; i++) {
01437 if ((i+step >= 0) &&
01438 (i+step < width_i)) {
01439 xcorr[step-STEP_MIN] += ((double)line_t[i] - mean_t) *
01440 ((double)line_i[i+step] - mean_i) *
01441 norm ;
01442 nval++ ;
01443 }
01444 }
01445 xcorr[step-STEP_MIN] /= (double)nval ;
01446 }
01447 *xcorr_max = xcorr[0] ;
01448 *maxpos = 0 ;
01449 for (i=0 ; i<nsteps ; i++) {
01450 if (xcorr[i]>*xcorr_max) {
01451 *maxpos = i ;
01452 *xcorr_max = xcorr[i];
01453 }
01454 }
01455 (*delta) = + (STEP_MIN + *maxpos);
01456 return xcorr ;
01457 }
01458
01459 /* FILE ELEMENT: nev_ille.c */
01460 /* */
01461 /**********************************************************************/
01462 /* */
01463 /* double nev_ille() */
01464 /* */
01465 /**********************************************************************/
01466 /* */
01467 /* DESCRIPTION: */
01468 /* For a given table (x , f(x )), i = 0(1)n and a given argument z */
01469 /* the function computes the interpolated value for the argument z */
01470 /* using Neville's interpolation/ extrapolation algorithm. */
01471 /* */
01472 /* FUNCTIONS CALLED: */
01473 /* System library: <stdio.h> printf(), fabs(); */
01474 /* Numlib library: None */
01475 /* Local functions: nevtable(); */
01476 /* User supplied: None */
01477 /* */
01478 /* PROGRAMMED BY: T.Haavie */
01479 /* DATE/VERSION: 88-07-06/1.0 */
01480 /* */
01481 /**********************************************************************/
01482
01483 float nev_ille(float x[], float f[], int n, float z, int* flag)
01484 /* PARAMETERS(input): */
01485 /* float x[]; Abscissae values in the table. */
01486 /* float f[]; Function values in the table. */
01487 /* int n; The number of elements in the table is n+1. */
01488 /* float z; Argument to be used in interpolation/extrapolation. */
01489
01490
01491 /* PARAMETERS(input/output): */
01492 /* int *flag; Flag parameter(output): */
01493 /* = 0, n < 0 and/or eps < 0, should be positive+. */
01494 /* = 1, required rel.err. is not obtained. */
01495 /* = 2, required rel. err. is obtained. */
01496
01497 /* the computed estimate for the interpolated/extrapolated value re- */
01498 /* turned through function name nev_ille. */
01499
01500 {
01501 float p[11]; /* Array used for storing the new row elements */
01502 /* in the interpolation/extrapolation table. */
01503 float q[11]; /* Array used for storing the old row elements */
01504 /* in the interpolation/extrapolation table */
01505
01506 float factor;
01507 int m, k;
01508
01509
01510
01511 if (n < 0 )
01512 {
01513 *flag = 0;
01514 return(0.);
01515 }
01516
01517
01518 q[0] = f[0]; /* Set initial value in the table. */
01519
01520 for (k = 1; k <= n; k++) /* k counts rows in the table. */
01521 {
01522 p[0] = f[k];
01523 for (m = 1; m <= k; m++) /* m counts element in row. */
01524 {
01525 factor = (z - x[k]) / (x[k] - x[k-m]);
01526 p[m] = p[m-1] + factor * (p[m-1] - q[m-1]);
01527 }
01528
01529
01530 for (m = 0; m <= k; m++) /* Shift old row to new row. */
01531 q[m] = p[m];
01532
01533 } /* End of k-loop. */
01534
01535 *flag = 1; /* Required rel.error is not obtained. */
01536 return(p[n]);
01537
01538 } /* End of nev_ille(). */
01539
01540 /*--------------------------------------------------------------------------*/
1.2.13.1 written by Dimitri van Heesch,
© 1997-2001