Initial commit

1 files changed, 93 insertions, 0 deletions

diff --git a/math/deboor/spline.x b/math/deboor/spline.x
new file mode 100644
index 00000000..d77c8122
--- /dev/null
+++ b/math/deboor/spline.x
@@ -0,0 +1,93 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+
+include	"bspln.h"
+
+.help spline 2 "math library"
+.ih
+NAME
+spline -- generalized spline interpolation with b-splines.
+.ih
+USAGE
+spline (x, y, n, bspln, q, k, ier)
+.ih
+PARAMETERS
+.ls x,y
+(real[n]).  Abcissae and ordinates of the N data points.
+.le
+.ls n
+The number of input data points, and the number of spline coefficients
+to be generated.
+.le
+.ls bspln
+(real[2*n+30]).  Output array of N b-spline coefficients, N+K knots,
+and the spline descriptor array.
+.le
+.ls q
+(real[(2*K-1)*N]).  On output, will contain the
+triangular factorization of the coefficient matrix of the linear
+system for the b-coefficients of the spline interpolant.  If a new
+data set differs only in the Y values, Q may be input to FSPLIN to
+efficiently solve for the b-spline coefficients of the new data set.
+.le
+.ls k
+The order of the spline (2-20).  Must be EVEN at present
+(k=4 is the cubic spline).
+.le
+.ls ier
+Zero if ok, nonzero if error (invalid input parameters).
+.le
+.ih
+DESCRIPTION
+General spline interpolation.  SPLINE is a thinly disguised version of SPLINT
+(DeBoor, pg. 204).  SPLINE differs from SPLINT in that the knot spacing T is
+calculated automatically, using the not-a-knot boundary conditions.  Only even
+K are permitted at present (K = 2, 4, 6,...).  The cubic spline interpolant
+corresponds to K = 4.  SPLINE saves a complete description of the spline in the
+BSPLN array, for use later by SEVAL (used to evaluate the fitted spline) and
+FSPLIN (used after a call to SPLINE to more efficiently fit subsequent data
+sets).
+.ih
+SOURCE
+See the listing of SPLINT, or Chapter XIII of DeBoors book.
+.ih
+SEE ALSO
+seval(2), fsplin(2).
+.endhelp ______________________________________________________________
+
+
+procedure spline (x, y, n, bspln, q, k, ier)
+
+real	x[n], y[n]
+int	n, k, ier
+real	q[ARB], bspln[ARB]		#q[(2*k-1)*n], bspln[2*n+30]
+int	m, i, knot
+
+begin
+	if (k < 2 || mod (k, 2) != 0) {
+	    ier = 1
+	    return
+	}
+
+	NCOEF = n			#set up spline descriptor
+	ORDER = k
+	XMIN = x[1]
+	XMAX = x[n]
+
+	knot = int (KNOT1) - 1		#offset to knot array in bspln
+	KINDEX = knot + k		#initial posn for SEVAL
+
+	m = knot + n
+	do i = 1, k {			#not-a-knot knot boundary conditions
+	    bspln[knot+i] = XMIN
+	    bspln[m+i] = XMAX
+	}
+
+	m = k / 2			#use data x-values inside
+	do i = k+1, n
+	    bspln[knot+i] = x[i+m-k]
+
+	call splint (x, y, bspln[knot+1], n, k, q, bspln[COEF1], ier)
+	if (ier == 1)			#1=success, 2==failure
+	    ier = 0
+end