recipes.c

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 /*--------------------------------------------------------------------------*/

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