1 files changed, 70 insertions, 0 deletions
diff --git a/pkg/xtools/extrema.x b/pkg/xtools/extrema.x
new file mode 100644
index 00000000..0a373aa5
--- /dev/null
+++ b/pkg/xtools/extrema.x
@@ -0,0 +1,70 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+define	ORDER		4	# The order of the spline
+
+# EXTREMA -- Find the extrema in an array of x and y points.
+# The input data points are fitted with a cubic interpolation spline.  The
+# spline is then searched for points where the first derivative changes sign.
+# The minimum step size of this search is controlled by the parameter dx.
+# The positions of these extrema are returned in the x array, the value of the
+# spline at the extrema are returned in the y array, and the curvature or
+# second derivative of the spline at the extrema are returned in the
+# curvature array.  The function returns the number of extrema found.
+
+int procedure extrema (x, y, curvature, npts, dx)
+
+real	x[npts], y[npts]	# Input data points and output extrema
+real	curvature[npts]		# 2nd deriv. of cubic spline at extrema
+int	npts			# Number of input data points
+real	dx			# Precision of extrema positions
+
+int	i, ier, nextrema
+real	xeval, left_deriv, right_deriv
+pointer	sp, bspln, q
+real	seval()
+errchk	salloc, seval
+
+begin
+	# Allocate working arrays for spline routines
+	call smark (sp)
+	call salloc (bspln, 2 * npts + 30, TY_REAL)
+	call salloc (q, (2 * ORDER - 1) * npts, TY_REAL)
+
+	# Calculate the spline coefficients
+	call spline (x, y, npts, Memr[bspln], Memr[q], ORDER, ier)
+	if (ier != 0) {
+	    call sfree (sp)
+	    return (0)
+	}
+
+	# Initialize the curvature array
+	call aclrr (curvature, npts)
+
+	# Find the extrema defined by a change in sign in the first derivative.
+	nextrema = 0
+	left_deriv = seval (x[1], 1, Memr[bspln])
+	do i = 2, npts {
+	    xeval = x[i]
+	    right_deriv = seval (xeval, 1, Memr[bspln])
+	    if (left_deriv * right_deriv <= 0.) {
+		for (xeval = x[i - 1] + dx; xeval <= x[i]; xeval = xeval + dx) {
+		    right_deriv = seval (xeval, 1, Memr[bspln])
+	    	    if (left_deriv * right_deriv <= 0.)
+			break
+		    left_deriv = right_deriv
+		}
+		nextrema = nextrema + 1
+		x[nextrema] = xeval
+		y[nextrema] = seval (xeval, 0, Memr[bspln])
+		curvature[nextrema] = seval (xeval, 2, Memr[bspln])
+		if (curvature[nextrema] == 0.)
+		    nextrema = nextrema - 1
+		if (nextrema == npts)
+		    break
+	    }
+	    left_deriv = right_deriv
+	}
+
+	call sfree (sp)
+	return (nextrema)
+end