Initial commit

1 files changed, 1080 insertions, 0 deletions

diff --git a/math/iminterp/ii_bieval.x b/math/iminterp/ii_bieval.x
new file mode 100644
index 00000000..3469128e
--- /dev/null
+++ b/math/iminterp/ii_bieval.x
@@ -0,0 +1,1080 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+include <math.h>
+
+# II_BINEAREST -- Procedure to evaluate the nearest neighbour interpolant.
+# The real array coeff contains the coefficients of the 2D interpolant.
+# The procedure assumes that 1 <= x <= nxpix and 1 <= y <= nypix and that
+# coeff[1+first_point] = datain[1,1].
+ 
+procedure ii_binearest (coeff, first_point, len_coeff, x, y, zfit, npts)
+
+real	coeff[ARB]	# 1D coefficient array
+int	first_point	# offset of first data point
+int	len_coeff	# row length of coeff
+real	x[npts]		# array of x values
+real	y[npts]		# array of y values
+real	zfit[npts]	# array of interpolated values
+int	npts		# number of points to be evaluated
+
+int	nx, ny
+int	index
+int	i
+
+begin
+	do i = 1, npts {
+
+	    nx = x[i] + 0.5
+	    ny = y[i] + 0.5
+
+	    # define pointer to data[nx,ny]
+	    index = first_point + (ny - 1) * len_coeff + nx
+
+	    zfit[i] = coeff[index]
+
+	}
+end
+
+
+# II_BILINEAR -- Procedure to evaluate the bilinear interpolant.
+# The real array coeff contains the coefficients of the 2D interpolant.
+# The procedure assumes that 1 <= x <= nxpix and 1 <= y <= nypix
+# and that coeff[1+first_point] = datain[1,1].
+
+procedure ii_bilinear (coeff, first_point, len_coeff, x, y, zfit, npts)
+
+real	coeff[ARB]	# 1D array of coefficients
+int	first_point	# offset of first data point
+int	len_coeff	# row length of coeff
+real	x[npts]		# array of x values
+real	y[npts]		# array of y values
+real	zfit[npts]	# array of interpolated values
+int	npts		# number of data points
+
+int	nx, ny
+int	index
+int	i
+real	sx, sy, tx, ty
+
+begin
+	do i = 1, npts {
+
+	    nx = x[i]
+	    ny = y[i]
+
+	    sx = x[i] - nx
+	    tx = 1. - sx
+	    sy = y[i] - ny
+	    ty = 1. - sy
+
+	    # define pointer to data[nx,ny]
+	    index = first_point + (ny - 1) * len_coeff + nx
+
+	    zfit[i] = tx * ty * coeff[index] + sx * ty * coeff[index + 1] +
+		      sy * tx * coeff[index+len_coeff] +
+		      sx * sy * coeff[index+len_coeff+1]
+	}
+end
+
+
+# II_BIPOLY3 -- Procedure to evaluate the bicubic polynomial interpolant.
+# The real array coeff contains the coefficients of the 2D interpolant.
+# The procedure assumes that 1 <= x <= nxpix and  1 <= y <= nypix
+# and that coeff[1+first_point] = datain[1,1]. The interpolant is
+# evaluated using Everett's central difference formula.
+
+procedure ii_bipoly3 (coeff, first_point, len_coeff, x, y, zfit, npts)
+
+real	coeff[ARB]	# 1D array of coefficients
+int	first_point	# offset first point
+int	len_coeff	# row length of the coefficient array
+real	x[npts]		# array of x values
+real	y[npts]		# array of y values
+real	zfit[npts]	# array of fitted values
+int	npts		# number of points to be evaluated
+
+int	nxold, nyold, nx, ny
+int	first_row, index
+int	i, j
+real	sx, tx, sx2m1, tx2m1, sy, ty
+real	cd20[4], cd21[4], ztemp[4], cd20y, cd21y
+
+begin
+	nxold = -1
+	nyold = -1
+
+	do i = 1, npts {
+
+	    nx = x[i]
+	    sx = x[i] - nx
+	    tx = 1. - sx
+	    sx2m1 = sx * sx - 1.
+	    tx2m1 = tx * tx - 1.
+
+	    ny = y[i]
+	    sy = y[i] - ny
+	    ty = 1. - sy
+
+	    # define pointer to datain[nx,ny-1]
+	    first_row = first_point + (ny - 2) * len_coeff + nx
+
+	    # loop over the 4 surrounding rows of data
+	    # calculate the  central differences at each value of y
+
+	    # if new data point caculate the central differnences in x
+	    # for each y
+
+	    index = first_row
+	    if (nx != nxold || ny != nyold) {
+	        do j = 1, 4 {
+		   cd20[j] = 1./6. * (coeff[index+1] - 2. * coeff[index] +
+		   	     coeff[index-1])
+		   cd21[j] = 1./6. * (coeff[index+2] - 2. * coeff[index+1] +
+			     coeff[index])
+		    index = index + len_coeff
+	        }
+	    }
+
+	    # interpolate in x at each value of y
+	    index = first_row
+	    do j = 1, 4 {
+		ztemp[j] = sx * (coeff[index+1] + sx2m1 * cd21[j]) +
+			   tx * (coeff[index] + tx2m1 * cd20[j])
+		index = index + len_coeff
+	    }
+
+	    # calculate y central differences
+	    cd20y = 1./6. * (ztemp[3] - 2. * ztemp[2] + ztemp[1])
+	    cd21y = 1./6. * (ztemp[4] - 2. * ztemp[3] + ztemp[2])
+
+	    # interpolate in y
+	    zfit[i] = sy * (ztemp[3] + (sy * sy - 1.) * cd21y) +
+		      ty * (ztemp[2] + (ty * ty - 1.) * cd20y)
+
+	    nxold = nx
+	    nyold = ny
+
+	}
+end
+
+
+# II_BIPOLY5 -- Procedure to evaluate a biquintic polynomial.
+# The real array coeff contains the coefficents of the 2D interpolant.
+# The routine assumes that 1 <= x <= nxpix and 1 <= y <= nypix
+# and that coeff[1+first_point] = datain[1,1]. The interpolant is evaluated
+# using Everett's central difference formula.
+
+procedure ii_bipoly5 (coeff, first_point, len_coeff, x, y, zfit, npts)
+
+real	coeff[ARB]	# 1D array of coefficients
+int	first_point	# offset to first data point
+int	len_coeff	# row length of coeff
+real	x[npts]		# array of x values
+real	y[npts]		# array of y values
+real	zfit[npts]	# array of fitted values
+int	npts		# number of points
+
+int	nxold, nyold, nx, ny
+int	first_row, index
+int	i, j
+real	sx, sx2, sx2m1, sx2m4, tx, tx2, tx2m1, tx2m4, sy, sy2, ty, ty2
+real	cd20[6], cd21[6], cd40[6], cd41[6], ztemp[6]
+real	cd20y, cd21y, cd40y, cd41y
+
+begin
+	nxold = -1
+	nyold = -1
+
+	do i = 1, npts {
+
+	    nx = x[i]
+	    sx = x[i] - nx
+	    sx2 = sx * sx
+	    sx2m1 = sx2 - 1.
+	    sx2m4 = sx2 - 4.
+	    tx = 1. - sx
+	    tx2 = tx * tx
+	    tx2m1 = tx2 - 1.
+	    tx2m4 = tx2 - 4.
+
+	    ny = y[i]
+	    sy = y[i] - ny
+	    sy2 = sy * sy
+	    ty = 1. - sy
+	    ty2 = ty * ty
+
+	    # calculate value of pointer to data[nx,ny-2]
+	    first_row = first_point + (ny - 3) * len_coeff + nx
+
+	    # calculate the central differences in x at each value of y
+	    index = first_row
+	    if (nx != nxold || ny != nyold) {
+		do j = 1, 6 {
+		    cd20[j] = 1./6. * (coeff[index+1] - 2. * coeff[index] +
+		    	      coeff[index-1])
+		    cd21[j] = 1./6. * (coeff[index+2] - 2. * coeff[index+1] +
+			      coeff[index])
+		    cd40[j] = 1./120. * (coeff[index-2] - 4. * coeff[index-1] +
+			      6. * coeff[index] - 4. * coeff[index+1] +
+			      coeff[index+2])
+		    cd41[j] = 1./120. * (coeff[index-1] - 4. * coeff[index] +
+			      6. * coeff[index+1] - 4. * coeff[index+2] +
+			      coeff[index+3])
+		    index = index + len_coeff
+		}
+	    }
+
+	    # interpolate in x at each value of y
+	    index = first_row
+	    do j = 1, 6 {
+		ztemp[j] = sx * (coeff[index+1] + sx2m1 * (cd21[j] + sx2m4 *
+			   cd41[j])) + tx * (coeff[index] + tx2m1 *
+			   (cd20[j] + tx2m4 * cd40[j]))
+		index = index + len_coeff
+	    }
+
+	    # central differences in y
+	    cd20y = 1./6. * (ztemp[4] - 2. * ztemp[3] + ztemp[2])
+	    cd21y = 1./6. * (ztemp[5] - 2. * ztemp[4] + ztemp[3])
+	    cd40y = 1./120. * (ztemp[1] - 4. * ztemp[2] + 6. * ztemp[3] -
+		    4. * ztemp[4] + ztemp[5])
+	    cd41y = 1./120. * (ztemp[2] - 4. * ztemp[3] + 6. * ztemp[4] -
+		    4. * ztemp[5] + ztemp[6])
+
+	    # interpolate in y
+	    zfit[i] = sy * (ztemp[4] + (sy2 - 1.) * (cd21y + (sy2 - 4.) *
+		      cd41y)) + ty * (ztemp[3] + (ty2 - 1.) * (cd20y +
+		      (ty2 - 4.) * cd40y))
+
+	    nxold = nx
+	    nyold = ny
+
+	}
+end
+
+
+# II_BISPLINE3 -- Procedure to evaluate a bicubic spline.
+# The real array coeff contains the B-spline coefficients.
+# The procedure assumes that 1 <= x <= nxpix and 1 <= y <= nypix
+# and that coeff[1+first_point] = B-spline[2].
+
+procedure ii_bispline3 (coeff, first_point, len_coeff, x, y, zfit, npts)
+
+real	coeff[ARB]	# 1D array of coefficients
+int	first_point	# offset to first data point
+int	len_coeff	# row length of coeff
+real	x[npts]		# array of x values
+real	y[npts]		# array of y values
+real	zfit[npts]	# array of interpolated values
+int	npts		# number of points to be evaluated
+
+int	nx, ny
+int	first_row, index
+int	i, j
+real	sx, tx, sy, ty
+real	bx[4], by[4], accum, sum
+
+begin
+	do i = 1, npts {
+
+	    nx = x[i]
+	    sx = x[i] - nx
+	    tx = 1. - sx
+
+	    ny = y[i]
+	    sy = y[i] - ny
+	    ty = 1. - sy
+
+
+	    # calculate the x B-splines
+	    bx[1] = tx ** 3 
+	    bx[2] = 1. + tx * (3. + tx * (3. - 3. * tx))
+	    bx[3] = 1. + sx * (3. + sx * (3. - 3. * sx))
+	    bx[4] = sx ** 3
+
+	    # calculate the y B-splines
+	    by[1] = ty ** 3
+	    by[2] = 1. + ty * (3. + ty * (3. - 3. * ty))
+	    by[3] = 1. + sy * (3. + sy * (3. - 3. * sy))
+	    by[4] = sy ** 3
+
+	    # calculate the pointer to data[nx,ny-1]
+	    first_row = first_point + (ny - 2) * len_coeff + nx
+
+	    # evaluate spline
+	    accum = 0.
+	    index = first_row
+	    do j = 1, 4 {
+		sum = coeff[index-1] * bx[1] + coeff[index] * bx[2] +
+		      coeff[index+1] * bx[3] + coeff[index+2] * bx[4]	
+		accum = accum + sum * by[j]
+		index = index + len_coeff
+	    }
+
+	    zfit[i] = accum
+	}
+end
+
+
+# II_BISINC -- Procedure to evaluate the 2D sinc function.  The real array
+# coeff contains the data. The procedure assumes that 1 <= x <= nxpix and
+# 1 <= y <= nypix and that coeff[1+first_point] = datain[1,1]. The since
+# truncation length is nsinc. The taper is a cosbell function 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 is returned. This
+# routine does not use precomputed arrays.
+
+procedure ii_bisinc (coeff, first_point, len_coeff, len_array, x, y, zfit,
+	npts, nsinc, mindx, mindy)
+
+real	coeff[ARB]	# 1D array of coefficients
+int	first_point	# offset to first data point
+int	len_coeff	# row length of coeff
+int	len_array	# column length of coeff
+real	x[npts]		# array of x values
+real	y[npts]		# array of y values
+real	zfit[npts]	# array of interpolated values
+int	npts		# the number of input points.
+int	nsinc		# sinc truncation length
+real	mindx		# interpolation mininmum in x
+real	mindy		# interpolation mininmum in y
+
+int	i, j, k, nconv, nx, ny, index, mink, maxk, offk, minj, maxj, offj
+int	last_point
+pointer	sp, taper, ac, ar
+real	sconst, a2, a4, sdx, dx, dy, dxn, dyn, ax, ay, px, py, sumx, sumy, sum
+real	dx2
+
+begin
+	# Compute the length of the convolution.
+	nconv = 2 * nsinc + 1
+
+	# Allocate working array space.
+	call smark (sp)
+	call salloc (taper, nconv, TY_REAL)
+	call salloc (ac, nconv, TY_REAL)
+	call salloc (ar, nconv, TY_REAL)
+
+	# Compute the constants for the cosine bell taper.
+        sconst = (HALFPI / nsinc) ** 2
+        a2 = -0.49670
+        a4 = 0.03705
+
+	# Precompute the taper array. Incorporate the sign change portion
+	# of the sinc interpolator into the taper array.
+	if (mod (nsinc, 2) == 0)
+	    sdx = 1.0
+	else
+	    sdx = -1.0
+        do j = -nsinc,  nsinc {
+	    dx2 = sconst * j * j
+            Memr[taper+j+nsinc] = sdx * (1.0 + a2 * dx2 + a4 * dx2 * dx2) ** 2
+            sdx = -sdx
+        }
+
+	do i = 1, npts {
+
+	    # define the fractional pixel interpolation.
+	    nx = nint (x[i])
+	    ny = nint (y[i])
+	    if (nx < 1 || nx > len_coeff || ny < 1 || ny > len_array) {
+		zfit[i] = 0.0
+		next
+	    }
+	    dx = x[i] - nx
+	    dy = y[i] - ny
+
+	    # define pointer to data[nx,ny]
+	    if (abs (dx) < mindx && abs (dy) < mindy) {
+	        index = first_point + (ny - 1) * len_coeff + nx
+		zfit[i] = coeff[index]
+		next
+	    }
+
+	    # initialize.
+	    #dxn = -1-nsinc-dx
+	    #dyn = -1-nsinc-dy
+	    dxn = 1 + nsinc + dx
+	    dyn = 1 + nsinc + dy
+
+	    # Compute the x and y sinc arrays using a cosbell taper.
+	    sumx = 0.0
+	    sumy = 0.0
+	    do j = 1, nconv {
+
+		#ax = j + dxn
+		#ay = j + dyn
+		ax = dxn - j
+		ay = dyn - j
+		if (ax == 0.0)
+		    px = 1.0
+		else if (dx == 0.0)
+		    px = 0.0
+		else
+		    px = Memr[taper+j-1] / ax 
+		if (ay == 0.0)
+		    py = 1.0
+		else if (dy == 0.0)
+		    py = 0.0
+		else
+		    py = Memr[taper+j-1] / ay 
+
+		Memr[ac+j-1] = px
+		Memr[ar+j-1] = py
+		sumx = sumx + px
+		sumy = sumy + py
+	    }
+
+	    # Compute the limits of the convolution.
+	    minj = max (1, ny - nsinc)
+	    maxj = min (len_array, ny + nsinc)
+	    offj = ar - ny + nsinc
+	    mink = max (1, nx - nsinc)
+	    maxk = min (len_coeff, nx + nsinc)
+	    offk = ac - nx + nsinc
+
+	    # Initialize
+	    zfit[i] = 0.0
+
+	    # Do the convolution.
+	    do j = ny - nsinc, minj - 1 {
+		sum = 0.0
+		do k = nx - nsinc, mink - 1
+		    sum = sum + Memr[k+offk] * coeff[first_point+1] 
+	        do k = mink, maxk
+		    sum = sum + Memr[k+offk] * coeff[first_point+k] 
+		do k = maxk + 1, nx + nsinc
+		    sum = sum + Memr[k+offk] * coeff[first_point+len_coeff] 
+
+		zfit[i] = zfit[i] + Memr[j+offj] * sum
+	    }
+
+	    do j = minj, maxj {
+	    	index = first_point + (j - 1) * len_coeff
+		sum = 0.0
+		do k = nx - nsinc, mink - 1
+		    sum = sum + Memr[k+offk] * coeff[index+1] 
+	        do k = mink, maxk
+		    sum = sum + Memr[k+offk] * coeff[index+k] 
+		do k = maxk + 1, nx + nsinc
+		    sum = sum + Memr[k+offk] * coeff[index+len_coeff] 
+
+		zfit[i] = zfit[i] + Memr[j+offj] * sum
+	    }
+
+	    do j = maxj + 1, ny + nsinc {
+	        last_point = first_point + (len_array - 1) * len_coeff
+		sum = 0.0
+		do k = nx - nsinc, mink - 1
+		    sum = sum + Memr[k+offk] * coeff[last_point+1] 
+	        do k = mink, maxk
+		    sum = sum + Memr[k+offk] * coeff[last_point+k] 
+		do k = maxk + 1, nx + nsinc
+		    sum = sum + Memr[k+offk] * coeff[last_point+len_coeff] 
+
+		zfit[i] = zfit[i] + Memr[j+offj] * sum
+	    }
+
+	    # Normalize.
+	    zfit[i] = zfit[i] / sumx / sumy
+	}
+
+	call sfree (sp)
+end
+
+
+# II_BILSINC -- Procedure to evaluate the 2D sinc function.  The real array
+# coeff contains the data. The procedure assumes that 1 <= x <= nxpix and
+# 1 <= y <= nypix and that coeff[1+first_point] = datain[1,1]. The since
+# truncation length is nsinc. The taper is a cosbell function 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 is returned. This
+# routine does use precomputed arrays.
+
+procedure ii_bilsinc (coeff, first_point, len_coeff, len_array, x, y, zfit,
+	npts, ltable, nconv, nxincr, nyincr, mindx, mindy)
+
+real	coeff[ARB]			   # 1D array of coefficients
+int	first_point			   # offset to first data point
+int	len_coeff			   # row length of coeff
+int	len_array			   # column length of coeff
+real	x[npts]				   # array of x values
+real	y[npts]				   # array of y values
+real	zfit[npts]			   # array of interpolated values
+int	npts				   # the number of input points.
+real	ltable[nconv,nconv,nxincr,nyincr]  # the pre-computed look-up table
+int	nconv				   # the sinc truncation full width
+int	nxincr				   # the interpolation resolution in x
+int	nyincr				   # the interpolation resolution in y
+real	mindx				   # interpolation mininmum in x
+real	mindy				   # interpolation mininmum in y
+
+int	i, j, k, nsinc, xc, yc, lutx, luty, minj, maxj, offj, mink, maxk, offk
+int	index, last_point
+real	dx, dy, sum
+
+begin
+	nsinc = (nconv - 1) / 2
+	do i = 1, npts {
+
+	    # Return zero outside of data.
+	    xc = nint (x[i])
+	    yc = nint (y[i])
+	    if (xc < 1 || xc > len_coeff || yc < 1 || yc > len_array) {
+		zfit[i] = 0.0
+		next
+	    }
+
+	    dx = x[i] - xc
+	    dy = y[i] - yc
+	    if (abs(dx) < mindx && abs(dy) < mindy) {
+		index = first_point + (yc - 1) * len_coeff + xc
+		zfit[i] = coeff[index]
+	    }
+
+	    # Find the correct look-up table entry.
+	    if (nxincr == 1)
+		lutx = 1
+	    else 
+		lutx = nint ((-dx + 0.5) * (nxincr - 1)) + 1
+		#lutx = int ((-dx + 0.5) * (nxincr - 1) + 0.5) + 1
+	    if (nyincr == 1)
+		luty = 1
+	    else
+		luty = nint ((-dy + 0.5) * (nyincr - 1)) + 1
+		#luty = int ((-dy + 0.5) * (nyincr - 1) + 0.5) + 1
+
+	    # Compute the convolution limits.
+	    minj = max (1, yc - nsinc)
+	    maxj = min (len_array, yc + nsinc)
+	    offj = 1 - yc + nsinc
+	    mink = max (1, xc - nsinc)
+	    maxk = min (len_coeff, xc + nsinc)
+	    offk = 1 - xc + nsinc
+
+	    # Initialize
+	    zfit[i] = 0.0
+
+	    # Do the convolution.
+	    do j = yc - nsinc, minj - 1 {
+		sum = 0.0
+		do k = xc - nsinc, mink - 1
+		    sum = sum + ltable[k+offk,j+offj,lutx,luty] *
+		        coeff[first_point+1] 
+	        do k = mink, maxk
+		    sum = sum + ltable[k+offk,j+offj,lutx,luty] *
+		        coeff[first_point+k] 
+		do k = maxk + 1, xc + nsinc
+		    sum = sum + ltable[k+offk,j+offj,lutx,luty] *
+		        coeff[first_point+len_coeff] 
+		zfit[i] = zfit[i] +  sum
+	    }
+
+	    do j = minj, maxj {
+	    	index = first_point + (j - 1) * len_coeff
+		sum = 0.0
+		do k = xc - nsinc, mink - 1
+		    sum = sum + ltable[k+offk,j+offj,lutx,luty] * coeff[index+1]
+	        do k = mink, maxk
+		    sum = sum + ltable[k+offk,j+offj,lutx,luty] * coeff[index+k]
+		do k = maxk + 1, xc + nsinc
+		    sum = sum + ltable[k+offk,j+offj,lutx,luty] *
+		        coeff[index+len_coeff] 
+		zfit[i] = zfit[i] + sum
+	    }
+
+	    do j = maxj + 1, yc + nsinc {
+	        last_point = first_point + (len_array - 1) * len_coeff
+		sum = 0.0
+		do k = xc - nsinc, mink - 1
+		    sum = sum + ltable[k+offk,j+offj,lutx,luty] *
+		        coeff[last_point+1] 
+	        do k = mink, maxk
+		    sum = sum + ltable[k+offk,j+offj,lutx,luty] *
+		        coeff[last_point+k] 
+		do k = maxk + 1, xc + nsinc
+		    sum = sum + ltable[k+offk,j+offj,lutx,luty] *
+		        coeff[last_point+len_coeff] 
+		zfit[i] = zfit[i] + sum
+	    }
+
+	}
+end
+
+
+# II_BIDRIZ -- Procedure to evaluate the drizzle interpolant.
+# The real array coeff contains the coefficients of the 2D interpolant.
+# The procedure assumes that 1 <= x <= nxpix and 1 <= y <= nypix and that
+# coeff[1+first_point] = datain[1,1]. Each x and y value is a set of 4
+# values describing the vertices of a quadrilateral in the input data. The
+# integration routine was adapted from the one developed by Bill Sparks at
+# ST and used the DITHER / DRIZZLE software. The 4 points describing the
+# corners of each quadrilateral integration region must be in order, e.g.
+# describe the vertices of a polygon in either CW or CCW order.
+ 
+procedure ii_bidriz (coeff, first_point, len_coeff, x, y, zfit, npts,
+	xfrac, yfrac, badval)
+
+real	coeff[ARB]	# 1D coefficient array
+int	first_point	# offset of first data point
+int	len_coeff	# row length of coeff
+real	x[ARB]		# array of x values
+real	y[ARB]		# array of y values
+real	zfit[npts]	# array of interpolated values
+int	npts		# number of points to be evaluated
+real	xfrac, yfrac	# the x and y drizzle pixel fractions
+real	badval		# undefined pixel value
+
+int	i, ii, jj, kk, index, nearax, nearbx, nearay, nearby
+real	px[5], py[5], dx, xmin, xmax, m, c, ymin, ymax, xtop
+real	ovlap, accum, waccum, dxfrac, dyfrac, hxfrac, hyfrac, dhxfrac, dhyfrac
+bool	negdx
+
+begin
+	dxfrac = max (0.0, min (1.0, 1.0 - xfrac))
+	hxfrac = max (0.0, min (0.5, dxfrac / 2.0))
+	dhxfrac = max (0.5, min (1.0, 1.0 - hxfrac))
+	dyfrac = max (0.0, min (1.0, 1.0 - yfrac))
+	hyfrac = max (0.0, min (0.5, dyfrac / 2.0))
+	dhyfrac = max (0.5, min (1.0, 1.0 - hyfrac))
+
+	do i = 1, npts {
+
+	    # Compute the limits of the integration in x and y.
+	    nearax = nint (min (x[4*i-3], x[4*i-2], x[4*i-1], x[4*i]))
+	    nearbx = nint (max (x[4*i-3], x[4*i-2], x[4*i-1], x[4*i]))
+	    nearay = nint (min (y[4*i-3], y[4*i-2], y[4*i-1], y[4*i]))
+	    nearby = nint (max (y[4*i-3], y[4*i-2], y[4*i-1], y[4*i]))
+
+	    # Initialize.
+	    accum = 0.0
+	    waccum = 0.0
+
+	    # Loop over all pixels which contribute to the integral.
+	    do jj = nearay, nearby {
+		index = first_point + (jj - 1) * len_coeff
+		do kk = 1, 4
+		    py[kk] = y[4*i+kk-4] - jj + 0.5
+		py[5] = py[1]
+		do ii = nearax, nearbx {
+
+		    # Transform the coordinates relative to a unit
+		    # square centered at the origin of the pixel. We
+		    # are going to approximate the new pixel area by
+		    # a quadilateral. Close the quadilateral.
+
+		    do kk = 1, 4
+			px[kk] = x[4*i+kk-4] - ii + 0.5 
+		    px[5] = px[1]
+
+		    # Compute the area overlap of the output pixel with
+		    # the input pixels. 
+		    ovlap = 0.0
+		    do kk = 1, 4 {
+
+			# Check for vertical line segment.
+			dx = px[kk+1] - px[kk]
+			if (dx == 0.0)
+			    next
+
+			# Order the x integration limits.
+			if (px[kk] < px[kk+1]) {
+			    xmin = px[kk]
+			    xmax = px[kk+1]
+			} else {
+			    xmin = px[kk+1]
+			    xmax = px[kk]
+			}
+
+			# Check the x limits ignoring y for now.
+			if ((xmin >= dhxfrac) || (xmax <= hxfrac))
+			    next
+			xmin = max (xmin, hxfrac)
+			xmax = min (xmax, dhxfrac)
+
+			# Get basic info about the line y = mx + c.
+			if (dx < 0.0)
+			    negdx = true
+			else
+			    negdx = false
+			m = (py[kk+1] - py[kk]) / dx
+			c = py[kk] - m * px[kk]
+			ymin = m * xmin + c
+			ymax = m * xmax + c
+
+			# Trap segment entirely below axis.
+			if (ymin <= hyfrac && ymax <= hyfrac)
+			    next
+
+			# Adjust bounds if segment crosses axis in order
+			# to exclude anything below the axis.
+			if (ymin < hyfrac) {
+			    ymin = hyfrac
+			    xmin = (hyfrac - c) / m
+			}
+			if (ymax < hyfrac) {
+			    ymax = hyfrac
+			    xmax = (hyfrac - c) / m
+			}
+
+			# There are four possibilities.
+
+			# Both y above 1.0 - hyfrac. Line segment is entirely
+			# above square.
+			if ((ymin >= dhyfrac) && (ymax >= dhyfrac)) {
+
+			    if (negdx)
+				ovlap = ovlap + (xmin - xmax) * yfrac
+			    else
+				ovlap = ovlap + (xmax - xmin) * yfrac
+
+			# Both y below 1.0 - hyfrac. Segment is entirely
+			# within square.
+			} else if ((ymin <= dhyfrac) && (ymax <= dhyfrac)) {
+
+			    if (negdx)
+				ovlap = ovlap + 0.5 * (xmin - xmax) *
+				    (ymax + ymin - dyfrac)
+			    else
+				ovlap = ovlap + 0.5 * (xmax - xmin) *
+				    (ymax + ymin - dyfrac)
+
+			# One of each. Segment must cross top of square.
+			} else {
+
+			    xtop = (dhyfrac - c) / m
+
+			    # insert precision check ?
+
+			    if (ymin < dhyfrac) {
+				if (negdx)
+				    ovlap = ovlap - (0.5 * (xtop - xmin) *
+				        (ymin + 1.0 - 3.0 * hyfrac) +
+				        (xmax - xtop) * yfrac)
+				else
+				    ovlap = ovlap + (0.5 * (xtop - xmin) *
+				        (ymin + 1.0 - 3.0 * hyfrac) +
+				        (xmax - xtop) * yfrac)
+			     } else {
+				if (negdx)
+				    ovlap = ovlap - (0.5 * (xmax - xtop) *
+				        (ymax + 1.0 - 3.0 * hyfrac) +
+					(xtop - xmin) * yfrac)
+				else
+				    ovlap = ovlap + (0.5 * (xmax - xtop) *
+				        (ymax + 1.0 - 3.0 * hyfrac) +
+					(xtop - xmin) * yfrac)
+			     }
+
+			}
+		    }
+
+		    accum = accum + coeff[index+ii] * ovlap
+		    waccum = waccum + ovlap
+		}
+	    }
+
+	    if (waccum == 0.0)
+		zfit[i] = badval
+	    else
+		zfit[i] = accum / waccum
+	}
+end
+
+
+# II_BIDRIZ1 -- Procedure to evaluate the drizzle interpolant when xfrac and
+# yfrac are 1.0. The real array coeff contains the coefficients of the 2D
+# interpolant. The procedure assumes that 1 <= x <= nxpix and 1 <= y <= nypix
+# and that coeff[1+first_point] = datain[1,1]. Each x and y point is a set of 4
+# values describing the vertices of a quadrilateral in the input data. The
+# integration routine was adapted from the one developed by Bill Sparks at
+# ST and used the DITHER / DRIZZLE software. The 4 points describing the
+# corners of each quadrilateral integration region must be in order, e.g.
+# describe the vertices of a polygon in either CW or CCW order.
+ 
+procedure ii_bidriz1 (coeff, first_point, len_coeff, x, y, zfit, npts, badval)
+
+real	coeff[ARB]	# 1D coefficient array
+int	first_point	# offset of first data point
+int	len_coeff	# row length of coeff
+real	x[ARB]		# array of x values
+real	y[ARB]		# array of y values
+real	zfit[npts]	# array of interpolated values
+int	npts		# number of points to be evaluated
+real	badval		# undefined pixel value
+
+int	i, ii, jj, kk, index, nearax, nearbx, nearay, nearby
+real	px[5], py[5], dx, xmin, xmax, m, c, ymin, ymax, xtop
+real	ovlap, accum, waccum
+bool	negdx
+
+begin
+	do i = 1, npts {
+
+	    # Compute the limits of the integration in x and y.
+	    nearax = nint (min (x[4*i-3], x[4*i-2], x[4*i-1], x[4*i]))
+	    nearbx = nint (max (x[4*i-3], x[4*i-2], x[4*i-1], x[4*i]))
+	    nearay = nint (min (y[4*i-3], y[4*i-2], y[4*i-1], y[4*i]))
+	    nearby = nint (max (y[4*i-3], y[4*i-2], y[4*i-1], y[4*i]))
+
+	    # Initialize.
+	    accum = 0.0
+	    waccum = 0.0
+
+	    # Loop over all pixels which contribute to the integral.
+	    do jj = nearay, nearby {
+		index = first_point + (jj - 1) * len_coeff
+		do kk = 1, 4
+		    py[kk] = y[4*i+kk-4] - jj + 0.5
+		py[5] = py[1]
+		do ii = nearax, nearbx {
+
+		    # Transform the coordinates relative to a unit
+		    # square centered at the origin of the pixel. We
+		    # are going to approximate the new pixel area by
+		    # a quadilateral. Close the polygon.
+
+		    do kk = 1, 4
+			px[kk] = x[4*i+kk-4] - ii + 0.5
+		    px[5] = px[1]
+
+		    # Compute the area overlap of the output pixel with
+		    # the input pixels. 
+		    ovlap = 0.0
+		    do kk = 1, 4 {
+
+			# Check for vertical line segment.
+			dx = px[kk+1] - px[kk]
+			if (dx == 0.0)
+			    next
+
+			# Order the x integration limits.
+			if (px[kk] < px[kk+1]) {
+			    xmin = px[kk]
+			    xmax = px[kk+1]
+			} else {
+			    xmin = px[kk+1]
+			    xmax = px[kk]
+			}
+
+			# Check the x limits ignoring y for now.
+			if (xmin >= 1.0 || xmax <= 0.0)
+			    next
+			xmin = max (xmin, 0.0)
+			xmax = min (xmax, 1.0)
+
+			# Get basic info about the line y = mx + c.
+			if (dx < 0.0)
+			    negdx = true
+			else
+			    negdx = false
+			m = (py[kk+1] - py[kk]) / dx
+			c = py[kk] - m * px[kk]
+			ymin = m * xmin + c
+			ymax = m * xmax + c
+
+			# Trap segment entirely below axis.
+			if (ymin <= 0.0 && ymax <= 0.0)
+			    next
+
+			# Adjust bounds if segment crosses axis in order
+			# to exclude anything below the axis.
+			if (ymin < 0.0) {
+			    ymin = 0.0
+			    xmin = - c / m
+			}
+			if (ymax < 0.0) {
+			    ymax = 0.0
+			    xmax = - c / m
+			}
+
+			# There are four possibilities.
+
+			# Both y above 1.0. Line segment is entirely above
+			# square.
+			if (ymin >= 1.0 && ymax >= 1.0) {
+
+			    if (negdx)
+				ovlap = ovlap + (xmin - xmax)
+			    else
+				ovlap = ovlap + (xmax - xmin)
+
+			# Both y below 1.0. Segment is entirely within square.
+			} else if (ymin <= 1.0 && ymax <= 1.0) {
+
+			    if (negdx)
+				ovlap = ovlap + 0.5 * (xmin - xmax) *
+				    (ymax + ymin)
+			    else
+				ovlap = ovlap + 0.5 * (xmax - xmin) *
+				    (ymax + ymin)
+
+			# One of each. Segment must cross top of square.
+			} else {
+
+			    xtop = (1.0 - c) / m
+
+			    # insert precision check, e.g. possible pixel
+			    # overlap ? need to decide what action to take ...
+
+			    if (ymin < 1.0) {
+				if (negdx)
+				    ovlap = ovlap - (0.5 * (xtop - xmin) *
+				        (1.0 + ymin) + (xmax - xtop))
+				else
+				    ovlap = ovlap + (0.5 * (xtop - xmin) *
+				        (1.0 + ymin) + (xmax - xtop))
+			     } else {
+				if (negdx)
+				    ovlap = ovlap - (0.5 * (xmax - xtop) *
+				        (1.0 + ymax) + (xtop - xmin))
+				else
+				    ovlap = ovlap + (0.5 * (xmax - xtop) *
+				        (1.0 + ymax) + (xtop - xmin))
+			     }
+
+			}
+		    }
+
+		    accum = accum + coeff[index+ii] * ovlap
+		    waccum = waccum + ovlap
+		}
+	    }
+
+	    if (waccum == 0.0)
+		zfit[i] = badval
+	    else
+		zfit[i] = accum / waccum
+	}
+end
+
+
+# II_BIDRIZ0-- Procedure to evaluate the drizzle interpolant when xfrac and
+# yfrac are 0.0. The real array coeff contains the coefficients of the 2D
+# interpolant. The procedure assumes that 1 <= x <= nxpix and 1 <= y <= nypix
+# and that coeff[1+first_point] = datain[1,1]. Each x and y point is a set of 4
+# values describing the vertices of a quadrilateral in the input data. The
+# integration routine determines whether a pixel is in, out, on the edge
+# of or at a vertex of a polygon. The 4 points describing the corners of
+# each quadrilateral integration region must be in order, e.g.  describe
+# the vertices of a polygon in either CW or CCW order.
+# THIS ROUTINE IS NOT CURRENTLY BEING USED.
+ 
+procedure ii_bidriz0 (coeff, first_point, len_coeff, x, y, zfit, npts, badval)
+
+real	coeff[ARB]	# 1D coefficient array
+int	first_point	# offset of first data point
+int	len_coeff	# row length of coeff
+real	x[ARB]		# array of x values
+real	y[ARB]		# array of y values
+real	zfit[npts]	# array of interpolated values
+int	npts		# number of points to be evaluated
+real	badval		# the undefined pixel value
+
+bool	boundary, vertex
+int	i, jj, ii, kk, nearax, nearbx, nearay, nearby, ninter
+real	accum, waccum, px[5], py[5], lx, ld, u1, u2, u1u2, dx, dy, dd
+real	xi, ovlap, xmin, xmax
+
+begin
+	do i = 1, npts {
+
+	    # Compute the limits of the integration in x and y.
+	    nearax = nint (min (x[4*i-3], x[4*i-2], x[4*i-1], x[4*i]))
+	    nearbx = nint (max (x[4*i-3], x[4*i-2], x[4*i-1], x[4*i]))
+	    nearay = nint (min (y[4*i-3], y[4*i-2], y[4*i-1], y[4*i]))
+	    nearby = nint (max (y[4*i-3], y[4*i-2], y[4*i-1], y[4*i]))
+
+	    # Initialize.
+	    accum = 0.0
+	    waccum = 0.0
+
+	    # Loop over all pixels which contribute to the integral.
+	    do jj = nearay, nearby {
+		do ii = nearax, nearbx {
+
+		    # Transform the coordinates relative to a unit
+		    # square centered at the origin of the pixel. We
+		    # are going to approximate the new pixel area by
+		    # a quadilateral. Close the quadrilateral.
+
+		    do kk = 1, 4 {
+			px[kk] = x[4*i+kk-4] - ii + 0.5
+			py[kk] = y[4*i+kk-4] - jj + 0.5
+		    }
+		    px[5] = px[1]
+		    py[5] = py[1]
+
+		    # Initialize the integration.
+		    ovlap = 0.0
+		    ninter = 0
+
+		    # Define a line segment which begins at the point x = 0.5
+		    # y = 0.5 and runs parallel to the y axis.
+		    call alimr (px, 5, xmin, xmax)
+		    lx = xmax - xmin
+		    ld = 0.5 * lx
+		    u1 = -lx * py[1] + ld
+		    boundary = false
+		    vertex = false
+		    do kk = 2, 5 {
+
+		        u2 = -lx * py[kk] + ld
+			u1u2 = u1 * u2
+
+			# No intersection.
+			if (u1*u2 > 0.0) {
+			    ;
+
+			# Intersection with polygon line segment.
+			} else if (u1 != 0.0 && u2 != 0.0) {
+			    dy = py[kk-1] - py[kk]
+			    dx = px[kk-1] - px[kk]
+			    dd = px[kk-1] * py[kk] - py[kk-1] * px[kk]
+			    xi = (dx * ld - lx * dd) / (dy * lx)
+			    if (xi > 0.5)
+			        ninter = ninter + 1
+			    if (xi == 0.5)
+				boundary = true
+
+			# Collinearity.
+			} else if (u1 == 0.0 && u2 == 0.0) {
+			    xmin = min (px[kk-1], px[kk])
+			    xmax = max (px[kk-1], px[kk])
+			    if (xmin == 0.5 || xmax == 0.5)
+				vertex = true
+			    else if (xmin < 0.5 && xmax > 0.5)
+				boundary = true
+
+			# Vertex.
+			} else if (u1 != 0.0) {
+			    if (px[kk] == 0.5)
+				vertex = true
+			}
+
+			u1 = u2
+		    }
+
+		    if (vertex)
+			ovlap = 0.25
+		    else if (boundary)
+			ovlap = 0.5
+		    else if (mod (ninter, 2) == 0)
+			ovlap = 0.0
+		    else
+		        ovlap = 1.0
+		    waccum = waccum + ovlap
+		    accum = accum + ovlap *
+		        coeff[first_point+(jj-1)*len_coeff+ii]
+		}
+	    }
+
+	    if (waccum == 0.0)
+		zfit[i] = badval
+	    else
+		zfit[i] = accum / waccum
+	}
+
+end