path: root/noao/onedspec/irsiids/t_bswitch.x

include	<error.h>
include	<imhdr.h>
include	<mach.h>
include	<time.h>
include	<smw.h>

define	MAX_NR_BEAMS	100	# Max number of instrument apertures
define	MIN_RANGES	100	# Minimum spectra per beam if not given

# T_BSWITCH -- Beam switch a series of spectra to produce a single
#	     sky subtracted spectrum.
#
# The spectra may be extinction corrected if not already done.
#
# The summation may include an optional statistical weighting
# based on the total countrate summed over a user definable
# piece of the spectrum. If the countrate is <= 0, the
# spectrum is given zero weight.
#
# The data may be organized as data from the IIDS/IRS are usually
# obtained - where the telescope is beam-switched so that the
# object is first in one aperture while sky is observed in the other,
# and then the process is reversed.
#
# If the instrument offers many apertures, "nebular" mode can be used
# to obtain the same effect. Here all apertures observe the object(s)
# at one time; then the telescope is moved so all apertures are observing
# sky.
#
# Both these methods are considered "idsmode". But if there are a different
# number of sky observations than object, an imbalance exists. 
# To account for this possibility, all summations are performed by computing
# an average countrate over all observations. Sky countrates can then be
# subtracted from the object. Later the differential countrate is returned
# to an "equivalent" count by multiplying by the exposure time.
#
# Spectra must be dispersion corrected to employ either
# weighting or extinction correction.
#
# The series of spectra may be accumulated in subsets rather than
# over the entire series by specifying a subset rate. (E.g. for
# IIDS data a subset rate of 4 would produce a summed pair for
# every quadruple.)

# Revisions made for WCS support and change from idsmtn.h structure to shdr.h
# structure.  Because this program is an awful mess the changes were made a
# small as possible without altering the structure.  (5/1/91, Valdes)

procedure t_bswitch ()

char	image[SZ_FNAME,MAX_NR_BEAMS+1]
char	rec_numbers[SZ_LINE], title[SZ_LINE,MAX_NR_BEAMS]
char	ofile[SZ_FNAME], stat_fil[SZ_FNAME]
int	sfd, nrecsx
int	i, infile, nrecs, def_beam, start_rec, nimage, sub_rate
int	records[300], beam_stat[MAX_NR_BEAMS], ncols[MAX_NR_BEAMS]
bool	idsmode, extinct, stat, weight, eof_test
pointer	ids[MAX_NR_BEAMS+1] 
pointer	imnames[MAX_NR_BEAMS]	# Hold pointers to pointers of image names
pointer	imin, sp, obs

# The following arrays are suffixed by either 'o' for object or 's' for sky

int	ico   [MAX_NR_BEAMS],   ics   [MAX_NR_BEAMS]	# nr obs in beam
real	expo  [MAX_NR_BEAMS],   exps  [MAX_NR_BEAMS]	# exposure times
pointer	accumo[MAX_NR_BEAMS+1], accums[MAX_NR_BEAMS+1]	# beam accumulators
pointer	counto[MAX_NR_BEAMS],   counts[MAX_NR_BEAMS]	# counts in each obs

int	clpopni(), clgeti(), get_next_image(), decode_ranges()
int	open(), mod()
pointer	immap()
bool	clgetb(), streq()

