diff --git a/sptable.par b/sptable.par index e5d9575..db0c637 100644 --- a/sptable.par +++ b/sptable.par @@ -1,3 +1,3 @@ # Package parameters. -version,s,h,"V1.0: 20140915" +version,s,h,"V1.0: 20180612" diff --git a/src/mwcs/mwlu.x b/src/mwcs/mwlu.x index f6a606f..38b8e4f 100644 --- a/src/mwcs/mwlu.x +++ b/src/mwcs/mwlu.x @@ -1,14 +1,11 @@ # Copyright(c) 1986 Association of Universities for Research in Astronomy Inc. -include - # MULU -- Matrix utilities for MWCS. # # mw_ludecompose performs LU decomposition of a square matrix # mw_lubacksub performs backsubstitution to solve a system # -# These routines are derived from routines in the book Numerical Recipes, -# Press et. al. 1986. +# These routines simply call the LU composition routines provided by LAPACK. # MW_LUDECOMPOSE -- Replace an NxN matrix A by the LU decomposition of a @@ -21,78 +18,10 @@ double a[ndim,ndim] #U matrix to be inverted; inverted matrix int ix[ndim] #O vector describing row permutation int ndim #I dimension of square matrix -pointer sp, vv -int d, i, j, k, imax -double aamax, sum, dum +int status begin - call smark (sp) - call salloc (vv, ndim, TY_DOUBLE) - - # Keep track of the number of row interchanges, odd or even (not used). - d = 1 - - # Loop over rows to get implicit scaling information. - do i = 1, ndim { - aamax = 0.0 - do j = 1, ndim - if (abs(a[i,j]) > aamax) - aamax = abs(a[i,j]) - if (aamax == 0.0) - call error (1, "singular matrix") - Memd[vv+i-1] = 1.0 / aamax - } - - # Loop over columns using Crout's method. - do j = 1, ndim { - 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 - } - - # Search for the largest pivot element. - aamax = 0.0 - do i = j, ndim { - sum = a[i,j] - do k = 1, j-1 - sum = sum - a[i,k] * a[k,j] - a[i,j] = sum - - # Figure of merit for the pivot. - dum = Memd[vv+i-1] * abs(sum) - if (dum >= aamax) { - imax = i - aamax = dum - } - } - - # Do we need to interchange rows? - if (j != imax) { - # Yes, do so... - do k = 1, ndim { - dum = a[imax,k] - a[imax,k] = a[j,k] - a[j,k] = dum - } - d = -d - Memd[vv+imax-1] = Memd[vv+j-1] - } - - ix[j] = imax - if (a[j,j] == 0.0) - a[j,j] = EPSILOND - - # Divide by the pivot element. - if (j != ndim) { - dum = 1.0 / a[j,j] - do i = j+1, ndim - a[i,j] = a[i,j] * dum - } - } - - call sfree (sp) + call dgetrf(ndim, ndim, a, ndim, ix, status) end @@ -109,35 +38,8 @@ int ix[ndim] #I permutation vector for A double b[ndim] #U rhs vector; solution vector int ndim #I dimension of system -int ii, ll, i, j -double sum +int status begin - # Do the forward substitution, unscrambling the permutation as we - # go. When II is set to a positive value, it will become the index - # of the first nonvanishing element of B. - - ii = 0 - do i = 1, ndim { - ll = ix[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 - } - - # Now do the backsubstitution. - do i = ndim, 1, -1 { - sum = b[i] - if (i < ndim) - do j = i+1, ndim - sum = sum - a[i,j] * b[j] - b[i] = sum / a[i,i] - } + call dgetrs('N', ndim, 1, a, ndim, ix, b, ndim, status) end diff --git a/src/xonedspec/splot/spdeblend.x b/src/xonedspec/splot/spdeblend.x index a07cd52..13a485b 100644 --- a/src/xonedspec/splot/spdeblend.x +++ b/src/xonedspec/splot/spdeblend.x @@ -1,3 +1,4 @@ +include include include @@ -789,31 +790,34 @@ 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. +# Copyright(c) 2017 Anastasia Galkin +# Reference: G. E. P. Box and Mervin E. Muller, A Note on the Generation of +# Random Normal Deviates, The Annals of Mathematical Statistics +# (1958), Vol. 29, No. 2 pp. 610–611 real procedure gasdev (seed) -long seed # Seed for random numbers +long seed + +int count +data count/0/ -real v1, v2, r, fac, urand() -int iset -data iset/0/ +real u1, u2, x +real urand() 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) + if (count == 0) { + repeat { + u1 = 2 * urand (seed) - 1. + } until (u1 > 0) + repeat { + u2 = 2 * urand (seed) - 1. + } until (u1 > 0) + x = sqrt(-2 * log(u1)) * cos(2*PI*u2); + count = 1 } else { - iset = 0 - return (v2 * fac) + x = sqrt(-2 * log(u1)) * sin(2*PI*u2); + count = 0 } + return (x) end diff --git a/src/xrv/numrep.x b/src/xrv/numrep.x index 75c6f9e..4be7106 100644 --- a/src/xrv/numrep.x +++ b/src/xrv/numrep.x @@ -1,75 +1,7 @@ include include -# NUMREP.X - A collection of routines recoded from Numerical Recipes by -# Press, Flannery, Teukolsky, and Vetterling. Used by permission of the -# authors. Copyright(c) 1986 Numerical Recipes Software. - - -# 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. - -procedure four1 (data, nn, isign) - -real data[ARB] #U Data array (returned as FFT) -int nn #I No. of points in data array -int isign #I 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 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 - +# NUMREP.X - A collection of routines API compatible to Numerical Recipes # 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 @@ -83,79 +15,45 @@ end procedure realft (data, N, isign) -real data[ARB] #U Input data array & output FFT +real data[ARB] #U Input data array & output FFT int N #I No. of points -int isign #I 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 +int isign #I Direction of transfer +pointer sp, wsave +real last +int j 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) - } - + call smark(sp) + call salloc(wsave, 4*N+15, TY_REAL) + call rffti(2*N, Memr[wsave]) + + if (isign == 1) { + call rfftf(2*N, data, Memr[wsave]) + last = data[2*N] + do j=2*N-1,3,-2 { + data[j+1] = -data[j] + data[j] = data[j-1] + } + data[2] = last + } else { + data[1] = data[1]/2.0 + last = data[2]/2.0 + do j=2,2*N-2,2 { + data[j] = data[j+1]/2.0 + data[j+1] = -data[j+2]/2.0 + } + data[2*N] = last + call rfftb(2*N, data, Memr[wsave]) + } + + 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. +# N, this routine 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. procedure twofft (data1, data2, fft1, fft2, N) @@ -163,37 +61,23 @@ real data1[ARB], data2[ARB] #I Input data arrays real fft1[ARB], fft2[ARB] #O Output FFT arrays int N #I No. of points -int nn3, nn2, jj, j -real rep, rem, aip, aim - +int j +pointer sp, wsave 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 - } - + call smark(sp) + call salloc(wsave, 4*N+15, TY_REAL) + call cffti(N, Memr[wsave]) + do j=1, N { + fft1[2*j-1] = data1[j] + fft1[2*j] = 0.0 + fft2[2*j-1] = data2[j] + fft2[2*j] = 0.0 + } + call cfftf(N, fft1, Memr[wsave]) + call cfftf(N, fft2, Memr[wsave]) + do j=1, N { + fft1[2*j] = -fft1[2*j] + fft2[2*j] = -fft2[2*j] + } + call sfree(sp) end