Repatch (from linux) of OSX IRAF

1 files changed, 92 insertions, 0 deletions

diff --git a/math/surfit/iseval.x b/math/surfit/iseval.x
new file mode 100644
index 00000000..08ef5f89
--- /dev/null
+++ b/math/surfit/iseval.x
@@ -0,0 +1,92 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <math/surfit.h>
+include "surfitdef.h"
+
+# ISEVAL -- Procedure to evaluate the fitted surface at a single point.
+# The SF_NXCOEFF(sf) by SF_NYCOEFF(sf) coefficients are stored in the
+# SF_NYCOEFF(sf) by SF_NXCOEFF(sf) matrix COEFF. The j-th element of the ith
+# row of COEFF contains the coefficient of the i-th basis function in x and
+# the j-th basis function in y.
+
+real procedure iseval (sf, x, y)
+
+pointer	sf		# pointer to surface descriptor structure
+real	x		# x value
+real	y		# y value
+
+real	sum, accum
+int	i, k, leftx, lefty, yorder
+pointer	sp, xb, xzb, yb, yzb, czptr
+
+begin
+	# allocate space for the basis functions
+	call smark (sp)
+	call salloc (xb, SF_XORDER(sf), MEM_TYPE)
+	xzb = xb - 1
+	call salloc (yb, SF_YORDER(sf), MEM_TYPE)
+	yzb = yb - 1
+
+	# calculate the basis functions
+	switch (SF_TYPE(sf)) {
+	case SF_CHEBYSHEV:
+	    leftx = 0
+	    lefty = 0
+	    czptr = SF_COEFF(sf) - 1
+	    call sf_b1cheb (x, SF_XORDER(sf), SF_XMAXMIN(sf), SF_XRANGE(sf),
+		XBS(xb))
+	    call sf_b1cheb (y, SF_YORDER(sf), SF_YMAXMIN(sf), SF_YRANGE(sf),
+		YBS(yb))
+
+	case SF_LEGENDRE:
+	    leftx = 0
+	    lefty = 0
+	    czptr = SF_COEFF(sf) - 1
+	    call sf_b1leg (x, SF_XORDER(sf), SF_XMAXMIN(sf), SF_XRANGE(sf),
+		XBS(xb))
+	    call sf_b1leg (y, SF_YORDER(sf), SF_YMAXMIN(sf), SF_YRANGE(sf),
+		YBS(yb))
+
+	case SF_SPLINE3:
+	    call sf_b1spline3 (x, SF_NXPIECES(sf), -SF_XMIN(sf),
+	        SF_XSPACING(sf), XBS(xb), leftx)
+	    call sf_b1spline3 (y, SF_NYPIECES(sf), -SF_YMIN(sf),
+	        SF_YSPACING(sf), YBS(yb), lefty)
+	    czptr = SF_COEFF(sf) - 1 + lefty + leftx * SF_NYCOEFF(sf)
+
+	case SF_SPLINE1:
+	    call sf_b1spline1 (x, SF_NXPIECES(sf), -SF_XMIN(sf),
+	        SF_XSPACING(sf), XBS(xb), leftx)
+	    call sf_b1spline1 (y, SF_NYPIECES(sf), -SF_YMIN(sf),
+	        SF_YSPACING(sf), YBS(yb), lefty)
+	    czptr = SF_COEFF(sf) - 1 + lefty + leftx * SF_NYCOEFF(sf)
+	}
+
+	# initialize accumulator
+	# basis functions
+	sum = 0.
+
+	# loop over y basis functions
+	yorder = SF_YORDER(sf)
+	do i = 1, SF_XORDER(sf) {
+
+	    # loop over the x basis functions
+	    accum = 0.
+	    do k = 1, yorder {
+		accum = accum + COEFF(czptr+k) * YBS(yzb+k) 
+	    }
+	    accum = accum * XBS(xzb+i)
+	    sum = sum + accum
+
+	    # elements of COEFF where neither k = 1 or i = 1
+	    # are not calculated if SF_XTERMS(sf) = NO
+	    if (SF_XTERMS(sf) == NO)
+		yorder = 1
+
+	    czptr = czptr + SF_NYCOEFF(sf)
+	}
+
+	call sfree (sp)
+
+	return (sum)
+end