/* @(#)mutil.c 8.9 (ESO-IPG) 12/01/94 10:38:21 */ /* @(#)mutil.c 1.2 (ESO-IPG) 2/11/94 09:31:25 */ /*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ .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 #include #include /****************************************************** * lfit(): general linear least squares solution. * ******************************************************/ #ifdef __STDC__ void lfit (double x[], double y[], int ndata, double a[], int mfit, void (*funcs) (double, double *afunc, int mfit)) #else void lfit (x, y, ndata, a, mfit, funcs) double x[], y[], a[]; int ndata, mfit; void (*funcs) (); #endif { int k, j, i; double **covar, **beta, *afunc; #ifdef __STDC__ void gaussj (double **, int, double **, int); #else void gaussj (), free_matrix (), free_dvector (); #endif 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 = 0; 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; } #ifdef __STDC__ void gaussj (double **a, int n, double **b, int m) #else void gaussj (a, n, b, m) double **a, **b; int m, n; #endif /* IN : a[1..n][1..n] */ /* IN/OUT: b[1..n][1..m] */ { int *indxc, *indxr, *ipiv; /* for bookeeping in the pivoting */ int i, icol, irow, j, k; double big, dum, pivinv; 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 ("mutil.c -> gaussj: Singular Matrix -- ipiv[k]>1"); 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 ("mutil.c -> gaussj: Singular Matrix -- a[icol][icol] == 0.0"); 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. * ********************************************************************/ #ifdef __STDC__ void polint (float xa[], float ya[], int n, float x, float *y, float *dy) #else void polint (xa, ya, n, x, y, dy) float xa[], ya[], x, *y, *dy; int n; #endif { int i, m, ns = 1; float den, dif, dift, ho, hp, w; float *c, *d; 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--])); } free_fvector (c, 1, n); free_fvector (d, 1, n); } /************************************************************ * eval_fpoly(), eval_dpoly(): fast polynomial evaluation. * ************************************************************/ #ifdef __STDC__ float eval_fpoly (float x, float *a, int n) #else float eval_fpoly (x, a, n) float x, *a; int n; #endif /* *a: coefficients (range: 1..n) */ { int i; float y = 0.0; for (i = n; i >= 1; i--) y = x * y + a[i]; return (y); } #ifdef __STDC__ double eval_dpoly (double x, double *a, int n) #else double eval_dpoly (x, a, n) double x, *a; int n; #endif /* *a: coefficients (range: 1..n) */ { int i; double y = 0.0; for (i = n; i >= 1; i--) y = x * y + a[i]; return (y); } #ifdef __STDC__ void dpoly (double x, double p[], int np) #else void dpoly (x, p, np) double x, p[]; int np; #endif { int j; p[1] = 1.0; for (j = 2; j <= np; j++) p[j] = p[j - 1] * x; } /************************************************************ * median(): efficient search of the median. * ************************************************************/ #ifdef __STDC__ float median (float x[], int n) #else float median (x, n) float x[]; int n; #endif { int n2, n2p; #ifdef __STDC__ void piksrt (int, float[]); #else void piksrt (); #endif piksrt (n, x); n2p = (n2 = n / 2) + 1; return ((n % 2) ? x[n2p] : 0.5 * (x[n2] + x[n2p])); } /************************************************************ * piksrt(): sorting by straight insertion. * ************************************************************/ #ifdef __STDC__ void piksrt (int n, float arr[]) #else void piksrt (n, arr) float arr[]; int n; #endif { 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; } } #ifdef __STDC__ void dsvdcmp (double **a, int m, int n, double w[], double **v) #else void dsvdcmp (a, m, n, w, v) double **a, w[], **v; int m, n; #endif { #ifdef __STDC__ double dpythag (double, double); #else double dpythag (); #endif int flag, i, its, j, jj, k, l, nm; double anorm, c, f, g, h, s, scale, x, y, z, *rv1; rv1 = dvector (1, n); g = scale = anorm = 0.0; for (i = 1; i <= n; i++) { l = i + 1; rv1[i] = scale * g; g = s = scale = 0.0; if (i <= m) { for (k = i; k <= m; k++) scale += fabs (a[k][i]); /* if (scale) */ if ((int) scale) { for (k = i; k <= m; k++) { a[k][i] /= scale; s += a[k][i] * a[k][i]; } f = a[i][i]; g = -SIGN (sqrt (s), f); h = f * g - s; a[i][i] = f - g; for (j = l; j <= n; j++) { for (s = 0.0, k = i; k <= m; k++) s += a[k][i] * a[k][j]; f = s / h; for (k = i; k <= m; k++) a[k][j] += f * a[k][i]; } for (k = i; k <= m; k++) a[k][i] *= scale; } } w[i] = scale * g; g = s = scale = 0.0; if (i <= m && i != n) { for (k = l; k <= n; k++) scale += fabs (a[i][k]); /* if (scale) */ if ((int) scale) { for (k = l; k <= n; k++) { a[i][k] /= scale; s += a[i][k] * a[i][k]; } f = a[i][l]; g = -SIGN (sqrt (s), f); h = f * g - s; a[i][l] = f - g; for (k = l; k <= n; k++) rv1[k] = a[i][k] / h; for (j = l; j <= m; j++) { for (s = 0.0, k = l; k <= n; k++) s += a[j][k] * a[i][k]; for (k = l; k <= n; k++) a[j][k] += s * rv1[k]; } for (k = l; k <= n; k++) a[i][k] *= scale; } } anorm = DMAX (anorm, (fabs (w[i]) + fabs (rv1[i]))); } for (i = n; i >= 1; i--) { if (i < n) { /* if (g) */ if ((int) g) { for (j = l; j <= n; j++) v[j][i] = (a[i][j] / a[i][l]) / g; for (j = l; j <= n; j++) { for (s = 0.0, k = l; k <= n; k++) s += a[i][k] * v[k][j]; for (k = l; k <= n; k++) v[k][j] += s * v[k][i]; } } for (j = l; j <= n; j++) v[i][j] = v[j][i] = 0.0; } v[i][i] = 1.0; g = rv1[i]; l = i; } for (i = IMIN (m, n); i >= 1; i--) { l = i + 1; g = w[i]; for (j = l; j <= n; j++) a[i][j] = 0.0; if ((int) g) { g = 1.0 / g; for (j = l; j <= n; j++) { for (s = 0.0, k = l; k <= m; k++) s += a[k][i] * a[k][j]; f = (s / a[i][i]) * g; for (k = i; k <= m; k++) a[k][j] += f * a[k][i]; } for (j = i; j <= m; j++) a[j][i] *= g; } else for (j = i; j <= m; j++) a[j][i] = 0.0; /* ++a[i][i]; */ a[i][i] = a[i][i]+1; } for (k = n; k >= 1; k--) { for (its = 1; its <= 30; its++) { flag = 1; for (l = k; l >= 1; l--) { nm = l - 1; if ((double) (fabs (rv1[l]) + anorm) == anorm) { flag = 0; break; } if ((double) (fabs (w[nm]) + anorm) == anorm) break; } if (flag) { c = 0.0; s = 1.0; for (i = l; i <= k; i++) { f = s * rv1[i]; rv1[i] = c * rv1[i]; if ((double) (fabs (f) + anorm) == anorm) break; g = w[i]; h = dpythag (f, g); w[i] = h; h = 1.0 / h; c = g * h; s = -f * h; for (j = 1; j <= m; j++) { y = a[j][nm]; z = a[j][i]; a[j][nm] = y * c + z * s; a[j][i] = z * c - y * s; } } } z = w[k]; if (l == k) { if (z < 0.0) { w[k] = -z; for (j = 1; j <= n; j++) v[j][k] = -v[j][k]; } break; } if (its == 30) nrerror ("no convergence in 30 dsvdcmp iterations"); x = w[l]; nm = k - 1; y = w[nm]; g = rv1[nm]; h = rv1[k]; f = ((y - z) * (y + z) + (g - h) * (g + h)) / (2.0 * h * y); g = dpythag (f, 1.0); f = ((x - z) * (x + z) + h * ((y / (f + SIGN (g, f))) - h)) / x; c = s = 1.0; for (j = l; j <= nm; j++) { i = j + 1; g = rv1[i]; y = w[i]; h = s * g; g = c * g; z = dpythag (f, h); rv1[j] = z; c = f / z; s = h / z; f = x * c + g * s; g = g * c - x * s; h = y * s; y *= c; for (jj = 1; jj <= n; jj++) { x = v[jj][j]; z = v[jj][i]; v[jj][j] = x * c + z * s; v[jj][i] = z * c - x * s; } z = dpythag (f, h); w[j] = z; if ((int) z) { z = 1.0 / z; c = f * z; s = h * z; } f = c * g + s * y; x = c * y - s * g; for (jj = 1; jj <= m; jj++) { y = a[jj][j]; z = a[jj][i]; a[jj][j] = y * c + z * s; a[jj][i] = z * c - y * s; } } rv1[l] = 0.0; rv1[k] = f; w[k] = x; } } free_dvector (rv1, 1, n); } #ifdef __STDC__ void dsvbksb (double **u, double w[], double **v, int m, int n, double b[], double x[]) #else void dsvbksb (u, w, v, m, n, b, x) double **u, w[], **v, b[], x[]; int m, n; #endif { int jj, j, i; double s, *tmp; tmp = dvector (1, n); for (j = 1; j <= n; j++) { s = 0.0; if ((int) w[j]) { for (i = 1; i <= m; i++) s += u[i][j] * b[i]; s /= w[j]; } tmp[j] = s; } for (j = 1; j <= n; j++) { s = 0.0; for (jj = 1; jj <= n; jj++) s += v[j][jj] * tmp[jj]; x[j] = s; } free_dvector (tmp, 1, n); } #ifdef __STDC__ double dpythag (double a, double b) #else double dpythag (a, b) double a, b; #endif { double absa, absb; absa = fabs (a); absb = fabs (b); if (absa > absb) return absa * sqrt (1.0 + DSQR (absb / absa)); else return (absb == 0.0 ? 0.0 : absb * sqrt (1.0 + DSQR (absa / absb))); } #ifdef __STDC__ double select_pos (unsigned long k, unsigned long n, double arr[]) #else double select_pos (k, n, arr) unsigned long k, n; double arr[]; #endif /* select the kth smallest value in arr[1...n], NR, pg. 342 */ #define SWAPD(a, b) {double temp = (a); (a) = (b); (b) = temp; } { unsigned long i, ir, j, l, mid; /* double a, temp; */ double a; l = 1; ir = n; for (;;) { if (ir <= l + 1) { if (ir == l + 1 && arr[ir] < arr[l]) { SWAPD (arr[l], arr[ir]) } return arr[k]; } else { mid = (l + ir) >> 1; SWAPD (arr[mid], arr[l + 1]) if (arr[l + 1] > arr[ir]) { SWAPD (arr[l + 1], arr[ir]) } if (arr[l] > arr[ir]) { SWAPD (arr[l], arr[ir]) } if (arr[l + 1] > arr[l]) { SWAPD (arr[l + 1], arr[l]) } i = l + 1; j = ir; a = arr[l]; for (;;) { do i++; while (arr[i] < a); do j--; while (arr[j] > a); if (j < i) break; SWAPD (arr[i], arr[j]) } arr[l] = arr[j]; arr[j] = a; if (j >= k) ir = j - 1; if (j <= k) l = i; } } } #undef SWAPD #undef SWAP /*---------------------------------------------------------------------------*/ #ifdef __STDC__ float pik_median (int n, float arr[]) #else float pik_median ( n, arr ) float arr[]; int n; #endif /* find median of distribution (use piksort) NR page 330 modified for mos by Sabine ?? input array arr[] is conserved by buffer b[] -- here we start with arr[0,..,n-1] both in contrary to select_pos TS-0896 */ { int i, j, index; float a, b[101]; for (j = 0; j < n; j++) b[j] = arr[j]; for (j = 1; j < n; j++) { a = b[j]; i = j - 1; while (i >= 0 && b[i] > a) { b[i + 1] = b[i]; i--; } b[i + 1] = a; } index = (n - 1) / 2; return(b[index]); } #ifdef __STDC__ float heap_median (int n, float arr[]) #else float heap_median (n, arr) float arr[]; int n; #endif /* find median of distribution (heapsort) */ { int i, j, index; float a, b[101]; for (j = 0; j < n; j++) b[j] = arr[j]; for (j = 1; j < n; j++) { a = b[j]; i = j - 1; while (i >= 0 && b[i] > a) { b[i + 1] = b[i]; i--; } b[i + 1] = a; } index = (n - 1) / 2; return(b[index]); } /*---------------------------------------------------------------------------*/ #ifdef __STDC__ void fit_poly (float inimage[], float outimage[], double start, double step, int npix, int fit_deg) #else void fit_poly ( inimage, outimage, start, step, npix, fit_deg ) double start,step; float inimage[],outimage[]; int npix,fit_deg; #endif /* fit a polynomial (Numerical Recipes) */ { float xout; double *af; double *ad, *xin, *yin; int i; xin = dvector (0, npix - 1); yin = dvector (0, npix - 1); for (i = 0; i < npix; i++) { xin[i] = start + i * step; yin[i] = (double) inimage[i]; } ad = dvector (1, fit_deg); af = dvector (1, fit_deg); lfit (xin, yin, npix, ad, fit_deg, dpoly); for (i = 1; i <= fit_deg; i++) { af[i] = ad[i]; } for (i = 0; i < npix; i++) { xout = start + step * i; outimage[i] = (float) eval_dpoly (xout, ad, fit_deg); } free_dvector(xin, 0, npix - 1), free_dvector(yin, 0, npix - 1); free_dvector(ad, 1, fit_deg), free_dvector(af, 1, fit_deg); }