begin
	call smark (sp)
	call aclri (ids, MAX_NR_BEAMS+1)

	# Open input filename template
	infile = clpopni ("input")

	# Get range specification
	call clgstr ("records", rec_numbers, SZ_LINE)
	if (decode_ranges (rec_numbers, records, 100, nrecs) == ERR)
	    call error (0, "Bad range specification")

	# If no ranges is given, filename expansion will occur, so
	# we must will need some indication of the number of spectra.
	if (nrecs == MAX_INT)
	    nrecsx = MIN_RANGES
	else
	    nrecsx = nrecs

	# Get root name for new records and starting record number
	call clgstr ("output", ofile, SZ_FNAME)
	start_rec = clgeti ("start_rec")

	# Get filename for statistics
	call clgstr ("stats", stat_fil, SZ_FNAME)

	# Assume spectra are in quadruples?
	idsmode = clgetb ("ids_mode")

	# Perform de-extinction?
	extinct = clgetb ("extinct")

	# Use weighting?
	weight = clgetb ("weighting")

	# Accumulate by subsets? - A very large number implies no subsetting
	sub_rate = clgeti ("subset")

	# Open statistics file if any
	if (streq (stat_fil, "")) {
	    sfd = NULL
	    stat = false
	} else {
	    stat = true
	    sfd = open (stat_fil, APPEND, TEXT_FILE)
	}

	# Initialize beam-switch status
	obs = NULL
	call init_file (extinct, def_beam, ico, ics, beam_stat)


	# Begin cycling through all images - accumulate if possible
	# by beam number

	# Initialize range decoder
	call reset_next_image ()

	# Set up for subsets
	nimage = 0
	eof_test = false

	repeat {

	while (get_next_image (infile, records, nrecs, image[1,def_beam],
	    SZ_FNAME) != EOF) {

	    # Attempt to open image with extended header -
	    iferr (imin = immap (image[1,def_beam], READ_ONLY, 0)) {
		call eprintf ("[%s]")
		    call pargstr (image[1,def_beam])
		call error (0, "Image not found or header info not available")
	    }

	    # Add in to accumlators
	    call accum_image (imin, ids, accumo, accums, counto, counts, 
		ico, ics, expo, exps, image, beam_stat, idsmode, extinct, 
		weight, nrecsx, ncols, title, imnames, sfd, obs)

	    call printf ("[%s] added\n")
		call pargstr (image[1,def_beam])
	    call flush (STDOUT)

	    # Close current image
	    call imunmap (imin)

	    # Test for subsets
	    nimage = nimage + 1
	    if (mod (nimage, sub_rate) == 0)
		go to 10
	}

	# Get here by running out of data
	eof_test = true

	# Must be careful not to write out the last sums if subsets are
	# in effect because the subset check would already have done so
	# We can check because "nimage" will not have been bumped
	# if EOF was encountered.

	if (mod (nimage, sub_rate) != 0) {
	
	    # All data has been summed - generate spectra of the accumlations
10	    call wrt_accum (ids, image, title, accumo, accums, ico, ics, 
		counto, counts, expo, exps, ncols, beam_stat, idsmode, weight, 
		extinct, ofile, start_rec, sub_rate)

	    # Generate statistics output for this beam
	    if (stat)
		call wrt_stats (sfd, accumo, accums, ico, ics, counto, counts, 
		    expo, exps, image, beam_stat, title, imnames, weight)

	    # Clear counters and accumulators
	    call reset_beams (accumo, accums, expo, exps, ico, ics, beam_stat, 
		ncols)
	}

	} until (eof_test)

	# Put current record counter back into the parameter file for
	# subsequent invocations
	call clputi ("next_rec", start_rec)

	# Close out inputs, outputs, and space
	do i = 1, MAX_NR_BEAMS+1
	    call shdr_close (ids[i])
	if (obs != NULL)
	    call obsclose (obs)
	call clpcls (infile)
	call close (sfd)
	call sfree (sp)
end

# ACCUM_IMAGE -- Opens current pixel file, loads header elements,
#                adds current spectrum to accumulator array(s),
#                and updates the accumulator status array.
#		 If not in IDSMODE, then returns both object and
#		 sky sums for further consideration.
#		 IDSMODE requires an equal number of each, object and sky,  in
#		 a sequence of OSSO-OSSO or OSSO-SOOS groups.

procedure accum_image (imin, ids, accumo, accums, counto, counts, ico, ics,
                expo, exps, image, beam_stat, idsmode, extinct, weight, nrecs,
		ncols, title, imnames, sfd, obs)

pointer	imin, ids[ARB]
pointer	imnames[ARB]			# Saved image names for stat printout
pointer	sfd				# Statistics file
pointer	obs				# Observatory

