Initial commit

1 files changed, 859 insertions, 0 deletions

diff --git a/math/iminterp/ii_greval.x b/math/iminterp/ii_greval.x
new file mode 100644
index 00000000..6b5d75b2
--- /dev/null
+++ b/math/iminterp/ii_greval.x
@@ -0,0 +1,859 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <math.h>
+include <math/iminterp.h>
+
+# II_GRNEAREST -- Procedure to evaluate the nearest neighbour interpolant on
+# a rectangular grid. The procedure assumes that 1 <= x <= nxpix and
+# that 1 <= y <= nypix. The x and y vectors must be sorted in increasing
+# value of x and y such that x[i] < x[i+1] and y[i] < y[i+1]. The routine
+# outputs a grid of nxpix by nypix points using the coeff array where
+# coeff[1+first_point] = datain[1,1]
+
+procedure ii_grnearest (coeff, first_point, len_coeff, x, y, zfit, nxpts,
+		        nypts, len_zfit)
+
+real	coeff[ARB]		# 1D coefficient array
+int	first_point		# offset of first data point
+int	len_coeff		# row length of coeff 
+real	x[nxpts]		# array of x values
+real	y[nypts]		# array of y values
+real	zfit[len_zfit,ARB]	# array of interpolatedvalues
+int	nxpts			# number of x values
+int	nypts			# number of y values
+int	len_zfit		# row length of zfit
+
+int	ny
+int	index
+int	i, j
+pointer	sp, nx
+
+errchk	smark, salloc, sfree
+
+begin
+	call smark (sp)
+	call salloc (nx, nxpts, TY_INT)
+
+	# calculate the nearest x
+	do i = 1, nxpts
+	    Memi[nx+i-1] = x[i] + 0.5
+
+	# loop over the rows
+	do j = 1, nypts {
+
+	    # calculate pointer to the ny-th row of data
+	    ny = y[j] + 0.5
+	    index = first_point + (ny - 1) * len_coeff
+
+	    # loop over the columns
+	    do i = 1, nxpts
+		zfit[i,j] = coeff[index + Memi[nx+i-1]]
+	}
+
+	call sfree (sp)
+end
+
+
+# II_GRLINEAR -- Procedure to evaluate the bilinear interpolant
+# on a rectangular grid. The procedure assumes that 1 <= x <= nxpix and that
+# 1 <= y <= nypix. The x and y vectors are assumed to be sorted in increasing
+# order of x and y such that x[i] < x[i+1] and y[i] < y[i+1]. The routine 
+# produces a grid of nxpix * nypix evaluated points using the coeff array
+# where coeff[1+first_point] = datain[1,1].
+
+procedure ii_grlinear (coeff, first_point, len_coeff, x, y, zfit, nxpts,
+		       nypts, len_zfit)
+
+real	coeff[ARB]		# 1D array of coefficients
+int	first_point		# offset of first data point
+int	len_coeff		# row length of coeff
+real	x[nxpts]		# array of x values
+real	y[nypts]		# array of y values
+real	zfit[len_zfit,ARB]	# array of interpolated values
+int	nxpts			# number of x values
+int	nypts			# number of y values
+int	len_zfit		# row length of zfit
+
+int	i, j, ny
+int	nymin, nymax, nylines
+int	row_index, xindex
+pointer	sp, nx, sx, tx, work, lbuf1, lbuf2
+real	sy, ty
+
+errchk	smark, salloc, sfree
+
+begin
+	# calculate the x and y limits
+	nymin = y[1]
+	nymax = int (y[nypts]) + 1
+	nylines = nymax - nymin + 1
+
+	# allocate storage for work array
+	call smark (sp)
+	call salloc (nx, nxpts, TY_INT)
+	call salloc (sx, nxpts, TY_REAL)
+	call salloc (tx, nxpts, TY_REAL)
+	call salloc (work, nxpts * nylines, TY_REAL)
+
+	# initialize
+	call achtri (x, Memi[nx], nxpts)
+	do i = 1, nxpts {
+	    Memr[sx+i-1] = x[i] - Memi[nx+i-1]
+	    Memr[tx+i-1] = 1. - Memr[sx+i-1]
+	}
+
+	# for each value of y interpolate in x and store in work array
+	lbuf1 = work
+	do j = 1, nylines {
+
+	    # define pointer to appropriate row
+	    row_index = first_point + (j + nymin - 2) * len_coeff
+
+	    # interpolate in x at each y
+	    do i = 1, nxpts {
+		xindex = row_index + Memi[nx+i-1]
+		Memr[lbuf1+i-1] = Memr[tx+i-1] * coeff[xindex] +
+		    Memr[sx+i-1] * coeff[xindex+1]
+	    }
+
+	    lbuf1 = lbuf1 + nxpts
+	}
+
+	# at each x interpolate in y and store in temporary work array
+	do j = 1, nypts {
+
+	    ny = y[j]
+	    sy = y[j] - ny
+	    ty = 1. - sy
+
+	    lbuf1 = work + nxpts * (ny - nymin)
+	    lbuf2 = lbuf1 + nxpts
+
+	    call awsur (Memr[lbuf1], Memr[lbuf2], zfit[1,j], nxpts,
+		ty, sy)
+
+	}
+
+	# deallocate work space
+	call sfree (sp)
+end
+
+
+# II_GRPOLY3 -- Procedure to evaluate the bicubic polynomial interpolant
+# on a rectangular grid. The points to be evaluated are assumed
+# to lie in the range 1 <= x <= nxpix and 1 <= y <= nypix. The x and y vectors
+# are assumed to be sorted in increasing order of x and y such that
+# x[i] < x[i+1] and y[i] < y[i+1]. The interpolation is done using
+# Everett's central difference formula and separation of variables
+# and assuming that coeff[1+first_point] = datain[1,1].
+
+procedure ii_grpoly3 (coeff, first_point, len_coeff, x, y, zfit, nxpts, nypts,
+		      len_zfit)
+
+real	coeff[ARB]		# 1D array of coefficients
+int	first_point		# offset of first data point
+int	len_coeff		# length of row of coeffcient
+real	x[nxpts]		# array of x values
+real	y[nypts]		# array of y values
+real	zfit[len_zfit,ARB]	# array of interpolatedvalues
+int	nxpts			# number of x points
+int	nypts			# number of y points
+int	len_zfit		# row length of zfit
+
+int	nymin, nymax, nylines
+int	nxold, nyold
+int	row_index, xindex
+int	i, j, ny
+pointer	sp, nx, sx, sx2m1, tx, tx2m1, work
+pointer	lbuf, lbufp1, lbufp2, lbufm1
+real	cd20x, cd21x, cd20y, cd21y
+real	sy, ty, sy2m1, ty2m1
+
+errchk smark, salloc, sfree
+
+begin
+	# find y limits
+	nymin = int (y[1]) - 1
+	nymax = int (y[nypts]) + 2
+	nylines = nymax - nymin + 1
+
+	# allocate work space
+	call smark (sp)
+	call salloc (nx, nxpts, TY_INT)
+	call salloc (sx, nxpts, TY_REAL)
+	call salloc (sx2m1, nxpts, TY_REAL)
+	call salloc (tx, nxpts, TY_REAL)
+	call salloc (tx2m1, nxpts, TY_REAL) 
+	call salloc (work, nxpts * nylines, TY_REAL)
+
+	# initialize
+	call achtri (x, Memi[nx], nxpts)
+	do i = 1, nxpts {
+	    Memr[sx+i-1] = x[i] - Memi[nx+i-1]
+	    Memr[sx2m1+i-1] = Memr[sx+i-1] * Memr[sx+i-1] - 1.
+	    Memr[tx+i-1] = 1. - Memr[sx+i-1]
+	    Memr[tx2m1+i-1] = Memr[tx+i-1] * Memr[tx+i-1] - 1.
+	}
+
+	# for each value of y interpolate in x
+	lbuf = work
+	do j = 1, nylines {
+
+	    # calculate pointer to a row
+	    row_index = first_point + (j + nymin - 2) * len_coeff
+
+	    # interpolate in x at each y
+	    nxold = -1
+	    do i = 1, nxpts {
+
+		xindex= row_index + Memi[nx+i-1]
+
+		if (Memi[nx+i-1] != nxold) {
+		    #cd20x = 1./6. * (coeff[xindex+1] - 2. * coeff[xindex] +
+		           #coeff[xindex-1])
+		    #cd21x = 1./6. * (coeff[xindex+2] - 2. * coeff[xindex+1] +
+			    #coeff[xindex])
+		    cd20x = (coeff[xindex+1] - 2. * coeff[xindex] +
+		           coeff[xindex-1]) / 6.
+		    cd21x = (coeff[xindex+2] - 2. * coeff[xindex+1] +
+			    coeff[xindex]) / 6.0
+		}
+
+		Memr[lbuf+i-1] = Memr[sx+i-1] * (coeff[xindex+1] +
+		                Memr[sx2m1+i-1] * cd21x) +
+			        Memr[tx+i-1] * (coeff[xindex] +
+				Memr[tx2m1+i-1] * cd20x)
+
+		nxold = Memi[nx+i-1]
+	    }
+
+	    lbuf = lbuf + nxpts
+	}
+
+	# interpolate in y at each x
+	nyold = -1
+	do j = 1, nypts {
+
+	    ny = y[j]
+	    sy = y[j] - ny
+	    ty = 1. - sy
+	    sy2m1 = sy ** 2 - 1.
+	    ty2m1 = ty ** 2 - 1.
+
+	    lbuf = work + nxpts * (ny - nymin) 
+	    lbufm1 = lbuf - nxpts
+	    lbufp1 = lbuf + nxpts
+	    lbufp2 = lbufp1 + nxpts
+
+	    do i = 1, nxpts {
+
+		# calculate central differences in y
+		#if (nyold != ny) {
+		    #cd20y = 1./6. * (Memr[lbufp1+i-1] - 2. * Memr[lbuf+i-1] +
+			    #Memr[lbufm1+i-1])
+		    #cd21y = 1./6. * (Memr[lbufp2+i-1] - 2. *
+		    	    #Memr[lbufp1+i-1] + Memr[lbuf+i-1])
+		    cd20y = (Memr[lbufp1+i-1] - 2. * Memr[lbuf+i-1] +
+			    Memr[lbufm1+i-1]) / 6.0
+		    cd21y = (Memr[lbufp2+i-1] - 2. * Memr[lbufp1+i-1] +
+		    	    Memr[lbuf+i-1]) / 6.0
+		#}
+
+		# interpolate in y
+		zfit[i,j] =  sy * (Memr[lbufp1+i-1] + sy2m1 * cd21y) +
+			     ty * (Memr[lbuf+i-1] + ty2m1 * cd20y)
+
+	    }
+
+	    #nyold = ny
+	}
+
+	# free work space
+	call sfree (sp)
+end
+
+
+# II_GRPOLY5 -- Procedure to evaluate the biquintic polynomial interpolant
+# on a rectangular grid. The routine assumes that 1 <= x <= nxpix and
+# 1 <= y <= nypix.  The vectors x and y are assumed to be sorted in
+# increasing order such that x[i] < x[i+1] and y[i] < y[i+1]. The
+# interpolation is done using Everett's interpolation formula and
+# separation of variables and assuming that coeff[1+first_point] =
+# datain[1,1].
+
+procedure ii_grpoly5 (coeff, first_point, len_coeff, x, y, zfit, nxpts,
+		      nypts, len_zfit)
+
+real	coeff[ARB]		# 1D array of coefficients
+int	first_point		# offset of first data point
+int	len_coeff		# row length of coeff
+real	x[nxpts]		# array of x values
+real	y[nypts]		# array of y values
+real	zfit[len_zfit,ARB]	# array of fitted values
+int	nxpts			# number of x points
+int	nypts			# number of y points
+int	len_zfit		# row length of zfit
+
+int	nymax, nymin, nylines, nxold, nyold
+int	row_index, xindex
+int	i, j, ny
+pointer	sp, nx, sx, tx, sx2m1, sx2m4, tx2m1, tx2m4, work
+pointer	lbuf, lbufp1, lbufp2, lbufp3, lbufm1, lbufm2
+real	cd20x, cd21x, cd40x, cd41x
+real	cd20y, cd21y, cd40y, cd41y
+real	sy, ty, sy2m1, sy2m4, ty2m1, ty2m4
+
+errchk	smark, salloc, sfree
+
+begin
+	# find the y limits
+	nymin = int (y[1]) - 2
+	nymax = int (y[nypts]) + 3
+	nylines = nymax - nymin + 1
+
+	# allocate space
+	call smark (sp)
+	call salloc (nx, nxpts, TY_INT)
+	call salloc (sx, nxpts, TY_REAL)
+	call salloc (sx2m1, nxpts, TY_REAL)
+	call salloc (sx2m4, nxpts, TY_REAL)
+	call salloc (tx, nxpts, TY_REAL)
+	call salloc (tx2m1, nxpts, TY_REAL)
+	call salloc (tx2m4, nxpts, TY_REAL)
+	call salloc (work, nxpts * nylines, TY_REAL)
+
+	# intialize
+	call achtri (x, Memi[nx], nxpts)
+	do i = 1, nxpts {
+	    Memr[sx+i-1] = x[i] - Memi[nx+i-1]
+	    Memr[sx2m1+i-1] = Memr[sx+i-1] ** 2 - 1.
+	    Memr[sx2m4+i-1] = Memr[sx2m1+i-1] - 3.
+	    Memr[tx+i-1] = 1. - Memr[sx+i-1]
+	    Memr[tx2m1+i-1] = Memr[tx+i-1] ** 2 - 1.
+	    Memr[tx2m4+i-1] = Memr[tx2m1+i-1] - 3.
+	}
+
+
+	# for each value of y interpolate in x
+	lbuf = work
+	do j = 1, nylines {
+
+	    # calculate pointer to a row
+	    row_index = first_point + (j + nymin - 2) * len_coeff
+
+	    # interpolate in x at each y
+	    nxold = -1
+	    do i = 1, nxpts {
+
+		xindex = row_index + Memi[nx+i-1]
+
+		if (Memi[nx+i-1] != nxold) {
+		    #cd20x = 1./6. * (coeff[xindex+1] - 2. * coeff[xindex] +
+		            #coeff[xindex-1])
+		    #cd21x = 1./6. * (coeff[xindex+2] - 2. * coeff[xindex+1] +
+			    #coeff[xindex])
+		    cd20x = (coeff[xindex+1] - 2. * coeff[xindex] +
+		            coeff[xindex-1]) / 6.0
+		    cd21x = (coeff[xindex+2] - 2. * coeff[xindex+1] +
+			    coeff[xindex]) / 6.0
+		    #cd40x = 1./120. * (coeff[xindex-2] - 4. * coeff[xindex-1] +
+			    #6. * coeff[xindex] - 4. * coeff[xindex+1] +
+			    #coeff[xindex+2])
+		    #cd41x = 1./120. * (coeff[xindex-1] - 4. * coeff[xindex] +
+			    #6. * coeff[xindex+1] - 4. * coeff[xindex+2] +
+			    #coeff[xindex+3])
+		    cd40x = (coeff[xindex-2] - 4. * coeff[xindex-1] +
+			    6. * coeff[xindex] - 4. * coeff[xindex+1] +
+			    coeff[xindex+2]) / 120.0
+		    cd41x = (coeff[xindex-1] - 4. * coeff[xindex] +
+			    6. * coeff[xindex+1] - 4. * coeff[xindex+2] +
+			    coeff[xindex+3]) / 120.0
+		}
+
+		Memr[lbuf+i-1] = Memr[sx+i-1] * (coeff[xindex+1] +
+				 Memr[sx2m1+i-1] * (cd21x +
+				 Memr[sx2m4+i-1] * cd41x)) +
+			         Memr[tx+i-1] * (coeff[xindex] +
+				 Memr[tx2m1+i-1] * (cd20x +
+				 Memr[tx2m4+i-1] * cd40x))
+
+
+		nxold = Memi[nx+i-1]
+
+	    }
+
+	    lbuf = lbuf + nxpts
+	}
+
+	# at each x interpolate in y
+	nyold = -1
+	do j = 1, nypts {
+
+	    ny = y[j]
+	    sy = y[j] - ny
+	    sy2m1 = sy ** 2 - 1.
+	    sy2m4 = sy2m1 - 3.
+	    ty = 1. - sy
+	    ty2m1 = ty ** 2 - 1.
+	    ty2m4 = ty2m1 - 3.
+
+	    lbuf = work + nxpts * (ny - nymin)
+	    lbufp1 = lbuf + nxpts
+	    lbufp2 = lbufp1 + nxpts
+	    lbufp3 = lbufp2 + nxpts
+	    lbufm1 = lbuf - nxpts
+	    lbufm2 = lbufm1 - nxpts
+
+	    do i = 1, nxpts {
+
+		# calculate central differences
+		#if (nyold != ny) {
+		    #cd20y = 1./6. * (Memr[lbufp1+i-1] - 2. * Memr[lbuf+i-1] +
+			    #Memr[lbufm1+i-1])
+		    #cd21y = 1./6. * (Memr[lbufp2+i-1] - 2. *
+		    	    #Memr[lbufp1+i-1] + Memr[lbuf+i-1])
+		    cd20y = (Memr[lbufp1+i-1] - 2. * Memr[lbuf+i-1] +
+			    Memr[lbufm1+i-1]) / 6.
+		    cd21y = (Memr[lbufp2+i-1] - 2. *
+		    	    Memr[lbufp1+i-1] + Memr[lbuf+i-1]) / 6.
+		    #cd40y = 1./120. * (Memr[lbufm2+i-1] -
+		    	    #4. * Memr[lbufm1+i-1] + 6. * Memr[lbuf+i-1] -
+			    #4. * Memr[lbufp1+i-1] + Memr[lbufp2+i-1])
+		    #cd41y = 1./120. * (Memr[lbufm1+i-1] - 4. * 
+		    	    #Memr[lbuf+i-1] + 6. * Memr[lbufp1+i-1] - 4. *
+			    #Memr[lbufp2+i-1] + Memr[lbufp3+i-1])
+		    cd40y = (Memr[lbufm2+i-1] -
+		    	    4. * Memr[lbufm1+i-1] + 6. * Memr[lbuf+i-1] -
+			    4. * Memr[lbufp1+i-1] + Memr[lbufp2+i-1]) / 120.
+		    cd41y = (Memr[lbufm1+i-1] - 4. * 
+		    	    Memr[lbuf+i-1] + 6. * Memr[lbufp1+i-1] - 4. *
+			    Memr[lbufp2+i-1] + Memr[lbufp3+i-1]) / 120.0
+		#}
+
+		# interpolate in y
+		zfit[i,j] = sy * (Memr[lbufp1+i-1] + sy2m1 *
+			    (cd21y + sy2m4 * cd41y)) +
+			    ty * (Memr[lbuf+i-1] + ty2m1 *
+			    (cd20y + ty2m4 * cd40y))
+
+	    }
+
+	    #nyold = ny
+	}
+
+	# release work space
+	call sfree (sp)
+
+end
+
+
+# II_GRSPLINE3 -- Procedure to evaluate the bicubic spline interpolant
+# on a rectangular grid.  The program assumes that 1 <= x <= nxpix and
+# 1 <= y <= nypix. The routine assumes that x and y vectors are sorted
+# such that x[i] < x[i+1] and y[i] < y[i+1]. The interpolant is evaluated
+# by calculating the polynomial coefficients in x and y.
+
+procedure ii_grspline3 (coeff, first_point, len_coeff, x, y, zfit, nxpts,
+			nypts, len_zfit)
+
+real	coeff[ARB]		# 1D array of coefficients
+int	first_point		# offset of first data point
+int	len_coeff		# row length of coeff
+real	x[nxpts]		# array of x values
+real	y[nypts]		# array of y values
+real	zfit[len_zfit,ARB]	# array of interpolated values
+int	nxpts			# number of x values
+int	nypts			# number of y values
+int	len_zfit		# row length of zfit
+
+int	ny, nymin, nymax, nylines
+int	row_index, xindex
+int	i, j
+pointer	sp, nx, sx, tx, sx3, tx3, work, lbuf, lbufp1, lbufp2, lbufm1
+real	sy, ty, ty3, sy3
+
+errchk	smark, salloc, sfree
+
+begin
+	# find the y limits
+	nymin = int (y[1]) - 1
+	nymax = int (y[nypts]) + 2
+	nylines = nymax - nymin + 1
+
+	# allocate space for work array
+	call smark (sp)
+	call salloc (nx, nxpts, TY_INT)
+	call salloc (sx, nxpts, TY_REAL)
+	call salloc (sx3, nxpts, TY_REAL)
+	call salloc (tx, nxpts, TY_REAL)
+	call salloc (tx3, nxpts, TY_REAL)
+	call salloc (work, nylines * nxpts, TY_REAL)
+
+	# intialize
+	call achtri (x, Memi[nx], nxpts)
+	do j = 1, nxpts {
+	    Memr[sx+j-1] = x[j] - Memi[nx+j-1]
+	    Memr[tx+j-1] = 1. - Memr[sx+j-1]
+	}
+	call apowkr (Memr[sx], 3, Memr[sx3], nxpts)
+	call apowkr (Memr[tx], 3, Memr[tx3], nxpts)
+	do j = 1, nxpts {
+	    Memr[sx+j-1] = 1. + Memr[sx+j-1] * (3. + Memr[sx+j-1] *
+			   (3. - 3. * Memr[sx+j-1])) 
+	    Memr[tx+j-1] = 1. + Memr[tx+j-1] * (3. + Memr[tx+j-1] *
+			   (3. -  3. * Memr[tx+j-1])) 
+	}
+
+	# interpolate in x for each y
+	lbuf = work
+	do i = 1, nylines {
+
+	    # find appropriate row
+	    row_index = first_point + (i + nymin - 2) * len_coeff
+
+	    # x interpolation
+	    do j = 1, nxpts {
+
+		xindex = row_index + Memi[nx+j-1]
+		Memr[lbuf+j-1] = Memr[tx3+j-1] * coeff[xindex-1] +
+				 Memr[tx+j-1] * coeff[xindex] +
+				 Memr[sx+j-1] * coeff[xindex+1] +
+				 Memr[sx3+j-1] * coeff[xindex+2]
+	    }
+	    lbuf = lbuf + nxpts
+	}
+
+	# interpolate in y
+	do i = 1, nypts {
+
+	    ny = y[i]
+	    sy = y[i] - ny
+	    ty = 1. - sy
+	    sy3 = sy ** 3
+	    ty3 = ty ** 3
+	    sy = 1. + sy * (3. + sy * (3. - 3. * sy))
+	    ty = 1. + ty * (3. + ty * (3. - 3. * ty))
+
+	    lbuf = work + nxpts * (ny - nymin)
+	    lbufp1 = lbuf + nxpts
+	    lbufp2 = lbufp1 + nxpts
+	    lbufm1 = lbuf - nxpts
+
+	    do j = 1, nxpts
+		zfit[j,i] = ty3 * Memr[lbufm1+j-1] + ty * Memr[lbuf+j-1] +
+			    sy * Memr[lbufp1+j-1] + sy3 * Memr[lbufp2+j-1]
+	}
+
+	# release working space
+	call sfree (sp)
+end
+
+# II_GRSINC -- Procedure to evaluate the sinc interpolant on a rectangular
+# grid. The procedure assumes that 1 <= x <= nxpix and  that 1 <= y <= nypix.
+# The x and y vectors must be sorted in increasing value of x and y such that
+# x[i] < x[i+1] and y[i] < y[i+1]. The routine outputs a grid of nxpix by
+# nypix points using the coeff array where coeff[1+first_point] = datain[1,1]
+# The sinc truncation length is nsinc. The taper is a cosine bell function
+# which is approximated by a quartic polynomial which is valid for 0  <= x
+# <= PI / 2 (Abramowitz and Stegun 1972, Dover Publications, p 76). If the
+# point to be interpolated is less than mindx and mindy from a data point
+# no interpolation is done and the data point itself is returned.
+
+procedure ii_grsinc (coeff, first_point, len_coeff, len_array, x, y, zfit,
+		nxpts, nypts, len_zfit, nsinc, mindx, mindy)
+
+real	coeff[ARB]		# 1D coefficient array
+int	first_point		# offset of first data point
+int	len_coeff		# row length of coeff 
+int	len_array		# column length of coeff 
+real	x[nxpts]		# array of x values
+real	y[nypts]		# array of y values
+real	zfit[len_zfit,ARB]	# array of interpolatedvalues
+int	nxpts			# number of x values
+int	nypts			# number of y values
+int	len_zfit		# row length of zfit
+int	nsinc			# sinc interpolant truncation length
+real	mindx, mindy		# the precision of the interpolant.
+
+int	i, j, k, nconv, nymin, nymax, nylines
+int	ixy, index, minj, maxj, offj
+pointer	sp, taper, ac, ixn, work, pac, pwork, ppwork
+real	sconst, a2, a4, dxy, dxyn, dx2, axy, pxy, sumxy, fdxy
+
+begin
+	# Compute the limits of the convolution in y.
+	nconv = 2 * nsinc + 1
+	nymin = max (1, nint (y[1]) - nsinc)
+	#nymin = max (1, int (y[1]) - nsinc)
+	nymax = min (len_array, nint (y[nypts]) + nsinc)
+	#nymax = min (len_array, int (y[nypts]) + nsinc)
+	nylines = nymax - nymin + 1
+
+	# Allocate working space.
+	call smark (sp)
+	call salloc (taper, nconv, TY_REAL)
+	call salloc (ac, nconv * max (nxpts, nypts), TY_REAL)
+	call salloc (ixn, max (nxpts, nypts), TY_INT)
+	call salloc (work, nxpts * nylines, TY_REAL)
+
+	# Compute the parameters of the cosine bell taper.
+	sconst = (HALFPI / nsinc) ** 2
+	a2 = -0.49670
+	a4 = 0.03705
+	if (mod (nsinc, 2) == 0)
+	    fdxy = 1.0
+	else
+	    fdxy = -1.0
+	do i = -nsinc, nsinc {
+	    dx2 = sconst * i * i
+	    Memr[taper+i+nsinc] = fdxy * (1.0 + a2 * dx2 + a4 * dx2 * dx2) ** 2
+	    fdxy = -fdxy
+	}
+
+	# Compute the x interpolants for each shift in x.
+	pac = ac
+	do i = 1, nxpts {
+	    ixy = nint (x[i])
+	    Memi[ixn+i-1] = ixy
+	    dxy = x[i] - ixy
+	    #dxyn = -1 - nsinc - dxy
+	    dxyn = 1 + nsinc + dxy
+	    sumxy = 0.0
+	    do j = 1, nconv {
+		#axy = j + dxyn
+		axy = dxyn - j
+		if (axy == 0.0)
+		    pxy = 1.0
+		else if (dxy == 0.0)
+		    pxy = 0.0
+		else
+		    pxy = Memr[taper+j-1] / axy 
+		Memr[pac+j-1] = pxy
+		sumxy = sumxy + pxy
+	    }
+	    call adivkr (Memr[pac], sumxy, Memr[pac], nconv)
+	    pac = pac + nconv
+	}
+
+	# Do the convolutions in the x direction.
+	pwork = work
+	do k = nymin, nymax {
+	    index = first_point + (k - 1) * len_coeff
+	    pac = ac
+	    do i = 1, nxpts {
+		sumxy = 0.0
+		ixy = Memi[ixn+i-1]
+		minj = max (1, ixy - nsinc)
+		maxj = min (len_coeff, ixy + nsinc)
+		offj = -ixy + nsinc 
+		do j = ixy - nsinc, minj - 1
+		    sumxy = sumxy + Memr[pac+j+offj] * coeff[index+1] 
+		do j = minj, maxj
+		    sumxy = sumxy + Memr[pac+j+offj] * coeff[index+j] 
+		do j = maxj + 1, ixy + nsinc
+		    sumxy = sumxy + Memr[pac+j+offj] * coeff[index+len_coeff] 
+	        Memr[pwork+i-1] = sumxy
+		pac = pac + nconv
+	    }
+	    pwork = pwork + nxpts
+	}
+
+	# Compute the y interpolants for each shift in y.
+	pac = ac
+	do i = 1, nypts {
+	    ixy = nint (y[i])
+	    dxy = y[i] - ixy
+	    Memi[ixn+i-1] = ixy - nsinc - nymin + 1
+	    #dxyn = -1 - nsinc - dxy
+	    dxyn = 1 + nsinc + dxy
+	    sumxy = 0.0
+	    do j = 1, nconv {
+		#axy = j + dxyn
+		axy = dxyn - j
+		if (axy == 0.0)
+		    pxy = 1.0
+		else if (dxy == 0.0)
+		    pxy = 0.0
+		else
+		    pxy = Memr[taper+j-1] / axy 
+		Memr[pac+j-1] = pxy
+		sumxy = sumxy + pxy
+	    }
+	    call adivkr (Memr[pac], sumxy, Memr[pac], nconv)
+	    pac = pac + nconv
+	}
+
+	# Do the interpolation in y.
+	do k = 1, nxpts {
+	    pwork = work + k - 1
+	    pac = ac
+	    do i = 1, nypts {
+		ixy = min (nylines, max (1, Memi[ixn+i-1]))
+		ppwork = pwork + (ixy - 1) * nxpts
+		sumxy = 0.0
+		do j = 1, nconv {
+		    sumxy = sumxy + Memr[pac+j-1] * Memr[ppwork]
+		    ppwork = ppwork + nxpts
+		}
+		pac = pac + nconv
+		zfit[k,i] = sumxy 
+	    }
+	}
+
+	call sfree (sp)
+end
+
+
+# II_GRLSINC -- Procedure to evaluate the sinc interpolant on a rectangular
+# grid. The procedure assumes that 1 <= x <= nxpix and  that 1 <= y <= nypix.
+# The x and y vectors must be sorted in increasing value of x and y such that
+# x[i] < x[i+1] and y[i] < y[i+1]. The routine outputs a grid of nxpix by
+# nypix points using the coeff array where coeff[1+first_point] = datain[1,1]
+# The sinc truncation length is nsinc. The taper is a cosine bell function
+# which is approximated by a quartic polynomial which is valid for 0  <= x
+# <= PI / 2 (Abramowitz and Stegun 1972, Dover Publications, p 76). If the
+# point to be interpolated is less than mindx and mindy from a data point
+# no interpolation is done and the data point itself is returned.
+
+procedure ii_grlsinc (coeff, first_point, len_coeff, len_array, x, y, zfit,
+		nxpts, nypts, len_zfit, ltable, nconv, nxincr, nyincr,
+		mindx, mindy)
+
+real	coeff[ARB]				# 1D coefficient array
+int	first_point				# offset of first data point
+int	len_coeff				# row length of coeff 
+int	len_array				# column length of coeff 
+real	x[nxpts]				# array of x values
+real	y[nypts]				# array of y values
+real	zfit[len_zfit,ARB]			# array of interpolated values
+int	nxpts					# number of x values
+int	nypts					# number of y values
+int	len_zfit				# row length of zfit
+real	ltable[nconv,nconv,nxincr,nyincr]	# pre-computed sinc lut
+int	nconv					# sinc trunction full-width
+int	nxincr, nyincr				# resolution of look-up table
+real	mindx, mindy				# the precision of interpolant
+
+int	j
+pointer	sp, ytmp
+
+begin
+	# Allocate working space.
+	call smark (sp)
+	call salloc (ytmp, nxpts, TY_REAL)
+
+	do j = 1, nypts {
+	    call amovkr (y[j], Memr[ytmp], nxpts)
+	    call ii_bilsinc (coeff, first_point, len_coeff, len_array, x,
+	        Memr[ytmp], zfit[1,j], nxpts, ltable, nconv, nxincr, nyincr,
+		mindx, mindy)
+	}
+
+	call sfree (sp)
+end
+
+
+# II_GRDRIZ -- Procedure to evaluate the drizzle interpolant on a rectangular
+# grid. The procedure assumes that the x and y intervals are ordered from
+# smallest to largest
+
+procedure ii_grdriz (coeff, first_point, len_coeff, len_array, x, y, zfit,
+	nxpts, nypts, len_zfit, xfrac, yfrac, badval)
+
+real	coeff[ARB]		# 1D coefficient array
+int	first_point		# offset of first data point
+int	len_coeff		# row length of coeff 
+int	len_array		# column length of coeff 
+real	x[ARB]			# array of x values
+real	y[ARB]			# array of y values
+real	zfit[len_zfit,ARB]	# array of interpolatedvalues
+int	nxpts			# number of x values
+int	nypts			# number of y values
+int	len_zfit		# row length of zfit
+real	xfrac, yfrac		# the x and y pixel fractions
+real	badval			# bad value
+
+int	i, j, jj, nylmin, nylmax, nylines
+int	cindex, neara, nearb
+pointer	sp, work, xindex
+real	ymin, ymax, dy, accum, waccum, hyfrac
+
+begin
+	ymin = min (y[1], y[2])
+	ymax = max (y[2*nypts-1], y[2*nypts])
+	nylmin = int (ymin + 0.5) 
+	nylmax = int (ymax + 0.5)
+	nylines = nylmax - nylmin + 1
+
+	call smark (sp)
+	call salloc (work, nxpts * nylines, TY_REAL)
+
+	# For each in range y integrate in x.
+	cindex = 1 + first_point + (nylmin - 1) * len_coeff
+	xindex = work
+	do j = nylmin, nylmax {
+	    if (xfrac >= 1.0)
+	        call ii_driz1 (x, Memr[xindex], nxpts, coeff[cindex], badval)
+	    else
+	        call ii_driz (x, Memr[xindex], nxpts, coeff[cindex], xfrac,
+		    badval)
+	    xindex = xindex + nxpts
+	    cindex = cindex + len_coeff
+	}
+
+	# For each range in x integrate in y. This may need to be vectorized?
+	hyfrac = yfrac / 2.0
+	do i = 1, nxpts {
+
+	    xindex = work + i - 1
+
+	    do j = 1, nypts {
+
+		ymin = min (y[2*j-1], y[2*j])
+		ymax = max (y[2*j-1], y[2*j])
+		neara = ymin + 0.5
+		nearb = ymax + 0.5
+
+		accum = 0.0
+		waccum = 0.0
+		if (neara == nearb) {
+
+		    dy = min (ymax, nearb + hyfrac) - max (ymin,
+		       neara - hyfrac)
+		    if (dy > 0.0) {
+			accum = accum + dy * Memr[xindex+(neara-nylmin)*nxpts]
+			waccum = waccum + dy
+		    }
+
+		} else {
+
+		    # First segment.
+		    dy = neara + hyfrac - max (ymin, neara - hyfrac)
+		    if (dy > 0.0) {
+			accum = accum + dy * Memr[xindex+(neara-nylmin)*nxpts]
+			waccum = waccum + dy
+		    }
+
+		    # Interior segments.
+		    do jj = neara + 1, nearb - 1 {
+			accum = accum + yfrac * Memr[xindex+(jj-nylmin)*nxpts]
+			waccum = waccum + yfrac
+		    }
+
+		    # Last segment.
+		    dy = min (ymax, nearb + hyfrac) - (nearb - hyfrac)
+		    if (dy > 0.0) {
+			accum = accum + dy * Memr[xindex+(nearb-nylmin)*nxpts]
+			waccum = waccum + dy
+		    }
+		}
+
+		if (waccum <= 0.0)
+		    zfit[i,j] = 0.0
+		else
+		    zfit[i,j] = accum / waccum
+	    }
+	}
+
+	call sfree (sp)
+end