Repatch (from linux) of OSX IRAF

1 files changed, 154 insertions, 0 deletions

diff --git a/math/iminterp/asider.x b/math/iminterp/asider.x
new file mode 100644
index 00000000..73e93b4d
--- /dev/null
+++ b/math/iminterp/asider.x
@@ -0,0 +1,154 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <math/iminterp.h>
+include "im1interpdef.h"
+
+# ASIDER -- Calculate nder derivatives assuming that x lands in the region
+# 1 <= x <= npts.
+
+procedure asider (asi, x, der, nder)
+
+pointer	asi		# interpolant descriptor
+real	x[ARB]		# x value
+real	der[ARB]	# derivatives, der[1] is value der[2] is f prime 
+int	nder		# number items returned = 1 + number of derivatives
+
+int	nearx, i, j, k, nterms, nd
+pointer	c0ptr, n0
+real	deltax, accum, tmpx[2], pcoeff[MAX_NDERIVS], diff[MAX_NDERIVS]
+
+begin
+	# Return zero for derivatives that are zero.
+	do i = 1, nder
+	    der[i] = 0.
+
+	# Nterms is number of terms in case polynomial type.
+	nterms = 0
+
+	# (c0ptr + 1) is the pointer to the first data point in COEFF.
+	c0ptr = ASI_COEFF(asi) - 1 + ASI_OFFSET(asi)
+
+	switch (ASI_TYPE(asi))	{
+
+	case II_NEAREST:
+	    der[1] = COEFF(c0ptr + int(x[1] + 0.5))
+	    return
+
+	case II_LINEAR:
+	    nearx = x[1]
+	    der[1] = (x[1] - nearx) * COEFF(c0ptr + nearx + 1) +
+		     (nearx + 1 - x[1]) * COEFF(c0ptr + nearx)
+	    if (nder > 1)
+		der[2] = COEFF(c0ptr + nearx + 1) - COEFF(c0ptr + nearx)
+	    return
+
+	case II_SINC, II_LSINC:
+	    call ii_sincder (x[1], der, nder,
+		COEFF(ASI_COEFF(asi) + ASI_OFFSET(asi)), ASI_NCOEFF(asi),
+		ASI_NSINC(asi), DX)
+	    return
+
+	case II_DRIZZLE:
+	    if (ASI_PIXFRAC(asi) >= 1.0)
+	        call ii_driz1 (x, der[1], 1, COEFF(ASI_COEFF(asi) +
+	            ASI_OFFSET(asi)), ASI_BADVAL(asi))
+	    else
+	        call ii_driz (x, der[1], 1, COEFF(ASI_COEFF(asi) +
+	            ASI_OFFSET(asi)), ASI_PIXFRAC(asi), ASI_BADVAL(asi))
+	    if (nder > 1) {
+		deltax = x[2] - x[1]
+		if (deltax == 0.0)
+		    der[2] = 0.0
+		else {
+		    tmpx[1] = x[1]
+		    tmpx[2] = (x[1] + x[2]) / 2.0
+	    	    if (ASI_PIXFRAC(asi) >= 1.0)
+	                call ii_driz1 (tmpx, accum, 1, COEFF(ASI_COEFF(asi) +
+	                    ASI_OFFSET(asi)), ASI_BADVAL(asi))
+		    else
+	                call ii_driz (tmpx, accum, 1, COEFF(ASI_COEFF(asi) +
+	                    ASI_OFFSET(asi)), ASI_PIXFRAC(asi), ASI_BADVAL(asi))
+		    tmpx[1] = tmpx[2]
+		    tmpx[2] = x[2]
+	    	    if (ASI_PIXFRAC(asi) >= 1.0)
+	                call ii_driz1 (tmpx, der[2], 1, COEFF(ASI_COEFF(asi) +
+	                    ASI_OFFSET(asi)), ASI_BADVAL(asi))
+		    else
+	                call ii_driz (tmpx, der[2], 1, COEFF(ASI_COEFF(asi) +
+	                    ASI_OFFSET(asi)), ASI_PIXFRAC(asi), ASI_BADVAL(asi))
+		    der[2] = 2.0 * (der[2] - accum) / deltax
+		}
+	    }
+	    return
+
+	case II_POLY3:
+	    nterms = 4
+
+	case II_POLY5:
+	    nterms = 6
+
+	case II_SPLINE3:
+	    nterms = 4
+
+	default:
+	    call error (0, "ASIDER: Unknown interpolant type")
+	}
+
+	# Routines falls through to this point if the interpolant is one of
+	# the higher order polynomial types or a third order spline.
+
+	nearx = x[1]
+	n0 = c0ptr + nearx
+	deltax = x[1] - nearx
+
+	# Compute the number of derivatives needed.
+	nd = nder
+	if (nder > nterms)
+	    nd = nterms
+
+	# Generate the polynomial coefficients.
+
+	if (ASI_TYPE(asi) == II_SPLINE3) {
+
+	    pcoeff[1] = COEFF(n0-1) + 4. * COEFF(n0) + COEFF(n0+1)
+	    pcoeff[2] = 3. * (COEFF(n0+1) - COEFF(n0-1))
+	    pcoeff[3] = 3. * (COEFF(n0-1) - 2. * COEFF(n0) + COEFF(n0+1))
+	    pcoeff[4] = -COEFF(n0-1) + 3. * COEFF(n0) - 3. * COEFF(n0+1) +
+				    COEFF(n0+2)
+
+	# Newton's form written in line to get polynomial from data
+	} else {
+
+	    # Load data.
+	    do i = 1, nterms
+		diff[i] = COEFF(n0 - nterms/2 + i)
+
+	    # Generate difference table.
+	    do k = 1, nterms - 1
+		do i = 1, nterms - k
+		    diff[i] = (diff[i+1] - diff[i]) / k
+
+	    # Shift to generate polynomial coefficients.
+	    do k = nterms, 2, -1
+		do i = 2,k
+		    diff[i] = diff[i] + diff[i-1] * (k - i - nterms/2)
+	    do i = 1,nterms
+		pcoeff[i] = diff[nterms + 1 - i]
+	}
+
+	# Compute the derivatives. As the loop progresses pcoeff contains
+	# coefficients of higher and higher derivatives.
+
+	do k = 1, nd {
+
+	    # Evaluate using nested multiplication.
+	    accum = pcoeff[nterms - k + 1]
+	    do j = nterms - k, 1, -1
+		accum = pcoeff[j] + deltax * accum
+	    der[k] = accum
+
+	    # Differentiate polynomial.
+	    do j = 1, nterms - k
+		pcoeff[j] = j * pcoeff[j + 1]
+	}
+end