Repatch (from linux) of OSX IRAF

1 files changed, 585 insertions, 0 deletions

diff --git a/noao/onedspec/t_sbands.x b/noao/onedspec/t_sbands.x
new file mode 100644
index 00000000..d3c3d1e2
--- /dev/null
+++ b/noao/onedspec/t_sbands.x
@@ -0,0 +1,585 @@
+include	<error.h>
+include	<smw.h>
+
+# Band structure
+define	LEN_BAND	9			# length of structure
+define	BAND_ID		Memi[$1]		# ptr to band id string
+define	BAND_FILTER	Memi[$1+1]		# ptr to filter string
+define	BAND_WC		Memd[P2D($1+2)]		# center wavelength
+define	BAND_DW		Memd[P2D($1+4)]		# wavelength width
+define	BAND_FN		Memi[$1+6]		# no. of filter points
+define	BAND_FW		Memi[$1+7]		# ptr to filter wavelengths
+define	BAND_FR		Memi[$1+8]		# ptr to filter responses
+
+# Multiple bands for indices and equivalent widths.
+define	NBANDS		3			# maximum number of bands
+define	BAND1		1
+define	BAND2		2
+define	BAND3		3
+define	BAND		Memi[$1+($2-1)*NBANDS+($3-1)]
+
+
+# T_SBANDS -- Compute band fluxes, indices, and equivalent widths.
+# A list of bandpasses is supplied in a text file, and all of them are applied
+# to each spectrum in the list.  The output is written to an output file
+# in multicolumn format.
+
+procedure t_sbands ()
+
+pointer	inlist			# Input list of spectra
+pointer	output			# Output file name
+pointer	fbands			# Band file name
+pointer	apertures		# Aperture list string
+bool	norm			# Normalize bands by response?
+bool	mag			# Output magnitudes instead of fluxes?
+double	magzero			# Magnitude zeropoint for magnitude output
+bool	verbose			# Verbose header?
+
+int	i, nbands, nsubbands, nimages, fd
+pointer	bands, aps, im, smw, sh
+pointer	sp, input
+
+int	open(), imtgetim()
+bool	clgetb(), rng_elementi()
+double	clgetd()
+pointer	imtopenp(), immap(), smw_openim(), rng_open()
+
+begin
+	call smark (sp)
+	call salloc (input, SZ_FNAME, TY_CHAR)
+	call salloc (output, SZ_FNAME, TY_CHAR)
+	call salloc (fbands, SZ_FNAME, TY_CHAR)
+	call salloc (apertures, SZ_LINE, TY_CHAR)
+
+	# Get task parameters.
+	inlist = imtopenp ("input")
+	call clgstr ("output", Memc[output], SZ_FNAME)
+	call clgstr ("bands", Memc[fbands], SZ_FNAME)
+	call clgstr ("apertures", Memc[apertures], SZ_LINE)
+	norm = clgetb ("normalize")
+	mag = clgetb ("mag")
+	magzero = clgetd ("magzero")
+	verbose = clgetb ("verbose")
+
+	# Read bands from the band file.
+	fd = open (Memc[fbands], READ_ONLY, TEXT_FILE)
+	call sb_bands (fd, bands, nbands, nsubbands)
+	call close (fd)
+
+	# Open the aperture list.
+	iferr (aps = rng_open (Memc[apertures], INDEF, INDEF, INDEF))
+	    call error (1, "Bad aperture list")
+
+	# Loop over the input spectra.
+	fd = 0
+	nimages = 0
+	while (imtgetim (inlist, Memc[input], SZ_FNAME) != EOF) {
+	    nimages = nimages + 1
+
+	    # Open the input image and get the WCS.
+	    iferr (im = immap (Memc[input], READ_ONLY, 0)) {
+		call erract (EA_WARN)
+		next
+	    }
+	    iferr (smw = smw_openim (im)) {
+		call imunmap (im)
+		call erract (EA_WARN)
+		next
+	    }
+
+	    # Open output file and write a verbose header if desired.
+	    # It is delayed until now to avoid output if an error occurs
+	    # such as image not found.
+
+	    if (nimages == 1) {
+		fd = open (Memc[output], APPEND, TEXT_FILE)
+		if (verbose)
+		    call sb_header (fd, norm, mag, magzero,
+			Memc[fbands], bands, nbands, nsubbands)
+	    }
+
+	    # Measure selected apertures.
+	    do i = 1, SMW_NSPEC(smw) {
+		call shdr_open (im, smw, i, 1, INDEFI, SHHDR, sh)
+		if (!rng_elementi (aps, AP(sh)))
+		    next
+		call shdr_open (im, smw, i, 1, INDEFI, SHDATA, sh)
+
+		call sb_proc (fd, sh, bands, nbands, norm, mag, magzero)
+	    }
+
+	    # Finish with image.
+	    call shdr_close (sh)
+	    call smw_close (smw)
+	    call imunmap (im)
+	}
+
+	# Finish up.
+	call sb_free (bands, nbands)
+	if (fd != 0)
+	    call close (fd)
+	call imtclose (inlist)
+	call sfree (sp)
+end
+
+
+# SB_BANDS - Read bands from the band file and put them into an array
+# of band pointers.
+
+procedure sb_bands (fd, bands, nbands, nsubbands)
+
+int	fd			#I Bandpass file descriptor
+pointer	bands			#O Bandpass table descriptor
+int	nbands			#O Number of bandpasses
+int	nsubbands		#O Number of individual bands
+
+bool	bandok
+int	ip
+double	center, width
+pointer	sp, line, id, filter
+
+int	getline(), ctowrd(), ctod()
+
+begin
+	call smark (sp)
+	call salloc (line, SZ_LINE, TY_CHAR)
+	call salloc (id, SZ_FNAME, TY_CHAR)
+	call salloc (filter, SZ_FNAME, TY_CHAR)
+
+	# Read the bands.  If the first band is not seen
+	# skip the line.  Check for 1, 2, or 3 bandpasses.
+	# Can't use fscan() because fscan() will be called later to
+	# read any filter file.
+
+	bands = NULL
+	nbands = 0
+	nsubbands = 0
+	while (getline (fd, Memc[line]) != EOF) {
+	    ip = 1
+	    bandok = (ctowrd (Memc[line], ip, Memc[id], SZ_FNAME) > 0)
+	    bandok = (bandok && ctod (Memc[line], ip, center) > 0)
+	    bandok = (bandok && ctod (Memc[line], ip, width) > 0)
+	    bandok = (bandok && ctowrd (Memc[line],ip,Memc[filter],SZ_FNAME)>0)
+	    if (!bandok || Memc[id] == '#')
+		next
+
+	    # Allocate and reallocate the array of band pointers.
+	    if (nbands == 0)
+		call malloc (bands, 10 * NBANDS, TY_POINTER)
+	    else if (mod (nbands, 10) == 0)
+		call realloc (bands, (nbands + 10) * NBANDS, TY_POINTER)
+	    nbands = nbands + 1
+
+	    call sb_alloc (BAND(bands,nbands,BAND1),
+		Memc[id], Memc[filter], center, width)
+	    nsubbands = nsubbands + 1
+
+	    bandok = (ctowrd (Memc[line], ip, Memc[id], SZ_FNAME) > 0)
+	    bandok = (bandok && ctod (Memc[line], ip, center) > 0)
+	    bandok = (bandok && ctod (Memc[line], ip, width) > 0)
+	    bandok = (bandok && ctowrd (Memc[line],ip,Memc[filter],SZ_FNAME)>0)
+	    if (bandok) {
+		call sb_alloc (BAND(bands,nbands,BAND2),
+		    Memc[id], Memc[filter], center, width)
+		nsubbands = nsubbands + 1
+	    } else
+		BAND(bands,nbands,BAND2) = NULL
+
+	    bandok = (ctowrd (Memc[line], ip, Memc[id], SZ_FNAME) > 0)
+	    bandok = (bandok && ctod (Memc[line], ip, center) > 0)
+	    bandok = (bandok && ctod (Memc[line], ip, width) > 0)
+	    bandok = (bandok && ctowrd (Memc[line],ip,Memc[filter],SZ_FNAME)>0)
+	    if (bandok) {
+		call sb_alloc (BAND(bands,nbands,BAND3),
+		    Memc[id], Memc[filter], center, width)
+		nsubbands = nsubbands + 1
+	    } else
+		BAND(bands,nbands,BAND3) = NULL
+	}
+
+	call sfree (sp)
+end
+
+
+# SB_ALLOC -- Allocate a band structure.
+
+procedure sb_alloc (band, id, filter, center, width)
+
+pointer	band			#O Band pointer
+char	id[ARB]			#I Band id
+char	filter[ARB]		#I Band filter
+double	center			#I Band wavelength
+double	width			#I Band width
+
+int	fn, fd, strlen(), open(), fscan(), nscan()
+double	w, r
+pointer	fw, fr
+bool	streq()
+errchk	open()
+
+begin
+	call calloc (band, LEN_BAND, TY_STRUCT)
+	call malloc (BAND_ID(band), strlen(id), TY_CHAR)
+	call malloc (BAND_FILTER(band), strlen(filter), TY_CHAR)
+	call strcpy (id, Memc[BAND_ID(band)], ARB)
+	call strcpy (filter, Memc[BAND_FILTER(band)], ARB)
+	BAND_WC(band) = center
+	BAND_DW(band) = width
+	BAND_FN(band) = 0
+	BAND_FW(band) = NULL
+	BAND_FR(band) = NULL
+
+	if (streq (filter, "none"))
+	    return
+
+	# Read the filter file.
+	fd = open (filter, READ_ONLY, TEXT_FILE)
+	fn = 0
+	while (fscan (fd) != EOF) {
+	    call gargd (w)
+	    call gargd (r)
+	    if (nscan() != 2)
+		next
+	    if (fn == 0) {
+		call malloc (fw, 100, TY_DOUBLE)
+		call malloc (fr, 100, TY_DOUBLE)
+	    } else if (mod (fn, 100) == 0) {
+		call realloc (fw, fn+100, TY_DOUBLE)
+		call realloc (fr, fn+100, TY_DOUBLE)
+	    }
+	    Memd[fw+fn] = w
+	    Memd[fr+fn] = r
+	    fn = fn + 1
+	}
+	call close (fd)
+
+	BAND_FN(band) = fn
+	BAND_FW(band) = fw
+	BAND_FR(band) = fr
+end
+
+
+# SB_FREE -- Free band structures.
+
+procedure sb_free (bands, nbands)
+
+pointer	bands			#I bands descriptor
+int	nbands			#I number of bands
+
+int	i, j
+pointer	band
+
+begin
+	do i = 1, nbands {
+	    do j = 1, NBANDS {
+		band = BAND(bands,i,j)
+		if (band != NULL) {
+		    call mfree (BAND_ID(band), TY_CHAR)
+		    call mfree (BAND_FILTER(band), TY_CHAR)
+		    call mfree (BAND_FW(band), TY_DOUBLE)
+		    call mfree (BAND_FR(band), TY_DOUBLE)
+		    call mfree (band, TY_STRUCT)
+		}
+	    }
+	}
+	call mfree (bands, TY_POINTER)
+end
+
+
+# SB_HEADER -- Print output header.
+
+procedure sb_header (fd, norm, mag, magzero, fbands, bands, nbands, nsubbands)
+
+pointer	fd			#I Output file descriptor
+bool	norm			#I Normalization flag
+bool	mag			#I Magnitude flag
+double	magzero			#I Magnitude zeropoint
+char	fbands[ARB]		#I Band file
+pointer	bands			#I Pointer to array of bands
+int	nbands			#I Number of bands
+int	nsubbands		#I Number of subbands
+
+int	i, j
+pointer	sp, str, band
+
+begin
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	# Output a banner and task parameters.
+	call sysid (Memc[str], SZ_LINE)
+	call fprintf (fd, "\n# SBANDS: %s\n#  ")
+	    call pargstr (Memc[str])
+	if (fbands[1] != EOS) {
+	    call fprintf (fd, " bands = %s,")
+		call pargstr (fbands)
+	}
+	call fprintf (fd, " norm = %b, mag = %b")
+	    call pargb (norm)
+	    call pargb (mag)
+	if (mag) {
+	    call fprintf (fd, ", magzero = %.2f")
+		call pargd (magzero)
+	    call strcpy ("mag", Memc[str], SZ_LINE)
+	} else
+	    call strcpy ("flux", Memc[str], SZ_LINE)
+
+	# Output the bands.
+	call fprintf (fd, "\n# %10s %10s %10s %10s\n")
+	    call pargstr ("band")
+	    call pargstr ("filter")
+	    call pargstr ("wavelength")
+	    call pargstr ("width")
+	do i = 1, nbands {
+	    do j = 1, NBANDS {
+		band = BAND(bands,i,j)
+		if (band == NULL)
+		    next
+		call fprintf (fd, "# %10s %10s %10g %10g\n")
+		    call pargstr (Memc[BAND_ID(band)])
+		    call pargstr (Memc[BAND_FILTER(band)])
+		    call pargd (BAND_WC(band))
+		    call pargd (BAND_DW(band))
+	    }
+	}
+
+	# Output column headings.
+	call fprintf (fd,
+	    "#\n# %24s %7.7s %11.11s")
+	    call pargstr ("spectrum")
+	    call pargstr ("band")
+	    call pargstr (Memc[str])
+	if (nsubbands > nbands) {
+	    call fprintf (fd, " %7.7s %11.11s %9.9s %9.9s")
+		call pargstr ("band")
+		call pargstr (Memc[str])
+		call pargstr ("index")
+		call pargstr ("eqwidth")
+	}
+	call fprintf (fd, "\n")
+
+	call sfree (sp)
+end
+
+
+# SB_PROC -- Measure the band fluxes and possibly a band index and eq. width.
+
+procedure sb_proc (fd, sh, bands, nbands, norm, mag, magzero)
+
+int	fd			#I Output file descriptor
+pointer	sh			#I Spectrum descriptor
+pointer	bands			#I Bandpass table pointer
+int	nbands			#I Number of bandpasses
+bool	norm			#I Normalize?
+bool	mag			#I Magnitude output?
+double	magzero			#I Magnitude zero point
+
+int	i
+double	flux, contval, index, eqwidth
+double	flux1, norm1, flux2, norm2, flux3, norm3, a, b
+pointer	sp, imname, band1, band2, band3
+
+begin
+	call smark (sp)
+	call salloc (imname, SZ_FNAME, TY_CHAR)
+
+	call sprintf (Memc[imname], SZ_FNAME, "%s%s(%d)")
+	    call pargstr (IMNAME(sh))
+	    call pargstr (IMSEC(sh))
+	    call pargi (AP(sh))
+
+	# Loop over all bandpasses
+	do i = 1, nbands {
+	    # Measure primary band flux, normalize, and print result.
+	    band1 = BAND(bands,i,BAND1)
+	    call sb_flux (sh, band1, flux1, norm1)
+	    if (IS_INDEFD(flux1))
+		next
+
+	    if (norm) {
+		flux1 = flux1 / norm1
+		norm1 = 1
+	    }
+	    if (mag && flux1 > 0.)
+		flux = magzero - 2.5 * log10 (flux1)
+	    else
+		flux = flux1
+
+	    call fprintf (fd, "%26s %7.7s %11.6g")
+		call pargstr (Memc[imname])
+		call pargstr (Memc[BAND_ID(band1)])
+		call pargd (flux)
+		
+	    # Measure the alternate band fluxes and compute and output
+	    # the band index and equivalent width.
+
+	    band2 = BAND(bands,i,BAND2)
+	    band3 = BAND(bands,i,BAND3)
+	    call sb_flux (sh, band2, flux2, norm2)
+	    call sb_flux (sh, band3, flux3, norm3)
+	    if (IS_INDEFD(flux2) && IS_INDEFD(flux3)) {
+		call fprintf (fd, "\n")
+		next
+	    }
+
+	    if (norm) {
+		if (!IS_INDEFD(flux2)) {
+		    flux2 = flux2 / norm2
+		    norm2 = 1
+		}
+		if (!IS_INDEFD(flux3)) {
+		    flux3 = flux3 / norm3
+		    norm3 = 1
+		}
+	    }
+
+	    contval = INDEFD
+	    index = INDEFD
+	    eqwidth = INDEFD
+	    if (!IS_INDEFD(flux2) && !IS_INDEFD(flux3)) {
+		# Interpolate to the center of the primary band.
+		a = (flux2 / norm2 - flux3 / norm3) /
+		    (BAND_WC(band2) - BAND_WC(band3))
+		b = flux2 / norm2 - a * BAND_WC(band2)
+		contval = (a * BAND_WC(band1) + b) * norm1
+		call fprintf (fd, " %7.7s")
+		    call pargstr ("cont")
+	    } else if (!IS_INDEFD(flux2)) {
+		contval = flux2
+		call fprintf (fd, " %7.7s")
+		    call pargstr (Memc[BAND_ID(band2)])
+	    } else if (!IS_INDEFD(flux3)) {
+		contval = flux3
+		call fprintf (fd, " %7.7s")
+		    call pargstr (Memc[BAND_ID(band3)])
+	    }
+
+	    if (mag && contval > 0.)
+		flux = magzero - 2.5 * log10 (contval)
+	    else
+		flux = contval
+	    call fprintf (fd, " %11.6g")
+		call pargd (flux)
+
+	    if (flux1 > 0. && contval > 0.) {
+		index = flux1 / contval
+		eqwidth = (1 - index) * BAND_DW(band1)
+	    }
+	    if (mag) {
+		if (!IS_INDEFD(contval) && contval > 0.)
+		    contval = magzero - 2.5 * log10 (contval)
+		if (!IS_INDEFD(index))
+		    index = -2.5 * log10 (index)
+	    }
+
+	    call fprintf (fd, " %9.6g %9.6g\n")
+		call pargd (index)
+		call pargd (eqwidth)
+	}
+
+	# Flush output and finish up.
+	call flush (fd)
+	call sfree (sp)
+end
+
+
+# SB_FLUX - Compute the flux and total response in a given band.
+# Return INDEF if the band is outside of the spectrum.
+
+procedure sb_flux (sh, band, flux, norm)
+
+pointer	sh			#I spectrum descriptor
+pointer	band			#I band descriptor
+double	flux			#O flux
+double	norm			#O normalization
+
+int	i, i1, i2
+double	a, b, w1, w2, x1, x2, wt
+pointer	x, y
+double	sb_filter(), shdr_wl()
+
+begin
+	# Return if no band is defined.
+	flux = INDEFD
+	norm = 1
+	if (band == NULL)
+	    return
+
+	# Check end points.
+	a = BAND_WC(band) - BAND_DW(band) / 2.
+	b = BAND_WC(band) + BAND_DW(band) / 2.
+	w1 = min (a, b)
+	w2 = max (a, b)
+	a = shdr_wl (sh, w1)
+	b = shdr_wl (sh, w2)
+	x1 = min (a, b)
+	x2 = max (a, b)
+	i1 = nint (x1)
+	i2 = nint (x2)
+	if (x1 == x2 || i1 < 1 || i2 > SN(sh))
+	    return
+
+	x = SX(sh) + i1 - 1
+	y = SY(sh) + i1 - 1
+
+	if (i1 == i2) {
+	    wt = sb_filter (double(Memr[x]), band) * (x2 - x1)
+	    flux = wt * Memr[y]
+	    norm = wt
+	} else {
+	    wt = sb_filter (double(Memr[x]), band) * (i1 + 0.5 - x1)
+	    flux = wt * Memr[y]
+	    norm = wt
+	    x = x + 1
+	    y = y + 1
+	    for (i = i1 + 1; i <= i2 - 1; i = i + 1) {
+		wt = sb_filter (double(Memr[x]), band)
+		flux = flux + wt * Memr[y]
+		norm = norm + wt
+		x = x + 1
+		y = y + 1
+	    }
+	    wt = sb_filter (double(Memr[x]), band) * (x2 - i2 + 0.5)
+	    flux = flux + wt * Memr[y]
+	    norm = norm + wt
+	}
+end
+
+
+# SB_FILTER -- Given a filter array interpolate to the specified wavelength.
+
+double procedure sb_filter (w, band)
+
+double	w		# Wavelength desired
+pointer	band		# Band pointer
+
+int	i, n
+double	x1, x2
+pointer	x, y
+
+begin
+	n = BAND_FN(band)
+	if (n == 0)
+	    return (1.)
+
+	x = BAND_FW(band)
+	y = BAND_FR(band)
+	x1 = Memd[x]
+	x2 = Memd[x+n-1]
+
+	if (w <= x1)
+	    return (Memd[y])
+	else if (w >= x2)
+	    return (Memd[y+n-1])
+	
+	if ((w - x1) < (x2 - w))
+	    for (i = 1; w > Memd[x+i]; i=i+1)
+		;
+	else
+	    for (i = n - 1; w < Memd[x+i-1]; i=i-1)
+		;
+		
+	x1 = Memd[x+i-1]
+	x2 = Memd[x+i]
+	return ((w - x1) / (x2 - x1) * (Memd[y+i] - Memd[y+i-1]) + Memd[y+i-1])
+end