pointer	accumo[ARB], accums[ARB]	# Object and sky accumlators
pointer	counto[ARB], counts[ARB]	#                counting stats
real	expo  [ARB], exps  [ARB]	#                total exposure times
int	ico   [ARB], ics   [ARB]	#                number of observations

char	image[SZ_FNAME, MAX_NR_BEAMS+1], title[SZ_LINE,MAX_NR_BEAMS]
char	observatory[SZ_FNAME]
int	beam_stat[ARB], ncols[ARB]
int	dum_beam
bool	idsmode, extinct, weight, exflag, newobs, obshead
real	latitude

int	last_len[MAX_NR_BEAMS], name_nr[MAX_NR_BEAMS]
int	ifile, nr_beams, i, j, def_beam, beam_nr
int	nwaves, ic, nrecs
real	airm, wave1, wave2, wt
pointer	wave_tbl, extn_tbl, ipacc, ipc, mw

real	clgetr(), obsgetr()
pointer	smw_openim()
errchk	smw_openim, shdr_open, obsimopen

begin
	# Bump image file counter
	ifile = ifile + 1

	# Load header area
	mw = smw_openim (imin)
	call shdr_open (imin, mw, 1, 1, INDEFI, SHDATA, ids[def_beam])
	call smw_close (MW(ids[def_beam]))

	accumo[def_beam] = SY(ids[def_beam])

	# Check for proper flags
	call flag_chk (ids[def_beam], exflag)

	if (ifile == 1) {

	    # Get region for statistics to operate over -
	    # Currently only one set of wavelengths is available, but
	    # at some point, it may be desirable to extend this to
	    # provide a start and ending wavelength for each aperture
	    # since an aperture must be considered as an independent
	    # instrument.

	    # Insert defaults --> entire spectrum
	    # Now ask user for start and end - if =0.0, use defaults
	    wave1 = clgetr ("wave1")
	    wave2 = clgetr ("wave2")

	    if (wave1 == 0.0)
		wave1 = W0(ids[def_beam])
	    if (wave2 == 0.0)
		wave2 = W0(ids[def_beam]) + (IM_LEN(imin,1)-1) * 
		    WP(ids[def_beam])

	}

	# Determine beam number and add/sub in pixels
	# Remember that IIDS/IRS "beams" are 0-indexed

	beam_nr = BEAM(ids[def_beam]) + 1
	if (beam_nr > MAX_NR_BEAMS || beam_nr < 1)
	    call error (0, "Illegal beam number")

	# Allocate space for this aperture if not already done
	# Space must be allocated for 2 lines of spectra for
	# each aperture - Line 1 is used to sum up the most
	# recent object-sky spectra to maintain the local
	# statistics. Line 2 is used for the net accumulation
	# over the entire sequence. The statistics from Line 1
	# may be used to weigh the observations as they are
	# added into the Line 2 accumulation.
	#
	# For non-IDSMODE the two lines are used for separate
	# object and sky sums

	if (IS_INDEFI (beam_stat[beam_nr])) {
	    beam_stat[beam_nr] = 0

	    # Allocate space for the accumulators for this beam nr
	    call salloc (accumo[beam_nr], IM_LEN(imin,1), TY_REAL)
	    call salloc (accums[beam_nr], IM_LEN(imin,1), TY_REAL)

	    # Zero object and sky accumulators
	    call amovkr (0.0, Memr[accumo[beam_nr]], IM_LEN(imin,1))
	    call amovkr (0.0, Memr[accums[beam_nr]], IM_LEN(imin,1))


	    # Allocate space for statistics array - For each beam,
	    # a series of up to 'nrecs' spectra may be read, and we
	    # want to keep track of the stats (=countrates) for each
	    # observation. For non-idsmode, need sky rates too.
	    call salloc (counto[beam_nr], nrecs, TY_REAL)
	    if (!idsmode)
		call salloc (counts[beam_nr], nrecs, TY_REAL)

	    # Allocate space for the image names
	    call salloc (imnames[beam_nr], nrecs, TY_INT)
	    name_nr[beam_nr] = 1
	    do j = 1, nrecs
		call salloc (Memi[imnames[beam_nr]+j-1], SZ_LINE, TY_CHAR)

	    # Save number of points for checking purposes
	    last_len[beam_nr] = IM_LEN(imin,1)
	    ncols[beam_nr] = last_len[beam_nr]

	    # Initialize exposure time
	    expo[beam_nr] = 0.0
	    exps[beam_nr] = 0.0

	    nr_beams = nr_beams + 1
	}

	# If this is an object observation, save the image name
	if (OFLAG(ids[def_beam]) == 1) {
	    call strcpy (image[1,def_beam], Memc[Memi[imnames[beam_nr]+
		name_nr[beam_nr]-1]], SZ_LINE)
	    name_nr[beam_nr] = name_nr[beam_nr] + 1
	}

	# If an object observation, save the header elements --
	# NOTE that if we get >1 objects before getting a sky, only
	# the last  observation header is saved!

	# The pixel data will be the sum of all objects until the
	# |object-sky| count = 0 -- Thus, beam switching does not
	# necessarily accumulate by pairs, but depends on how the
	# sequence of observations are presented to the program.

	# The following test has been deleted so that headers
	# will be saved for sky frames as well. This is necessary
	# if BSWITCH is to perform the function of EXTINCTION
	# only when sky frames are to be written out as well.

	if (OFLAG(ids[def_beam]) == 1 || !idsmode) {
	    # Save headers - could probably be done faster by AMOV
	    call shdr_copy (ids[def_beam], ids[beam_nr], NO)

	    # Fix airmass if necessary
	    if (extinct && IS_INDEF (AM(ids[beam_nr]))) {
		call clgstr ("observatory", observatory, SZ_FNAME)
		call obsimopen (obs, imin, observatory, NO, newobs, obshead)
		if (newobs) {
		    call obslog (obs, "BSWITCH", "latitude", STDOUT)
		    if (sfd != NULL)
			call obslog (obs, "BSWITCH", "latitude", sfd)
		}
		latitude = obsgetr (obs, "latitude")
		call get_airm (RA(ids[beam_nr]), DEC(ids[beam_nr]),
		    HA(ids[beam_nr]), ST(ids[beam_nr]), latitude,
		    AM(ids[beam_nr]))
	    }

	    call strcpy (image[1,def_beam], image[1,beam_nr], SZ_FNAME)

	    # Save length - Each beam may be independent sizes
	    ncols[beam_nr] = IM_LEN(imin,1)

	    # Save title, too for same reason
	    call strcpy (IM_TITLE(imin), title[1,beam_nr], SZ_LINE)
	}

	# Verify length
	if (last_len[beam_nr] != ncols[beam_nr]) {
	    call eprintf ("[%s] -- Length not consistent %d\n")
		call pargstr (image[1,beam_nr])
		call pargi   (ncols[beam_nr])
	    ncols[beam_nr] = min (ncols[beam_nr], last_len[beam_nr])
	}
	last_len[beam_nr] = ncols[beam_nr]


	# Check to see if a pair is obtained - then perform statistics
	# and add into global accumulator

	if (idsmode) {

	    # Add spectrum to local accumulation buffer --> Use SKY buffer
	    # At this point of deriving a sequentially local sum, weighting
	    # is not used. 

	    call add_spec (Memr[accumo[def_beam]], Memr[accums[beam_nr]],
		beam_stat[beam_nr], OFLAG(ids[def_beam]), last_len[beam_nr])

	    # IDSMODE requires that every 2N observations produce an
	    # OBJECT-SKY pair
	    if (mod (ifile, 2*nr_beams) == 0)

		# Review all beams in use for non-zero pairings
		do i = 1, MAX_NR_BEAMS
		    if (!IS_INDEFI (beam_stat[i]) && beam_stat[i] != 0)
			call error (0, "Spectra are not in quadruples")


	    # Object and sky exposure times must be equal.
	    if (OFLAG(ids[def_beam]) == 1) {
		expo[beam_nr] = expo[beam_nr] + IT(ids[def_beam])

		# Increment number of object observations for this beam
		ico[beam_nr] = ico[beam_nr] + 1
	    }


	    if (beam_stat[beam_nr] == 0) {
		# Add up all counts within a region for statistics of objects
		# This must be kept separately for each beam number and for
		# each observation

		# First convert to counts per second (CPS)
		call adivkr (Memr[accums[beam_nr]], IT(ids[def_beam]),
		    Memr[accums[beam_nr]], last_len[beam_nr])

		# Sum CPS in statistics region
		call sum_spec (Memr[accums[beam_nr]], wave1, wave2, 
		    W0(ids[def_beam]), WP(ids[def_beam]), Memr[counto[beam_nr]+
		    ico[beam_nr]-1], last_len[beam_nr])

		# De-extinct spectrum
		if (extinct && !exflag) {
		    airm = AM(ids[beam_nr])
		    call de_ext_spec (Memr[accums[beam_nr]], airm, 
			W0(ids[def_beam]), WP(ids[def_beam]), Memr[wave_tbl], 
			Memr[extn_tbl], nwaves, last_len[beam_nr])
		}

		# Add to global accumulator
		# Use weights which are proportional to countrate, if desired
		if (weight) {
		    wt = Memr[counto[beam_nr]+ico[beam_nr]-1]
		    call amulkr (Memr[accums[beam_nr]], wt,
			Memr[accums[beam_nr]], last_len[beam_nr])
		}

		# And add into global sum
		call aaddr (Memr[accums[beam_nr]], Memr[accumo[beam_nr]],
		    Memr[accumo[beam_nr]], last_len[beam_nr])
	    }

	} else {
	    # Non IDSMODE -accumulate separate object and sky CPS sums

	    # Set pointers and update obj-sky parameters
	    if (OFLAG(ids[def_beam]) == 1) {
		beam_stat[beam_nr] = beam_stat[beam_nr] + 1
		ipacc = accumo[beam_nr]
		ipc   = counto[beam_nr]
		ico[beam_nr] = ico[beam_nr] + 1
		ic = ico[beam_nr]
		expo[beam_nr] = expo[beam_nr] + IT(ids[def_beam])
	    } else {
		beam_stat[beam_nr] = beam_stat[beam_nr] - 1
		ipacc = accums[beam_nr]
		ipc   = counts[beam_nr]
		ics[beam_nr] = ics[beam_nr] + 1
		ic = ics[beam_nr]
		exps[beam_nr] = exps[beam_nr] + IT(ids[def_beam])
	    }

	    # First convert to counts per second (CPS)
	    call adivkr (Memr[accumo[def_beam]], IT(ids[def_beam]),
		    Memr[accumo[def_beam]], last_len[beam_nr])

	    # Get counting stats
	    call sum_spec (Memr[accumo[def_beam]], wave1, wave2, 
		W0(ids[def_beam]), WP(ids[def_beam]), Memr[ipc+ic-1], 
		last_len[beam_nr])

	    # De-extinct spectrum
	    if (extinct && !exflag) {
		airm = AM(ids[beam_nr])
		call de_ext_spec (Memr[accumo[def_beam]], airm, 
		    W0(ids[def_beam]), WP(ids[def_beam]), Memr[wave_tbl], 
		    Memr[extn_tbl], nwaves, last_len[beam_nr])
	    }

	    if (weight) {
		wt = Memr[ipc+ic-1]
		call amulkr (Memr[accumo[def_beam]], wt, Memr[accumo[def_beam]],
		    last_len[beam_nr])
	    }

	    # Add into appropriate accumulator
	    call aaddr (Memr[accumo[def_beam]], Memr[ipacc], Memr[ipacc],
		last_len[beam_nr])
	}
	
	return

