Repatch (from linux) of OSX IRAF

1 files changed, 150 insertions, 0 deletions

diff --git a/math/surfit/islfit.x b/math/surfit/islfit.x
new file mode 100644
index 00000000..7b60c361
--- /dev/null
+++ b/math/surfit/islfit.x
@@ -0,0 +1,150 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <math/surfit.h>
+include "surfitdef.h"
+
+# ISLFIT -- Procedure to fit a single line of a surface. The inner products
+# of the x basis functions are calculated and accumulated into the
+# SF_XORDER(sf) by SF_NXCOEFF(sf) matrix XMATRIX. The main diagonal is stored
+# in the first row of XMATRIX followed by the non-zero lower diagonals. This
+# method of storage is very efficient for the large symmetric banded matrices
+# generated by spline fits. The Cholesky factorization of XMATRIX is calculated
+# and stored in XMATRIX destroying the original data. The inner products
+# of the x basis functions and the data ordinates are stored in the lineno-th
+# row of the SF_NXCOEFF(sf) by SF_NLINES(sf) matrix XCOEFF. The x coefficients
+# for each line are calculated by forward and back substitution and
+# stored in the lineno-th line of XCOEFF destroying the original data.
+
+procedure islfit (sf, cols, lineno, z, w, ncols, wtflag, ier)
+
+pointer	sf		# pointer to the surface descriptor
+int	cols[ncols]	# array of columns
+int	lineno    	# lineno
+real	z[ncols]	# the surface values
+real	w[ncols]	# array of weights
+int	ncols		# the number of columns
+int	wtflag		# type of weighting
+int	ier		# error codes
+
+int	i, ii, j, k
+pointer	xbzptr, xbptr
+pointer	xmzptr, xmindex
+pointer	xczptr, xcindex
+pointer	xlzptr
+pointer	left, rows, bw
+pointer	sp
+real	adotr()
+
+begin
+	# index the pointers
+	xbzptr = SF_XBASIS(sf) - 1
+	xmzptr = SF_XMATRIX(sf)
+	xczptr = SF_XCOEFF(sf) + (lineno - 1) * SF_NXCOEFF(sf) - 1
+
+	# zero the accumulators
+	call aclrr (XMATRIX(SF_XMATRIX(sf)), SF_NXCOEFF(sf) * SF_XORDER(sf))
+	call aclrr (XCOEFF(xczptr + 1), SF_NXCOEFF(sf))
+
+	# free space used by previous islrefit calls
+	if (SF_WZ(sf) != NULL)
+	    call mfree (SF_WZ(sf), MEM_TYPE)
+	if (SF_TLEFT(sf) != NULL)
+	    call mfree (SF_TLEFT(sf), TY_INT)
+
+	# reset number of points
+	SF_NXPTS(sf) = ncols
+
+	# calculate the weights, default is uniform weighting
+	switch (wtflag) {
+	case SF_UNIFORM:
+	    call amovkr (1.0, w, ncols)
+	case SF_USER:
+	    # do not alter weights
+	default:
+	    call amovkr (1.0, w, ncols)
+	}
+
+	# allocate temporary space
+	call smark (sp)
+	call salloc (bw, ncols, TY_REAL)
+
+	switch (SF_TYPE(sf)) {
+	case SF_LEGENDRE, SF_CHEBYSHEV:
+
+	    do i = 1, SF_XORDER(sf) {
+		do j = 1, ncols
+		    Memr[bw+j-1] = w[j] * XBASIS(xbzptr+cols[j])
+
+		xcindex = xczptr + i
+		XCOEFF(xcindex) = XCOEFF(xcindex) + adotr (Memr[bw], z, ncols)
+
+		xbptr = xbzptr
+		ii = 0
+		do k = i, SF_XORDER(sf) {
+		    xmindex = xmzptr + ii
+		    do j = 1, ncols
+		        XMATRIX(xmindex) = XMATRIX(xmindex) + Memr[bw+j-1] *
+			    XBASIS(xbptr+cols[j])
+		    ii = ii + 1
+		    xbptr = xbptr + SF_NCOLS(sf)
+		}
+
+		xbzptr = xbzptr + SF_NCOLS(sf)
+		xmzptr = xmzptr + SF_XORDER(sf)
+	    }
+
+	case SF_SPLINE3, SF_SPLINE1:
+	    xlzptr = SF_XLEFT(sf) - 1
+
+	    call salloc (left, ncols, TY_INT)
+	    call salloc (rows, ncols, TY_INT)
+
+	    do j = 1, ncols
+		Memi[left+j-1] = XLEFT(xlzptr+cols[j])
+	    call amulki (Memi[left], SF_XORDER(sf), Memi[rows], ncols)
+	    call aaddki (Memi[rows], SF_XMATRIX(sf), Memi[rows], ncols)
+	    call aaddki (Memi[left], xczptr, Memi[left], ncols)
+
+	    do i = 1, SF_XORDER(sf) {
+		do j = 1, ncols {
+		    Memr[bw+j-1] = w[j] * XBASIS(xbzptr+cols[j])
+		    xcindex = Memi[left+j-1] + i
+		    XCOEFF(xcindex) = XCOEFF(xcindex) + Memr[bw+j-1] * z[j]
+		}
+
+		xbptr = xbzptr
+		ii = 0
+		do k = i, SF_XORDER(sf) {
+		    do j = 1, ncols {
+			xmindex = Memi[rows+j-1] + ii
+			XMATRIX(xmindex) = XMATRIX(xmindex) + Memr[bw+j-1] *
+			    XBASIS(xbptr+cols[j])
+		    }
+		    ii = ii + 1
+		    xbptr = xbptr + SF_NCOLS(sf)
+		}
+
+		xbzptr = xbzptr + SF_NCOLS(sf)
+		call aaddki (Memi[rows], SF_XORDER(sf), Memi[rows], ncols)
+	    }
+	}
+
+	# release space
+	call sfree (sp)
+
+	# return if not enough data points
+	ier = OK
+	if ((SF_NXPTS(sf) - SF_NXCOEFF(sf)) < 0) {
+	    ier = NO_DEG_FREEDOM
+	    return
+	}
+
+	# calculate the Cholesky factorization of XMATRIX
+	call sfchofac (XMATRIX(SF_XMATRIX(sf)), SF_XORDER(sf), SF_NXCOEFF(sf),
+	    XMATRIX(SF_XMATRIX(sf)), ier)
+
+	# calculate the x coefficients for lineno-th image line and place in the
+	# lineno-th row of XCOEFF
+	call sfchoslv (XMATRIX(SF_XMATRIX(sf)), SF_XORDER(sf), SF_NXCOEFF(sf),
+	    XCOEFF(xczptr + 1), XCOEFF(xczptr + 1))
+end