STScI Logo

t_sums



include	<error.h>
include	<imhdr.h>

define	MAX_NR_BEAMS	100	# Max number of instrument apertures

# T_SUMS -- Compute sums of strings of spectra according to
#	    Aperture number and object/sky flag. So for IIDS/IRS
#           type spectra, 4 sums will be generated.
#           In general, there will be 2N sums where N is the number
#           apertures.

procedure t_sums ()

pointer	image					# Image name to be added
pointer	images					# Image name to be added
pointer	ofile					# Output image file name
pointer	recstr					# Record number string
int	recs					# Spectral record numbers
int	root, nrecs				# CL and ranges flags
real	expo					# Exposure time
pointer	bstat[2]				# Status of each aperture
pointer	npts[2]					# Length of spectrum
pointer	esum[2]					# Accumulated exposure time
pointer	accum[2]				# Pointers to beam accumulators
pointer	title[2]
int	beam, object
int	start_rec

int	i, j
pointer	sp, work, im

real	imgetr()
int	clgeti(), clpopni(), imgeti()
int	get_next_image(), decode_ranges()
pointer	immap()

begin
	call smark (sp)
	call salloc (image, SZ_FNAME, TY_CHAR)
	call salloc (images, MAX_NR_BEAMS, TY_POINTER)
	call salloc (ofile, SZ_FNAME, TY_CHAR)
	call salloc (recstr, SZ_LINE, TY_CHAR)
	call salloc (recs, 300, TY_INT)
	call salloc (accum, MAX_NR_BEAMS, TY_POINTER)
	call salloc (title, MAX_NR_BEAMS, TY_POINTER)
	call amovki (NULL, Memi[images], MAX_NR_BEAMS)
	call salloc (work, 2*5*MAX_NR_BEAMS, TY_STRUCT)
	bstat[1] = work
	bstat[2] = work + MAX_NR_BEAMS
	npts[1] = work + 2 * MAX_NR_BEAMS
	npts[2] = work + 3 * MAX_NR_BEAMS
	esum[1] = work + 4 * MAX_NR_BEAMS
	esum[2] = work + 5 * MAX_NR_BEAMS
	accum[1] = work + 6 * MAX_NR_BEAMS
	accum[2] = work + 7 * MAX_NR_BEAMS
	title[1] = work + 8 * MAX_NR_BEAMS
	title[2] = work + 9 * MAX_NR_BEAMS

	# Get task parameters.
	root = clpopni ("input")

	# Get input record numbers
	call clgstr ("records", Memc[recstr], SZ_LINE)
	if (decode_ranges (Memc[recstr], Memi[recs], 100, nrecs) == ERR)
	    call error (0, "Bad range specification")

	call clgstr ("output", Memc[ofile], SZ_LINE)

	start_rec = clgeti ("start_rec")

	call reset_next_image ()

	# Clear all beam status flags
	call amovki (INDEFI, Memi[bstat[1]], MAX_NR_BEAMS*2)
	call aclrr (Memr[esum[1]], MAX_NR_BEAMS*2)

	call printf ("Accumulating spectra --\n")
	call flush (STDOUT)

	while (get_next_image (root, Memi[recs], nrecs, Memc[image],
	    SZ_FNAME) != EOF) {
	    iferr (im = immap (Memc[image], READ_ONLY, 0)) {
		call erract (EA_WARN)
		next
	    }

	    # Load header
	    iferr (beam = imgeti (im, "BEAM-NUM"))
		beam = 0
	    if (beam < 0  || beam > MAX_NR_BEAMS-1)
		call error (0, "Invalid aperture number")

	    # Select array: Object = array 2; sky = array 1
	    iferr (object = imgeti (im, "OFLAG"))
		object = 1
	    if (object == 1)
		object = 2
	    else
		object = 1

	    iferr (expo = imgetr (im, "EXPOSURE"))
		iferr (expo = imgetr (im, "ITIME"))
		    iferr (expo = imgetr (im, "EXPTIME"))
			expo = 1

	    # Add spectrum into accumulator
	    if (IS_INDEFI (Memi[bstat[object]+beam])) {
	        Memi[npts[object]+beam] = IM_LEN (im,1)
	        call salloc (Memi[accum[object]+beam], IM_LEN(im,1), TY_REAL)
	        call aclrr (Memr[Memi[accum[object]+beam]], IM_LEN(im,1))
	        Memi[bstat[object]+beam] = 0

	        call salloc (Memi[title[object]+beam], SZ_LINE, TY_CHAR)
	        call strcpy (IM_TITLE(im), Memc[Memi[title[object]+beam]],
		    SZ_LINE)
	    }

	    call su_accum_spec (im, Memi[npts[1]], expo, Memi[bstat[1]],
		beam+1, Memi[accum[1]], Memr[esum[1]], Memi[title[1]], object)

	    call printf ("[%s] %s spectrum added to aperture %1d\n")
		call pargstr (Memc[image])
		if (object == 2)
		    call pargstr ("object")
		else
		    call pargstr ("sky   ")
		call pargi (beam)
	    call flush (STDOUT)

	    if (Memi[images+beam] == NULL)
		call salloc (Memi[images+beam], SZ_FNAME, TY_CHAR)
	    call strcpy (Memc[image], Memc[Memi[images+beam]], SZ_FNAME)
	    call imunmap (im)
	}

	# Review all apertures containing data and write sums
	do i = 0, MAX_NR_BEAMS-1
	    do j = 1, 2
	    if (!IS_INDEFI (Memi[bstat[j]+i])) {
		call wrt_spec (Memc[Memi[images+i]], Memr[Memi[accum[j]+i]],
		    Memr[esum[j]+i], Memc[ofile], start_rec,
		    Memc[Memi[title[j]+i]], Memi[npts[j]+i], i, j)

		start_rec = start_rec + 1
	    }

	call clputi ("next_rec", start_rec)
	call sfree (sp)
	call clpcls (root)
