Merge pull request #2 from olebole/replace-numerical-recipes

Replace numerical recipes code

sptable.par

@@ -1,3 +1,3 @@
 # Package parameters.
 
-version,s,h,"V1.0: 20140915"
+version,s,h,"V1.0: 20180612"

src/mwcs/mwlu.x

@@ -1,14 +1,11 @@
 # Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
 
-include	<mach.h>
-
 # 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

src/xonedspec/splot/spdeblend.x

@@ -1,3 +1,4 @@
+include	<math.h>
 include	<error.h>
 include	<gset.h>
 
@@ -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