1 files changed, 127 insertions, 0 deletions
diff --git a/math/surfit/isresolve.x b/math/surfit/isresolve.x
new file mode 100644
index 00000000..7af6a8da
--- /dev/null
+++ b/math/surfit/isresolve.x
@@ -0,0 +1,127 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <math/surfit.h>
+include "surfitdef.h"
+
+# ISRESOLVE -- Procedure to solve the x coefficient matrix for the surface
+# coefficients assuming that the lines array is unchanged since the last
+# call to SIFSOLVE. The Cholesky factorization of YMATRIX is assumed to be
+# stored in YMATRIX. The inner product of the y basis functions and i-th
+# column of XCOEFF containing the i-th x coefficients for each line are
+# calculated and stored in the i-th row of the SF_NYCOEFF(sf) by
+# SF_NXCOEFF(sf) array COEFF. Each of the SF_NXCOEFF(sf) rows of COEFF
+# is solved to determine the SF_NYCOEFF(sf) by SF_NXCOEFF(sf) surface
+# coefficients. After a call to SIFSOLVE the coefficient for the i-th
+# y basis function and the j-th x coefficient is found in the j-th row
+# and i-th column of COEFF.
+
+procedure isresolve (sf, lines, ier)
+
+pointer	sf		# pointer to the surface descriptor structure
+int	lines[ARB]	# line numbers included in the fit
+int	ier		# error code
+
+int	i, j, k, nxcoeff
+pointer	ybzptr
+pointer	ylzptr
+pointer	xczptr, xcptr, xcindex
+pointer	czptr, cptr
+pointer	left, tleft
+pointer	sp
+
+begin
+	# define pointers
+	ybzptr = SF_YBASIS(sf) - 1
+	xczptr = SF_XCOEFF(sf) - SF_NXCOEFF(sf) - 1
+	czptr = SF_COEFF(sf) - 1
+
+	# zero out coefficient matrix
+	call aclrr (COEFF(SF_COEFF(sf)), SF_NXCOEFF(sf) * SF_NYCOEFF(sf))
+
+	switch (SF_TYPE(sf)) {
+	case SF_LEGENDRE, SF_CHEBYSHEV:
+
+	    # loop over the y basis functions
+	    nxcoeff = SF_NXCOEFF(sf)
+	    do i = 1, SF_YORDER(sf) {
+		cptr = czptr + i
+		do k = 1, nxcoeff {
+		    xcptr = xczptr + k
+		    do j = 1, SF_NYPTS(sf) {
+			xcindex = xcptr + lines[j] * SF_NXCOEFF(sf)
+		        COEFF(cptr) = COEFF(cptr) +
+		    		  YBASIS(ybzptr+lines[j]) * XCOEFF(xcindex)
+		    }
+		    cptr = cptr + SF_NYCOEFF(sf)
+		}
+
+		ybzptr = ybzptr + SF_NYPTS(sf)
+
+		# if SF_XTERMS(sf) = NO do not accumulate elements
+		# of COEFF where i != 1 and k != 1
+		if (SF_XTERMS(sf) == NO)
+		    nxcoeff = 1
+	    }
+
+	case SF_SPLINE3, SF_SPLINE1:
+	    call smark (sp)
+	    call salloc (left, SF_NYPTS(sf), TY_INT)
+	    call salloc (tleft, SF_NYPTS(sf), TY_INT)
+
+	    ylzptr = SF_YLEFT(sf) - 1
+	    do j = 1, SF_NYPTS(sf)
+		Memi[left+j-1] = YLEFT(ylzptr+lines[j])
+	    call aaddki (Memi[left], czptr, Memi[left], SF_NYPTS(sf))
+
+	    nxcoeff = SF_NXCOEFF(sf)
+	    do i = 1, SF_YORDER(sf) {
+		call aaddki (Memi[left], i, Memi[tleft], SF_NYPTS(sf))
+		do k = 1, nxcoeff {
+		    xcptr = xczptr + k
+		    do j = 1, SF_NYPTS(sf) {
+			cptr = Memi[tleft+j-1]
+			xcindex = xcptr + lines[j] * SF_NXCOEFF(sf)
+			COEFF(cptr) = COEFF(cptr) + YBASIS(ybzptr+lines[j]) *
+			    XCOEFF(xcindex)
+		    }
+		    call aaddki (Memi[tleft], SF_NYCOEFF(sf), Memi[tleft],
+			SF_NYPTS(sf))
+		}
+
+		ybzptr = ybzptr + SF_NYPTS(sf)
+	    }
+
+	    call sfree (sp)
+	}
+
+	# return if not enough points
+	ier = OK
+	if ((SF_NYPTS(sf) - SF_NYCOEFF(sf)) < 0) {
+	    ier = NO_DEG_FREEDOM
+	    return
+	}
+
+	if (SF_XTERMS(sf) == YES) {
+
+	    # solve for the nxcoeff right sides
+	    cptr = SF_COEFF(sf)
+	    do i = 1, SF_NXCOEFF(sf) {
+	        call sfchoslv (YMATRIX(SF_YMATRIX(sf)), SF_YORDER(sf),
+	    		       SF_NYCOEFF(sf), COEFF(cptr), COEFF(cptr))
+	        cptr = cptr + SF_NYCOEFF(sf)
+	    }
+
+	} else {
+
+	    # solve for the y coefficients
+	    cptr = SF_NYCOEFF(sf)
+	    call sfchoslv (YMATRIX(SF_YMATRIX(sf)), SF_YORDER(sf),
+			      SF_NYCOEFF(sf), COEFF(cptr), COEFF(cptr))
+
+	    # solve for the x coefficients
+	    do i = 2, SF_NXCOEFF(sf) {
+		cptr = czptr + SF_NYCOEFF(sf)
+		COEFF(cptr) = COEFF(cptr) / SF_NYPTS(sf)
+	    }
+	}
+end