1 files changed, 74 insertions, 0 deletions
diff --git a/math/surfit/islrefit.x b/math/surfit/islrefit.x
new file mode 100644
index 00000000..b0430e5b
--- /dev/null
+++ b/math/surfit/islrefit.x
@@ -0,0 +1,74 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <math/surfit.h>
+include "surfitdef.h"
+
+# ISLREFIT -- Procedure to refit the data assuming that the cols and w
+# arrays do not change. SIFLREFIT assumes that the Cholesky factorization
+# of XMATRIX is stored in XMATRIX. The inner products of the x basis
+# functions and the data ordinates are accumulated into the lineno-th
+# row of the SF_NXCOEFF(sf) by SF_NLINES(sf) matrix XCOEFF. The coefficients
+# for line number lineno are calculated by forward and back
+# substitution and placed in the lineno-th row of XCOEFF replacing the
+# original data.
+
+procedure islrefit (sf, cols, lineno, z, w)
+
+pointer	sf		# pointer to surface descriptor
+int	cols[ARB]	# columns to be fit 
+int	lineno		# line number
+real	z[ARB]		# surface values
+real	w[ARB]		# weight values
+
+int	i, j
+pointer	xbzptr, xczptr, xcindex, xlzptr
+
+begin
+	# set pointers
+	xbzptr = SF_XBASIS(sf) - 1
+	xczptr = SF_XCOEFF(sf) + (lineno - 1) * SF_NXCOEFF(sf) - 1
+	xlzptr = SF_XLEFT(sf) - 1
+
+	# reset lineno-th row of the x coefficient matrix
+	call aclrr (XCOEFF(xczptr+1), SF_NXCOEFF(sf))
+
+	if (SF_WZ(sf) == NULL)
+	    call malloc (SF_WZ(sf), SF_NXPTS(sf), MEM_TYPE)
+
+	# calculate new right sides
+	call amulr (w, z, Memr[SF_WZ(sf)], SF_NXPTS(sf))
+
+	switch (SF_TYPE(sf)) {
+	case SF_LEGENDRE, SF_CHEBYSHEV:
+
+	    do i = 1, SF_XORDER(sf) {
+		xcindex = xczptr + i
+		do j = 1, SF_NXPTS(sf)
+		    XCOEFF(xcindex) = XCOEFF(xcindex) + Memr[SF_WZ(sf)+j-1] *
+			XBASIS(xbzptr+cols[j])
+		xbzptr = xbzptr + SF_NCOLS(sf)
+	    }
+
+	case SF_SPLINE3, SF_SPLINE1:
+
+	    if (SF_TLEFT(sf) == NULL)
+	        call malloc (SF_TLEFT(sf), SF_NXPTS(sf), TY_INT)
+
+	    do i = 1, SF_NXPTS(sf)
+		Memi[SF_TLEFT(sf)+i-1] = XLEFT(xlzptr+cols[i]) + xczptr
+
+	    do i = 1, SF_XORDER(sf) {
+		do j = 1, SF_NXPTS(sf) {
+		    xcindex = Memi[SF_TLEFT(sf)+j-1] + i
+		    XCOEFF(xcindex) = XCOEFF(xcindex) + Memr[SF_WZ(sf)+j-1] *
+			XBASIS(xbzptr+cols[j])
+		}
+		xbzptr = xbzptr + SF_NCOLS(sf)
+	    }
+
+	}
+
+	# solve for the new x coefficients for line number lineno
+	call sfchoslv (XMATRIX(SF_XMATRIX(sf)), SF_XORDER(sf), SF_NXCOEFF(sf),
+	    XCOEFF(xczptr+1), XCOEFF(xczptr+1))
+end