# INIT_FILE -- Zero the file initializer, the beam counter, beam stats
#	       and read the extinction data if necessary

entry init_file (extinct, dum_beam, ico, ics, beam_stat)

	ifile = 0
	nr_beams = 0
	def_beam = MAX_NR_BEAMS + 1
	dum_beam = def_beam

	do i = 1, MAX_NR_BEAMS {
	    beam_stat[i] = INDEFI
	    ico[i] = 0
	    ics[i] = 0
	}

	# If extinction required, read in extinction file, and sensitivity file
	if (extinct)
	    call get_extn (wave_tbl, extn_tbl, nwaves)

	return

# INIT_NAME -- Reset name index counter for a beam number

entry init_name (dum_beam)

	name_nr[dum_beam] = 1
	return
end

# ACCUM_OUT -- Checks accumulator flags and writes out a new summed
#              image if the count is zero

procedure accum_out (accum, image, ncols, title, root, rec, beam_nr,
	bsflag, itm, exflag)

real	accum[ARB], itm
char	image[SZ_FNAME], title[SZ_LINE], root[SZ_FNAME]
int	ncols, rec, beam_nr
int	bsflag, exflag

pointer	imin, imout, spec
char	bs_image[SZ_FNAME]

pointer	immap(), impl1r()

begin
	# Create new image with user area
	# Use ROOT for spectrum name and increment starting record number

	call sprintf (bs_image, SZ_FNAME, "%s.%04d")
	    call pargstr (root)
	    call pargi (rec)

	rec = rec + 1

	# Provide user info
	call printf ("writing: [%s] %s\n")
	    call pargstr (bs_image)
	    call pargstr (title)
	call flush (STDOUT)

	imin = immap (image, READ_ONLY, 0)
	imout = immap (bs_image, NEW_COPY, imin)

	# Add standard image header
	IM_NDIM(imout) = 1
	IM_LEN(imout,1) = ncols
	IM_PIXTYPE(imout) = TY_REAL
	call strcpy (title, IM_TITLE(imout), SZ_LINE)

	# Write out pixels
	spec = impl1r (imout)
	call amovr (accum, Memr[spec], ncols)

	# Update changed parameters
	if(bsflag == 1)
	    call imaddi (imout, "BS-FLAG", bsflag)
	call imaddr (imout, "EXPTIME", itm)
	call imaddi (imout, "EX-FLAG", exflag)

	call imunmap (imin)
	call imunmap (imout)

	# Store new image name back into image
	call strcpy (bs_image, image, SZ_FNAME)
