1 files changed, 137 insertions, 0 deletions
diff --git a/noao/twodspec/multispec/t_findpeaks.x b/noao/twodspec/multispec/t_findpeaks.x
new file mode 100644
index 00000000..2e4cf79e
--- /dev/null
+++ b/noao/twodspec/multispec/t_findpeaks.x
@@ -0,0 +1,137 @@
+include <imhdr.h>
+include <fset.h>
+include "ms.h"
+
+# T_FIND_PEAKS -- Find the spectra peaks in a MULTISPEC image and record
+# their positions in the database.
+#
+# An average of naverage lines from the MULTISPEC image is searched
+# for peaks satisfying constraints on the minimum and maximum number,
+# columns, peak values, and separation between peaks.  The positions
+# of the peaks satisfying these constraints is entered in the database.
+# It is an error if fewer than the minimum number of peaks is found
+# or if the number of peaks differs from a previously determined number.
+# The peak finding is done by the function FIND_PEAKS which is numerical
+# and may be used outside the MULTISPEC package.
+
+procedure t_find_peaks ()
+
+# CL parameters:
+char	image[SZ_FNAME]		# Image to be searched
+int	lines[3, MAX_RANGES]	# Image lines in which to find spectra
+int	min_npeaks		# Minimum number of spectra to be found
+int	max_npeaks		# Maximum number of spectra to be accepted
+int	separation		# Minimum pixel separation between spectra
+int	edge			# Minimum distance to edge of image
+real	threshold		# Minimum peak value
+real	contrast		# Max contrast between strongest and weakest
+int	columns[3, MAX_RANGES]	# Spectra positions limited to these columns
+int	naverage		# Number of image lines to average
+bool	debug			# Print debugging information
+
+char	comment[SZ_LINE]
+int	i, j, k, line, sample, nsamples, npoints, nspectra
+pointer	ms, im
+pointer	sp, data, x, samples
+
+int	find_peaks(), get_sample_lines()
+int	clgeti(), clgranges()
+real	clgetr()
+bool	clgetb(), is_in_range()
+pointer	msmap(), immap()
+
+begin
+	# Get task parameters and access files.
+	call clgstr ("image", image, SZ_FNAME)
+	ms = msmap (image, READ_WRITE, 0)
+	im = immap (image, READ_ONLY, 0)
+	i = clgranges ("lines", 1, IM_LEN(im, 2), lines, MAX_RANGES)
+	min_npeaks = clgeti ("min_npeaks")
+	max_npeaks = clgeti ("max_npeaks")
+	separation = clgeti ("separation")
+	edge = clgeti ("edge")
+	threshold = clgetr ("threshold")
+	contrast = clgetr ("contrast")
+	i = clgranges ("columns", 1, IM_LEN(im, 1), columns, MAX_RANGES)
+	naverage = clgeti ("naverage")
+	debug = clgetb ("debug")
+
+	call fseti (STDOUT, F_FLUSHNL, YES)
+
+	# Allocate working memory.
+	npoints = IM_LEN(im, 1)
+	call smark (sp)
+	call salloc (samples, MS_NSAMPLES(ms), TY_INT)
+	call salloc (data, npoints, TY_REAL)
+	call salloc (x, npoints, TY_REAL)
+
+	# Get the sample lines.
+	nsamples = get_sample_lines (ms, lines, Memi[samples])
+
+	# Loop through each sample line.
+	do i = 1, nsamples {
+	    sample = Memi[samples + i - 1]
+	    line = LINE(ms, sample)
+
+	    # Get the image data with averaging.
+	    call  msgimage (im, line, naverage, Memr[data])
+
+	    # Mark columns which are to be ignored with INDEFR.
+	    do j = 1, npoints
+	        if (!is_in_range (columns, j))
+		    Memr[data + j - 1] = INDEFR
+
+	    # Find the peaks.
+	    nspectra = find_peaks (Memr[data], Memr[x], npoints,
+		contrast, separation, edge, max_npeaks, threshold, debug)
+
+	    if (debug) {
+	        call printf ("  Number of spectra found in line %d = %d.\n")
+		    call pargi (line)
+	            call pargi (nspectra)
+	    }
+	    if (nspectra < min_npeaks)
+		call error (MS_ERROR, "Too few spectra found")
+
+	    # Enter the spectra found in the database.  If the number of
+	    # spectra has not been previously set in the database then
+	    # enter the number of spectra and make entries in the
+	    # database.  Otherwise check that the number of spectra found
+	    # agrees with that already in the database.
+
+	    if (MS_NSPECTRA(ms) == 0) {
+	        if (nspectra == 0)
+		    next
+		MS_NSPECTRA(ms) = nspectra
+	        call dbenter (MS_DB(ms), NAME(ms, I0), nspectra * SZ_REAL,
+		    MS_NSAMPLES(ms))
+	        call dbenter (MS_DB(ms), NAME(ms, X0), nspectra * SZ_REAL,
+		    MS_NSAMPLES(ms))
+	    } else if (MS_NSPECTRA(ms) != nspectra)
+		call error (MS_ERROR, "Attempt to change the number of spectra")
+
+	    call msgparam (ms, X0, sample)
+	    call amovr (Memr[x], PARAMETER(ms, X0, 1), nspectra)
+	    call mspparam (ms, X0, sample)
+
+	    # The peak scale is taken and the pixel value at the peak.
+	    call msgparam (ms, I0, sample)
+	    do j = 1, nspectra {
+		k = PARAMETER(ms, X0, j)
+		PARAMETER(ms, I0, j) = Memr[data + k - 1]
+	    }
+	    call mspparam (ms, I0, sample)
+
+	    # Enter a comment in the database.
+	    call sprintf (comment, SZ_LINE,
+		"Spectra located in sample line %d.")
+		call pargi (sample)
+	    call history (ms, comment)
+	}
+
+	# Update the database and close the database and image.
+	call msphdr (ms)
+	call msunmap (ms)
+	call imunmap (im)
+	call sfree (sp)
+end