From 40e5a5811c6ffce9b0974e93cdd927cbcf60c157 Mon Sep 17 00:00:00 2001
From: Joe Hunkeler <jhunkeler@gmail.com>
Date: Tue, 11 Aug 2015 16:51:37 -0400
Subject: Repatch (from linux) of OSX IRAF

---
 math/iminterp/ii_sinctable.x | 123 +++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 123 insertions(+)
 create mode 100644 math/iminterp/ii_sinctable.x

(limited to 'math/iminterp/ii_sinctable.x')

diff --git a/math/iminterp/ii_sinctable.x b/math/iminterp/ii_sinctable.x
new file mode 100644
index 00000000..e062e9d0
--- /dev/null
+++ b/math/iminterp/ii_sinctable.x
@@ -0,0 +1,123 @@
+include <math.h>
+
+# II_SINCTABLE -- Compute the 1D sinc function look-up tables given the
+# width of the sinc function and the number of increments.
+
+procedure ii_sinctable (table, nconv, nincr, xshift)
+
+real	table[nconv,nincr]		#O the computed look-up table
+int	nconv				#I the sinc truncation length
+int	nincr				#I the number of look-up tables
+real	xshift				#I the shift of the look up table
+
+int	i, j, nsinc
+real	sconst, a2, a4, fsign, xsign, sum, dx, dx2, x, f
+
+begin
+	# Set up some constants.
+	nsinc = (nconv - 1) / 2
+	sconst = (HALFPI / nsinc) ** 2
+	a2 = -0.49670
+	a4 = 0.03705
+	if (mod (nsinc, 2) == 0)
+	    fsign = 1.0
+        else
+	    fsign = -1.0
+
+	# Create a one entry look-up table.
+	if (! IS_INDEFR(xshift)) {
+
+	    dx = xshift
+	    x = -nsinc - dx
+	    xsign = fsign
+	    sum = 0.0
+	    do j = 1, nconv {
+		if (x == 0.0)
+		    f = 1.0
+		else if (dx == 0.0)
+		    f = 0.0
+		else {
+		    dx2 = sconst * (j - nsinc - 1) ** 2
+		    f = xsign / x * (1.0 + a2 * dx2 + a4 * dx2 * dx2) ** 2
+		}
+		table[j,1] = f
+		sum = sum + f
+		x = x + 1.0
+		xsign = -xsign
+	    }
+	    do j = 1, nconv
+	        table[j,1] = table[j,1] / sum
+
+	# Create a multi-entry evenly spaced look-up table.
+	} else {
+
+	    do i = 1, nincr {
+	        dx = -0.5 + real (i - 1) / real (nincr - 1)
+	        x = -nsinc + dx
+		xsign = fsign
+	        sum = 0.0
+	        do j = 1, nconv {
+		    if ((x >= - 0.5 / (nincr - 1)) && (x < 0.5 / (nincr - 1)))
+		        f = 1.0
+		    else if ((dx >= -0.5 / (nincr - 1)) &&
+		        (dx < 0.5 / (nincr - 1)))
+		        f = 0.0
+		    else {
+		        dx2 = sconst * (j - nsinc - 1) ** 2
+		        f = xsign / x * (1.0 + a2 * dx2 + a4 * dx2 * dx2) ** 2
+		    }
+		    table[j,i] = f
+		    sum = sum + f
+		    x = x + 1.0
+		    xsign = -xsign
+	        }
+	        do j = 1, nconv
+		    table[j,i] = table[j,i] / sum
+	    }
+	}
+end
+
+
+# II_BISINCTABLE -- Compute the 2D sinc function look-up tables given the
+# width of the sinc function and the number of increments.
+
+procedure ii_bisinctable (table, nconv, nxincr, nyincr, xshift, yshift)
+
+real	table[nconv,nconv,nxincr,nyincr] #O the computed look-up table
+int	nconv				 #I the sinc truncation length
+int	nxincr, nyincr			 #I the number of look-up tables
+real	xshift, yshift			 #I the shift of the look up table
+
+int	j, ii, jj
+pointer	sp, fx, fy
+
+begin
+	# Allocate some working memory.
+	call smark (sp)
+	call salloc (fx, nconv * nxincr, TY_REAL)
+	call salloc (fy, nconv * nyincr, TY_REAL)
+
+	# Create a one entry look-up table.
+	if (! IS_INDEFR(xshift) && ! IS_INDEFR(yshift)) {
+
+	    call ii_sinctable (Memr[fx], nconv, 1, xshift)
+	    call ii_sinctable (Memr[fy], nconv, 1, yshift)
+	    do j = 1, nconv {
+		call amulkr (Memr[fx], Memr[fy+j-1], table[1,j,1,1], nconv)
+	    }
+
+	} else {
+
+	    call ii_sinctable (Memr[fx], nconv, nxincr, xshift)
+	    call ii_sinctable (Memr[fy], nconv, nyincr, yshift)
+	    do jj = 1, nyincr {
+		do ii = 1, nxincr {
+		    do j = 1, nconv
+			call amulkr (Memr[fx+(ii-1)*nconv],
+			    Memr[fy+(jj-1)*nconv+j-1], table[1,j,ii,jj], nconv)
+		}
+	    }
+	}
+
+	call sfree (sp)
+end
-- 
cgit