end

# ACCUM_NORM - Normalize weighted rate and convert to counts

procedure accum_norm (accum, nr, counts, exp, ncols, weight)

real	accum[ARB], counts[ARB], exp
int	nr, ncols
bool	weight

real	sum_wt
int	i

begin
	# The accumulation is an array weighted by non-normalized weights
	# Normalize to total weight to produce a true weighted average
	# and multiply by the total exposure to produce
	# an equivalent sum

	# Add up all weighting factors
	if (weight) {
	    sum_wt = 0.0
	    do i = 1, nr
		sum_wt = sum_wt + counts[i]
	} else
	    sum_wt = real (nr)

	if (sum_wt == 0.0)
	    sum_wt = 1.0

	# Correct for exposure time
	sum_wt = exp / sum_wt

	call amulkr (accum, sum_wt, accum, ncols)
end

# WRT_ACCUM -- Write out accumulations as spectra

procedure wrt_accum (ids, image, title, accumo, accums, ico, ics, 
	counto, counts, expo, exps, ncols, beam_stat, idsmode, weight, 
	extinct, ofile, start_rec, sub_rate)

pointer	ids[ARB]
char	image[SZ_FNAME,MAX_NR_BEAMS+1], title[SZ_LINE,MAX_NR_BEAMS]

pointer	accumo[ARB], accums[ARB]
pointer	counto[ARB], counts[ARB]
int	ico   [ARB], ics   [ARB]
real	expo  [ARB], exps  [ARB]
int	ncols[ARB]
int	beam_stat[ARB]
bool	idsmode, weight, extinct
char	ofile[SZ_FNAME]
int	start_rec, sub_rate, bsflag

