Repatch (from linux) of OSX IRAF

1 files changed, 169 insertions, 0 deletions

diff --git a/math/surfit/issolve.x b/math/surfit/issolve.x
new file mode 100644
index 00000000..7790ef60
--- /dev/null
+++ b/math/surfit/issolve.x
@@ -0,0 +1,169 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <math/surfit.h>
+include "surfitdef.h"
+
+# ISSOLVE -- Procedure to solve the x coefficient matrix for the surface
+# coefficients. The inner products of the y basis functions are accumulated
+# and stored in the SF_YORDER(sf) by SF_NYCOEFF(sf) array YMATRIX. The
+# main diagonal of YMATRIX is stored in the first row of YMATRIX followed
+# by the remaining non-zero lower diagonals. The Cholesky factorization
+# of YMATRIX is calculated and stored on top of YMATRIX destroying the
+# original data. The inner products of the y basis functions and 
+# and the i-th column of the SF_NXCOEFF(sf) by SF_NLINES(sf) matrix XCOEFF
+# containing the i-th x coefficients for each line are calculated and
+# placed in the i-th row of the SF_NYCOEFF(sf) by SF_NXCOEFF(sf) matrix
+# COEFF. Each of the SF_NXCOEFF(sf) rows of COEFF is solved to determine
+# the SF_NXCOEFF(sf) by SF_NYCOEFF(sf) surface coefficients. After a
+# call to SIFSOLVE the coefficient of the i-th x basis function and the
+# j-th y basis function will be found in the j-th column and i-th row
+# of COEFF.
+
+procedure issolve (sf, lines, nlines, ier)
+
+pointer	sf		# pointer to the curve descriptor structure
+int	lines[ARB]	# line numbers included in the fit
+int	nlines		# number of lines fit
+int	ier		# error code
+
+int	i, ii, j, k, nxcoeff
+pointer	ybzptr, ybptr
+pointer	ylzptr
+pointer	ymzptr, ymindex
+pointer	xczptr, xcptr, xcindex
+pointer	czptr, cptr
+pointer	left, tleft, rows
+pointer	sp
+
+begin
+	# define pointers
+	ybzptr = SF_YBASIS(sf) - 1
+	ymzptr = SF_YMATRIX(sf)
+	xczptr = SF_XCOEFF(sf) - SF_NXCOEFF(sf) - 1
+	czptr = SF_COEFF(sf) - 1
+
+	# zero out coefficient matrix and the y coefficient matrix
+	call aclrr (YMATRIX(SF_YMATRIX(sf)), SF_YORDER(sf) * SF_NYCOEFF(sf))
+	call aclrr (COEFF(SF_COEFF(sf)), SF_NXCOEFF(sf) * SF_NYCOEFF(sf))
+
+	# increment the number of points
+	SF_NYPTS(sf) =  nlines
+
+	switch (SF_TYPE(sf)) {
+	case SF_LEGENDRE, SF_CHEBYSHEV:
+
+	    # accumulate the y value in the y matrix 
+	    nxcoeff = SF_NXCOEFF(sf)
+	    do i = 1, SF_YORDER(sf) {
+		cptr = czptr + i
+		do k = 1, nxcoeff {
+		    xcptr = xczptr + k
+		    do j = 1, nlines {
+		        xcindex = xcptr + lines[j] * SF_NXCOEFF(sf)
+		        COEFF(cptr) = COEFF(cptr) +
+			    YBASIS(ybzptr+lines[j]) * XCOEFF(xcindex)
+		    }
+		    cptr = cptr + SF_NYCOEFF(sf)
+		}
+		ii = 0
+		ybptr = ybzptr
+		do k = i, SF_YORDER(sf) {
+		    ymindex = ymzptr + ii
+		    do j = 1, nlines  
+		        YMATRIX(ymindex) = YMATRIX(ymindex) +
+				     YBASIS(ybzptr+lines[j]) *
+				     YBASIS(ybptr+lines[j])
+		    ii = ii + 1
+		    ybptr = ybptr + SF_NLINES(sf)
+		}
+
+		if (SF_XTERMS(sf) == NO)
+		    nxcoeff = 1
+
+		ybzptr = ybzptr + SF_NLINES(sf)
+		ymzptr = ymzptr + SF_YORDER(sf)
+	    }
+
+	case SF_SPLINE3, SF_SPLINE1:
+
+	    call smark (sp)
+	    call salloc (left, nlines, TY_INT)
+	    call salloc (tleft, nlines, TY_INT)
+	    call salloc (rows, nlines, TY_INT)
+
+	    ylzptr = SF_YLEFT(sf) - 1
+	    do j = 1, nlines
+		Memi[left+j-1] = YLEFT(ylzptr+lines[j])
+	    call amulki (Memi[left], SF_YORDER(sf), Memi[rows], nlines)
+	    call aaddki (Memi[rows], SF_YMATRIX(sf), Memi[rows], nlines)
+	    call aaddki (Memi[left], czptr, Memi[left], nlines)
+
+	    # accumulate the y value in the y matrix 
+	    nxcoeff = SF_NXCOEFF(sf)
+	    do i = 1, SF_YORDER(sf) {
+		call aaddki (Memi[left], i, Memi[tleft], nlines)
+		do k = 1, nxcoeff {
+		    xcptr = xczptr + k
+		    do j = 1, nlines {
+			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],
+		        nlines)
+		}
+		ii = 0
+		ybptr = ybzptr
+		do k = i, SF_YORDER(sf) {
+		    do j = 1, nlines {  
+		        ymindex = Memi[rows+j-1] + ii
+		        YMATRIX(ymindex) = YMATRIX(ymindex) +
+				     YBASIS(ybzptr+lines[j]) *
+				     YBASIS(ybptr+lines[j])
+		    }
+		    ii = ii + 1
+		    ybptr = ybptr + SF_NLINES(sf)
+		}
+
+		ybzptr = ybzptr + SF_NLINES(sf)
+		call aaddki (Memi[rows], SF_YORDER(sf), Memi[rows], nlines)
+	    }
+
+	    call sfree (sp)
+
+	}
+
+	# return if not enough points
+	ier = OK
+	if ((SF_NYPTS(sf) - SF_NYCOEFF(sf)) < 0) {
+	    ier = NO_DEG_FREEDOM
+	    return
+	}
+
+	# calculate the Cholesky factorization of the y matrix
+	call sfchofac (YMATRIX(SF_YMATRIX(sf)), SF_YORDER(sf), SF_NYCOEFF(sf),
+		       YMATRIX(SF_YMATRIX(sf)), ier)
+
+	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 {
+
+	    cptr = SF_COEFF(sf)
+	    call sfchoslv (YMATRIX(SF_YMATRIX(sf)), SF_YORDER(sf),
+			      SF_NYCOEFF(sf), COEFF(cptr), COEFF(cptr))
+
+	    do i = 2, SF_NXCOEFF(sf) {
+		cptr = cptr + SF_NYCOEFF(sf)
+		COEFF(cptr) = COEFF(cptr) / SF_NYPTS(sf)
+	    }
+	}
+end