Initial commit

1 files changed, 160 insertions, 0 deletions

diff --git a/math/deboor/spllsq.x b/math/deboor/spllsq.x
new file mode 100644
index 00000000..8760c8a0
--- /dev/null
+++ b/math/deboor/spllsq.x
@@ -0,0 +1,160 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include	"bspln.h"
+
+.help spllsq 2 "math library"
+.ih ___________________________________________________________________________
+NAME
+spllsq -- general smoothing b-spline with uniform knots.
+.ih
+USAGE
+spllsq (x, y, w, npts, bspln, q, work, k, n, wflg, ier)
+.ih
+PARAMETERS
+See the listing of SPLSQV if a more detailed discussion of the parameters
+than that given here is desired.
+.ls x,y,w
+(real[npts]).  Abscissae, ordinates, and weights of the data points.
+X[1] and X[NPTS] determine the range of the spline (the interval over which
+it can be evaluated).  Dummy data points with zero weight can be supplied
+to fix the range if desired.
+.le
+.ls npts
+The number of data points.
+.le
+.ls bspln
+(real[2*n+30]).  Output array of N b-spline coefficients, N+K knots,
+and the spline descriptor structure.
+.le
+.ls q
+(real[n,k]).  Work array.
+.le
+.ls work
+(real[n]).  Work array.
+.le
+.ls k
+The order of the desired spline (1-20) (Cubic == 4).
+.le
+.ls n
+N is the number of b-spline coefficients to be output.  The relationship
+between N and the number of polynomial pieces in the spline (NPP)
+is given by NPP = N-K-1.  Note that NPP is the important parameter.
+N is used in the parameter list rather than NPP because it is used in
+dimensioning the arrays.
+.le
+.ls wflg
+Weight generation flag.  Options:
+
+.nf
+    0  pass W on to SPLSQV as is
+    1  calculate default weights
+    2  same as 1, but set W[i] only
+         if W[i] already NONZERO
+.fi
+
+If WFLG==2, data points may be rejected from the fit by setting their
+corresponding weight to zero, before calling SPLLSQ (good points must be
+assigned nonzero ws before calling SPLLSQ).  Use WFLG==1 if no points
+are to be rejected, to avoid having to initialize the W array on every
+call.
+.le
+.ls ier
+Error return.  Zero if no error.  Nonzero indicates invalid combination
+of N, K, and NPTS.
+.le
+.ih
+DESCRIPTION
+A general algorithm for smoothing with uniform B-splines of arbitrary order.
+Calls SPLSQV, which is adapted from L2APPR, as described in "A Practical
+Guide To Splines", by C. DeBoor.
+ 
+SPLLSQ is identical to SPLSQV, except that (1): the array T[n+k], giving
+the x-positions of the knots of the b-spline, is generated internally,
+rather than by the calling program, and (2): the weights associated with
+the data points may be calculated automatically.  The knots are generated
+with a uniform spacing.  Note that the x-values of the data points do not
+have to be uniformly spaced, but they must be monotonically increasing, and
+(1) both must span the same range, and (2) if N spline coefficients are
+desired, there must be at least N+K data points.
+.ih
+BUGS
+A routine FSPLSQ needs to be written to make use of the Cholesky factorization
+returned by SPLSQV in W1, for more efficient fitting of subsequent data sets
+that differ only in the y-values of the data points.
+.ih
+SEE ALSO
+seval(2)
+.endhelp _______________________________________________________________________
+
+
+procedure spllsq (x, y, w, npts, bspln, q, work, k, n, wflg, ier)
+
+real	x[npts], y[npts], w[npts]
+real	q[n,k], work[n]
+real	bspln[ARB]
+int	wflg, npts, n, k, ier
+int	i, npp, km1, knot
+real	dx
+
+begin
+	ier = 0				#successful return value
+	npp = n - k + 1			#number of polynomial pieces
+	if (npp <= 0) {		
+	    ier = 1
+	    return
+	}				
+
+	# Set up spline descriptor in the array BCOEF, for later evaluation
+	# of the spline by BSPLN (see "bcoef.h" for definitions of the fields).
+
+	NCOEF = n				#number of b-spline coeff
+	ORDER = k				#order of the spline
+	XMIN = x[1]				#x at left endpoint
+	XMAX = x[npts]				#x at right endpoint
+	KINDEX = KNOT1 - 1 + k			#for SEVAL
+
+
+	# Set values of knots.  First K knots take on the value of the first
+	# breakpoint, next npp-1 knots have spacing DX, and the last K knots
+	# take on the value of the last breakpoint.
+
+	dx = (XMAX - XMIN ) / npp		#span(x) = span(t)
+
+	km1 = k - 1
+	knot = KNOT1 - 1
+
+	do i = 1, km1
+	    bspln[knot+i] = XMIN
+
+	knot = knot + km1
+	do i = 1, npp
+	    bspln[knot+i] = XMIN + (i-1) * dx
+
+	knot = knot + npp
+	do i = 1, k
+	    bspln[knot+i] = XMAX
+
+
+	# Calculate default weights.  The default weight of a datapoint is
+	# proportional to the spacing.
+
+	if (wflg > 0) {
+	    if (w[1] > 0  ||  wflg == 1)
+		w[1] = x[2] - x[1]
+
+	    do i = 2, npts - 1 {
+		if (w[i] > 0  ||  wflg == 1)
+		    w[i] = (x[i+1] - x[i-1]) / 2.
+	    }
+	    if (w[npts] > 0  ||  wflg == 1)
+		w[npts] = x[npts] - x[npts-1]
+	}
+
+
+	# Call SPLSQV to do the actual spline fit.  The N b-spline coefficients
+	# of order K, N+K knots, and the spline descriptor are returned in
+	# BSPLN, for evaluation of the spline via SEVAL.
+
+	call splsqv (x, y, w, npts, bspln[nint(KNOT1)], n, k, q,
+	    work, bspln[COEF1], ier)
+end