int	i, nr_beams
real	exp_ratio

begin
	# First compute number of beams
	nr_beams = 0
	do i = 1, MAX_NR_BEAMS
	    if (!IS_INDEFI (beam_stat[i]) && ((ico[i] > 0) || (ics[i] > 0)))
		nr_beams = nr_beams + 1

	# For all present apertures, write out a spectrum
	do i = 1, MAX_NR_BEAMS {

	    if (!IS_INDEFI (beam_stat[i]) && ((ico[i] > 0) || (ics[i] > 0))) {
		if (beam_stat[i] != 0 && idsmode) {
		    call eprintf ("Non-equal number of obj-sky observations")
		    call eprintf (" beam: %d - residual: %d\n")
			call pargi (i-1)
		        call pargi (beam_stat[i])

		    # Reset to 0 and force output
		    beam_stat[i] = 0
		}

		# The accumulator has a total CPS using non-normalized
		# weights - apply normalization and exposure time to
		# generate an equivalent COUNT sum.
		call accum_norm (Memr[accumo[i]], ico[i], Memr[counto[i]], 
			expo[i], ncols[i], weight)

		if (!idsmode) {
		    # Separate object and sky sums require sky info
		    call accum_norm (Memr[accums[i]], ics[i], Memr[counts[i]],
			exps[i], ncols[i], weight)

		    # Then normalize sky exposure time to that of object
		    if (exps[i] != 0.0)
			exp_ratio = expo[i]/exps[i]
		    else
			exp_ratio = 1.0

		    # Check that some object observtion was made
		    # If not, then we only have sky data so multiply by -1
		    # so that the subsequent subtraction will produce a
		    # positive sky
		    if (expo[i] == 0.0)
			exp_ratio = -1.0

		    if (exp_ratio != 1.0)
			call amulkr (Memr[accums[i]], exp_ratio,
			    Memr[accums[i]], ncols[i])

		    # Finally subtract sky from object equivalent counts
		    call asubr (Memr[accumo[i]], Memr[accums[i]], 
			Memr[accumo[i]], ncols[i])

		}
		# Set header flags
		# BS flag is not set if the subset rate equals the
		# number of apertures since each record in is copied out
		if (sub_rate > nr_beams)
		    bsflag = 1
		else
		    bsflag = -1

		if (OFLAG(ids[i]) == 1)
		    IT(ids[i]) = expo[i]
		else
		    IT(ids[i]) = exps[i]

		if (extinct)
		    EC(ids[i]) = 0

		# And write out spectrum, at last
		call accum_out (Memr[accumo[i]], image[1,i], 
		    ncols[i], title[1,i], ofile, start_rec, i,
		    bsflag, IT(ids[i]), EC(ids[i]))

		# Reset name entry counter
		call init_name (i)
	    }

	}
