1 files changed, 70 insertions, 0 deletions
diff --git a/pkg/xtools/peaks.x b/pkg/xtools/peaks.x
new file mode 100644
index 00000000..dfad9abb
--- /dev/null
+++ b/pkg/xtools/peaks.x
@@ -0,0 +1,70 @@
+# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
+
+# PEAKS -- find the peaks in an array of x and y points.
+# The extrema in the input data points are found using extrema(xtools).
+# The extrema are located to a precision of dx.
+# The extrema with negative curvature (peaks) are selected and returned
+# in the x array.  The spline value is returned in the y array.  The
+# background is estimated by linear interpolation of the neighboring
+# minima (extrema of positive curvature) to the position of the peak.
+# The background is returned in the background array.  The number of
+# peaks found is returned as the function value.
+
+int procedure peaks (x, y, background, npts, dx)
+
+real	x[npts], y[npts]	# Input data points and output peaks
+real	background[npts]	# Background estimate
+int	npts			# Number of input data points
+real	dx			# Precision of peak positions
+
+int	i, j, k, npeaks, nextrema
+pointer	sp, a, b, c
+
+int	extrema()
+errchk	salloc
+
+begin
+	nextrema = extrema (x, y, background, npts, dx)
+
+	if (nextrema == 0)
+	    return (0)
+
+	# Allocate working storage
+	call smark (sp)
+	call salloc (a, nextrema, TY_REAL)
+	call salloc (b, nextrema, TY_REAL)
+	call salloc (c, nextrema, TY_REAL)
+
+	npeaks = 0
+	do i = 1, nextrema {
+	    if (background[i] < 0.) {
+		Memr[a + npeaks] = x[i]
+		Memr[b + npeaks] = y[i]
+		for (j = i - 1; j > 0; j = j - 1)
+		    if (background[j] > 0.)
+			break;
+		for (k = i + 1; k <= npts; k = k + 1)
+		    if (background[k] > 0.)
+			break;
+		if ((j >= 1) && (k <= npts))
+		    Memr[c + npeaks] =
+			(y[k] - y[j]) / (x[k] - x[j]) * (x[i] - x[j]) + y[j]
+		else if (j >= 1)
+		    Memr[c + npeaks] = y[j]
+		else if (k <= npts)
+		    Memr[c + npeaks] = y[k]
+		else {
+		    call sfree (sp)
+		    call error (1, "No background points")
+		}
+		npeaks = npeaks + 1
+	    }
+	}
+
+	call amovr (Memr[a], x, npeaks)
+	call amovr (Memr[b], y, npeaks)
+	call amovr (Memr[c], background, npeaks)
+
+	call sfree (sp)
+	return (npeaks)
+end