/* @(#)mutil.c 17.1.1.1 (ES0-DMD) 01/25/02 17:55:51 */ /*=========================================================================== Copyright (C) 1995 European Southern Observatory (ESO) This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, MA 02139, USA. Corresponding concerning ESO-MIDAS should be addressed as follows: Internet e-mail: midas@eso.org Postal address: European Southern Observatory Data Management Division Karl-Schwarzschild-Strasse 2 D 85748 Garching bei Muenchen GERMANY ===========================================================================*/ /* @(#)mutil.c 17.1.1.1 (ESO) 01/25/02 17:55:51 */ /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ /* .COPYRIGHT (C) 1993 European Southern Observatory */ /* .IDENT .prg */ /* .AUTHORS Pascal Ballester (ESO/Garching) */ /* Cristian Levin (ESO/La Silla) */ /* .KEYWORDS Spectroscopy, Long-Slit */ /* .PURPOSE */ /* .VERSION 1.0 Package Creation 17-MAR-1993 */ /* ------------------------------------------------------- */ /* @(#)mutil.c 1.0.0.0 (ESO-La Silla) 10/08/91 12:00:00 */ /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ .IDENT mutil.c .MODULE subroutines .LANGUAGE C .AUTHOR Cristian Levin - ESO La Silla .PURPOSE Miscelanous mathematical algorithms: - linear square fitting. - fast polynomial evaluation. - search of the median. - polynomial interpolation. .KEYWORDS interpolation, fitting, polynomial evaluation. .COMMENTS Most of this routines were taken as they are from the book "Numerical Recipes in C" -- 1st edition. .VERSION 1.0 1-Apr-1991 Implementation .ENVIRONMENT UNIX ------------------------------------------------------------*/ #include #include /****************************************************** * * lfit(): general linear least squares solution. * ******************************************************/ void lfit( x, y, ndata, a, mfit, funcs ) int ndata, mfit; /* IN: no. of values, no. of coefficients */ double x[], y[], a[]; /* IN: (Xi,Yi) OUT: (Ai) */ void (*funcs)(); { int k, j, i; double **covar, **beta, *afunc; void gaussj(), free_dmatrix(), free_dvector(); double **dmatrix(), *dvector(); beta = dmatrix( 1, mfit, 1, 1 ); covar = dmatrix( 1, mfit, 1, mfit ); afunc = dvector( 1, mfit ); for ( j = 1; j <= mfit; j++ ) { for ( k = 1; k <= mfit; k++ ) covar[j][k] = 0.0; beta[j][1] = 0.0; } for ( i = 1; i <= ndata; i++ ) { (*funcs)( x[i], afunc, mfit ); for ( j = 1; j <= mfit; j++ ) { for ( k = 1; k <= j; k++ ) covar[j][k] += afunc[j] * afunc[k]; beta[j][1] += y[i] * afunc[j]; } } for ( j = 2; j <= mfit; j++ ) for ( k = 1; k <= j-1; k++ ) covar[k][j] = covar[j][k]; gaussj( covar, mfit, beta, 1 ); /* Matrix solution */ for ( j = 1; j <= mfit; j++ ) /* solution to coefficients a */ a[j] = beta[j][1]; free_dvector( afunc, 1, mfit ); free_dmatrix( beta, 1, mfit, 1, 1 ); free_dmatrix( covar, 1, mfit, 1, mfit ); } /****************************************************************** * * gaussj(): linear equation solution by Gauss-Jordan elimination * ******************************************************************/ #define SWAP(a, b) {float temp = (a); (a) = (b); (b) = temp; } void gaussj( a, n, b, m ) double **a, /* IN : a[1..n][1..n] */ **b; /* IN/OUT: b[1..n][1..m] */ int n, m; { int *indxc, *indxr, *ipiv; /* for bookeeping in the pivoting */ int i, icol, irow, j, k; int *ivector(); double big, dum, pivinv; void nrerror(), free_ivector(); indxc = ivector( 1, n ); indxr = ivector( 1, n ); ipiv = ivector( 1, n ); for ( i = 1; i <= n; i++ ) ipiv[i] = 0.0; for ( i = 1; i <= n; i++ ) { /* main loop over columns to be reduced */ big = 0.0; for ( j = 1; j <= n; j++ ) if ( ipiv[j] != 1 ) for ( k = 1; k <= n; k++ ) { if ( ipiv[k] == 0 ) { if ( fabs(a[j][k]) >= big ) { big = fabs(a[j][k]); irow = j; icol = k; } } else if ( ipiv[k] > 1 ) { nrerror( "Tolerance too small for fitting" ); return; } } ++(ipiv[icol]); if ( irow != icol ) { for ( j = 1; j <= n; j++ ) SWAP( a[irow][j], a[icol][j] ); for ( j = 1; j <= m; j++ ) SWAP( b[irow][j], b[icol][j] ); } indxr[i] = irow; indxc[i] = icol; if ( a[icol][icol] == 0.0 ) { nrerror( "Tolerance too small for fitting" ); return; } pivinv = 1.0 / a[icol][icol]; a[icol][icol] = 1.0; for ( j = 1; j <= n; j++ ) a[icol][j] *= pivinv; for ( j = 1; j <= m; j++ ) b[icol][j] *= pivinv; for ( j = 1; j <= n; j++ ) if ( j != icol ) { dum = a[j][icol]; a[j][icol] = 0.0; for ( k = 1; k <= n; k++ ) a[j][k] -= a[icol][k] * dum; for ( k = 1; k <= m; k++ ) b[j][k] -= b[icol][k] * dum; } } /* main loop */ for ( i = n; i >= 1; i-- ) if ( indxr[i] != indxc[i] ) for ( k = 1; k <= n; k++ ) SWAP( a[k][indxr[i]], a[k][indxc[i]] ); free_ivector( ipiv, 1, n ); free_ivector( indxc, 1, n ); free_ivector( indxr, 1, n ); } /******************************************************************** * * polint(): polynomial interpolation or extrapolation from N input * points, using the Neville's algorithm. * ********************************************************************/ void polint( xa, ya, n, x, y, dy ) float xa[], ya[], x, *y, *dy; int n; { int i, m, ns = 1; float den, dif, dift, ho, hp, w; float *c, *d, *fvector(); void nrerror(), free_vector(); dif = fabs(x - xa[1]); c = fvector(1, n); d = fvector(1, n); for ( i = 1; i <= n; i++ ) { /* find the index ns of the closest */ if ( (dift = fabs(x - xa[i])) < dif ) { /* table entry. */ ns = i; dif = dift; } c[i] = ya[i]; /* initialize the tableau of c's and d's */ d[i] = ya[i]; } *y = ya[ns--]; /* initial aprox. of y */ for ( m = 1; m < n; m++ ) { for ( i = 1; i <= n-m; i++ ) { ho = xa[i] - x; hp = xa[i+m] - x; w = c[i+1] - d[i]; if ( (den = ho - hp) == 0.0 ) { nrerror( "Error in polynomial interpolation" ); return; } den = w / den; d[i] = hp * den; c[i] = ho * den; } *y += (*dy = (2*ns < (n-m) ? c[ns+1] : d[ns--])); } fvector(c, 1, n); fvector(d, 1, n); } /************************************************************ * * eval_fpoly(), eval_dpoly(): fast polynomial evaluation. * ************************************************************/ float eval_fpoly(x, a, n) float x; float *a; /* coefficients (range: 1..n) */ int n; { int i; float y = 0.0; for ( i = n; i >= 1; i-- ) y = x*y + a[i]; return(y); } double eval_dpoly(x, a, n) double x; double *a; /* coefficients (range: 1..n) */ int n; { int i; double y = 0.0; for ( i = n; i >= 1; i-- ) y = x*y + a[i]; return(y); } /************************************************************ * * median(): efficient search of the median. * ************************************************************/ float median(x, n) float x[]; int n; { int n2 ,n2p; void piksrt(); piksrt(n, x); n2p = (n2 = n/2) + 1; return((n % 2) ? x[n2p] : 0.5 * (x[n2] + x[n2p])); } /************************************************************ * * piksrt(): sorting by straight insertion. * ************************************************************/ void piksrt(n,arr) int n; float arr[]; { int i,j; float a; for (j=2;j<=n;j++) { a=arr[j]; i=j-1; while (i > 0 && arr[i] > a) { arr[i+1]=arr[i]; i--; } arr[i+1]=a; } }