diff options
Diffstat (limited to 'pkg/xtools/numrecipes.x')
| -rw-r--r-- | pkg/xtools/numrecipes.x | 689 |
1 files changed, 689 insertions, 0 deletions
diff --git a/pkg/xtools/numrecipes.x b/pkg/xtools/numrecipes.x new file mode 100644 index 00000000..ae437b6d --- /dev/null +++ b/pkg/xtools/numrecipes.x @@ -0,0 +1,689 @@ +# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. + +include <math.h> +include <mach.h> + +# GAMMLN -- Return natural log of gamma function. +# POIDEV -- Returns Poisson deviates for a given mean. +# GASDEV -- Return a normally distributed deviate of zero mean and unit var. +# MR_SOLVE -- Levenberg-Marquardt nonlinear chi square minimization. +# MR_EVAL -- Evaluate curvature matrix. +# MR_INVERT -- Solve a set of linear equations using Householder transforms. +# TWOFFT -- Returns the complex FFTs of two input real arrays. +# REALFT -- Calculates the FFT of a set of 2N real valued data points. +# FOUR1 -- Computes the forward or inverse FFT of the input array. + + +# GAMMLN -- Return natural log of gamma function. +# Argument must greater than 0. Full accuracy is obtained for values +# greater than 1. For 0<xx<1, the reflection formula can be used first. +# +# Based on Numerical Recipes by Press, Flannery, Teukolsky, and Vetterling. +# Used by permission of the authors. +# Copyright(c) 1986 Numerical Recipes Software. + +real procedure gammln (xx) + +real xx # Value to be evaluated + +int j +double cof[6], stp, x, tmp, ser +data cof, stp / 76.18009173D0, -86.50532033D0, 24.01409822D0, + -1.231739516D0,.120858003D-2,-.536382D-5,2.50662827465D0/ + +begin + x = xx - 1.0D0 + tmp = x + 5.5D0 + tmp = (x + 0.5D0) * log (tmp) - tmp + ser = 1.0D0 + do j = 1, 6 { + x = x + 1.0D0 + ser = ser + cof[j] / x + } + return (tmp + log (stp * ser)) +end + + +# POIDEV -- Returns Poisson deviates for a given mean. +# The real value returned is an integer. +# +# Based on Numerical Recipes by Press, Flannery, Teukolsky, and Vetterling. +# Used by permission of the authors. +# Copyright(c) 1986 Numerical Recipes Software. +# +# Modified to return zero for input values less than or equal to zero. + +real procedure poidev (xm, seed) + +real xm # Poisson mean +long seed # Random number seed + +real oldm, g, em, t, y, ymin, ymax, sq, alxm, gammln(), urand(), gasdev() +data oldm /-1./ + +begin + if (xm <= 0) + em = 0 + else if (xm < 12) { + if (xm != oldm) { + oldm = xm + g = exp (-xm) + } + em = 0 + for (t = urand (seed); t > g; t = t * urand (seed)) + em = em + 1 + } else if (xm < 100) { + if (xm != oldm) { + oldm = xm + sq = sqrt (2. * xm) + ymin = -xm / sq + ymax = (1000 - xm) / sq + alxm = log (xm) + g = xm * alxm - gammln (xm+1.) + } + repeat { + repeat { + y = tan (PI * urand(seed)) + } until (y >= ymin) + em = int (sq * min (y, ymax) + xm) + t = 0.9 * (1 + y**2) * exp (em * alxm - gammln (em+1) - g) + } until (urand(seed) <= t) + } else + em = xm + sqrt (xm) * gasdev (seed) + return (em) +end + + +# GASDEV -- Return a normally distributed deviate with zero mean and unit +# variance. The method computes two deviates simultaneously. +# +# Based on Numerical Recipes by Press, Flannery, Teukolsky, and Vetterling. +# Used by permission of the authors. +# Copyright(c) 1986 Numerical Recipes Software. + +real procedure gasdev (seed) + +long seed # Seed for random numbers + +real v1, v2, r, fac, urand() +int iset +data iset/0/ + +begin + if (iset == 0) { + repeat { + v1 = 2 * urand (seed) - 1. + v2 = 2 * urand (seed) - 1. + r = v1 ** 2 + v2 ** 2 + } until ((r > 0) && (r < 1)) + fac = sqrt (-2. * log (r) / r) + + iset = 1 + return (v1 * fac) + } else { + iset = 0 + return (v2 * fac) + } +end + + +# MR_SOLVE -- Levenberg-Marquardt nonlinear chi square minimization. +# +# Use the Levenberg-Marquardt method to minimize the chi squared of a set +# of paraemters. The parameters being fit are indexed by the flag array. +# To initialize the Marquardt parameter, MR, is less than zero. After that +# the parameter is adjusted as needed. To finish set the parameter to zero +# to free memory. This procedure requires a subroutine, DERIVS, which +# takes the derivatives of the function being fit with respect to the +# parameters. There is no limitation on the number of parameters or +# data points. For a description of the method see NUMERICAL RECIPES +# by Press, Flannery, Teukolsky, and Vetterling, p523. +# +# These routines have their origin in Numerical Recipes, MRQMIN, MRQCOF, +# but have been completely redesigned. + +procedure mr_solve (x, y, npts, params, flags, np, nfit, mr, chisq) + +real x[npts] # X data array +real y[npts] # Y data array +int npts # Number of data points +real params[np] # Parameter array +int flags[np] # Flag array indexing parameters to fit +int np # Number of parameters +int nfit # Number of parameters to fit +real mr # MR parameter +real chisq # Chi square of fit + +int i +real chisq1 +pointer new, a1, a2, delta1, delta2 + +errchk mr_invert + +begin + # Allocate memory and initialize. + if (mr < 0.) { + call mfree (new, TY_REAL) + call mfree (a1, TY_REAL) + call mfree (a2, TY_REAL) + call mfree (delta1, TY_REAL) + call mfree (delta2, TY_REAL) + + call malloc (new, np, TY_REAL) + call malloc (a1, nfit*nfit, TY_REAL) + call malloc (a2, nfit*nfit, TY_REAL) + call malloc (delta1, nfit, TY_REAL) + call malloc (delta2, nfit, TY_REAL) + + call amovr (params, Memr[new], np) + call mr_eval (x, y, npts, Memr[new], flags, np, Memr[a2], + Memr[delta2], nfit, chisq) + mr = 0.001 + } + + # Restore last good fit and apply the Marquardt parameter. + call amovr (Memr[a2], Memr[a1], nfit * nfit) + call amovr (Memr[delta2], Memr[delta1], nfit) + do i = 1, nfit + Memr[a1+(i-1)*(nfit+1)] = Memr[a2+(i-1)*(nfit+1)] * (1. + mr) + + # Matrix solution. + call mr_invert (Memr[a1], Memr[delta1], nfit) + + # Compute the new values and curvature matrix. + do i = 1, nfit + Memr[new+flags[i]-1] = params[flags[i]] + Memr[delta1+i-1] + call mr_eval (x, y, npts, Memr[new], flags, np, Memr[a1], + Memr[delta1], nfit, chisq1) + + # Check if chisq has improved. + if (chisq1 < chisq) { + mr = 0.1 * mr + chisq = chisq1 + call amovr (Memr[a1], Memr[a2], nfit * nfit) + call amovr (Memr[delta1], Memr[delta2], nfit) + call amovr (Memr[new], params, np) + } else + mr = 10. * mr + + if (mr == 0.) { + call mfree (new, TY_REAL) + call mfree (a1, TY_REAL) + call mfree (a2, TY_REAL) + call mfree (delta1, TY_REAL) + call mfree (delta2, TY_REAL) + } +end + + +# MR_EVAL -- Evaluate curvature matrix. This calls procedure DERIVS. + +procedure mr_eval (x, y, npts, params, flags, np, a, delta, nfit, chisq) + +real x[npts] # X data array +real y[npts] # Y data array +int npts # Number of data points +real params[np] # Parameter array +int flags[np] # Flag array indexing parameters to fit +int np # Number of parameters +real a[nfit,nfit] # Curvature matrix +real delta[nfit] # Delta array +int nfit # Number of parameters to fit +real chisq # Chi square of fit + +int i, j, k +real ymod, dy, dydpj, dydpk +pointer sp, dydp + +begin + call smark (sp) + call salloc (dydp, np, TY_REAL) + + do j = 1, nfit { + do k = 1, j + a[j,k] = 0. + delta[j] = 0. + } + + chisq = 0. + do i = 1, npts { + call derivs (x[i], params, ymod, Memr[dydp], np) + dy = y[i] - ymod + do j = 1, nfit { + dydpj = Memr[dydp+flags[j]-1] + delta[j] = delta[j] + dy * dydpj + do k = 1, j { + dydpk = Memr[dydp+flags[k]-1] + a[j,k] = a[j,k] + dydpj * dydpk + } + } + chisq = chisq + dy * dy + } + + do j = 2, nfit + do k = 1, j-1 + a[k,j] = a[j,k] + + call sfree (sp) +end + + +# MR_INVERT -- Solve a set of linear equations using Householder transforms. +# This calls a routine published in in "Solving Least Squares Problems", +# by Charles L. Lawson and Richard J. Hanson, Prentice Hall, 1974. + +procedure mr_invert (a, b, n) + +real a[n,n] # Input matrix and returned inverse +real b[n] # Input RHS vector and returned solution +int n # Dimension of input matrices + +int krank +real rnorm +pointer sp, h, g, ip + +begin + call smark (sp) + call salloc (h, n, TY_REAL) + call salloc (g, n, TY_REAL) + call salloc (ip, n, TY_INT) + + call hfti (a, n, n, n, b, n, 1, 0.001, krank, rnorm, + Memr[h], Memr[g], Memi[ip]) + + call sfree (sp) +end + + +# TWOFFT - Given two real input arrays DATA1 and DATA2, each of length +# N, this routine calls cc_four1() and returns two complex output arrays, +# FFT1 and FFT2, each of complex length N (i.e. real length 2*N), which +# contain the discrete Fourier transforms of the respective DATAs. As +# always, N must be an integer power of 2. +# +# Based on Numerical Recipes by Press, Flannery, Teukolsky, and Vetterling. +# Used by permission of the authors. +# Copyright(c) 1986 Numerical Recipes Software. + +procedure twofft (data1, data2, fft1, fft2, N) + +real data1[ARB], data2[ARB] # Input data arrays +real fft1[ARB], fft2[ARB] # Output FFT arrays +int N # No. of points + +int nn3, nn2, jj, j +real rep, rem, aip, aim + +begin + nn2 = 2 + N + N + nn3 = nn2 + 1 + + jj = 2 + for (j=1; j <= N; j = j + 1) { + fft1[jj-1] = data1[j] # Pack 'em into one complex array + fft1[jj] = data2[j] + jj = jj + 2 + } + + call four1 (fft1, N, 1) # Transform the complex array + fft2[1] = fft1[2] + fft2[2] = 0.0 + fft1[2] = 0.0 + for (j=3; j <= N + 1; j = j + 2) { + rep = 0.5 * (fft1[j] + fft1[nn2-j]) + rem = 0.5 * (fft1[j] - fft1[nn2-j]) + aip = 0.5 * (fft1[j + 1] + fft1[nn3-j]) + aim = 0.5 * (fft1[j + 1] - fft1[nn3-j]) + fft1[j] = rep + fft1[j+1] = aim + fft1[nn2-j] = rep + fft1[nn3-j] = -aim + fft2[j] = aip + fft2[j+1] = -rem + fft2[nn2-j] = aip + fft2[nn3-j] = rem + } + +end + + +# REALFT - Calculates the Fourier Transform of a set of 2N real valued +# data points. Replaces this data (which is stored in the array DATA) by +# the positive frequency half of it's complex Fourier Transform. The real +# valued first and last components of the complex transform are returned +# as elements DATA(1) and DATA(2) respectively. N must be an integer power +# of 2. This routine also calculates the inverse transform of a complex +# array if it is the transform of real data. (Result in this case must be +# multiplied by 1/N). A forward transform is perform for isign == 1, other- +# wise the inverse transform is computed. +# +# Based on Numerical Recipes by Press, Flannery, Teukolsky, and Vetterling. +# Used by permission of the authors. +# Copyright(c) 1986 Numerical Recipes Software. + +procedure realft (data, N, isign) + +real data[ARB] # Input data array & output FFT +int N # No. of points +int isign # Direction of transfer + +double wr, wi, wpr, wpi, wtemp, theta # Local variables +real c1, c2, h1r, h1i, h2r, h2i +real wrs, wis +int i, i1, i2, i3, i4 +int N2P3 + +begin + # Initialize + theta = PI/double(N) + c1 = 0.5 + + if (isign == 1) { + c2 = -0.5 + call four1 (data,n,1) # Forward transform is here + } else { + c2 = 0.5 + theta = -theta + } + + wtemp = sin (0.5 * theta) + wpr = -2.0d0 * wtemp * wtemp + wpi = dsin (theta) + wr = 1.0D0 + wpr + wi = wpi + n2p3 = 2*n + 3 + + for (i=2; i<=n/2; i = i + 1) { + i1 = 2 * i - 1 + i2 = i1 + 1 + i3 = n2p3 - i2 + i4 = i3 + 1 + wrs = sngl (wr) + wis = sngl (wi) + # The 2 transforms are separated out of Z + h1r = c1 * (data[i1] + data[i3]) + h1i = c1 * (data[i2] - data[i4]) + h2r = -c2 * (data[i2] + data[i4]) + h2i = c2 * (data[i1] - data[i3]) + # Here they are recombined to form the true + # transform of the original real data. + data[i1] = h1r + wr*h2r - wi*h2i + data[i2] = h1i + wr*h2i + wi*h2r + data[i3] = h1r - wr*h2r + wi*h2i + data[i4] = -h1i + wr*h2i + wi*h2r + + wtemp = wr # The reccurrence + wr = wr * wpr - wi * wpi + wr + wi = wi * wpr + wtemp * wpi + wi + } + + if (isign == 1) { + h1r = data[1] + data[1] = h1r + data[2] + data[2] = h1r - data[2] + } else { + h1r = data[1] + data[1] = c1 * (h1r + data[2]) + data[2] = c1 * (h1r - data[2]) + call four1 (data,n,-1) + } + +end + + +# FOUR1 - Replaces DATA by it's discrete transform, if ISIGN is input +# as 1; or replaces DATA by NN times it's inverse discrete Fourier transform +# if ISIGN is input as -1. Data is a complex array of length NN or, equiv- +# alently, a real array of length 2*NN. NN *must* be an integer power of +# two. +# +# Based on Numerical Recipes by Press, Flannery, Teukolsky, and Vetterling. +# Used by permission of the authors. +# Copyright(c) 1986 Numerical Recipes Software. + +procedure four1 (data, nn, isign) + +real data[ARB] # Data array (returned as FFT) +int nn # No. of points in data array +int isign # Direction of transform + +double wr, wi, wpr, wpi # Local variables +double wtemp, theta +real tempr, tempi +int i, j, istep +int n, mmax, m + +begin + n = 2 * nn + j = 1 + for (i=1; i<n; i = i + 2) { + if (j > i) { # Swap 'em + tempr = data[j] + tempi = data[j+1] + data[j] = data[i] + data[j+1] = data[i+1] + data[i] = tempr + data[i+1] = tempi + } + m = n / 2 + while (m >= 2 && j > m) { + j = j - m + m = m / 2 + } + j = j + m + } + mmax = 2 + while (n > mmax) { + istep = 2 * mmax + theta = TWOPI / double (isign*mmax) + wtemp = dsin (0.5*theta) + wpr = -2.d0 * wtemp * wtemp + wpi = dsin (theta) + wr = 1.d0 + wi = 0.d0 + for (m=1; m < mmax; m = m + 2) { + for (i=m; i<=n; i = i + istep) { + j = i + mmax + tempr = real (wr) * data[j] - real (wi) * data[j+1] + tempi = real (wr) * data[j + 1] + real (wi) * data[j] + data[j] = data[i] - tempr + data[j+1] = data[i+1] - tempi + data[i] = data[i] + tempr + data[i+1] = data[i+1] + tempi + } + wtemp = wr + wr = wr * wpr - wi * wpi + wr + wi = wi * wpr + wtemp * wpi + wi + } + mmax = istep + } +end + + +################################################################################ +# LU Decomosition +################################################################################ +define TINY (1E-20) # Number of numerical limit + +# Given an N x N matrix A, with physical dimension N, this routine +# replaces it by the LU decomposition of a rowwise permutation of +# itself. A and N are input. A is output, arranged as in equation +# (2.3.14) above; INDX is an output vector which records the row +# permutation effected by the partial pivioting; D is output as +/-1 +# depending on whether the number of row interchanges was even or odd, +# respectively. This routine is used in combination with LUBKSB to +# solve linear equations or invert a matrix. +# +# Based on Numerical Recipes by Press, Flannery, Teukolsky, and Vetterling. +# Used by permission of the authors. +# Copyright(c) 1986 Numerical Recipes Software. + +procedure ludcmp (a, n, np, indx, d) + +real a[np,np] +int n +int np +int indx[n] +real d + +int i, j, k, imax +real aamax, sum, dum +pointer vv + +begin + # Allocate memory. + call malloc (vv, n, TY_REAL) + + # Loop over rows to get the implict scaling information. + d = 1. + do i = 1, n { + aamax = 0. + do j = 1, n { + if (abs (a[i,j]) > aamax) + aamax = abs (a[i,j]) + } + if (aamax == 0.) { + call mfree (vv, TY_REAL) + call error (1, "Singular matrix") + } + Memr[vv+i-1] = 1. / aamax + } + + # This is the loop over columns of Crout's method. + do j = 1, n { + do i = 1, j-1 { + sum = a[i,j] + do k = 1, i-1 + sum = sum - a[i,k] * a[k,j] + a[i,j] = sum + } + + aamax = 0. + do i = j, n { + sum = a[i,j] + do k = 1, j-1 + sum = sum - a[i,k] * a[k,j] + a[i,j] = sum + dum = Memr[vv+i-1] * abs (sum) + if (dum >= aamax) { + imax = i + aamax = dum + } + } + + if (j != imax) { + do k = 1, n { + dum = a[imax,k] + a[imax,k] = a[j,k] + a[j,k] = dum + } + d = -d + Memr[vv+imax-1] = Memr[vv+j-1] + } + indx[j] = imax + + # Now, finally, divide by the pivot element. + # If the pivot element is zero the matrix is signular (at + # least to the precission of the algorithm. For some + # applications on singular matrices, it is desirable to + # substitute TINY for zero. + + if (a[j,j] == 0.) + a[j,j] = TINY + if (j != n) { + dum = 1. / a[j,j] + do i = j+1, n + a[i,j] = a[i,j] * dum + } + } + + call mfree (vv, TY_REAL) +end + + +# Solves the set of N linear equations AX = B. Here A is input, not +# as the matrix of A but rather as its LU decomposition, determined by +# the routine LUDCMP. INDX is input as the permuation vector returned +# by LUDCMP. B is input as the right-hand side vector B, and returns +# with the solution vector X. A, N, NP and INDX are not modified by +# this routine and can be left in place for successive calls with +# different right-hand sides B. This routine takes into account the +# possiblity that B will begin with many zero elements, so it is +# efficient for use in matrix inversion. +# +# Based on Numerical Recipes by Press, Flannery, Teukolsky, and Vetterling. +# Used by permission of the authors. +# Copyright(c) 1986 Numerical Recipes Software. + +procedure lubksb (a, n, np, indx, b) + +real a[np,np] +int n +int np +int indx[n] +real b[n] + +int i, j, ii, ll +real sum + +begin + ii = 0 + do i = 1, n { + ll = indx[i] + sum = b[ll] + b[ll] = b[i] + if (ii != 0) { + do j = ii, i-1 + sum = sum - a[i,j] * b[j] + } else if (sum != 0.) + ii = i + b[i] = sum + } + + do i = n, 1, -1 { + sum = b[i] + if (i < n) { + do j = i+1, n + sum = sum - a[i,j] * b[j] + } + b[i] = sum / a[i,i] + } +end + + +# Invert a matrix using LU decomposition using A as both input and output. + +procedure luminv (a, n, np) + +real a[np,np] +int n +int np + +int i, j +real d +pointer y, indx + +begin + # Allocate working memory. + call calloc (y, n*n, TY_REAL) + call malloc (indx, n, TY_INT) + + # Setup identify matrix. + do i = 0, n-1 + Memr[y+(n+1)*i] = 1. + + # Do LU decomposition. + call ludcmp (a, n, np, Memi[indx], d) + + # Find inverse by columns. + do j = 0, n-1 + call lubksb (a, n, np, Memi[indx], Memr[y+n*j]) + + # Return inverse in a. + do i = 1, n + do j = 1, n + a[i,j] = Memr[y+n*(j-1)+(i-1)] + + call mfree (y, TY_REAL) +end +################################################################################ |