1 files changed, 47 insertions, 0 deletions
diff --git a/noao/digiphot/daophot/daolib/bicubic.x b/noao/digiphot/daophot/daolib/bicubic.x
new file mode 100644
index 00000000..50bec983
--- /dev/null
+++ b/noao/digiphot/daophot/daolib/bicubic.x
@@ -0,0 +1,47 @@
+# BICUBIC -- Perform bicubic interpolation on an array. The point at which
+# the function is to be estimated lies between the second and third columns
+# and second and third rows of f, at a distance of (dx, dy), from the (2,2)
+# element of f.
+
+real procedure bicubic (f, nbox, dx, dy, dfdx, dfdy)
+
+real	f[nbox,nbox]		# input real array to be interpolated
+int	nbox			# size of the input array
+real	dx, dy			# point at which array is to be interpolated
+real	dfdx, dfdy		# output derivative of the interpolant
+
+int	jy
+real	c1, c2, c3, c4, temp[4], dfdxt[4], interp
+
+begin
+	# Interpolate first in x.
+	do jy = 1, 4 {
+	    c1 = 0.5 * (f[3,jy] - f[1,jy])
+	    c4 = f[3,jy] - f[2,jy] - c1
+	    c2 = 3.0 * c4 - 0.5 * (f[4,jy] - f[2,jy]) + c1
+	    c3 = c4 - c2
+	    c4 = dx * c3
+	    temp[jy] = dx * (dx * (c4 + c2) + c1) + f[2,jy]
+	    dfdxt[jy] = dx * (c4 * 3.0 + 2.0 * c2) + c1
+	}
+
+	# Interpolate next in y.
+	c1 = 0.5 * (temp[3] - temp[1])
+	c4 = temp[3] - temp[2] - c1
+	c2 = 3.0 * c4 - 0.5 * (temp[4] - temp[2]) + c1
+	c3 = c4 - c2
+	c4 = dy * c3
+
+	# Get the result.
+	interp = dy * (dy * (c4 + c2) + c1) + temp[2]
+
+	# Compute the derivatives.
+	dfdy = dy * (c4 * 3.0 + 2.0 * c2) + c1
+	c1 = 0.5 * (dfdxt[3] - dfdxt[1])
+	c4 = dfdxt[3] - dfdxt[2] - c1
+	c2 = 3.0 * c4 - 0.5 * (dfdxt[4] - dfdxt[2]) + c1
+	c3 = c4 - c2
+	dfdx = dy * (dy * (dy * c3 + c2) + c1) + dfdxt[2]
+
+	return (interp)
+end