end

# ACCUM_SPEC -- Accumulate spectra by beams

procedure su_accum_spec (im, len, expo, beam_stat, beam, accum, expo_sum,
	title, object)

pointer	im, accum[MAX_NR_BEAMS,2], title[MAX_NR_BEAMS,2]
real	expo, expo_sum[MAX_NR_BEAMS,2]
int	beam_stat[MAX_NR_BEAMS,2], beam, len[MAX_NR_BEAMS,2]
int	object

int	npts
pointer	pix

pointer	imgl1r()

begin
	npts = IM_LEN (im, 1)

	# Map pixels and optionally correct for coincidence
	pix = imgl1r (im)

	# Add in the current data
	npts = min (npts, len[beam, object])

	call aaddr (Memr[pix], Memr[accum[beam, object]], 
	    Memr[accum[beam, object]], npts)

	beam_stat[beam, object] = beam_stat[beam, object] + 1
	expo_sum [beam, object] = expo_sum [beam, object] + expo
end

# WRT_SPEC -- Write out normalized spectrum

procedure wrt_spec (image, accum, expo_sum, ofile, start, title, npts, object,
	beam)

char	image[SZ_FNAME]
real	accum[ARB], expo_sum
int	start, npts
char	ofile[SZ_FNAME]
char	title[SZ_LINE]
int	object, beam

char	output[SZ_FNAME], temp[SZ_LINE]
pointer	im, imnew, newpix

pointer	immap(), impl1r()
int	strlen()

begin
	im = immap (image, READ_ONLY, 0)
10	call strcpy (ofile, output, SZ_FNAME)
	call sprintf (output[strlen (output) + 1], SZ_FNAME, ".%04d")
	    call pargi (start)

	# Create new image with a user area
	# If an error occurs, ask user for another name to try
	# since many open errors result from trying to overwrite an
	# existing image.

	iferr (imnew = immap (output, NEW_COPY, im)) {
	    call eprintf ("Cannot create [%s] -- Already exists??\07\n")
		call pargstr (output)
	    call clgstr ("newoutput", ofile, SZ_FNAME)
	    go to 10
	}

	call strcpy ("Summation:", temp, SZ_LINE)
	call sprintf (temp[strlen (temp) + 1], SZ_LINE, "%s")
	    call pargstr (title)
	call strcpy (temp, IM_TITLE (imnew), SZ_LINE)

	newpix = impl1r (imnew)
	call amovr (accum, Memr[newpix], npts)

	call imaddr (imnew, "EXPOSURE", expo_sum)
	call imunmap (im)
	call imunmap (imnew)

	call printf ("%s sum for aperture %1d --> [%s]\n")
	    if (object == 1)
		call pargstr ("Object")
	    else
		call pargstr ("Sky   ")
	    call pargi (beam)
	    call pargstr (output)
	call flush (STDOUT)
end

Source Code · Search Form · STSDAS