1 files changed, 195 insertions, 0 deletions
diff --git a/noao/onedspec/irsiids/t_addsets.x b/noao/onedspec/irsiids/t_addsets.x
new file mode 100644
index 00000000..cd145b4a
--- /dev/null
+++ b/noao/onedspec/irsiids/t_addsets.x
@@ -0,0 +1,195 @@
+include	<error.h>
+include	<imhdr.h>
+
+
+# T_ADDSETS -- Add a series of spectra by subsets. A single spectrum
+#  is produced for every "subset" number of input spectra. The input
+#  list is accumulated until "subset" number of spectra have been
+#  encountered. The result is then written out.
+#
+#  If the input data are calibrated (CA_FLAG = 0) then the result
+#  is an average over the subset size, but the header exposure
+#  time is updated.
+#
+#  If the data is uncalibrated then the resulting spectrum is a sum
+#  of the total counts.
+
+procedure t_addsets ()
+
+pointer	image
+pointer	recstr, ofile
+int	root, start_rec, subset
+int	nrecs
+int	nrem, ifile, ca_flag
+real	itm, expo, wt, wtsum
+bool	weight
+pointer	sp, recs, im, cur_pix, sp_sum
+
+real	imgetr()
+int	clpopni(), clgeti(), imgeti()
+int	get_next_image(), decode_ranges()
+bool	clgetb()
+pointer	immap(), imgl1r()
+
+begin
+	call smark (sp)
+	call salloc (image, SZ_FNAME, TY_CHAR)
+	call salloc (ofile, SZ_FNAME, TY_CHAR)
+	call salloc (recstr, SZ_LINE, TY_CHAR)
+	call salloc (recs, 300, TY_INT)
+
+	# Open input file name template
+	root   = clpopni ("input")
+
+	# Get range specification if any
+	call clgstr ("records", Memc[recstr], SZ_LINE)
+	if (decode_ranges (Memc[recstr], Memi[recs], 100, nrecs) == ERR)
+	    call error (0, "Bad range specification")
+
+	# Get rootname for output files and starting record
+	call clgstr ("output", Memc[ofile], SZ_FNAME)
+	start_rec = clgeti ("start_rec")
+
+	# Get subset size
+	subset = clgeti ("subset")
+
+	# Apply integration time weighting?
+	weight = clgetb ("weighting")
+
+	# Initialize range decoder
+	call reset_next_image ()
+
+	#Initialize file counter
+	ifile = 0
+	wtsum = 0.0
+
+	# Loop over all input images by subsets
+	while (get_next_image (root, Memi[recs], nrecs, Memc[image], 
+	    SZ_FNAME) != EOF) {
+
+	    # Open image
+	    iferr (im = immap (Memc[image], READ_ONLY, 0)) {
+		call erract (EA_WARN)
+		next
+	    }
+
+	    # Allocate space for current subset
+	    if (mod (ifile, subset) == 0) {
+		call calloc (sp_sum, IM_LEN (im,1), TY_REAL)
+
+		# Zero exposure counter
+		expo = 0.0
+	    }
+
+	    # Add in current spectrum
+	    iferr (itm = imgetr (im, "EXPOSURE"))
+		iferr (itm = imgetr (im, "ITIME"))
+		    iferr (itm = imgetr (im, "EXPTIME"))
+			itm = 1
+	    iferr (ca_flag = imgeti (im, "CA-FLAG"))
+		ca_flag = -1
+	    cur_pix = imgl1r (im)
+
+	    # Apply integration time weighting
+	    if (weight)
+		wt = itm
+	    else
+		wt = 1.0
+
+	    if (ca_flag != 0)
+		wt = 1.0
+
+	    wtsum = wtsum + wt
+	    call amulkr (Memr[cur_pix], wt, Memr[cur_pix], IM_LEN(im,1))
+	    call aaddr (Memr[cur_pix], Memr[sp_sum], Memr[sp_sum], IM_LEN(im,1))
+	    expo = expo + itm
+
+	    # Issue status report
+	    call printf ("[%s] added\n")
+		call pargstr (Memc[image])
+
+	    ifile = ifile + 1
+	    if (mod (ifile, subset) == 0) {
+		call wrt_set (Memr[sp_sum], subset, im, Memc[ofile], start_rec,
+		    expo, wtsum, ca_flag)
+		wtsum = 0.0
+		call mfree (sp_sum, TY_REAL)
+	    } else
+		call imunmap (im)
+
+	}
+	# Check that there are no remaining spectra in an unfulfilled subset
+	nrem = mod (ifile, subset)
+	if (nrem != 0) {
+	    call wrt_set (Memr[sp_sum], nrem, im, Memc[ofile], start_rec, 
+		expo, wtsum, ca_flag)
+	    wtsum = 0.0
+	    call mfree (sp_sum, TY_REAL)
+
+	    call eprintf ("Unfulfilled subset accumulation written - ")
+	    call eprintf ("missing %d spectra\n")
+		call pargi (subset - nrem)
+	}
+
+	# Update record number
+	call clputi ("next_rec", start_rec)
+
+	# Free space
+	call sfree (sp)
+	call clpcls (root)
+end
+
+# WRT_SET -- Write spectra ccumulated from the set
+
+procedure wrt_set (sp_sum, subset, im, ofile, start_rec, expo, wtsum, ca_flag)
+
+real	sp_sum[ARB]
+int	subset, start_rec, ca_flag
+pointer	im
+char	ofile[SZ_FNAME]
+real	expo, wtsum
+
+char	newfile[SZ_FNAME]
+pointer imnew, newpix
+
+pointer	impl1r(), immap()
+int	strlen()
+
+begin
+	# Create new spectrum - first make up a name
+	call strcpy (ofile, newfile, SZ_FNAME)
+	call sprintf (newfile[strlen (newfile) + 1], SZ_FNAME, ".%04d")
+	    call pargi (start_rec)
+
+	imnew = immap (newfile, NEW_COPY, im)
+
+	IM_NDIM (imnew) = 1
+	IM_LEN (imnew,1) = IM_LEN (im,1)
+	IM_PIXTYPE (imnew) = TY_REAL
+	call strcpy (IM_TITLE(im), IM_TITLE(imnew), SZ_LINE)
+
+	call imunmap (im)
+
+	newpix = impl1r (imnew)
+
+	# If this spectrum is calibrated, perform an average
+	# weighted by integration time and copy new pixels into image
+	if (ca_flag == 0)
+	    call adivkr (sp_sum, real (wtsum), Memr[newpix], IM_LEN(imnew,1))
+	else
+	    call amovr (sp_sum, Memr[newpix], IM_LEN(imnew,1))
+
+	# Update keyword
+	call imaddr (imnew, "EXPOSURE", expo)
+
+	# Send user report
+	call printf ("writing [%s]: %s\n")
+	    call pargstr (newfile)
+	    call pargstr (IM_TITLE(imnew))
+	call flush (STDOUT)
+
+	call imunmap (imnew)
+
+	# Update record counter
+	start_rec = start_rec + 1
+end