1 files changed, 41 insertions, 0 deletions
diff --git a/sys/vops/ak/advzd.x b/sys/vops/ak/advzd.x
new file mode 100644
index 00000000..ca5bb0da
--- /dev/null
+++ b/sys/vops/ak/advzd.x
@@ -0,0 +1,41 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+# ADVZ -- Vector divide with checking for zero divisors.  If the result of a
+# divide would be undefined a user supplied function is called to get the
+# output pixel value.
+#
+# NOTE: in the interests of simplicity a somewhat arbitrary tolerance is used
+# to check for an undefined divide, i.e., a divide by zero or a divide by a
+# number small enough to cause floating point overflow.  A better way to do
+# this would be to provide a machine dependent version of this operator in
+# host$as which catches the hardware exception rather than using a comparison.
+
+procedure advzd (a, b, c, npix, errfcn)
+
+double	a[ARB], b[ARB], c[ARB]	# numerator, divisor, and output arrays
+int	npix			# number of pixels
+double	errfcn()		# user function, called on divide by zero
+
+int	i
+double	divisor
+double	tol
+extern	errfcn()
+errchk	errfcn
+
+begin
+	    tol = 1.0D-20
+
+	do i = 1, npix {
+	    divisor = b[i]
+		# The following is most efficient when the data tends to be
+		# positive.
+
+		if (divisor < tol)
+		    if (divisor > -tol) {
+			c[i] = errfcn (a[i])
+			next
+		    }
+		c[i] = a[i] / divisor
+
+	}
+end