end

# RESET_BEAMS -- Zeroes the counters and accumulators for additional
#                cases

procedure reset_beams (accumo, accums, expo, exps, ico, ics, beam_stat, ncols)

pointer	accumo[ARB], accums[ARB]
real	expo  [ARB], exps  [ARB]
int	ico   [ARB], ics   [ARB]
int	beam_stat[ARB], ncols[ARB]

int	i

begin
	do i = 1, MAX_NR_BEAMS
	    if (!IS_INDEFI (beam_stat[i])) {

		expo[i] = 0.0
		exps[i] = 0.0
		ico[i]  = 0
		ics[i]  = 0
		call amovkr (0.0, Memr[accumo[i]], ncols[i])
		call amovkr (0.0, Memr[accums[i]], ncols[i])
	    }
end

# WRT_STATS -- Write out statistics file

procedure wrt_stats (fd, accumo, accums, ico, ics, counto, counts,
     expo, exps, image, beam_stat, title, imnames, weight)

int	fd
pointer	accumo[ARB], accums[ARB], counto[ARB], counts[ARB]
real	expo[ARB], exps[ARB]
int	ico[ARB], ics[ARB], beam_stat[ARB]
char	image[SZ_FNAME,MAX_NR_BEAMS+1]
char	title[ARB]
pointer	imnames[ARB]
bool	weight

int	i, j
real	cmaxo, cmaxs
char	ctime[SZ_TIME]

long	clktime()

begin
	# Issue time stamp
	call cnvtime (clktime (long(0)), ctime, SZ_TIME)
	call fprintf (fd, "%s\n\n")
	    call pargstr (ctime)

	# Issue message if weighted sums in effect
	if (weight)
	    call fprintf (fd, "--> Using weighted averages <--\n\n")

	# Cycle over beams
	do i = 1, MAX_NR_BEAMS {
	    if (!IS_INDEFI (beam_stat[i])) {

		# Write out Object stats if any
		if (ico[i] > 0) {
		    call fprintf (fd, "Object statistics for beam %d -->[%s]\n")
			call pargi (i-1)
			call pargstr (image[1,i])
		    call fprintf (fd, "Title: %s\n")
			call pargstr (title)

		    # Find maximum count value for this beam
		    cmaxo = Memr[counto[i]]

		    do j = 1, ico[i]
		        cmaxo = max (cmaxo, Memr[counto[i]+j-1])

		    call fprintf (fd,
			"Obs  Relative CPS   Image%12wPeak CPS = %10.3g\n")
			call pargr (cmaxo)

		    if (cmaxo == 0.0)
			cmaxo = 1.0

		    do j = 1, ico[i] {
			call fprintf (fd, "%3d    %8.2f    [%s]\n")
			    call pargi (j)
			    call pargr (Memr[counto[i]+j-1] / cmaxo)
			    call pargstr (Memc[Memi[imnames[i]+j-1]])
		    }
		}

		# Write out sky stats if any
		if (ics[i] > 0) {
		    call fprintf (fd, "Sky statistics for beam %d\n")
			call pargi (i-1)

		    cmaxs = Memr[counts[i]]

		    do j = 1, ics[i]
		        cmaxs = max (cmaxs, Memr[counts[i]+j-1])

		    call fprintf (fd, "Obs  Relative CPS   Peak CPS = %10.3g\n")
			call pargr (cmaxs)

		    if (cmaxs == 0.0)
			cmaxs = 1.0

		    do j = 1, ics[i] {
			call fprintf (fd, "%3d    %8.2f\n")
			    call pargi (j)
			    call pargr (Memr[counts[i]+j-1] / cmaxs)
		    }
		}

		call fprintf (fd, "\n\n")
	    }
	}
end


# ADD_SPEC -- Accumulate spectrum into array - either add or subtract
#             Returns status = net number of object - sky apectra
#                            = 0 for equal numbers to indicate further
#                              processing may take place

procedure add_spec (inspec, accum, stat, flag, len)

real	inspec[ARB], accum[ARB]
int	stat, flag, len

int	i, add_sub

begin
	add_sub = 0

	# Is this an Object or Sky?
	# If flag is neither 0 or 1, spectrum is ignored
	if (flag == 0)
	    add_sub = -1
	if (flag == 1)
	    add_sub = +1

	if (add_sub == 0) {
	    stat = INDEFI
	    return
	}

	# Is accumulator to be cleared?
	if (IS_INDEFI (stat) || stat == 0) {
	    call amulkr (inspec, real (add_sub), accum, len)
	    stat = add_sub

	} else {
	    # Add into accumulator
	    do i = 1, len
		accum[i] = accum[i] + add_sub * inspec[i]

	    stat = stat + add_sub
	}
end

# FLAG_CHK -- Check header flags prior to beam switching

procedure flag_chk (ids, exflag)

pointer	ids
bool	exflag

int	bsflag, imgeti()

begin
	# BS requires 
	# 1. dispersion corrected spectra
	# 2. non-beam switched
	# 3. may be either extinction corrected or not

	if (DC(ids) != DCLINEAR)
	    call error (0, "Spectrum not dispersion corrected")

	iferr (bsflag = imgeti (IM(ids), "BS-FLAG"))
	    bsflag = -1
	if (bsflag == 1)
	    call error (0, "Spectrum already beam-switched")

	if (EC(ids) == ECYES)
	    exflag = true
	else
	    exflag = false
end