Repatch (from linux) of OSX IRAF

32 files changed, 3675 insertions, 0 deletions

diff --git a/noao/onedspec/sensfunc/mkpkg b/noao/onedspec/sensfunc/mkpkg
new file mode 100644
index 00000000..4fdd11c0
--- /dev/null
+++ b/noao/onedspec/sensfunc/mkpkg
@@ -0,0 +1,38 @@
+# SENSFUNC
+
+$checkout libpkg.a ../
+$update   libpkg.a
+$checkin  libpkg.a ../
+$exit
+
+libpkg.a:
+	sfadd.x		sensfunc.h <gset.h>
+	sfapertures.x	sensfunc.h
+	sfcgraph.x	sensfunc.h <gset.h>
+	sfcolon.x	sensfunc.h <error.h> <gset.h>
+	sfcolors.x	sensfunc.h <gset.h>
+	sfcomposite.x	sensfunc.h
+	sfdata.x	sensfunc.h
+	sfdelete.x	sensfunc.h <gset.h>
+	sfeout.x	<ctype.h> <error.h> <mach.h> <math/curfit.h>
+	sfextinct.x	sensfunc.h <pkg/gtools.h>
+	sffit.x		sensfunc.h <math/curfit.h>
+	sfginit.x	sensfunc.h <gset.h>
+	sfgraph.x	sensfunc.h <error.h> <gset.h> <math/curfit.h>
+	sfimage.x	sensfunc.h <smw.h> <gset.h> <math/curfit.h>
+	sfmarks.x	sensfunc.h <gset.h>
+	sfmove.x	sensfunc.h <gset.h>
+	sfnearest.x	sensfunc.h <gset.h> <mach.h>
+	sfoutput.x	sensfunc.h <imhdr.h> <mach.h>
+	sfreset.x	sensfunc.h
+	sfrms.x		sensfunc.h
+	sfsensfunc.x	sensfunc.h <error.h> <gset.h> <mach.h>
+	sfshift.x	sensfunc.h
+	sfstats.x	sensfunc.h
+	sfstds.x	sensfunc.h
+	sftitle.x	sensfunc.h
+	sfundelete.x	sensfunc.h <gset.h>
+	sfvstats.x	sensfunc.h
+	sfweights.x	sensfunc.h
+	t_sensfunc.x	sensfunc.h
+	;
diff --git a/noao/onedspec/sensfunc/sensfunc.h b/noao/onedspec/sensfunc/sensfunc.h
new file mode 100644
index 00000000..9e3afb77
--- /dev/null
+++ b/noao/onedspec/sensfunc/sensfunc.h
@@ -0,0 +1,64 @@
+# SENSFUNC definitions.
+
+define	SF_NGRAPHS	4		# Number of graphs per frame
+define	SF_INCLUDE	1		# Include observation
+define	SF_EXCLUDE	2		# Exclude observation
+define	SF_DELETE	3		# Delete observation
+
+# SENSFUNC Standard Star Data Structure.
+
+define	SZ_STDIMAGE	63		# Length of standard image name
+define	SZ_STDTITLE	63		# Length of standard title
+
+define	LEN_STD		115		# Length of standard obs. structure
+
+define	STD_IMAGE	Memc[P2C($1)]	# Standard image name
+define	STD_SKY		Memc[P2C($1+32)] # Standard image sky
+define	STD_TITLE	Memc[P2C($1+64)] # Standard title
+define	STD_FLAG	Memi[$1+96]	# Flag
+define	STD_BEAM	Memi[$1+97]	# Beam number of spectrum
+define	STD_NPTS	Memi[$1+98]	# Number of points in spectrum
+define	STD_EXPTIME	Memr[P2R($1+99)]  # Exposure time
+define	STD_AIRMASS	Memr[P2R($1+100)] # Airmass of spectrum
+define	STD_WSTART	Memr[P2R($1+101)] # Starting wavelength of spectrum
+define	STD_WEND	Memr[P2R($1+102)] # Ending wavelength of spectrum
+define	STD_SHIFT	Memr[P2R($1+103)] # Added shift
+define	STD_NWAVES	Memi[$1+104]	# Number of calibration wavelengths
+define	STD_WAVES	Memi[$1+105]	# Pointer to wavelengths
+define	STD_FLUXES	Memi[$1+106]	# Pointer to standard flux values
+define	STD_DWAVES	Memi[$1+107]	# Pointer to flux bandwidths
+define	STD_COUNTS	Memi[$1+108]	# Pointer to counts
+define	STD_SENS	Memi[$1+109]	# Pointer to sensitivities
+define	STD_FIT		Memi[$1+110]	# Pointer to fitted sensitivities
+define	STD_WTS		Memi[$1+111]	# Pointer to weights
+define	STD_IWTS	Memi[$1+112]	# Pointer to weights
+define	STD_X		Memi[$1+114]	# Pointer to plotted x values
+define	STD_Y		Memi[$1+115]	# Pointer to plotted y values
+
+# Graphics structure
+
+define	GP_SZTITLE	79			# Size of title string
+
+define	LEN_GP		75			# Length of structure
+
+define	GP_GIO		Memi[$1]		# GIO pointer
+define	GP_TITLE	Memc[P2C($1+1)]		# General title
+define	GP_GRAPHS	Memc[P2C($1+41)+$2-1]	# Graphs
+define	GP_IMAGES	Memi[$1+44+$2-1]	# Pointer to image names
+define	GP_SKYS		Memi[$1+48+$2-1]	# Pointer to sky names
+define	GP_MARK		Memi[$1+52]		# Mark type
+define	GP_SZMARK	Memr[P2R($1+53)]	# Mark size
+define	GP_CMARK	Memi[$1+54]		# Mark color
+define	GP_MDEL		Memi[$1+55]		# Deleted mark
+define	GP_SZMDEL	Memr[P2R($1+56)]	# Size of deleted mark
+define	GP_CDEL		Memi[$1+57]		# Color of deleted mark
+define	GP_MADD		Memi[$1+58]		# Mark type
+define	GP_CADD		Memi[$1+59]		# Mark color
+define	GP_PLCOLOR	Memi[$1+60]		# Line color
+define  GP_WSTART       Memr[P2R($1+61)]        # Starting wavelength for plots
+define  GP_WEND         Memr[P2R($1+62)]        # Ending wavelength for plots
+define	GP_LOG		Memi[$1+63]		# Log flux plots?
+define	GP_FMIN		Memr[P2R($1+64)]	# Minimum flux plot limit
+define	GP_FMAX		Memr[P2R($1+65)]	# Maximum flux plot limit
+define	GP_SHDR		Memi[$1+65+$2]		# SHDR pointer
+define	GP_AIRMASS	Memr[P2R($1+69+$2)]	# Airmass range of plots
diff --git a/noao/onedspec/sensfunc/sensfunc.key b/noao/onedspec/sensfunc/sensfunc.key
new file mode 100644
index 00000000..5f09739f
--- /dev/null
+++ b/noao/onedspec/sensfunc/sensfunc.key
@@ -0,0 +1,81 @@
+		SENSFUNC: Determine Sensitivity Function
+
+SUMMARY:
+
+? Help            a Add data         c Composite data      d Delete data
+e Extinction      f Fit (overplot)   g Fit (redraw)        i Info 
+m Move data       o Original data    q Quit                r Redraw
+s Shift data      u Undelete data    w Change weights      I Interrupt
+
+:function [type]    :graphs [types]     :images [images]   :marks [types]
+:order [value]      :skys [images]      :stats [file]      :vstats [file]
+:colors [colors]
+
+Graph types: a=(resid,airmass), c=(composite,lambda), e=(extinction,lambda)
+	     i=(Fluxed image,lambda), l=(Log of fluxed image, lambda),
+	     r=(resid, lambda), s=(Sensitivity,lambda)
+
+
+CURSOR KEYS:
+
+?	Print help
+a	Add a point at the cursor position
+c	Toggle composite points
+d	Delete point, star, or wavelength nearest the cursor
+e	Toggle residual extinction correction
+f	Fit data with a sensitivity function and overplot the fit
+g	Fit data with a sensitivity function and redraw the graphs
+i	Print information about point nearest the cursor
+m	Move point, star, wavelength nearest the cursor to new sensitivity
+o	Reset to original data
+q	Quit and write sensitivity function for current aperture
+r	Redraw graph(s)
+s	Toggle shift of standard stars to eliminate mean deviations
+u	Undelete point, star, or wavelength nearest the cursor
+w	Change weights of point, star, or wavelength nearest the cursor
+I	Interrupt task immediately
+
+
+COLON COMMANDS AND ARGUMENTS:
+
+:flux [min] [max]  Limits for flux calibrated graphs (INDEF for autoscale)
+:function [type]   Function to be fit to sensitivity data.  The types are:
+			chebyshev - Chebyshev polynomial
+			legendre  - Legendre polynomial
+			spline1   - Linear spline
+			spline3   - Cubic spline
+:graphs [types]    Graphs to be displayed (up to four).  The types are:
+			a - Residual sensitivity vs airmass
+			c - Composite residuals and error bars vs wavelength
+			e - Extinction (and revised extinction) vs wavelength
+			i - Flux calibrated image vs wavelength
+			l - Log of flux calibrated image vs wavelength
+			r - Residual sensitivity vs wavelength
+			s - Sensitivity vs wavelength
+:images [images]   Images to flux calibrate and plot (up to four images)
+:colors [colors]   Line and mark colors to use for line and included, deleted,
+			and added points.  The colors are specified as four
+			whitespace separated integers between 1 and 9.
+:marks [marks]     Mark types to use for included, deleted, and added points:
+			point, box, plus, cross, diamond, hline, vline,
+			hebar, vebar, circle
+:order [order]     Order of function (polynomial terms or spline pieces)
+:skys [images]     Sky images for flux calibration (optional, up to four images)
+:stats [file]      Statistics about standard stars and sensitivity fit
+:vstats [file]     Verbose statistics about standard stars and sensitivity fit
+
+
+EXAMPLES:
+
+:func spline3		Select cubic spline function
+:g srae			Graph sensitivity, residuals, airmass, and extinction
+:g sii			Graph sensitivity and two images
+:i n1.0004 n1.0008	Set first two images to graph (the defaults are
+			    taken from the standard star list)
+:skys n1.0005		Subtract this sky image from first image for calibration
+:colors 2		Change only the line color to 2
+:colors	2 5 4 3		Change the line and mark colors
+:m plus			Change the mark type for included points and don't
+			    change the deleted or added point mark type
+:stats			Print statistics to terminal
+:vstats stdstats	Print verbose statistics to file
diff --git a/noao/onedspec/sensfunc/sfadd.x b/noao/onedspec/sensfunc/sfadd.x
new file mode 100644
index 00000000..c77694d7
--- /dev/null
+++ b/noao/onedspec/sensfunc/sfadd.x
@@ -0,0 +1,105 @@
+include	<gset.h>
+include	"sensfunc.h"
+
+# SF_ADD -- Add a point to the added point observation structure.
+# The added star is the next to last of the standard stars.
+
+procedure sf_add (gp, stds, nstds, cv, wx, wy, wc)
+
+pointer	gp			# Graphics structure
+pointer	stds[nstds]		# Standard star structures
+int	nstds			# Number of standard stars
+pointer	cv			# Sensitivity function curve
+real	wx			# Cursor X value
+real	wy			# Cursor Y value
+int	wc			# Cursor WCS
+
+int	nwaves
+real	wave, sen, fit, cveval()
+pointer	std, waves, sens, fits, wts, iwts, x, y
+errchk	malloc, realloc
+
+begin
+	# Convert from particular WCS to wavelength and sensitivity.  In
+	# order to add a point the graph must be either sensitivity or
+	# residual verses wavelength.  If not then return without adding
+	# a point.
+
+	switch (GP_GRAPHS(gp,wc)) {
+	case 's':
+	    wave = wx
+	    fit = cveval (cv, wx)
+	    sen = wy
+	case 'r':
+	    wave = wx
+	    fit = cveval (cv, wx)
+	    sen = wy + fit
+	default:
+	    return
+	}
+
+	# Add the point to the next to last standard star.  Allocate
+	# or reallocate memory as needed.  Turn the added star on by
+	# setting INCLUDE flag.
+
+	std = stds[nstds-1]
+	nwaves = STD_NWAVES(std) + 1
+	waves = STD_WAVES(std)
+	if (waves == NULL) {
+	    call malloc (waves, nwaves, TY_REAL)
+	    call malloc (sens, nwaves, TY_REAL)
+	    call malloc (fits, nwaves, TY_REAL)
+	    call malloc (wts, nwaves, TY_REAL)
+	    call malloc (iwts, nwaves, TY_REAL)
+	    call malloc (x, nwaves, TY_REAL)
+	    call malloc (y, nwaves, TY_REAL)
+	} else {
+	    sens = STD_SENS(std)
+	    fits = STD_FIT(std)
+	    wts = STD_WTS(std)
+	    iwts = STD_IWTS(std)
+	    x = STD_X(std)
+	    y = STD_Y(std)
+	    call realloc (waves, nwaves, TY_REAL)
+	    call realloc (sens, nwaves, TY_REAL)
+	    call realloc (fits, nwaves, TY_REAL)
+	    call realloc (wts, nwaves, TY_REAL)
+	    call realloc (iwts, nwaves, TY_REAL)
+	    call realloc (x, nwaves, TY_REAL)
+	    call realloc (y, nwaves, TY_REAL)
+	}
+	STD_FLAG(std) = SF_INCLUDE
+	STD_NWAVES(std) = nwaves
+	STD_WAVES(std) = waves
+	STD_SENS(std) = sens
+	STD_FIT(std) = fits
+	STD_WTS(std) = wts
+	STD_IWTS(std) = iwts
+	STD_X(std) = x
+	STD_Y(std) = y
+
+	Memr[waves+nwaves-1] = wave
+	Memr[sens+nwaves-1] = sen
+	Memr[fits+nwaves-1] = fit
+	Memr[wts+nwaves-1] = 1
+	Memr[iwts+nwaves-1] = 1
+
+	# Mark the added point on all graphs.
+	for (wc = 1; GP_GRAPHS(gp,wc) != EOS; wc=wc+1) {
+	    call gseti (GP_GIO(gp), G_WCS, wc)
+	    call gseti (GP_GIO(gp), G_PLCOLOR, GP_CADD(gp))
+	    switch (GP_GRAPHS(gp,wc)) {
+	    case 's':
+ 		call gmark (GP_GIO(gp), wave, sen, GP_MADD(gp), GP_SZMARK(gp),
+		    GP_SZMARK(gp))
+	    case 'r':
+		wy = sen - cveval (cv, wave)
+ 		call gmark (GP_GIO(gp), wave, wy, GP_MADD(gp), GP_SZMARK(gp),
+		    GP_SZMARK(gp))
+	    case 'a':
+		wy = sen - cveval (cv, wave)
+ 		call gmark (GP_GIO(gp), STD_AIRMASS(std), wy, GP_MADD(gp),
+		    GP_SZMARK(gp), GP_SZMARK(gp))
+	    }
+	}
+end
diff --git a/noao/onedspec/sensfunc/sfapertures.x b/noao/onedspec/sensfunc/sfapertures.x
new file mode 100644
index 00000000..7eb2b6f8
--- /dev/null
+++ b/noao/onedspec/sensfunc/sfapertures.x
@@ -0,0 +1,27 @@
+include "sensfunc.h"
+
+# SF_APERTURES -- Determine the apertures in use.
+
+procedure sf_apertures (stds, nstds, apertures, napertures)
+
+pointer	stds[nstds]		# Standard star data
+int	nstds			# Number of standard stars
+pointer	apertures		# Pointer to apertures (returned)
+int	napertures		# Number of apertures (returned)
+
+int	i, j, aperture
+
+errchk	malloc, realloc
+
+begin
+	call malloc (apertures, nstds, TY_INT)
+	napertures = 0
+	do i = 1, nstds {
+	    aperture = STD_BEAM(stds[i])
+	    for (j=1; (j<=napertures)&&(aperture!=Memi[apertures+j-1]); j=j+1)
+		;
+	    napertures = max (napertures, j)
+	    Memi[apertures+j-1] = aperture
+	}
+	call realloc (apertures, napertures, TY_INT)
+end
diff --git a/noao/onedspec/sensfunc/sfcgraph.x b/noao/onedspec/sensfunc/sfcgraph.x
new file mode 100644
index 00000000..2f689d47
--- /dev/null
+++ b/noao/onedspec/sensfunc/sfcgraph.x
@@ -0,0 +1,104 @@
+include	<gset.h>
+include	"sensfunc.h"
+
+# SF_CGRAPH -- Graph of composite points and errors
+
+procedure sf_cgraph (gp, stds, nstds, cv)
+
+pointer	gp		# Graphics structure
+pointer	stds[nstds]	# Standard star data
+int	nstds		# Number of standard stars
+pointer	cv		# Sensitivity function curve
+
+int	i, j, n, nwaves
+real	w, s, ymin, ymax, cveval()
+double	sum, sum2
+pointer	sp, waves, sens, errors, xp, yp, zp, gio
+
+begin
+	nwaves = 0
+	do i = 1, nstds-2
+	    if (STD_FLAG(stds[i]) != SF_EXCLUDE)
+	        nwaves = nwaves + STD_NWAVES(stds[i])
+
+	call smark (sp)
+	call salloc (waves, nwaves, TY_REAL)
+	call salloc (sens, nwaves, TY_REAL)
+	call salloc (errors, nwaves, TY_REAL)
+	    
+	nwaves = 0
+	do i = 1, nstds-2 {
+	    if (STD_FLAG(stds[i]) == SF_EXCLUDE)
+		next
+	    n = STD_NWAVES(stds[i])
+	    xp = STD_WAVES(stds[i])
+	    yp = STD_SENS(stds[i])
+	    zp = STD_WTS(stds[i])
+	    do j = 1, n {
+		if (Memr[zp] != 0.) {
+		    Memr[waves+nwaves] = Memr[xp]
+		    Memr[sens+nwaves] = Memr[yp]
+		    nwaves = nwaves + 1
+		}
+		xp = xp + 1
+		yp = yp + 1
+		zp = zp + 1
+	    }
+	}
+	call xt_sort2 (Memr[waves], Memr[sens], nwaves)
+
+	n = 0
+	sum = 0.
+	sum2 = 0.
+	ymin = 0.
+	ymax = 0.
+	j = 0
+	do i = 1, nwaves {
+	    w = Memr[waves+i-1]
+	    s = Memr[sens+i-1]
+	    n = n + 1
+	    sum = sum + s
+	    sum2 = sum2 + s * s
+
+	    if ((i < nwaves) && (w == Memr[waves+i]))
+		next
+
+	    if (n > 2) {
+	        sum = sum / n
+	        sum2 = sum2 / n - sum * sum
+	        if (sum2 > 0)
+	            sum2 = sqrt (sum2 / n)
+	        else
+		    sum2 = 0.
+	        sum = sum - cveval (cv, w)
+
+		Memr[waves+j] = w
+		Memr[sens+j] = sum
+		Memr[errors+j] = sum2
+		j = j + 1
+
+		if (sum + sum2 > ymax)
+		    ymax = sum + sum2
+		if (sum - sum2 < ymin)
+		    ymin = sum - sum2
+	    }
+	    n = 0
+	    sum = 0.
+	    sum2 = 0.
+	}
+	nwaves = j
+
+	if (j == 0) {
+	    call printf ("No wavelength overlap for composite points")
+	} else {
+	    gio = GP_GIO(gp)
+	    call gswind (gio, GP_WSTART(gp), GP_WEND(gp), ymin, ymax)
+	    call glabax (gio, "Composite Points vs Wavelength", "", "")
+	    call gseti (gio, G_PLCOLOR, GP_CMARK(gp))
+	    do i = 0, nwaves-1
+	        call gmark (gio, Memr[waves+i], Memr[sens+i], GM_VEBAR,
+		    1., -Memr[errors+i])
+	}
+
+	call sfree (sp)
+end
diff --git a/noao/onedspec/sensfunc/sfcolon.x b/noao/onedspec/sensfunc/sfcolon.x
new file mode 100644
index 00000000..43a056e7
--- /dev/null
+++ b/noao/onedspec/sensfunc/sfcolon.x
@@ -0,0 +1,193 @@
+include	<error.h>
+include	<gset.h>
+include	"sensfunc.h"
+
+# SENSFUNC colon commands
+define	CMDS	"|stats|vstats|function|order|graphs|images|skys|marks\
+		 |fluxlimits|colors|"
+define	STATS		1	# Show results
+define	VSTATS		2	# Show verbose results
+define	FUNCTION	3	# Sensitivity function type
+define	ORDER		4	# Function order
+define	GRAPHS		5	# Select graphs
+define	IMAGES		6	# Select images
+define	SKYS		7	# Select skys
+define	MARKS		8	# Set graph mark types
+define	FLIMITS		9	# Flux graph limits
+define	COLORS		10	# Flux graph limits
+
+# SF_COLON -- Process SENSFUNC colon commands.
+# This procedure has so many arguments because of the STATS option.
+
+procedure sf_colon (cmd, gp, stds, nstds, cv, wextn, extn, nextn, ecv, function,
+    order, npts, rms, newfit, newgraph)
+
+char	cmd[ARB]		# Colon command
+pointer	gp			# Graphics structure
+pointer	stds[nstds]		# Standard star data
+int	nstds			# Number of standard stars
+pointer	cv			# Sensitivity function curve
+real	wextn[nextn]		# Extinction table wavelengths
+real	extn[nextn]		# Extinction table values
+int	nextn			# Number of extinction table values
+pointer	ecv			# Residual extinction curve
+char	function[ARB]		# Function type
+int	order			# Function order
+int	npts			# Number of points in fit
+real	rms			# RMS in fit
+int	newfit			# New function?
+int	newgraph		# New graphs?
+
+int	i, j, ncmd, ival, fd, nscan(), strdic(), open(), stridx()
+real	rval1, rval2
+bool	streq()
+pointer	sp, str
+errchk	open
+
+begin
+	# Match the command against a dictionary.
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+	call sscan (cmd)
+	call gargwrd (Memc[str], SZ_LINE)
+	ncmd = strdic (Memc[str], Memc[str], SZ_LINE, CMDS)
+
+	# Switch on the command.
+	switch (ncmd) {
+	case STATS:
+	    call gargwrd (Memc[str], SZ_LINE)
+	    iferr {
+		# If no argument write to STDOUT otherwise append to file.
+	        if (nscan() == 1) {
+		    call gdeactivate (GP_GIO(gp), AW_CLEAR)
+		    call sf_stats (STDOUT, stds, nstds, function, order, npts,
+			rms)
+		    call greactivate (GP_GIO(gp), AW_PAUSE)
+	        } else {
+		    fd = open (Memc[str], APPEND, TEXT_FILE)
+		    call sf_stats (fd, stds, nstds, function, order, npts, rms)
+		    call close (fd)
+	        }
+	    } then
+		call erract (EA_WARN)
+	case VSTATS:
+	    call gargwrd (Memc[str], SZ_LINE)
+	    iferr {
+	        if (nscan() == 1) {
+		    # If no argument page on STDOUT otherwise append to file.
+		    # A temp file is used in order to page output.
+
+		    call mktemp ("tmp$sf", Memc[str], SZ_LINE)
+		    fd = open (Memc[str], NEW_FILE, TEXT_FILE)
+		    call sf_stats (fd, stds, nstds, function, order, npts, rms)
+		    call sf_vstats (fd, stds, nstds, cv, wextn, extn, nextn,
+			ecv)
+		    call close (fd)
+		    call gpagefile (GP_GIO(gp), Memc[str], "sensfunc")
+		    call delete (Memc[str])
+	        } else {
+		    fd = open (Memc[str], APPEND, TEXT_FILE)
+		    call sf_stats (fd, stds, nstds, function, order, npts, rms)
+		    call sf_vstats (fd, stds, nstds, cv, wextn, extn, nextn,
+			ecv)
+		    call close (fd)
+		}
+	    } then
+		call erract (EA_WARN)
+	case FUNCTION:
+	    call gargwrd (Memc[str], SZ_LINE)
+	    if (nscan() == 2) {
+	        call strcpy (Memc[str], function, SZ_FNAME)
+		newfit = NO
+	    } else {
+		call printf ("function %s")
+		    call pargstr (function)
+	    }
+	case ORDER:
+	    call gargi (ival)
+	    if (nscan() == 2) {
+		order = ival
+		newfit = NO
+	    } else {
+		call printf ("order %d")
+		    call pargi (order)
+	    }
+	case GRAPHS:
+	    call gargstr (Memc[str], SZ_LINE)
+	    j = str
+	    for (i=str; Memc[i] != EOS; i=i+1) {
+		switch (Memc[i]) {
+		case 'a','c','e','i','l','r','s':
+		    Memc[j] = Memc[i]
+		    j = j + 1
+		}
+	    }
+	    Memc[j] = EOS
+	    if (Memc[str] != EOS) {
+		call strcpy (Memc[str], GP_GRAPHS(gp,1), SF_NGRAPHS)
+		newgraph = YES
+	    } else {
+		call printf ("graphs %s")
+		    call pargstr (GP_GRAPHS(gp,1))
+	    }
+	case IMAGES:
+	    # Note that changing the image automatically clears the sky.
+	    do i = 1, SF_NGRAPHS {
+	        call gargwrd (Memc[str], SZ_LINE)
+		if (nscan() == i + 1) {
+		    call strcpy (Memc[str], Memc[GP_IMAGES(gp,i)], SZ_FNAME)
+		    Memc[GP_SKYS(gp,i)] = EOS
+		    do j = 1, nstds
+	    		if (streq (Memc[str], STD_IMAGE(stds[j])))
+			    call strcpy (STD_SKY(stds[j]), Memc[GP_SKYS(gp,i)],
+				SZ_FNAME)
+		} else
+		    break
+	    }
+	    if (nscan() == 1) {
+	        call printf ("images %s %s %s %s")
+		    call pargstr (Memc[GP_IMAGES(gp,1)])
+		    call pargstr (Memc[GP_IMAGES(gp,2)])
+		    call pargstr (Memc[GP_IMAGES(gp,3)])
+		    call pargstr (Memc[GP_IMAGES(gp,4)])
+	    }
+	case SKYS:
+	    do i = 1, SF_NGRAPHS {
+	        call gargwrd (Memc[str], SZ_LINE)
+		if (nscan() == i + 1)
+		    call strcpy (Memc[str], Memc[GP_SKYS(gp,i)], SZ_FNAME)
+		else
+		    break
+	    }
+	    if (nscan() == 1) {
+	        call printf ("skys %s %s %s %s")
+		    call pargstr (Memc[GP_SKYS(gp,1)])
+		    call pargstr (Memc[GP_SKYS(gp,2)])
+		    call pargstr (Memc[GP_SKYS(gp,3)])
+		    call pargstr (Memc[GP_SKYS(gp,4)])
+	    }
+	case MARKS:
+	    call gargstr (Memc[str], SZ_LINE)
+	    call sf_marks (gp, Memc[str])
+	case FLIMITS:
+	    call gargr (rval1)
+	    call gargr (rval2)
+	    if (nscan() == 3) {
+		GP_FMIN(gp) = rval1
+		GP_FMAX(gp) = rval2
+		if (stridx (GP_GRAPHS(gp,1), "il") != 0)
+		    newgraph = YES
+	    } else {
+		call printf ("fluxlimits %g %g")
+		    call pargr (GP_FMIN(gp))
+		    call pargr (GP_FMAX(gp))
+	    }
+	case COLORS:
+	    call gargstr (Memc[str], SZ_LINE)
+	    call sf_colors (gp, Memc[str])
+	default:
+	    call printf ("Unrecognized or ambiguous command\007")
+	}
+
+	call sfree (sp)
+end
diff --git a/noao/onedspec/sensfunc/sfcolors.x b/noao/onedspec/sensfunc/sfcolors.x
new file mode 100644
index 00000000..db3f69af
--- /dev/null
+++ b/noao/onedspec/sensfunc/sfcolors.x
@@ -0,0 +1,28 @@
+include	<gset.h>
+include	"sensfunc.h"
+
+
+# SF_COLORS -- Set colors.
+
+procedure sf_colors (gp, colors)
+
+pointer	gp
+char	colors[ARB]
+
+int	i, nscan()
+
+begin
+	call sscan (colors)
+	call gargi (i)
+	if (nscan() == 1)
+	    GP_PLCOLOR(gp) = i
+	call gargi (i)
+	if (nscan() == 2)
+	    GP_CMARK(gp) = i
+	call gargi (i)
+	if (nscan() == 3)
+	    GP_CDEL(gp) = i
+	call gargi (i)
+	if (nscan() == 4)
+	    GP_CADD(gp) = i
+end
diff --git a/noao/onedspec/sensfunc/sfcomposite.x b/noao/onedspec/sensfunc/sfcomposite.x
new file mode 100644
index 00000000..506416e3
--- /dev/null
+++ b/noao/onedspec/sensfunc/sfcomposite.x
@@ -0,0 +1,147 @@
+include	"sensfunc.h"
+
+# SF_COMPOSITE -- Create a composite standard structure.
+# The composite star is the last of the standard stars.
+# When the composite star is created the other stars are turned off.
+# The function toggles.
+
+procedure sf_composite (stds, nstds, cv)
+
+pointer	stds[nstds]		# Standard star data
+int	nstds			# Number of standard stars
+pointer	cv			# Sensitivity pointer
+
+int	i, j, k, n, nwaves
+pointer	std, waves, sens, fit, wts, iwts, x, y, z
+errchk	malloc, realloc, xt_sort3
+
+begin
+	# If data is already composite toggle back to original data.
+	# Delete data points if composite point is deleted.
+	std = stds[nstds]
+	if (STD_FLAG(std) == SF_INCLUDE) {
+            do i = 1, nstds - 2 {
+	        if (STD_FLAG(stds[i]) == SF_EXCLUDE)
+		    next
+	        STD_FLAG(stds[i]) = SF_INCLUDE
+	    }
+	    STD_FLAG(std) = SF_EXCLUDE
+
+	    n = STD_NWAVES(std)
+	    x = STD_WAVES(std)
+	    z = STD_WTS(std)
+	    do i = 1, n {
+		if (Memr[z] == 0.) {
+		    do j = 1, nstds - 2 {
+		        if (STD_FLAG(stds[j]) == SF_EXCLUDE)
+			    next
+		        nwaves = STD_NWAVES(stds[j])
+			waves = STD_WAVES(stds[j])
+		        wts = STD_WTS(stds[j])
+			do k = 1, nwaves {
+			    if (Memr[waves] == Memr[x])
+		                Memr[wts] = 0.
+			    waves = waves + 1
+			    wts = wts + 1
+			}
+		    }
+		}
+		x = x + 1
+		z = z + 1
+	    }
+	    call printf ("Individual star data")
+	    return
+	}
+
+	# Initialize
+	if (STD_WAVES(std) != NULL) {
+	    call mfree (STD_WAVES(std), TY_REAL)
+	    call mfree (STD_SENS(std), TY_REAL)
+	    call mfree (STD_WTS(std), TY_REAL)
+	    call mfree (STD_IWTS(std), TY_REAL)
+	    call mfree (STD_X(std), TY_REAL)
+	    call mfree (STD_Y(std), TY_REAL)
+	}
+
+	# To bin the data we collect all the data and then sort by wavelength.
+	nwaves = 0
+	do i = 1, nstds - 2
+	    if (STD_FLAG(stds[i]) == SF_INCLUDE)
+	        nwaves = nwaves + STD_NWAVES(stds[i])
+
+	call malloc (waves, nwaves, TY_REAL)
+	call malloc (sens, nwaves, TY_REAL)
+	call malloc (wts, nwaves, TY_REAL)
+	    
+	nwaves = 0
+	do i = 1, nstds - 2 {
+	    if (STD_FLAG(stds[i]) != SF_INCLUDE)
+		next
+	    n = STD_NWAVES(stds[i])
+	    x = STD_WAVES(stds[i])
+	    y = STD_SENS(stds[i])
+	    z = STD_WTS(stds[i])
+	    do j = 1, n {
+		if (Memr[z] != 0.) {
+		    Memr[waves+nwaves] = Memr[x]
+		    Memr[sens+nwaves] = Memr[y]
+		    Memr[wts+nwaves] = Memr[z]
+		    nwaves = nwaves + 1
+		}
+		x = x + 1
+		y = y + 1
+		z = z + 1
+	    }
+	    STD_FLAG(stds[i]) = SF_DELETE
+	    STD_BEAM(std) = STD_BEAM(stds[i])
+	    STD_WSTART(std) = STD_WSTART(stds[i])
+	    STD_WEND(std) = STD_WEND(stds[i])
+	}
+#	STD_NWAVES(stds[nstds-1]) = 0
+
+	call xt_sort3 (Memr[waves], Memr[sens], Memr[wts], nwaves)
+
+	# Go through the wavelength sorted data and composite all points
+	# with the same wavelength.
+
+	n = 0
+	Memr[sens] = Memr[wts] * Memr[sens]
+	do i = 1, nwaves-1 {
+	    if (Memr[waves+i] == Memr[waves+n]) {
+		Memr[sens+n] = Memr[sens+n] + Memr[wts+i] * Memr[sens+i]
+		Memr[wts+n] = Memr[wts+n] + Memr[wts+i]
+	    } else {
+		n = n + 1
+		Memr[waves+n] = Memr[waves+i]
+		Memr[sens+n] = Memr[wts+i] * Memr[sens+i]
+		Memr[wts+n] = Memr[wts+i]
+	    }
+	}
+
+	nwaves = n + 1
+	do i = 0, nwaves-1
+	    Memr[sens+i] = Memr[sens+i] / Memr[wts+i]
+
+	# Store the composite data in the standard star structure.
+	call realloc (waves, nwaves, TY_REAL)
+	call realloc (sens, nwaves, TY_REAL)
+	call realloc (wts, nwaves, TY_REAL)
+	call malloc (iwts, nwaves, TY_REAL)
+	call malloc (fit, nwaves, TY_REAL)
+	call malloc (x, nwaves, TY_REAL)
+	call malloc (y, nwaves, TY_REAL)
+	call amovr (Memr[wts], Memr[iwts], nwaves)
+	call cvvector (cv, Memr[waves], Memr[fit], nwaves)
+
+	STD_FLAG(std) = SF_INCLUDE
+	STD_NWAVES(std) = nwaves
+	STD_WAVES(std) = waves
+	STD_SENS(std) = sens
+	STD_FIT(std) = fit
+	STD_WTS(std) = wts
+	STD_IWTS(std) = iwts
+	STD_X(std) = x
+	STD_Y(std) = y
+
+	call printf ("Composite star data")
+end
diff --git a/noao/onedspec/sensfunc/sfdata.x b/noao/onedspec/sensfunc/sfdata.x
new file mode 100644
index 00000000..94140049
--- /dev/null
+++ b/noao/onedspec/sensfunc/sfdata.x
@@ -0,0 +1,59 @@
+include	"sensfunc.h"
+
+# SF_DATA -- Compute the X and Y values for the particular graph.
+
+procedure sf_data (stds, nstds, graph)
+
+pointer	stds[nstds]		# Standard star structures
+int	nstds			# Number of standard stars
+char	graph			# Graph type
+
+real	a
+int	i, n
+pointer	wp, sp, fp, xp, yp
+
+begin
+	switch (graph) {
+	case 's':	# Sensitivity vs. Wavelength
+	    do i = 1, nstds {
+	        if (STD_FLAG(stds[i]) != SF_INCLUDE)
+		    next
+		n = STD_NWAVES(stds[i])
+		a = STD_AIRMASS(stds[i])
+		wp = STD_WAVES(stds[i])
+		sp = STD_SENS(stds[i])
+		xp = STD_X(stds[i])
+		yp = STD_Y(stds[i])
+		call amovr (Memr[wp], Memr[xp], n)
+		call amovr (Memr[sp], Memr[yp], n)
+	    }
+	case 'a':	# Residuals vs. Airmass
+	    do i = 1, nstds {
+	        if (STD_FLAG(stds[i]) != SF_INCLUDE)
+		    next
+		n = STD_NWAVES(stds[i])
+		a = STD_AIRMASS(stds[i])
+		wp = STD_WAVES(stds[i])
+		sp = STD_SENS(stds[i])
+		fp = STD_FIT(stds[i])
+		xp = STD_X(stds[i])
+		yp = STD_Y(stds[i])
+		call amovkr (a, Memr[xp], n)
+		call asubr (Memr[sp], Memr[fp], Memr[yp], n)
+	    }
+	case 'r':	# Residuals vs. Wavelength
+	    do i = 1, nstds {
+	        if (STD_FLAG(stds[i]) != SF_INCLUDE)
+		    next
+		n = STD_NWAVES(stds[i])
+		a = STD_AIRMASS(stds[i])
+		wp = STD_WAVES(stds[i])
+		sp = STD_SENS(stds[i])
+		fp = STD_FIT(stds[i])
+		xp = STD_X(stds[i])
+		yp = STD_Y(stds[i])
+		call amovr (Memr[wp], Memr[xp], n)
+		call asubr (Memr[sp], Memr[fp], Memr[yp], n)
+	    }
+	}
+end
diff --git a/noao/onedspec/sensfunc/sfdelete.x b/noao/onedspec/sensfunc/sfdelete.x
new file mode 100644
index 00000000..ff2d267c
--- /dev/null
+++ b/noao/onedspec/sensfunc/sfdelete.x
@@ -0,0 +1,127 @@
+include	<gset.h>
+include "sensfunc.h"
+
+# SF_DELETE -- Delete point, star, or wavelength identified by the
+# star index and index within the array of values.
+
+procedure sf_delete (gp, stds, nstds, key, istd, ipt)
+
+pointer	gp			# GIO pointer
+pointer	stds[nstds]		# Standard star data
+int	nstds			# Number of standard stars
+int	key			# Delete point, star, or wavelength
+int	istd			# Index of standard star
+int	ipt			# Index of point
+
+int	i, j, n, wcs, mark, mdel, cdel, stridx()
+real	wave, szmark, szmdel
+pointer	x, y, z, w, gio
+
+begin
+	gio = GP_GIO(gp)
+	mdel = GP_MDEL(gp)
+	cdel = GP_CDEL(gp)
+	szmdel = GP_SZMDEL(gp)
+	szmark = GP_SZMARK(gp)
+
+	# Delete the point or points from each displayed graph.
+	# When deleting multiple points check if point already deleted.
+	for (wcs = 1; GP_GRAPHS(gp,wcs) != EOS; wcs = wcs + 1) {
+	    if (stridx (GP_GRAPHS(gp,wcs), "ars") == 0)
+	        next
+
+	    call gseti (gio, G_WCS, wcs)
+	    call sf_data (stds, nstds, GP_GRAPHS(gp,wcs))
+	    switch (key) {
+	    case 'p':
+		if (istd != nstds-1)
+		    mark = GP_MARK(gp)
+		else
+		    mark = GP_MADD(gp)
+	        x = STD_X(stds[istd])+ipt-1
+	        y = STD_Y(stds[istd],1)+ipt-1
+		call gseti (gio, G_PMLTYPE, 0)
+	        call gmark (gio, Memr[x], Memr[y], mark, szmark, szmark)
+		call gseti (gio, G_PMLTYPE, 1)
+		call gseti (gio, G_PLCOLOR, cdel)
+	        call gmark (gio, Memr[x], Memr[y], mdel, szmdel, szmdel)
+	    case 's':
+		if (istd != nstds-1)
+		    mark = GP_MARK(gp)
+		else
+		    mark = GP_MADD(gp)
+	        n = STD_NWAVES(stds[istd])
+	        x = STD_X(stds[istd])
+	        y = STD_Y(stds[istd])
+	        w = STD_WTS(stds[istd])
+		do i = 1, n {
+		    if (Memr[w] != 0.) {
+		        call gseti (gio, G_PMLTYPE, 0)
+	                call gmark (gio, Memr[x], Memr[y], mark, szmark, szmark)
+		        call gseti (gio, G_PMLTYPE, 1)
+			call gseti (gio, G_PLCOLOR, cdel)
+	                call gmark (gio, Memr[x], Memr[y], mdel, szmdel, szmdel)
+		    }
+		    x = x + 1
+		    y = y + 1
+		    w = w + 1
+		}
+	    case 'w':
+	        wave = Memr[STD_WAVES(stds[istd])+ipt-1]
+	        do i = 1, nstds {
+	            if (STD_FLAG(stds[i]) != SF_INCLUDE)
+		        next
+		    if (i != nstds-1)
+			mark = GP_MARK(gp)
+		    else
+			mark = GP_MADD(gp)
+	            n = STD_NWAVES(stds[i])
+	            x = STD_X(stds[i])
+	            y = STD_Y(stds[i])
+		    z = STD_WAVES(stds[i])
+	            w = STD_WTS(stds[i])
+		    do j = 1, n {
+		        if ((Memr[z] == wave) && (Memr[w] != 0.)) {
+		            call gseti (gio, G_PMLTYPE, 0)
+	                    call gmark (gio, Memr[x], Memr[y], mark, szmark,
+				    szmark)
+		            call gseti (gio, G_PMLTYPE, 1)
+			    call gseti (gio, G_PLCOLOR, cdel)
+	                    call gmark (gio, Memr[x], Memr[y], mdel, szmdel,
+				szmdel)
+		        }
+			x = x + 1
+			y = y + 1
+			z = z + 1
+			w = w + 1
+		    }
+	        }
+	    }
+	}
+
+	# Mark the points as deleted by setting their weights to zero.
+	switch (key) {
+	case 'p':
+	    w = STD_WTS(stds[istd])+ipt-1
+	    Memr[w] = 0.
+	case 's':
+	    n = STD_NWAVES(stds[istd])
+	    w = STD_WTS(stds[istd])
+	    call aclrr (Memr[w], n)
+	case 'w':
+	    wave = Memr[STD_WAVES(stds[istd])+ipt-1]
+	    do i = 1, nstds {
+	        if (STD_FLAG(stds[i]) != SF_INCLUDE)
+		    next
+	        n = STD_NWAVES(stds[i])
+		z = STD_WAVES(stds[i])
+	        w = STD_WTS(stds[i])
+		do j = 1, n {
+		    if (Memr[z] == wave)
+	                Memr[w] = 0.
+		    w = w + 1
+		    z = z + 1
+		}
+	    }
+	}
+end
diff --git a/noao/onedspec/sensfunc/sfeout.x b/noao/onedspec/sensfunc/sfeout.x
new file mode 100644
index 00000000..8dae6301
--- /dev/null
+++ b/noao/onedspec/sensfunc/sfeout.x
@@ -0,0 +1,114 @@
+include	<error.h>
+include	<ctype.h>
+include	<mach.h>
+include	<math/curfit.h>
+
+define	NEPOINTS	100		# Number of points in extinction table
+
+# SF_EOUT -- Output a revised extinction table.  This is only done if there
+# is at least one residual extinction curve.  No assumption is made about
+# overlapping extinction curves.  In the overlap the extinction corrections
+# are averaged.
+
+procedure sf_eout (wextn, extn, nextn, ecvs, necvs)
+
+real	wextn[nextn]		# Standard extinction wavelengths
+real	extn[nextn]		# Standard extinction values
+int	nextn			# Number of standard extinction points
+pointer	ecvs[necvs]		# Residual extinction curves (one for each beam)
+int	necvs			# Number of residual extinction curves
+
+int	i, j, fd, open(), scan(), nscan()
+real	w, ext, wmin, wmax, dw, xmin, xmax, cvstatr(), cveval()
+pointer	sp, fname, waves, extns, navg, cv
+
+define	newfile_  99
+
+begin
+	# If there are no residual extinction values then do nothing.
+	for (i=1; (i<=necvs) && (ecvs[i]==NULL); i=i+1)
+	    ;
+	if (i > necvs)
+	    return
+
+	# The output table consists of NEPOINTS.
+	call smark (sp)
+	call salloc (fname, SZ_FNAME, TY_CHAR)
+	call salloc (waves, NEPOINTS, TY_REAL)
+	call salloc (extns, NEPOINTS, TY_REAL)
+	call salloc (navg, NEPOINTS, TY_INT)
+	call aclrr (Memr[extns], NEPOINTS)
+	call aclri (Memi[navg], NEPOINTS)
+
+	# Open the extinction table.  If it fails allow the user to
+	# enter a new name.
+
+	call clgstr ("newextinction", Memc[fname], SZ_FNAME)
+	for (i=fname; (Memc[i]!=EOS) && IS_WHITE(Memc[i]); i=i+1)
+	if (Memc[i] == EOS) {
+	    call sfree (sp)
+	    return
+	}
+
+newfile_
+	iferr (fd = open (Memc[i], NEW_FILE, TEXT_FILE)) {
+	    call printf ("Cannot create %s -- Enter new extinction filename: ")
+		call pargstr (Memc[fname])
+	    call flush (STDOUT)
+	    if (scan() != EOF) {
+	        call gargwrd (Memc[fname], SZ_FNAME)
+	        if (nscan() == 1)
+	    	    goto newfile_
+	    }
+	    call sfree (sp)
+	    call printf ("No revised extinction file created\n")
+	    return
+	}
+
+	# Determine the range of the extinction table.
+	wmin = MAX_REAL
+	wmax = -MAX_REAL
+	do i = 1, necvs {
+	    if (ecvs[i] == NULL)
+		next
+	    wmin = min (wmin, cvstatr (ecvs[i], CVXMIN))
+	    wmax = max (wmax, cvstatr (ecvs[i], CVXMAX))
+	}
+	dw = (wmax - wmin) / (NEPOINTS - 1)
+	do i = 1, NEPOINTS
+	    Memr[waves+i-1] = wmin + (i - 1) * dw
+
+	# Average the residual extinctions and add the original extinction.
+	do j = 1, necvs {
+	    if (ecvs[j] == NULL)
+		next
+	    cv = ecvs[j]
+	    xmin = cvstatr (cv, CVXMIN)
+	    xmax = cvstatr (cv, CVXMAX)
+	    do i = 1, NEPOINTS {
+		w = Memr[waves+i-1]
+		if ((w < xmin) || (w > xmax))
+		    next
+		Memr[extns+i-1] = Memr[extns+i-1] + cveval (cv, w)
+		Memi[navg+i-1] = Memi[navg+i-1] + 1
+	    }
+	}
+	do i = 1, NEPOINTS {
+	    if (Memi[navg+i-1] > 0)
+		Memr[extns+i-1] = Memr[extns+i-1] / Memi[navg+i-1]
+	    w = Memr[waves+i-1]
+	    call intrp (1, wextn, extn, nextn, w, ext, j)
+	    Memr[extns+i-1] = Memr[extns+i-1] + ext
+	}
+
+	# Output extinction table.
+	call fprintf (fd, "# Revised extinction table.\n")
+	do i = 1, NEPOINTS {
+	    call fprintf (fd, "%7.2f  %6.3f\n")
+		call pargr (Memr[waves+i-1])
+		call pargr (Memr[extns+i-1])
+	}
+	call close (fd)
+
+	call sfree (sp)
+end
diff --git a/noao/onedspec/sensfunc/sfextinct.x b/noao/onedspec/sensfunc/sfextinct.x
new file mode 100644
index 00000000..f7b95326
--- /dev/null
+++ b/noao/onedspec/sensfunc/sfextinct.x
@@ -0,0 +1,226 @@
+include	<pkg/gtools.h>
+include	"sensfunc.h"
+
+define	RANGE_AIRMASS	0.1	# Minimum airmass range
+define	SIGMA_AIRMASS	0.05	# Minimum sigma in airmass
+
+# SF_EXINCT -- Determine a residual extinction curve.  At each wavelength
+# for which there are multiple observations or neighboring wavelengths
+# such that there is a sufficient airmass range determine the slope
+# of the sensitivity vs airmass.  Residual sensitivity is used to minimize
+# wavelength scatter when multiple wavelengths are needed because of
+# nonoverlapping standard star data.  Each such slope is a measure of the
+# residual extinction at that wavelength.  To make the residual extinction
+# curve fit the extinction vs. wavelength using ICFIT.
+
+procedure sf_extinct (gp, stds, nstds, cv, ecv, function, order)
+
+pointer	gp			# Graphics structure
+pointer	stds[nstds]		# Standard star data
+int	nstds			# Number of standard stars
+pointer	cv			# Sensitivity function curve
+pointer	ecv			# Residual extinction curve
+char	function[ARB]		# Fitting function
+int	order			# Function order
+
+bool	ans
+int	i, j, n, nwaves, sum, npts, scan()
+real	a, amin, amax, rms, rms1, r2, sig, cveval()
+double	x, y, z, sumx, sumy, sumz, sumx2, sumxy
+pointer	sp, waves, sens, airm, xp, yp, fp, wp, ic
+pointer	gt, gt_init()
+errchk	salloc, xt_sort3, icg_fit, ic_open
+
+define	cancel_	99
+
+begin
+	# Cancel previous extinction if defined.
+	if (ecv != NULL)
+	    goto cancel_
+	
+	# Check for minimum airmass range and determine number of points.
+	# Ignore added objects and composite data.
+	amin = 100.
+	amax = 0.
+	nwaves = 0
+	do i = 1, nstds - 2 {
+	    if (STD_FLAG(stds[i]) != SF_INCLUDE)
+		next
+	    nwaves = nwaves + STD_NWAVES(stds[i])
+	    a = STD_AIRMASS(stds[i])
+	    amin = min (amin, a)
+	    amax = max (amax, a)
+	}
+	if (amax - amin < RANGE_AIRMASS) {
+	    call printf (
+		"Insufficient airmass coverage for extinction determination")
+	    return
+	}
+
+	# Extract data to be used and sort by wavelength.
+	# The data is wavelength, airmass, and residual sensitivity.
+	call smark (sp)
+	call salloc (waves, nwaves, TY_REAL)
+	call salloc (sens, nwaves, TY_REAL)
+	call salloc (airm, nwaves, TY_REAL)
+	    
+	nwaves = 0
+	do i = 1, nstds-2 {
+	    if (STD_FLAG(stds[i]) != SF_INCLUDE)
+		next
+	    n = STD_NWAVES(stds[i])
+	    a = STD_AIRMASS(stds[i])
+	    xp = STD_WAVES(stds[i])
+	    yp = STD_SENS(stds[i])
+	    fp = STD_FIT(stds[i])
+	    wp = STD_WTS(stds[i])
+	    do j = 1, n {
+		if (Memr[wp] != 0.) {
+		    Memr[airm+nwaves] = a
+		    Memr[waves+nwaves] = Memr[xp]
+		    Memr[sens+nwaves] = Memr[yp] - Memr[fp]
+		    nwaves = nwaves + 1
+		}
+		xp = xp + 1
+		yp = yp + 1
+		fp = fp + 1
+		wp = wp + 1
+	    }
+	}
+
+	call xt_sort3 (Memr[waves], Memr[sens], Memr[airm], nwaves)
+
+	# Bin points with common wavelengths or at least two points.
+	sum = 0
+	sumx = 0.
+	sumy = 0.
+	sumz = 0.
+	sumx2 = 0.
+	sumxy = 0.
+	n = 0
+	do i = 0, nwaves-2 {
+	    x = Memr[airm+i]
+	    y = Memr[sens+i]
+	    z = Memr[waves+i]
+	    sum = sum + 1
+	    sumx = sumx + x
+	    sumy = sumy + y
+	    sumx2 = sumx2 + x * x
+	    sumxy = sumxy + x * y
+	    sumz = sumz + z
+
+	    if ((z == Memr[waves+i+1]) || (sum < 2))
+		next
+
+	    x = sumx2 - sumx * sumx / sum
+	    if (x > SIGMA_AIRMASS) {
+	        Memr[waves+n] = sumz / sum
+	        Memr[sens+n] = (sumx * sumy / sum - sumxy) / x
+	        Memr[airm+n] = 1.
+	        n = n + 1
+		sum = 0
+		sumx = 0.
+		sumy = 0.
+		sumz = 0.
+		sumx2 = 0.
+		sumxy = 0.
+	    }
+	}
+	if (sum > 1) {
+	    x = sumx2 - sumx * sumx / sum
+	    if (x > SIGMA_AIRMASS) {
+	        Memr[waves+n] = sumz / sum
+	        Memr[sens+n] = (sumx * sumy / sum - sumxy) / x
+	        Memr[airm+n] = 1.
+		n = n + 1
+	    }
+	}
+
+	if (n < 2) {
+	    call printf ("Cannot determine residual extinction")
+	    call sfree (sp)
+	    return
+	}
+
+	# Fit residual extinction curve using ICFIT.
+	gt = gt_init()
+	call gt_sets (gt, GTTYPE, "mark")
+	call gt_seti (gt, GTCOLOR, GP_PLCOLOR(gp))
+	call ic_open (ic)
+	call ic_putr (ic, "xmin", min (STD_WSTART(stds[1]), STD_WEND(stds[1])))
+	call ic_putr (ic, "xmax", max (STD_WSTART(stds[1]), STD_WEND(stds[1])))
+	call ic_pstr (ic, "function", "chebyshev")
+	call ic_puti (ic, "order", 1)
+	call ic_pstr (ic, "xlabel", "wavelength")
+	call ic_pstr (ic, "ylabel", "residual extinction")
+	call ic_pstr (ic, "yunits", "mag")
+	call icg_fit (ic, GP_GIO(gp), "cursor", gt, ecv, Memr[waves],
+	    Memr[sens], Memr[airm], n)
+	call ic_closer (ic)
+	call gt_free (gt)
+
+	# Determine significance of the fit.
+	call sf_fit (stds, nstds, cv, function, order,
+	    min (GP_WSTART(gp), GP_WEND(gp)), max (GP_WSTART(gp), GP_WEND(gp)))
+	call sf_rms (stds, nstds, rms1, npts)
+	do i = 1, nstds - 2 {
+	    if (STD_FLAG(stds[i]) != SF_INCLUDE)
+		next
+	    n = STD_NWAVES(stds[i])
+	    a = STD_AIRMASS(stds[i])
+	    xp = STD_WAVES(stds[i])
+	    yp = STD_SENS(stds[i])
+	    call cvvector (ecv, Memr[xp], Memr[sens], n)
+	    call amulkr (Memr[sens], a, Memr[sens], n)
+	    call aaddr (Memr[yp], Memr[sens], Memr[yp], n)
+	}
+	call sf_fit (stds, nstds, cv, function, order,
+	    min (GP_WSTART(gp), GP_WEND(gp)), max (GP_WSTART(gp), GP_WEND(gp)))
+	call sf_rms (stds, nstds, rms, npts)
+	do i = 1, SF_NGRAPHS
+	    if (GP_SHDR(gp,i) != NULL)
+		call shdr_close (GP_SHDR(gp,i))
+
+	r2 = 1 - rms ** 2 / rms1 ** 2
+	sig = r2 * (nwaves - 2) / max (0.01, 1. - r2)
+	if (sig <= 0.0)
+	    sig = 0.
+	else
+	    sig = sqrt (sig)
+
+	# Apply to data if desired.
+	call printf (
+	"Significance = %4.1f sigma:  Apply residual extinction correction? ")
+	    call pargr (sig)
+	call flush (STDOUT)
+
+	ans = false
+	if (scan() != EOF)
+	    call gargb (ans)
+
+	# Undo last fit if not applying correction.
+	if (!ans)
+	    goto cancel_
+
+	call printf ("Residual extinction correction applied")
+	call sfree (sp)
+	return
+
+cancel_
+	do i = 1, nstds - 2 {
+	    if (STD_FLAG(stds[i]) == SF_EXCLUDE)
+		next
+	    n = STD_NWAVES(stds[i])
+	    a = STD_AIRMASS(stds[i])
+	    xp = STD_WAVES(stds[i])
+	    yp = STD_SENS(stds[i])
+	    do j = 1, n {
+		Memr[yp] = Memr[yp] - a * cveval (ecv, Memr[xp])
+		xp = xp + 1
+		yp = yp + 1
+	    }
+	}
+	call cvfree (ecv)
+	call printf ("Residual extinction correction canceled")
+	call sfree (sp)
+end
diff --git a/noao/onedspec/sensfunc/sffit.x b/noao/onedspec/sensfunc/sffit.x
new file mode 100644
index 00000000..3be306ad
--- /dev/null
+++ b/noao/onedspec/sensfunc/sffit.x
@@ -0,0 +1,78 @@
+include	<math/curfit.h>
+include	"sensfunc.h"
+
+define	FUNCTIONS	"|chebyshev|legendre|spline3|spline1|"
+
+procedure sf_fit (stds, nstds, cv, function, order, xmin, xmax)
+
+pointer	stds[nstds]
+int	nstds
+pointer	cv
+char	function[ARB]
+int	order
+real	xmin
+real	xmax
+
+int	i, n, func, ord, strdic()
+pointer	x, y, w
+
+int	functions[4]
+data	functions/CHEBYSHEV,LEGENDRE,SPLINE3,SPLINE1/
+
+begin
+	func = strdic (function, function, SZ_FNAME, FUNCTIONS)
+	func = functions[max (1, func)]
+	ord = order
+
+	while (ord > 0) {
+	    call cvfree (cv)
+	    call cvinit (cv, func, ord, xmin, xmax)
+	    do i = 1, nstds {
+		if (STD_FLAG(stds[i]) == SF_INCLUDE) {
+	            n = STD_NWAVES(stds[i])
+	            x = STD_WAVES(stds[i])
+	            y = STD_SENS(stds[i])
+	            w = STD_WTS(stds[i])
+	            call cvacpts (cv, Memr[x], Memr[y], Memr[w], n, WTS_USER)
+		}
+	    }
+	    call cvsolve (cv, i)
+	    if (i == OK)
+		break
+
+	    switch (func) {
+	    case CHEBYSHEV, LEGENDRE:
+		ord = ord - 1
+	    case SPLINE3:
+		ord = ord - 1
+		if (ord == 0) {
+		    func = CHEBYSHEV
+		    ord = 2
+		}
+	    case SPLINE1:
+		ord = ord - 1
+		if (ord == 0) {
+		    func = CHEBYSHEV
+		    ord = 1
+		}
+	    }
+	}
+
+	switch (i) {
+	case SINGULAR:
+	    call error (0, "Singular solution")
+	case NO_DEG_FREEDOM:
+	    call error (0, "No degrees of freedom")
+	}
+
+	# Set fitted values
+	do i = 1, nstds
+	    if (STD_FLAG(stds[i]) != SF_EXCLUDE) {
+	        n = STD_NWAVES(stds[i])
+		if (n < 1)
+		    next
+	        x = STD_WAVES(stds[i])
+		y = STD_FIT(stds[i])
+		call cvvector (cv, Memr[x], Memr[y], n)
+	    }
+end
diff --git a/noao/onedspec/sensfunc/sfginit.x b/noao/onedspec/sensfunc/sfginit.x
new file mode 100644
index 00000000..0214c7a7
--- /dev/null
+++ b/noao/onedspec/sensfunc/sfginit.x
@@ -0,0 +1,89 @@
+include	<gset.h>
+include	"sensfunc.h"
+
+# SF_GINIT -- Initialize graphics structure and open the graphics device.
+# This includes CL requests for the starting graphs (default is "sr"),
+# the mark types (default is "plus box"), and graphics device.
+
+procedure sf_ginit (gp)
+
+pointer	gp		# Graphics structure (returned)
+
+int	i, j
+pointer	sp, str, gopen()
+errchk	malloc, gopen
+
+begin
+	call smark (sp)
+	call salloc (str, SZ_FNAME, TY_CHAR)
+
+	call calloc (gp, LEN_GP, TY_STRUCT)
+	do i = 1, SF_NGRAPHS {
+	    call malloc (GP_IMAGES(gp,i), SZ_FNAME, TY_CHAR)
+	    Memc[GP_IMAGES(gp,i)] = EOS
+	    call malloc (GP_SKYS(gp,i), SZ_FNAME, TY_CHAR)
+	    Memc[GP_SKYS(gp,i)] = EOS
+	}
+
+	# Set the starting graph types.
+	call clgstr ("graphs", Memc[str], SZ_FNAME)
+	j = str
+	for (i=str; Memc[i] != EOS; i=i+1) {
+	    switch (Memc[i]) {
+	    case 'a','c','e','i','l','r','s':
+		Memc[j] = Memc[i]
+		j = j + 1
+	    }
+	}
+	Memc[j] = EOS
+	if (Memc[str] != EOS)
+	    call strcpy (Memc[str], GP_GRAPHS(gp,1), SF_NGRAPHS)
+	else
+	    call strcpy ("sr", GP_GRAPHS(gp,1), SF_NGRAPHS)
+
+	# Set the starting mark types and colors.
+	GP_MARK(gp) = GM_PLUS
+	GP_MDEL(gp) = GM_CROSS
+	GP_MADD(gp) = GM_BOX
+	GP_PLCOLOR(gp) = 2
+	GP_CMARK(gp) = 1
+	GP_CDEL(gp) = 3
+	GP_CADD(gp) = 4
+	call clgstr ("marks", Memc[str], SZ_FNAME)
+	call sf_marks (gp, Memc[str])
+	call clgstr ("colors", Memc[str], SZ_FNAME)
+	call sf_colors (gp, Memc[str])
+
+	# Set flux limits
+	GP_FMIN(gp) = INDEF
+	GP_FMAX(gp) = INDEF
+
+	# Open the graphics device.
+	call clgstr ("device", Memc[str], SZ_FNAME)
+	GP_GIO(gp) = gopen (Memc[str], NEW_FILE, STDGRAPH)
+
+	call sfree (sp)
+end
+
+
+# SF_GFREE -- Free the graphics structure and close the graphics device.
+
+procedure sf_gfree (gp)
+
+pointer	gp		# Graphics structure
+
+int	i
+
+begin
+	if (gp == NULL)
+	    return
+
+	call gclose (GP_GIO(gp))
+	do i = 1, SF_NGRAPHS {
+	    call mfree (GP_IMAGES(gp,i), TY_CHAR)
+	    call mfree (GP_SKYS(gp,i), TY_CHAR)
+	    if (GP_SHDR(gp,i) != NULL)
+		call shdr_close (GP_SHDR(gp,i))
+	}
+	call mfree (gp, TY_STRUCT)
+end
diff --git a/noao/onedspec/sensfunc/sfgraph.x b/noao/onedspec/sensfunc/sfgraph.x
new file mode 100644
index 00000000..bb6fb26f
--- /dev/null
+++ b/noao/onedspec/sensfunc/sfgraph.x
@@ -0,0 +1,289 @@
+include	<gset.h>
+include	<error.h>
+include	<math/curfit.h>
+include	"sensfunc.h"
+
+define	NCURVE		50		# Number of vectors in curve
+
+
+# SF_GRAPH -- Graph the desired data on the output graphics device.
+# This entry procedure determines the graph types, sets the device viewports
+# for each graph, and calls a procedure to make the graph.
+
+procedure sf_graph (gp, stds, nstds, cv, wextn, extn, nextn, ecv)
+
+pointer	gp			# Graphics structure
+pointer	stds[nstds]		# Standard star data
+int	nstds			# Number of standard stars
+pointer	cv			# Sensitivity function curve
+real	wextn[nextn]		# Extinction table wavelengths
+real	extn[nextn]		# Extinction table values
+int	nextn			# Number of values in extinction table
+pointer	ecv			# Residual extinction curve.
+
+int	i, image, ngraphs, strlen()
+real	fa[8]
+pointer	sp, id, gio
+
+data	fa/0.,1.,1.,0.,0.,0.,1.,1./
+
+begin
+	# Clear the graphs, write the title, and set the viewports based on
+	# the number of graphs.
+
+	call smark (sp)
+	call salloc (id, SZ_LINE, TY_CHAR)
+	call sysid (Memc[id], SZ_LINE)
+
+	gio = GP_GIO(gp)
+	call gclear (gio)
+	call gseti (gio, G_FACOLOR, 0)
+	call gseti (gio, G_NMINOR, 0)
+	ngraphs = strlen (GP_GRAPHS(gp,1))
+	switch (ngraphs) {
+	case 1:
+	    call gseti (gio, G_WCS, 1)
+	    call gsview (gio, .10, .97, .1, .9) 
+	    GP_SZMARK(gp) = 2.
+	    GP_SZMDEL(gp) = 2.
+	case 2:
+	    call gseti (gio, G_WCS, 1)
+	    call gsview (gio, .10, .97, .55, .9) 
+	    call gseti (gio, G_WCS, 2)
+	    call gsview (gio, .10, .97, .10, .45) 
+	    GP_SZMARK(gp) = 2.
+	    GP_SZMDEL(gp) = 2.
+	case 3:
+	    call gseti (gio, G_WCS, 1)
+	    call gsview (gio, .10, .97, .55, .9) 
+	    call gseti (gio, G_WCS, 2)
+	    call gsview (gio, .10, .50, .1, .45) 
+	    call gseti (gio, G_WCS, 3)
+	    call gsview (gio, .57, .97, .1, .45) 
+	    GP_SZMARK(gp) = 2.
+	    GP_SZMDEL(gp) = 2.
+	default:
+	    call gseti (gio, G_WCS, 1)
+	    call gsview (gio, .10, .50, .55, .9) 
+	    call gseti (gio, G_WCS, 2)
+	    call gsview (gio, .57, .97, .55, .9) 
+	    call gseti (gio, G_WCS, 3)
+	    call gsview (gio, .10, .50, .1, .45) 
+	    call gseti (gio, G_WCS, 4)
+	    call gsview (gio, .57, .97, .1, .45) 
+	    GP_SZMARK(gp) = .01
+	    GP_SZMDEL(gp) = .01
+	}
+
+	# For each graph select the viewport and call a procedure to make
+	# the graph.
+
+	image = 0
+	for (i = 1; GP_GRAPHS(gp,i) != EOS; i = i + 1) {
+	    call gseti (gio, G_WCS, i)
+	    if (i > 1)
+		call gfill (gio, fa, fa[5], 4, GF_SOLID)
+	    switch (GP_GRAPHS(gp,i)) {
+	    case 'a', 's', 'r':
+		call sf_data (stds, nstds, GP_GRAPHS(gp,i))
+		call sf_graph1 (gp, stds, nstds, GP_GRAPHS(gp,i))
+	    case 'e':
+		call sf_egraph (gp, wextn, extn, nextn, ecv)
+	    case 'c':
+		call sf_cgraph (gp, stds, nstds, cv)
+	    case 'i', 'l':
+		if (GP_GRAPHS(gp,i) == 'i')
+		    GP_LOG(gp) = NO
+		else
+		    GP_LOG(gp) = YES
+		image = image + 1
+		iferr (call sf_image (gp, image, stds, nstds, cv, wextn, extn,
+		    nextn, ecv))
+		    call erract (EA_WARN)
+	    }
+	}
+
+	call gseti (gio, G_WCS, 0)
+	call gtext (gio, 0.5, 1., Memc[id], "h=c,v=t,f=b")
+	call gtext (gio, 0.5, 0.97, GP_TITLE(gp), "h=c,v=t,f=b")
+
+	call sfree (sp)
+end
+
+
+# SF_GRAPH1 -- Make graph of sensitivity or residual sensitivity vs wavelength.
+
+procedure sf_graph1 (gp, stds, nstds, graph)
+
+pointer	gp			# Graphics structure
+pointer	stds[nstds]		# Standard star data
+int	nstds			# Number of standard stars
+char	graph			# Graph type
+
+int	i, j, n, mark, mdel, cdel, color
+real	szmark, szmdel, ymin, ymax, y1, y2
+pointer	x, y, w, gio
+
+begin
+	gio = GP_GIO(gp)
+
+	# Autoscale the included data in Y and set wavelength range.
+	j = 0
+	do i = 1, nstds {
+	    if (STD_FLAG(stds[i]) != SF_INCLUDE)
+		next
+	    j = j + 1
+	    n = STD_NWAVES(stds[i])
+	    x = STD_X(stds[i])
+	    y = STD_Y(stds[i])
+	    if (j == 1)
+		call alimr (Memr[y], n, ymin, ymax)
+	    else {
+		call alimr (Memr[y], n, y1, y2)
+		ymin = min (ymin, y1)
+		ymax = max (ymax, y2)
+	    }
+	}
+	y1 = 0.05 * (ymax - ymin)
+	ymin = ymin - y1
+	ymax = ymax + y1
+
+	# Draw axes and title based on type of graph.
+	switch (graph) {
+	case 'a':
+	    call gswind (gio, GP_AIRMASS(gp,1), GP_AIRMASS(gp,2), ymin, ymax)
+	    call glabax (gio, "Sensitivity Residuals vs Airmass", "", "")
+	case 's':
+	    call gswind (gio, GP_WSTART(gp), GP_WEND(gp), ymin, ymax)
+	    call glabax (gio, "Sensitivity vs Wavelength", "", "")
+	case 'r':
+	    call gswind (gio, GP_WSTART(gp), GP_WEND(gp), ymin, ymax)
+	    call glabax (gio, "Sensitivity Residuals vs Wavelength", "", "")
+	}
+
+	# Plot the data with appropriate mark types and sizes.
+	mdel = GP_MDEL(gp)
+	cdel = GP_CDEL(gp)
+	szmdel = GP_SZMDEL(gp)
+	szmark = GP_SZMARK(gp)
+	do i = 1, nstds {
+	    if (STD_FLAG(stds[i]) != SF_INCLUDE)
+		next
+	    if (i != nstds-1) {
+		mark = GP_MARK(gp)
+		color = GP_CMARK(gp)
+	    } else {
+		mark = GP_MADD(gp)
+		color = GP_CADD(gp)
+	    }
+	    n = STD_NWAVES(stds[i])
+	    x = STD_X(stds[i]) - 1
+	    y = STD_Y(stds[i]) - 1
+	    w = STD_WTS(stds[i]) - 1
+	    do j = 1, n {
+		if (Memr[w+j] == 0.) {
+		    call gseti (gio, G_PLCOLOR, cdel)
+		    call gmark (gio, Memr[x+j], Memr[y+j], mdel, szmdel, szmdel)
+		} else {
+		    call gseti (gio, G_PLCOLOR, color)
+	            call gmark (gio, Memr[x+j], Memr[y+j], mark, szmark, szmark)
+		}
+	    }
+	}
+end
+
+
+# SF_EGRAPH -- Graph extinction curves with and without residual extinction
+# correction.
+
+procedure sf_egraph (gp, wextn, extn, nextn, ecv)
+
+pointer	gp			# Graphics structure
+real	wextn[nextn]		# Extinction table wavelengths
+real	extn[nextn]		# Extinction table values
+int	nextn			# Number of extinction table values
+pointer	ecv			# Residual extinction curve
+
+int	i, j
+real	xmin, xmax, dx, x, cveval()
+pointer	sp, ext, extnew, gio
+
+begin
+	call smark (sp)
+	call salloc (ext, NCURVE, TY_REAL)
+
+	# Interpolate extinction table to a grid of wavelengths within
+	# the range of the sensitivity data.
+
+	gio = GP_GIO(gp)
+	xmin = GP_WSTART(gp)
+	xmax = GP_WEND(gp)
+	dx = (xmax - xmin) / (NCURVE - 1)
+	x = xmin
+	do i = 0, NCURVE-1 {
+	    call intrp (1, wextn, extn, nextn, x, Memr[ext+i], j)
+	    x = x + dx
+	}
+	call gswind (gio, xmin, xmax, INDEF, INDEF)
+	call gascale (gio, Memr[ext], NCURVE, 2)
+
+	# If there is a residual extinction curve determine a new extinction
+	# curve.
+
+	if (ecv != NULL) {
+	    call salloc (extnew, NCURVE, TY_REAL)
+	    call amovr (Memr[ext], Memr[extnew], NCURVE)
+	    x = xmin
+	    do i = 0, NCURVE-1 {
+	        Memr[extnew+i] = Memr[extnew+i] + cveval (ecv, x)
+		x = x + dx
+	    }
+	    call grscale (gio, Memr[extnew], NCURVE, 2)
+	}
+
+	# Draw the axes and title and extinction curves.
+	call glabax (gio, "Extinction vs Wavelength", "", "")
+	call gseti (gio, G_PLCOLOR, GP_PLCOLOR(gp))
+	call gvline (gio, Memr[ext], NCURVE, xmin, xmax)
+	if (ecv != NULL) {
+	    call gseti (gio, G_PLTYPE, 2)
+	    call gseti (gio, G_PLCOLOR, GP_PLCOLOR(gp)+1)
+	    call gvline (gio, Memr[extnew], NCURVE, xmin, xmax)
+	    call gseti (gio, G_PLTYPE, 1)
+	}
+
+	call sfree (sp)
+end
+
+
+# SF_FITGRAPH -- Overplot the fitted sensitivity curve.
+
+procedure sf_fitgraph (gp, cv)
+
+pointer	gp			# Graphics structure
+pointer	cv			# Sensitivity function curve
+
+int	i, j
+real	x1, x2, y1, y2, cveval()
+pointer	gio
+
+begin
+	gio = GP_GIO(gp)
+	call gseti (gio, G_PLCOLOR, GP_PLCOLOR(gp))
+
+	# Only plot on sensitivity curve graph types.
+	for (j = 1; GP_GRAPHS(gp,j) != EOS; j = j + 1) {
+	    if (GP_GRAPHS(gp,j) != 's')
+		next
+	    call gseti (gio, G_WCS, j)
+	    call ggwind (gio, x1, x2, y1, y2)
+	    x2 = (x2 - x1) / NCURVE
+	    y1 = cveval (cv, x1)
+	    call gamove (gio, x1, y1)
+	    do i = 1, NCURVE {
+	        x1 = x1 + x2
+	        y1 = cveval (cv, x1)
+	        call gadraw (gio, x1, y1)
+	    }
+	}
+end
diff --git a/noao/onedspec/sensfunc/sfimage.x b/noao/onedspec/sensfunc/sfimage.x
new file mode 100644
index 00000000..71edc213
--- /dev/null
+++ b/noao/onedspec/sensfunc/sfimage.x
@@ -0,0 +1,234 @@
+include	<gset.h>
+include	<math/curfit.h>
+include	"sensfunc.h"
+include	<smw.h>
+
+
+# SF_IMAGE -- Graph fluxed image data and possible standard flux points.
+# For efficiency the IMIO pointer, buffer, and associated data are kept
+# since a redraw is a common occurence and generating the data is slow.
+
+procedure sf_image (gp, wc, stds, nstds, cv, wextn, extn, nextn, ecv)
+
+pointer	gp			# Graphics structure
+int	wc			# WC of graph
+pointer	stds[nstds]		# Standard star data for flux points
+int	nstds			# Number of standard stars
+pointer	cv			# Sensitivity function curve
+real	wextn[nextn]		# Extinction table wavelengths
+real	extn[nextn]		# Extinction table values
+int	nextn			# Number of extinction table values
+pointer	ecv			# Residual extinction curve
+
+int	scale[SF_NGRAPHS], log[SF_NGRAPHS]
+
+bool	newobs, obshead
+int	i, j, n, err
+real	a, t, w, dw, e, sens, latitude, smin, smax, xmin, xmax
+pointer	im, mw, sh, skyim, skymw, skysh, std, gio, sp, str, x, y, z, obs
+pointer	immap(), smw_openim()
+real	cveval(), obsgetr(), cvstatr()
+double	shdr_lw()
+bool	streq(), strne()
+errchk	immap, smw_openim, obsimopen
+
+define	plot_	99
+
+begin
+	# Return if no image name.
+	if (Memc[GP_IMAGES(gp,wc)] == EOS)
+	    return
+
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	# Get the spectrum and sky subtract if necessary.
+	sh = GP_SHDR(gp,wc)
+	if (sh != NULL) {
+	    if (streq (Memc[GP_IMAGES(gp,wc)], IMNAME(sh))) {
+		if (GP_LOG(gp) == log[wc])
+		    goto plot_
+		else
+		    call shdr_close (sh)
+	    }
+	}
+
+	# Determine a valid standard star to get aperture number.
+	do i = 1, nstds
+	    if (STD_FLAG(stds[i]) != SF_EXCLUDE) {
+		std = stds[i]
+		break
+	    }
+
+	im = immap (Memc[GP_IMAGES(gp,wc)], READ_ONLY, 0)
+	mw = smw_openim (im)
+	call shdr_open (im, mw, 1, 1, STD_BEAM(std), SHDATA, sh)
+
+	# Check for dispersion correction
+	if (DC(sh) == DCNO) {
+	    call shdr_close (sh)
+	    call smw_close (mw)
+	    call imunmap (im)
+	    GP_SHDR(gp,wc) = NULL
+	    call sfree (sp)
+	    call printf ("-%s must be dispersion corrected-")
+		call pargstr (Memc[GP_IMAGES(gp,wc)])
+	    return
+	}
+
+	# Sky subtract if necessary
+	if (Memc[GP_SKYS(gp,wc)] != EOS) {
+	    skyim = immap (Memc[GP_SKYS(gp,wc)], READ_ONLY, 0)
+	    skymw = smw_openim (skyim)
+	    call shdr_open (skyim, skymw, 1, 1, STD_BEAM(std), SHDATA, skysh)
+	    call shdr_rebin (skysh, sh)
+	    call asubr (Memr[SY(sh)], Memr[SY(skysh)], Memr[SY(sh)], SN(sh))
+	    call shdr_close (skysh)
+	    call smw_close (skymw)
+	    call imunmap (skyim)
+	}
+
+	# Set airmass and exposure time
+	if (IS_INDEF (AM(sh))) {
+	    obs = NULL
+	    call clgstr ("observatory", Memc[str], SZ_LINE)
+	    call obsimopen (obs, im, Memc[str], NO, newobs, obshead)
+	    latitude = obsgetr (obs, "latitude")
+	    call obsclose (obs)
+	    call get_airm (RA(sh), DEC(sh), HA(sh), ST(sh), latitude,
+		AM(sh))
+	}
+	a = AM(sh)
+	if (IS_INDEF (IT(sh)))
+	    t = 1.
+	else
+	    t = IT(sh)
+
+	# Apply extinction correction if needed
+	if (EC(sh) == ECNO) {
+	    if (ecv != NULL) {
+	        xmin = cvstatr (ecv, CVXMIN)
+	        xmax = cvstatr (ecv, CVXMAX)
+	    }
+	    do i = 1, SN(sh) {
+		w = Memr[SX(sh)+i-1]
+	        call intrp (1, wextn, extn, nextn, w, e, err)
+	        if (ecv != NULL)
+		    e = e + cveval (ecv, min (xmax, max (w, xmin)))
+	        Memr[SY(sh)+i-1] = Memr[SY(sh)+i-1] * 10. ** (0.4 * a * e)
+	    }
+	} else {
+	    call printf ("-%s already extinction corrected-")
+	       call pargstr (Memc[GP_IMAGES(gp,wc)])
+	}
+
+	# Apply flux calibration if needed
+	if (FC(sh) == FCNO) {
+	    do i = 1, SN(sh) {
+	        w = Memr[SX(sh)+i-1]
+		dw = abs (shdr_lw (sh, double (i+0.5)) -
+		    shdr_lw (sh, double (i-0.5)))
+	        sens = cveval (cv, w)
+	        Memr[SY(sh)+i-1] = Memr[SY(sh)+i-1] / t / dw / 10.**(0.4*sens)
+	    }
+	} else {
+	    call printf ("-%s already flux calibrated-")
+		call pargstr (Memc[GP_IMAGES(gp,wc)])
+	}
+
+	# Set flux scaling
+	call alimr (Memr[SY(sh)], SN(sh), smin, smax)
+	if (smax < 0.)
+	    scale[wc] = 0.
+	else if (GP_LOG(gp) == NO) {
+	    scale[wc] = -log10 (smax) + 1
+	    w = 10. ** scale[wc]
+	    call amulkr (Memr[SY(sh)], w, Memr[SY(sh)], SN(sh))
+	} else {
+	    scale[wc] = INDEFI
+	    smin = smax / 1000.
+	    w = smax
+	    y = SY(sh)
+	    do i = 1, SN(sh) {
+		if (Memr[y] > smin)
+		    w = min (w, Memr[y])
+		y = y + 1
+	    }
+	    y = SY(sh)
+	    do i = 1, SN(sh) {
+		Memr[y] = log10 (max (Memr[y], w))
+		y = y + 1
+	    }
+	}
+	log[wc] = GP_LOG(gp)
+
+	# Save the spectrum for future redraw.
+	call smw_close (MW(sh))
+	call imunmap (im)
+	GP_SHDR(gp,wc) = sh
+
+plot_
+	# Plot scaled graph.
+	smin = GP_FMIN(gp)
+	smax = GP_FMAX(gp)
+	if (IS_INDEFI(scale[wc])) {
+	    call sprintf (Memc[str], SZ_LINE, "%s: Log Flux")
+	        call pargstr (Memc[GP_IMAGES(gp,wc)])
+	    if (!IS_INDEF(smin)) {
+		if (smin > 0.)
+		    smin = log10 (smin)
+		else
+		    smin = INDEF
+	    }
+	    if (!IS_INDEF(smax)) {
+		if (smax > 0.)
+		    smax = log10 (smax)
+		else
+		    smax = INDEF
+	    }
+	} else if (scale[wc] != 0) {
+	    call sprintf (Memc[str], SZ_LINE, "%s: Flux x 1E%d")
+	        call pargstr (Memc[GP_IMAGES(gp,wc)])
+	        call pargi (scale[wc])
+	    w = 10. ** scale[wc]
+	    if (!IS_INDEF(smin))
+		smin = w * smin
+	    if (!IS_INDEF(smax))
+		smax = w * smax
+	} else {
+	    call sprintf (Memc[str], SZ_LINE, "%s: Flux")
+	        call pargstr (Memc[GP_IMAGES(gp,wc)])
+	    w = 1.
+	}
+
+	gio = GP_GIO(gp)
+	call gascale (gio, Memr[SX(sh)], SN(sh), 1)
+	call gascale (gio, Memr[SY(sh)], SN(sh), 2)
+	call gswind (gio, INDEF, INDEF, smin, smax)
+	call glabax (gio, Memc[str], "", "")
+	call gseti (gio, G_PLCOLOR, GP_PLCOLOR(gp))
+	call gpline (gio, Memr[SX(sh)], Memr[SY(sh)], SN(sh))
+
+	call sfree (sp)
+
+	# Check if image is one of the standard stars and plot flux points.
+	do i = 1, nstds {
+	    if (strne (Memc[GP_IMAGES(gp,wc)], STD_IMAGE(stds[i])))
+		next
+	    n = STD_NWAVES(stds[i])
+	    x = STD_WAVES(stds[i])
+	    y = STD_FLUXES(stds[i])
+	    z = STD_DWAVES(stds[i])
+	    call gseti (gio, G_PMLTYPE, 1)
+	    call gseti (gio, G_PLCOLOR, GP_CMARK(gp))
+	    if (IS_INDEFI(scale[wc])) {
+	        do j = 0, n-1
+		    call gmark (gio, Memr[x+j], log10 (Memr[y+j]), GM_HEBAR,
+			-Memr[z+j], 1.)
+	    } else {
+	        do j = 0, n-1
+		    call gmark (gio, Memr[x+j], w * Memr[y+j], GM_HEBAR,
+			-Memr[z+j], 1.)
+	    }
+	}
+end
diff --git a/noao/onedspec/sensfunc/sfmarks.x b/noao/onedspec/sensfunc/sfmarks.x
new file mode 100644
index 00000000..39d85af6
--- /dev/null
+++ b/noao/onedspec/sensfunc/sfmarks.x
@@ -0,0 +1,46 @@
+include	<gset.h>
+include	"sensfunc.h"
+
+define	GMTYPES	"|point|box|plus|cross|diamond|hline|vline|hebar|vebar|circle|"
+
+
+# SF_MARKS -- Decode user mark types into GIO mark types.  The input string
+# consists of two whitespace separated mark types.
+
+procedure sf_marks (gp, marks)
+
+pointer	gp
+char	marks[ARB]
+
+int	i, nscan(), strdic()
+pointer	sp, str
+
+int	gmtypes[10]
+data	gmtypes /GM_POINT,GM_BOX,GM_PLUS,GM_CROSS,GM_DIAMOND,GM_HLINE,GM_VLINE,
+		GM_HEBAR,GM_VEBAR,GM_CIRCLE/
+
+begin
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+	call sscan (marks)
+	call gargwrd (Memc[str], SZ_LINE)
+	if (nscan() == 1) {
+	    i = strdic (Memc[str], Memc[str], SZ_LINE, GMTYPES)
+	    if (i != 0)
+		GP_MARK(gp) = gmtypes[i]
+	}
+	call gargwrd (Memc[str], SZ_LINE)
+	if (nscan() == 2) {
+	    i = strdic (Memc[str], Memc[str], SZ_LINE, GMTYPES)
+	    if (i != 0)
+		GP_MDEL(gp) = gmtypes[i]
+	}
+	call gargwrd (Memc[str], SZ_LINE)
+	if (nscan() == 3) {
+	    i = strdic (Memc[str], Memc[str], SZ_LINE, GMTYPES)
+	    if (i != 0)
+		GP_MADD(gp) = gmtypes[i]
+	}
+
+	call sfree (sp)
+end
diff --git a/noao/onedspec/sensfunc/sfmove.x b/noao/onedspec/sensfunc/sfmove.x
new file mode 100644
index 00000000..4451938e
--- /dev/null
+++ b/noao/onedspec/sensfunc/sfmove.x
@@ -0,0 +1,166 @@
+include	<gset.h>
+include "sensfunc.h"
+
+
+# SF_MOVE -- Move point, star, or wavelength.
+
+procedure sf_move (gp, stds, nstds, key, istd, ipt, shift)
+
+pointer	gp			# GIO pointer
+pointer	stds[nstds]		# Standard star data
+int	nstds			# Number of standard stars
+int	key			# Delete point, star, or wavelength
+int	istd			# Index of standard star
+int	ipt			# Index of point
+real	shift
+
+int	i, j, n, wcs, mark, mdel, cdel, color, stridx()
+real	wave, szmark, szmdel
+pointer	x, y, z, w, gio
+
+begin
+	gio = GP_GIO(gp)
+	mdel = GP_MDEL(gp)
+	cdel = GP_CDEL(gp)
+	szmdel = GP_SZMDEL(gp)
+	szmark = GP_SZMARK(gp)
+
+	# Move points in each displayed graph.
+	for (wcs = 1; GP_GRAPHS(gp,wcs) != EOS; wcs = wcs + 1) {
+	    if (stridx (GP_GRAPHS(gp,wcs), "ars") == 0)
+	        next
+
+	    call gseti (gio, G_WCS, wcs)
+	    call sf_data (stds, nstds, GP_GRAPHS(gp,wcs))
+	    switch (key) {
+	    case 'p':
+		if (istd != nstds-1) {
+		    mark = GP_MARK(gp)
+		    color = GP_CMARK(gp)
+		} else {
+		    mark = GP_MADD(gp)
+		    color = GP_CADD(gp)
+		}
+	        x = STD_X(stds[istd])+ipt-1
+	        y = STD_Y(stds[istd],1)+ipt-1
+	        w = STD_WTS(stds[istd])+ipt-1
+		if (Memr[w] != 0.) {
+		    call gseti (gio, G_PMLTYPE, 0)
+	            call gmark (gio, Memr[x], Memr[y], mark, szmark, szmark)
+		    call gseti (gio, G_PMLTYPE, 1)
+		    call gseti (gio, G_PLCOLOR, color)
+	            call gmark (gio, Memr[x], Memr[y]+shift, mark, szmark,
+			szmark)
+		} else {
+		    call gseti (gio, G_PMLTYPE, 0)
+	            call gmark (gio, Memr[x], Memr[y], mdel, szmdel, szmdel)
+		    call gseti (gio, G_PMLTYPE, 1)
+		    call gseti (gio, G_PLCOLOR, cdel)
+	            call gmark (gio, Memr[x], Memr[y]+shift, mdel, szmdel,
+			szmdel)
+		}
+	    case 's':
+		if (istd != nstds-1) {
+		    mark = GP_MARK(gp)
+		    color = GP_CMARK(gp)
+		} else {
+		    mark = GP_MADD(gp)
+		    color = GP_CADD(gp)
+		}
+	        n = STD_NWAVES(stds[istd])
+	        x = STD_X(stds[istd])
+	        y = STD_Y(stds[istd])
+	        w = STD_WTS(stds[istd])
+		do i = 1, n {
+		    if (Memr[w] != 0.) {
+		        call gseti (gio, G_PMLTYPE, 0)
+	                call gmark (gio, Memr[x], Memr[y], mark, szmark, szmark)
+		        call gseti (gio, G_PMLTYPE, 1)
+		        call gseti (gio, G_PLCOLOR, color)
+	                call gmark (gio, Memr[x], Memr[y]+shift, mark, szmark,
+			    szmark)
+		    } else {
+		        call gseti (gio, G_PMLTYPE, 0)
+	                call gmark (gio, Memr[x], Memr[y], mdel, szmdel, szmdel)
+		        call gseti (gio, G_PMLTYPE, 1)
+		        call gseti (gio, G_PLCOLOR, cdel)
+	                call gmark (gio, Memr[x], Memr[y]+shift, mdel, szmdel,
+			    szmdel)
+		    }
+		    x = x + 1
+		    y = y + 1
+		    w = w + 1
+		}
+	    case 'w':
+	        wave = Memr[STD_WAVES(stds[istd])+ipt-1]
+	        do i = 1, nstds {
+	            if (STD_FLAG(stds[i]) != SF_INCLUDE)
+		        next
+		    if (i != nstds-1) {
+			mark = GP_MARK(gp)
+			color = GP_CMARK(gp)
+		    } else {
+			mark = GP_MADD(gp)
+			color = GP_CADD(gp)
+		    }
+	            n = STD_NWAVES(stds[i])
+	            x = STD_X(stds[i])
+	            y = STD_Y(stds[i])
+		    z = STD_WAVES(stds[i])
+	            w = STD_WTS(stds[i])
+		    do j = 1, n {
+		        if (Memr[z] == wave) {
+			    if (Memr[w] != 0.) {
+		                call gseti (gio, G_PMLTYPE, 0)
+	                        call gmark (gio, Memr[x], Memr[y], mark, szmark,
+				    szmark)
+		                call gseti (gio, G_PMLTYPE, 1)
+		                call gseti (gio, G_PLCOLOR, color)
+	                        call gmark (gio, Memr[x], Memr[y]+shift, mark,
+				    szmark, szmark)
+			    } else {
+		                call gseti (gio, G_PMLTYPE, 0)
+	                        call gmark (gio, Memr[x], Memr[y], mdel, szmdel,
+				    szmdel)
+		                call gseti (gio, G_PMLTYPE, 1)
+		                call gseti (gio, G_PLCOLOR, cdel)
+	                        call gmark (gio, Memr[x], Memr[y]+shift, mdel,
+				    szmdel, szmdel)
+			    }
+		        }
+			x = x + 1
+			y = y + 1
+			z = z + 1
+			w = w + 1
+		    }
+	        }
+	    }
+	}
+
+	# Now add the shift to the data.
+	switch (key) {
+	case 'p':
+	    y = STD_SENS(stds[istd])+ipt-1
+	    Memr[y] = Memr[y] + shift
+	case 's':
+	    n = STD_NWAVES(stds[istd])
+	    y = STD_SENS(stds[istd])
+	    call aaddkr (Memr[y], shift, Memr[y], n)
+	    STD_SHIFT(stds[istd]) = STD_SHIFT(stds[istd]) + shift
+	case 'w':
+	    wave = Memr[STD_WAVES(stds[istd])+ipt-1]
+	    do i = 1, nstds {
+	        if (STD_FLAG(stds[i]) != SF_INCLUDE)
+		    next
+	        n = STD_NWAVES(stds[i])
+		z = STD_WAVES(stds[i])
+	        y = STD_SENS(stds[i])
+		do j = 1, n {
+		    if (Memr[z] == wave)
+	                Memr[y] = Memr[y] + shift
+		    w = w + 1
+		    y = y + 1
+		}
+	    }
+	}
+end
diff --git a/noao/onedspec/sensfunc/sfnearest.x b/noao/onedspec/sensfunc/sfnearest.x
new file mode 100644
index 00000000..540faad2
--- /dev/null
+++ b/noao/onedspec/sensfunc/sfnearest.x
@@ -0,0 +1,69 @@
+include	<gset.h>
+include	<mach.h>
+include "sensfunc.h"
+
+# SF_NEAREST -- Find the nearest point to the cursor.  Return the standard
+# star index and the wavelength point index.  The metric is in NDC.
+# The cursor is moved to the nearest point selected.  Return zero for
+# the standard star index if valid point not found.
+
+procedure sf_nearest (gp, stds, nstds, wx, wy, wcs, type, istd, ipt)
+
+pointer	gp			# Graphics pointer
+pointer	stds[nstds]		# Standard star data
+int	nstds			# Number of standard stars
+real	wx, wy			# Cursor position
+int	wcs			# WCS
+int	type			# Type of points (0=not del, 1=del, 2=both)
+int	istd			# Index of standard star (returned)
+int	ipt			# Index of point (returned)
+
+int	i, j, n, stridx()
+real	x0, y0, x1, y1, x2, y2, r2, r2min
+pointer	x, y, w, gio
+
+begin
+	# Check for valid wc.
+	istd = 0
+	if (stridx (GP_GRAPHS(gp,wcs), "ars") == 0)
+	    return
+
+	# Transform world cursor coordinates to NDC.
+	gio = GP_GIO(gp)
+	call gctran (gio, wx, wy, wx, wy, wcs, 0)
+
+	# Search for nearest point.
+	r2min = MAX_REAL
+	do i = 1, nstds {
+	    if (STD_FLAG(stds[i]) != SF_INCLUDE)
+		next
+	    n = STD_NWAVES(stds[i])
+	    x = STD_X(stds[i]) - 1
+	    y = STD_Y(stds[i]) - 1
+	    w = STD_WTS(stds[i]) - 1
+	    do j = 1, n {
+		if (type == 0) {
+		    if (Memr[w+j] == 0.)
+			next
+		} else if (type == 1) {
+		    if (Memr[w+j] != 0.)
+			next
+		}
+		x1 = Memr[x+j]
+		y1 = Memr[y+j]
+	        call gctran (gio, x1, y1, x0, y0, wcs, 0)
+	        r2 = (x0 - wx) ** 2 + (y0 - wy) ** 2
+	        if (r2 < r2min) {
+		    r2min = r2
+		    istd = i
+		    ipt = j
+		    x2 = x1
+		    y2 = y1
+	        }
+	    }
+	}
+
+	# Move the cursor to the selected point.
+	call gseti (gio, G_WCS, wcs)
+	call gscur (gio, x2, y2)
+end
diff --git a/noao/onedspec/sensfunc/sfoutput.x b/noao/onedspec/sensfunc/sfoutput.x
new file mode 100644
index 00000000..e21df280
--- /dev/null
+++ b/noao/onedspec/sensfunc/sfoutput.x
@@ -0,0 +1,114 @@
+include	<mach.h>
+include	<imhdr.h>
+include	"sensfunc.h"
+
+
+# SF_OUTPUT -- Write the sensitivity function image.
+
+procedure sf_output (stds, nstds, cv, output, ignoreaps)
+
+pointer	stds[nstds]		# Standard star data
+int	nstds			# Number of standard stars
+pointer	cv			# Sensitivity function curve
+char	output[SZ_FNAME]	# Output root image name (must be SZ_FNAME)
+bool	ignoreaps		# Ignore apertures?
+
+int	i, ap, nw, scan(), nscan()
+real	w1, w2, dw, dw1, aplow[2], aphigh[2], cveval()
+pointer	sp, fname, std, im, mw, buf, immap(), mw_open(), impl1r()
+errchk	imaddi, imaddr
+
+define	makeim_	99
+
+begin
+	# Return if no output root sensitivity imagename is specified.
+	if (output[1] == EOS)
+	    return
+
+	call smark (sp)
+	call salloc (fname, SZ_FNAME, TY_CHAR)
+
+	# Determine wavelength range and reference standard star.
+	w1 = MAX_REAL
+	w2 = -MAX_REAL
+	dw = MAX_REAL
+	do i = 1, nstds-2 {
+	    if (STD_FLAG(stds[i]) == SF_EXCLUDE)
+		next
+	    std = stds[i]
+	    dw1 = abs ((STD_WEND(std) - STD_WSTART(std)) / (STD_NPTS(std) - 1))
+	    w1 = min (w1, STD_WSTART(std), STD_WEND(std))
+	    w2 = max (w2, STD_WSTART(std), STD_WEND(std))
+	    dw = min (dw, dw1)
+	}
+	nw = (w2 - w1) / dw + 1.5
+
+	# Make output image name with aperture number appended.  If the
+	# image exists allow the user to change root name.
+makeim_
+	if (ignoreaps) {
+	    call strcpy (output, Memc[fname], SZ_FNAME)
+	} else {
+	    call sprintf (Memc[fname], SZ_FNAME, "%s.%04d")
+		call pargstr (output)
+		call pargi (STD_BEAM(std))
+	}
+
+	iferr (im = immap (Memc[fname], NEW_IMAGE, 0)) {
+	    call printf ("Cannot create %s --  Enter new name: ")
+		call pargstr (Memc[fname])
+	    call flush (STDOUT)
+	    if (scan() != EOF) {
+	        call gargwrd (Memc[fname], SZ_FNAME)
+	        if (nscan() == 1) {
+	    	    call strcpy (Memc[fname], output, SZ_FNAME)
+	    	    goto makeim_
+		}
+	    }
+	    call printf ("No sensitivity function created for aperture %2d\n")
+	        call pargi (STD_BEAM(std))
+	    call flush (STDOUT)
+	    return
+	}
+
+	# Define the image header.
+	IM_NDIM(im) = 1
+	IM_LEN(im,1) = nw
+	IM_PIXTYPE(im) = TY_REAL
+	if (ignoreaps) {
+	    call sprintf (IM_TITLE(im), SZ_FNAME,
+		"Sensitivity function for all apertures")
+	} else {
+	    call sprintf (IM_TITLE(im), SZ_FNAME,
+		"Sensitivity function for aperture %d")
+		call pargi (STD_BEAM(std))
+	}
+
+	mw = mw_open (NULL, 1)
+	call mw_newsystem (mw, "equispec", 1)
+	call mw_swtype (mw, 1, 1, "linear", "label=Wavelength units=Angstroms")
+	call smw_open (mw, NULL, im)
+	ap = STD_BEAM(std)
+	aplow[1] = INDEF
+	aphigh[1] = INDEF
+	aplow[2] = INDEF
+	aphigh[2] = INDEF
+	call smw_swattrs (mw, 1, 1, ap, STD_BEAM(std), 0,
+	    double(w1), double(dw), nw, 0D0, aplow, aphigh, "")
+	call smw_saveim (mw, im)
+	call smw_close (mw)
+
+	# Write sensitivity data.
+	buf = impl1r (im)
+	do i = 0, nw-1
+	    Memr[buf+i] = cveval (cv, w1 + i * dw)
+
+	# Notify user.
+	call printf ("%s --> %s\n")
+		call pargstr (IM_TITLE(im))
+		call pargstr (Memc[fname])
+	call flush (STDOUT)
+
+	call imunmap (im)
+	call sfree (sp)
+end
diff --git a/noao/onedspec/sensfunc/sfreset.x b/noao/onedspec/sensfunc/sfreset.x
new file mode 100644
index 00000000..fc4d974e
--- /dev/null
+++ b/noao/onedspec/sensfunc/sfreset.x
@@ -0,0 +1,62 @@
+include	"sensfunc.h"
+
+# SF_RESET -- Reset the standard star data to the original input.
+# This is called cancel changes and start over.
+
+procedure sf_reset (stds, nstds, wextn, extn, nextn, ecv, shift)
+
+pointer	stds[nstds]		# Standard star data
+int	nstds			# Number of standard stars
+real	wextn[nextn]		# Extinction table wavelengths
+real	extn[nextn]		# Extinction table values
+int	nextn			# Number of extinction values
+pointer	ecv			# Residual extinction curve
+int	shift			# Shift flag
+
+int	i, j, n, err
+real	exptime, airmass, ext
+pointer	waves, fluxes, dwaves, counts, sens, iwts, wts
+
+begin
+	# Reset the flags, sensitivity, and weight values.
+	shift = NO
+        do i = 1, nstds - 2 {
+	    if (STD_FLAG(stds[i]) == SF_EXCLUDE)
+		next
+	    STD_FLAG(stds[i]) = SF_INCLUDE
+	    STD_SHIFT(stds[i]) = 0.
+	    n = STD_NWAVES(stds[i])
+	    exptime = STD_EXPTIME(stds[i])
+	    airmass = STD_AIRMASS(stds[i])
+	    waves = STD_WAVES(stds[i])
+	    fluxes = STD_FLUXES(stds[i])
+	    dwaves = STD_DWAVES(stds[i])
+	    counts = STD_COUNTS(stds[i])
+	    sens = STD_SENS(stds[i])
+	    iwts = STD_IWTS(stds[i])
+	    wts = STD_WTS(stds[i])
+	    do j = 1, n {
+	        call intrp (1, wextn, extn, nextn, Memr[waves], ext, err)
+		Memr[sens] = Memr[counts] /
+		    (Memr[fluxes] * Memr[dwaves] * exptime)
+		Memr[sens] = 2.5 * log10 (Memr[sens]) + airmass * ext
+		Memr[wts] = Memr[iwts]
+
+		waves = waves + 1
+		fluxes = fluxes + 1
+		dwaves = dwaves + 1
+		counts = counts + 1
+		sens = sens + 1
+		iwts = iwts + 1
+		wts = wts + 1
+	    }
+	}
+
+	# Reset the added and composite stars.
+	STD_NWAVES(stds[nstds-1]) = 0
+	STD_FLAG(stds[nstds-1]) = SF_DELETE
+	STD_FLAG(stds[nstds]) = SF_EXCLUDE
+
+	# Clear the residual extinction curve.
+	call cvfree (ecv)
+end
diff --git a/noao/onedspec/sensfunc/sfrms.x b/noao/onedspec/sensfunc/sfrms.x
new file mode 100644
index 00000000..72b8ea98
--- /dev/null
+++ b/noao/onedspec/sensfunc/sfrms.x
@@ -0,0 +1,43 @@
+include	"sensfunc.h"
+
+
+# SF_RMS -- Compute the RMS of the sensitivity function fit.
+
+procedure sf_rms (stds, nstds, rms, npts)
+
+pointer	stds[nstds]		# Standard star data
+int	nstds			# Number of standard stars
+real	rms			# RMS about fit (returned)
+int	npts			# Number of points in fit (excluding zero wts.)
+
+int	i, j, f, n
+pointer	x, y, w
+
+begin
+	npts = 0
+	rms = 0.
+	do i = 1, nstds {
+	    if (STD_FLAG(stds[i]) != SF_INCLUDE)
+		next
+	    n = STD_NWAVES(stds[i])
+	    x = STD_WAVES(stds[i])
+	    y = STD_SENS(stds[i])
+	    f = STD_FIT(stds[i])
+	    w = STD_WTS(stds[i])
+	    do j = 1, n {
+		if (Memr[w] != 0.) {
+		    rms = rms + (Memr[y] - Memr[f]) ** 2
+		    npts = npts + 1
+		}
+		x = x + 1
+		y = y + 1
+		f = f + 1
+		w = w + 1
+	    }
+	}
+
+	if (npts > 1)
+	    rms = sqrt (rms / (npts - 1))
+	else
+	    rms = INDEF
+end
diff --git a/noao/onedspec/sensfunc/sfsensfunc.x b/noao/onedspec/sensfunc/sfsensfunc.x
new file mode 100644
index 00000000..ee2f1b2a
--- /dev/null
+++ b/noao/onedspec/sensfunc/sfsensfunc.x
@@ -0,0 +1,255 @@
+include	<error.h>
+include	<gset.h>
+include	<mach.h>
+include	"sensfunc.h"
+
+define	KEY		"noao$onedspec/sensfunc/sensfunc.key"
+define	PROMPT		"sensfunc options"
+
+
+# SF_SENSFUNC -- Interactive sensitivity function determination.
+
+procedure sf_sensfunc (gp, stds, nstds, wextn, extn, nextn, sensimage, logfile,
+    ecv, function, order, ignoreaps, interactive)
+
+pointer	gp			# Graphics structure
+pointer	stds[nstds]		# Pointer to standard observations
+int	nstds			# Number of standards
+real	wextn[nextn]		# Extinction table wavelengths
+real	extn[nextn]		# Extinction table values
+int	nextn			# Number of extinction table values
+char	sensimage[ARB]		# Output rootname
+char	logfile[ARB]		# Statistics filename
+pointer	ecv			# Residual extinction curve
+char	function[ARB]		# Fitting function type
+int	order			# Function order
+bool	ignoreaps		# Ignore apertures?
+int	interactive		# Interactive?
+
+char	cmd[SZ_FNAME]
+int	wc, key, newgraph, newfit
+real	wx, wy
+
+int	i, j, aperture, shift, npts, fd, open()
+real	xmin, xmax, rms, delta
+pointer	cv
+
+int	clgcur(), scan(), nscan(), clgwrd()
+errchk	open
+
+define	output_	99
+
+begin
+	# Initialize data and do the initial fit.
+	call sf_reset (stds, nstds, wextn, extn, nextn, ecv, shift)
+
+	xmin = MAX_REAL
+	xmax = -MAX_REAL
+	do i = 1, nstds - 2 {
+	    if (STD_FLAG(stds[i]) == SF_EXCLUDE)
+		next
+	    aperture = STD_BEAM(stds[i])
+	    xmin = min (xmin, STD_WSTART(stds[i]), STD_WEND(stds[i]))
+	    xmax = max (xmax, STD_WSTART(stds[i]), STD_WEND(stds[i]))
+	}
+	cv = NULL
+	call sf_fit (stds, nstds, cv, function, order, xmin, xmax)
+	call sf_rms (stds, nstds, rms, npts)
+
+	# If not interactive go to the output.
+	if (interactive == 3)
+	    goto output_
+	if (interactive != 4) {
+	    call printf ("Fit aperture %d interactively? ")
+		call pargi (aperture)
+	    interactive = clgwrd ("answer", cmd, SZ_FNAME, "|no|yes|NO|YES")
+	    switch (interactive) {
+	    case 1:
+		goto output_
+	    case 3:
+		call sf_gfree (gp)
+		goto output_
+	    }
+	}
+
+	# Initialize graphics structure parameters: airmass and wavelength
+	# limits and default images to plot.
+
+	if (gp == NULL)
+	    call sf_ginit (gp)
+	GP_AIRMASS(gp,1) = MAX_REAL
+	GP_AIRMASS(gp,2) = -MAX_REAL
+	j = 0
+	do i = 1, nstds - 2 {
+	    if (STD_FLAG(stds[i]) == SF_EXCLUDE)
+		next
+	    GP_AIRMASS(gp,1) = min (GP_AIRMASS(gp,1), STD_AIRMASS(stds[i]))
+	    GP_AIRMASS(gp,2) = max (GP_AIRMASS(gp,2), STD_AIRMASS(stds[i]))
+	    if (j < SF_NGRAPHS) {
+		j = j + 1
+	        call strcpy (STD_IMAGE(stds[i]), Memc[GP_IMAGES(gp,j)],
+		    SZ_FNAME)
+	        call strcpy (STD_SKY(stds[i]), Memc[GP_SKYS(gp,j)], SZ_FNAME)
+	    }
+	}
+	delta = GP_AIRMASS(gp,2) - GP_AIRMASS(gp,1)
+	GP_AIRMASS(gp,1) = GP_AIRMASS(gp,1) - 0.05 * delta
+	GP_AIRMASS(gp,2) = GP_AIRMASS(gp,2) + 0.05 * delta
+	GP_WSTART(gp) = xmin
+	GP_WEND(gp) = xmax
+	call sf_title (gp, aperture, function, order, npts, rms)
+
+	# Enter cursor loop by drawing the graphs.
+	key = 'r'
+	repeat {
+	    switch (key) {
+	    case '?':
+		call gpagefile (GP_GIO(gp), KEY, PROMPT)
+	    case ':':
+		call sf_colon (cmd, gp, stds, nstds, cv, wextn, extn, nextn,
+		    ecv, function, order, npts, rms, newfit, newgraph)
+	    case 'a':
+		call sf_add (gp, stds, nstds, cv, wx, wy, wc)
+	    case 'c':
+		call sf_composite (stds, nstds, cv)
+		newfit = YES
+		newgraph = YES
+	    case 'd':
+		call sf_data (stds, nstds, GP_GRAPHS(gp,wc))
+		call sf_nearest (gp, stds, nstds, wx, wy, wc, 0, i, j)
+		if (i > 0) {
+		    call printf (
+		        "%s - Delete p(oint), s(tar), or w(avelength):")
+		        call pargstr (STD_IMAGE(stds[i]))
+		    if (clgcur ("cursor", wx, wy, wc, key, cmd, SZ_FNAME)==EOF)
+		        break
+		    call printf ("\n")
+		    call sf_delete (gp, stds, nstds, key, i, j)
+		}
+	    case 'e':
+		call sf_extinct (gp, stds, nstds, cv, ecv, function, order)
+		newfit = YES
+		newgraph = YES
+	    case 'f':
+		newfit = YES
+	    case 'g':
+		newgraph = YES
+		newfit = YES
+	    case 'i':
+		call sf_data (stds, nstds, GP_GRAPHS(gp,wc))
+		call sf_nearest (gp, stds, nstds, wx, wy, wc, 2, i, j)
+		if (i > 0) {
+		    call printf (
+	    "%s: airmass=%6.3f wavelen=%6.3f sens=%6.3f fit=%6.3f weight=%3f")
+		        call pargstr (STD_IMAGE(stds[i]))
+			call pargr (STD_AIRMASS(stds[i]))
+		        call pargr (Memr[STD_WAVES(stds[i])+j-1])
+		        call pargr (Memr[STD_SENS(stds[i])+j-1])
+		        call pargr (Memr[STD_FIT(stds[i])+j-1])
+		        call pargr (Memr[STD_WTS(stds[i])+j-1])
+		}
+	    case 'm':
+		call sf_data (stds, nstds, GP_GRAPHS(gp,wc))
+		call sf_nearest (gp, stds, nstds, wx, wy, wc, 2, i, j)
+		if (i > 0) {
+		    call printf (
+		   "%s - Move p(oint), s(tar), or w(avelength) to cursor:")
+		        call pargstr (STD_IMAGE(stds[i]))
+		    if (clgcur ("cursor", wx, wy, wc, key, cmd, SZ_FNAME)==EOF)
+		        break
+		    call printf ("\n")
+		    delta = wy - Memr[STD_Y(stds[i])+j-1]
+		    call sf_move (gp, stds, nstds, key, i, j, delta)
+		}
+	    case 'o':
+		call sf_reset (stds, nstds, wextn, extn, nextn, ecv, shift)
+		newfit = YES
+		newgraph = YES
+	    case 'q':
+		break
+	    case 'I':
+		call fatal (0, "Interrupt")
+	    case 'r':
+		newgraph = YES
+	    case 's':
+		call sf_shift (stds, nstds, shift)
+		newfit=YES
+		newgraph=YES
+	    case 'u':
+		call sf_data (stds, nstds, GP_GRAPHS(gp,wc))
+		call sf_nearest (gp, stds, nstds, wx, wy, wc, 1, i, j)
+		if (i > 0) {
+		    call printf (
+		        "%s - Undelete p(oint), s(tar), or w(avelength):")
+		        call pargstr (STD_IMAGE(stds[i]))
+		    if (clgcur ("cursor", wx, wy, wc, key, cmd, SZ_FNAME)==EOF)
+		        break
+		    call printf ("\n")
+		    call sf_undelete (gp, stds, nstds, key, i, j)
+		}
+	    case 'w':
+		call sf_data (stds, nstds, GP_GRAPHS(gp,wc))
+		call sf_nearest (gp, stds, nstds, wx, wy, wc, 0, i, j)
+		if (i > 0) {
+		    call printf (
+		        "%s - Reweight p(oint), s(tar), or w(avelength):")
+		        call pargstr (STD_IMAGE(stds[i]))
+		    if (clgcur ("cursor", wx, wy, wc, key, cmd, SZ_FNAME)==EOF)
+		        break
+		    call printf ("New weight (%g):")
+			call pargr (Memr[STD_IWTS(stds[i])+j-1])
+		    call flush (STDOUT)
+		    if (scan() != EOF) {
+			call gargr (delta)
+			if (nscan() == 1)
+		    	    call sf_weights (stds, nstds, key, i, j, delta)
+		    }
+		    call printf ("\n")
+		}
+	    default:
+		call printf ("\007")
+	    }
+
+	    # Do a new fit and recompute the RMS, and title string.
+	    if (newfit == YES) {
+		call sf_fit (stds, nstds, cv, function, order, xmin, xmax)
+		call sf_rms (stds, nstds, rms, npts)
+		call sf_title (gp, aperture, function, order, npts, rms)
+		do i = 1, SF_NGRAPHS
+		    if (GP_SHDR(gp,i) != NULL)
+			call shdr_close (GP_SHDR(gp,i))
+	    }
+
+	    # Draw new graphs.
+	    if (newgraph == YES) {
+		call sf_graph (gp, stds, nstds, cv, wextn, extn, nextn, ecv)
+		newgraph = NO
+		newfit = YES
+	    }
+
+	    # Overplot new fit.
+	    if (newfit == YES) {
+		call sf_fitgraph (gp, cv)
+		newfit = NO
+	    }
+	} until (clgcur ("cursor", wx, wy, wc, key, cmd, SZ_FNAME) == EOF)
+
+	# Close any open images.
+	do i = 1, SF_NGRAPHS
+	    if (GP_SHDR(gp,i) != NULL)
+	        call shdr_close (GP_SHDR(gp,i))
+
+output_
+	# Output the sensitivity function and logfile statistics.
+	call sf_output (stds, nstds, cv, sensimage, ignoreaps)
+	if (logfile[1] != EOS) {
+	    iferr {
+	        fd = open (logfile, APPEND, TEXT_FILE)
+	        call sf_stats (fd, stds, nstds, function, order, npts, rms)
+	        call sf_vstats (fd, stds, nstds, cv, wextn, extn, nextn, ecv)
+	        call close (fd)
+	    } then
+		call erract (EA_WARN)
+	}
+	call cvfree (cv)
+end
diff --git a/noao/onedspec/sensfunc/sfshift.x b/noao/onedspec/sensfunc/sfshift.x
new file mode 100644
index 00000000..07b204f3
--- /dev/null
+++ b/noao/onedspec/sensfunc/sfshift.x
@@ -0,0 +1,81 @@
+include	"sensfunc.h"
+
+
+# SF_SHIFT -- Shift or unshift all standard stars to have zero mean residual.
+
+procedure sf_shift (stds, nstds, flag)
+
+pointer	stds[nstds]		# Standard star data
+int	nstds			# Number of standard stars
+int	flag			# Shift flag
+
+pointer	x, y, w, f
+int	i, j, n, nshift
+real	shift, shift1, minshift
+
+begin
+	# If flag is YES then unshift the data.
+	if (flag == YES) {
+	    do i = 1, nstds {
+	        if (STD_FLAG(stds[i]) == SF_EXCLUDE)
+		    next
+	        n = STD_NWAVES(stds[i])
+	        if (n == 0)
+		    next
+	        y = STD_SENS(stds[i])
+	        call asubkr (Memr[y], STD_SHIFT(stds[i]), Memr[y], n)
+	        STD_SHIFT(stds[i]) = 0.
+	    }
+	    flag = NO
+	    call printf ("Data unshifted")
+	    return
+	}
+
+	# Determine the shifts needed to make the mean residual zero.
+	# Also determine the minimum shift.
+
+	minshift = 0.
+	do i = 1, nstds {
+	    if (STD_FLAG(stds[i]) == SF_EXCLUDE)
+		next
+	    n = STD_NWAVES(stds[i])
+	    if (n == 0)
+		next
+	    x = STD_WAVES(stds[i])
+	    y = STD_SENS(stds[i])
+	    w = STD_WTS(stds[i])
+	    f = STD_FIT(stds[i])
+	    nshift = 0
+	    shift = 0.
+	    shift1 = 0.
+	    do j = 1, n {
+		shift1 = shift1 + Memr[f+j-1] - Memr[y+j-1]
+		if (Memr[w+j-1] > 0.) {
+		    shift = shift + Memr[f+j-1] - Memr[y+j-1]
+		    nshift = nshift + 1
+		}
+	    }
+	    if (nshift > 0) {
+		shift = STD_SHIFT(stds[i]) + shift / nshift
+		if (shift < minshift)
+		    minshift = shift
+	    } else
+		shift = STD_SHIFT(stds[i]) + shift1 / n
+	    STD_SHIFT(stds[i]) = shift
+	}
+
+	# Adjust the shifts to be upwards.
+	do i = 1, nstds {
+	    if (STD_FLAG(stds[i]) == SF_EXCLUDE)
+		next
+	    n = STD_NWAVES(stds[i])
+	    if (n == 0)
+		next
+	    y = STD_SENS(stds[i])
+	    shift = STD_SHIFT(stds[i]) - minshift
+	    call aaddkr (Memr[y], shift, Memr[y], n)
+	    STD_SHIFT(stds[i]) = shift
+	}
+	flag = YES
+	call printf ("Data shifted")
+end
diff --git a/noao/onedspec/sensfunc/sfstats.x b/noao/onedspec/sensfunc/sfstats.x
new file mode 100644
index 00000000..a94691a4
--- /dev/null
+++ b/noao/onedspec/sensfunc/sfstats.x
@@ -0,0 +1,152 @@
+include	"sensfunc.h"
+
+
+# SF_STATS -- Print basic statistics about the stars and the fit.
+
+procedure sf_stats (fd, stds, nstds, function, order, npts, rms)
+
+int	fd			# Output file descriptor (may be STDOUT)
+pointer	stds[nstds]		# Standard star data
+int	nstds			# Number of standard stars
+char	function[ARB]		# Fitted function
+int	order			# Order of function
+int	npts			# Number of points in fit
+real	rms			# RMS of fit
+
+int	i, j, n
+real	rms1, dev1, dev2, dev3
+pointer	sp, str, wts
+
+begin
+	# Start with system ID.
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+	call sysid (Memc[str], SZ_LINE)
+
+	# Determine beam from first standard star not excluded.
+	for (i=1; (i<nstds) && (STD_FLAG(stds[i])==SF_EXCLUDE); i=i+1)
+	    ;
+	call fprintf (fd, "%s\n")
+	    call pargstr (Memc[str])
+	call fprintf (fd, "Sensitivity function for aperture %d:\n")
+	    call pargi (STD_BEAM(stds[i]))
+	call fprintf (fd,
+        "Fitting function is %s of order %d with %d points and RMS of %6.4f.\n")
+	    call pargstr (function)
+	    call pargi (order)
+	    call pargi (npts)
+	    call pargr (rms)
+
+	call fprintf (fd, "%12s %7s %7s %7s %7s %7s %7s %7s\n")
+	    call pargstr ("Image")
+	    call pargstr ("Airmass")
+	    call pargstr ("Points")
+	    call pargstr ("Shift")
+	    call pargstr ("RMS Fit")
+	    call pargstr ("Dev 1")
+	    call pargstr ("Dev 2")
+	    call pargstr ("Dev 3")
+
+	do i = 1, nstds {
+	    if (STD_FLAG(stds[i]) == SF_EXCLUDE)
+		next
+
+	    n = 0
+	    wts = STD_WTS(stds[i]) - 1
+	    for (j=1; j<=STD_NWAVES(stds[i]); j=j+1)
+		if (Memr[wts+j] != 0.)
+		    n = n + 1
+	    if ((i == nstds-1) && (n == 0))
+		next
+
+	    call sf_devs (stds[i], rms1, dev1, dev2, dev3)
+
+	    call fprintf (fd, "%12s %7.3f %7d %7.4f %7.4f %7.4f %7.4f %7.4f")
+		call pargstr (STD_IMAGE(stds[i]))
+		call pargr (STD_AIRMASS(stds[i]))
+		call pargi (n)
+		call pargr (STD_SHIFT(stds[i]))
+		call pargr (rms1)
+		call pargr (dev1)
+		call pargr (dev2)
+		call pargr (dev3)
+
+	    if (n == 0) {
+	        call fprintf (fd, "%s")
+		    call pargstr (" <-- deleted")
+	    }
+	    call fprintf (fd, "\n")
+	}
+
+	# Trailing spacer
+	call fprintf (fd, "\n")
+end
+
+
+# SF_DEVS - Compute rms and mean deviations from the fit.
+# The deviations are computed in three segments.
+
+procedure  sf_devs (std, rms, dev1, dev2, dev3)
+
+pointer	std		# Standard star data
+real	rms		# RMS about fit
+real	dev1		# Average deviation in first third of data
+real	dev2		# Average deviation in second third of data
+real	dev3		# Average deviation in last third of data
+
+int	i, ndev1, ndev2, ndev3, nrms, nbin, nwaves
+real	dev
+pointer	sens, fit, wts
+
+begin
+	# Get elements froms standard star structure.
+	nwaves = STD_NWAVES(std)
+	sens = STD_SENS(std)
+	fit = STD_FIT(std)
+	wts = STD_WTS(std)
+
+	# Divide into thirds.
+	rms = 0.
+	ndev1 = 0
+	dev1 = 0.
+	nbin = nwaves / 3
+	for (i=1; i<= nbin; i=i+1)
+	    if (Memr[wts+i-1] != 0.) {
+		dev = Memr[sens+i-1] - Memr[fit+i-1]
+		dev1 = dev1 + dev
+		rms = rms + dev ** 2
+		ndev1 = ndev1 + 1
+	    }
+	if (ndev1 > 0)
+	    dev1 = dev1 / ndev1
+
+	ndev2 = 0
+	dev2 = 0.
+	nbin = 2 * nwaves / 3
+	for (; i<=nbin; i=i+1)
+	    if (Memr[wts+i-1] != 0.) {
+		dev = Memr[sens+i-1] - Memr[fit+i-1]
+		dev2 = dev2 + dev
+		rms = rms + dev ** 2
+		ndev2 = ndev2 + 1
+	    }
+	if (ndev2 > 0)
+	    dev2 = dev2 / ndev2
+
+	ndev3 = 0
+	dev3 = 0.
+	nbin = nwaves
+	for (; i<=nbin; i=i+1)
+	    if (Memr[wts+i-1] != 0.) {
+		dev = Memr[sens+i-1] - Memr[fit+i-1]
+		dev3 = dev3 + dev
+		rms = rms + dev ** 2
+		ndev3 = ndev3 + 1
+	    }
+	if (ndev3 > 0)
+	    dev3 = dev3 / ndev3
+
+	nrms = ndev1 + ndev2 + ndev3
+	if (nrms > 0)
+	    rms = sqrt (rms / nrms)
+end
diff --git a/noao/onedspec/sensfunc/sfstds.x b/noao/onedspec/sensfunc/sfstds.x
new file mode 100644
index 00000000..07219729
--- /dev/null
+++ b/noao/onedspec/sensfunc/sfstds.x
@@ -0,0 +1,266 @@
+include "sensfunc.h"
+
+
+# SF_STDS -- Get the standard observations for the specified apertures.
+# If ignoring aperture set all apertures to 1.
+# This routine knows the output of the task STANDARD.
+
+procedure sf_stds (standards, aps, ignoreaps, stds, nstds)
+
+char	standards	# Standard star data file
+pointer	aps		# Aperture list
+bool	ignoreaps	# Ignore apertures?
+pointer	stds		# Pointer to standard observations (returned)
+int	nstds		# Number of standard observations (returned)
+
+int	i, j, fd, beam, npts, nwaves, nalloc
+real	exptime, airmass, wstart, wend
+real	wavelength, flux, dwave, count
+pointer	sp, image, title, std
+pointer	waves, fluxes, dwaves, counts, sens, fit, wts, iwts, x, y
+
+bool	rng_elementi()
+int	open(), fscan(), nscan(), stridxs()
+errchk	open, malloc, realloc
+
+begin
+	call smark (sp)
+	call salloc (image, SZ_STDIMAGE, TY_CHAR)
+	call salloc (title, SZ_STDTITLE, TY_CHAR)
+
+	# Open the standard observation data file.
+	fd = open (standards, READ_ONLY, TEXT_FILE)
+
+	# Read the standard observations and create a structure for each one.
+	# The beginning of a new star is found by a line whose first word
+	# begins with the character '['.  Otherwise the line is interpreted
+	# as a data line.  All unrecognized formats are skipped.
+
+	nwaves = 0
+	nstds = 0
+	while (fscan (fd) != EOF) {
+	    call gargwrd (Memc[image], SZ_STDIMAGE)
+	    if (Memc[image] == '[') {
+	        call gargi (beam)
+	        call gargi (npts)
+		call gargr (exptime)
+	        call gargr (airmass)
+	        call gargr (wstart)
+	        call gargr (wend)
+	        call gargstr (Memc[title], SZ_STDTITLE)
+	        if (nscan() < 7)
+		    next
+		if (!rng_elementi (aps, beam))
+		    next
+		if (IS_INDEF (exptime) || exptime <= 0.) {
+		    call eprintf (
+		"%s: Warning - exposure time missing or zero, using 1 second\n")
+			call pargstr (Memc[image])
+		    exptime = 1.
+		}
+
+		# For the first one create the pointer to the array of
+		# structures.  For the following stars increase the size
+		# of the pointer array and finish up the previous standard
+		# star.
+
+	        if (nstds == 0) {
+		    nstds = nstds + 1
+		    call calloc (stds, nstds, TY_INT)
+		    call calloc (std, LEN_STD, TY_STRUCT)
+		    Memi[stds+nstds-1] = std
+	        } else {
+	    	    if (nwaves > 0) {
+			call realloc (waves, nwaves, TY_REAL)
+			call realloc (fluxes, nwaves, TY_REAL)
+			call realloc (dwaves, nwaves, TY_REAL)
+			call realloc (counts, nwaves, TY_REAL)
+			call realloc (wts, nwaves, TY_REAL)
+			call malloc (sens, nwaves, TY_REAL)
+			call malloc (fit, nwaves, TY_REAL)
+			call malloc (iwts, nwaves, TY_REAL)
+			call malloc (x, nwaves, TY_REAL)
+			call malloc (y, nwaves, TY_REAL)
+			call amovr (Memr[wts], Memr[iwts], nwaves)
+			STD_NWAVES(std) = nwaves
+			STD_WAVES(std) = waves
+			STD_FLUXES(std) = fluxes
+			STD_DWAVES(std) = dwaves
+			STD_COUNTS(std) = counts
+			STD_SENS(std) = sens
+			STD_FIT(std) = fit
+			STD_WTS(std) = wts
+			STD_IWTS(std) = iwts
+			STD_X(std) = x
+			STD_Y(std) = y
+
+			nstds = nstds + 1
+			call realloc (stds, nstds, TY_INT)
+			call calloc (std, LEN_STD, TY_STRUCT)
+			Memi[stds+nstds-1] = std
+		    }
+	        }
+
+		# Start a new standard star.
+	        std = Memi[stds+nstds-1]
+		if (ignoreaps)
+	            STD_BEAM(std) = 1
+		else
+	            STD_BEAM(std) = beam
+		STD_NPTS(std) = npts
+	        STD_EXPTIME(std) = exptime
+	        STD_AIRMASS(std) = airmass
+		STD_WSTART(std) = wstart
+		STD_WEND(std) = wend
+		STD_SHIFT(std) = 0.
+		STD_NWAVES(std) = 0
+
+		# Decode the image and sky strings.
+	        call strcpy (Memc[title], STD_TITLE(std), SZ_STDTITLE)
+		i = stridxs ("]", Memc[image])
+		if (Memc[image+i] == ']')
+		    i = i + 1
+		Memc[image+i-1] = EOS
+	        call strcpy (Memc[image+1], STD_IMAGE(std), SZ_STDIMAGE)
+		if (Memc[image+i] == '-') {
+		    i = i + 2
+		    j = stridxs ("]", Memc[image+i]) + i
+		    Memc[image+j-1] = EOS
+	            call strcpy (Memc[image+i], STD_SKY(std), SZ_STDIMAGE)
+		} else
+		    STD_SKY(std) = EOS
+	        nwaves = 0
+
+	    # Interprete the line as standard star wavelength point.
+	    } else if (nstds > 0) {
+		call reset_scan()
+		call gargr (wavelength)
+		call gargr (flux)
+		call gargr (dwave)
+		call gargr (count)
+		if (nscan() < 3)
+		    next
+		if (wavelength < min (wstart, wend) ||
+		    wavelength > max (wstart, wend) ||
+		    flux<=0. || dwave<=0. || count<=0.)
+		    next
+		if (!rng_elementi (aps, beam))
+		    next
+		nwaves = nwaves + 1
+
+		# Allocate in blocks to minimize the number of reallocations.
+		if (nwaves == 1) {
+		    nalloc = 100
+		    call malloc (waves, nalloc, TY_REAL)
+		    call malloc (fluxes, nalloc, TY_REAL)
+		    call malloc (dwaves, nalloc, TY_REAL)
+		    call malloc (counts, nalloc, TY_REAL)
+		    call malloc (wts, nalloc, TY_REAL)
+		} else if (nwaves > nalloc) {
+		    nalloc = nalloc + 100
+		    call realloc (waves, nalloc, TY_REAL)
+		    call realloc (fluxes, nalloc, TY_REAL)
+		    call realloc (dwaves, nalloc, TY_REAL)
+		    call realloc (counts, nalloc, TY_REAL)
+		    call realloc (wts, nalloc, TY_REAL)
+		}
+
+		# Record the data and compute the sensitivity.
+		Memr[waves+nwaves-1] = wavelength
+		Memr[fluxes+nwaves-1] = flux
+		Memr[dwaves+nwaves-1] = dwave
+		Memr[counts+nwaves-1] = count
+		Memr[wts+nwaves-1] = 1.
+	    }
+	}
+
+	# Finish up the last standard star and close the file.
+	if (nstds > 0) {
+	    STD_NWAVES(std) = nwaves
+    	    if (nwaves > 0) {
+		call realloc (waves, nwaves, TY_REAL)
+		call realloc (fluxes, nwaves, TY_REAL)
+		call realloc (dwaves, nwaves, TY_REAL)
+		call realloc (counts, nwaves, TY_REAL)
+		call realloc (wts, nwaves, TY_REAL)
+		call malloc (sens, nwaves, TY_REAL)
+		call malloc (fit, nwaves, TY_REAL)
+		call malloc (iwts, nwaves, TY_REAL)
+		call malloc (x, nwaves, TY_REAL)
+		call malloc (y, nwaves, TY_REAL)
+		call amovr (Memr[wts], Memr[iwts], nwaves)
+		STD_WAVES(std) = waves
+		STD_FLUXES(std) = fluxes
+		STD_DWAVES(std) = dwaves
+		STD_COUNTS(std) = counts
+		STD_SENS(std) = sens
+		STD_FIT(std) = fit
+		STD_WTS(std) = wts
+		STD_IWTS(std) = iwts
+		STD_X(std) = x
+		STD_Y(std) = y
+	    }
+	}
+	call close (fd)
+	call sfree (sp)
+
+	# Add standard stars for any added points and composite points.
+	nstds = nstds + 2
+	call realloc (stds, nstds, TY_INT)
+	call calloc (std, LEN_STD, TY_STRUCT)
+	Memi[stds+nstds-2] = std
+	call strcpy ("Added", STD_IMAGE(std), SZ_STDIMAGE)
+	STD_BEAM(std) = STD_BEAM(Memi[stds])
+	STD_NPTS(std) = STD_NPTS(Memi[stds])
+	STD_EXPTIME(std) = 1.
+	STD_AIRMASS(std) = 1.
+	STD_WSTART(std) = STD_WSTART(Memi[stds])
+	STD_WEND(std) = STD_WEND(Memi[stds])
+	STD_SHIFT(std) = 0.
+	STD_NWAVES(std) = 0
+	call calloc (std, LEN_STD, TY_STRUCT)
+	Memi[stds+nstds-1] = std
+	call strcpy ("Composite", STD_IMAGE(std), SZ_STDIMAGE)
+	STD_BEAM(std) = STD_BEAM(Memi[stds])
+	STD_NPTS(std) = STD_NPTS(Memi[stds])
+	STD_EXPTIME(std) = 1.
+	STD_AIRMASS(std) = 1.
+	STD_WSTART(std) = STD_WSTART(Memi[stds])
+	STD_WEND(std) = STD_WEND(Memi[stds])
+	STD_SHIFT(std) = 0.
+	STD_NWAVES(std) = 0
+end
+
+
+# SF_FREE -- Free the standard observations and aperture array.
+
+procedure sf_free (stds, nstds, apertures, napertures)
+
+pointer	stds		# Pointer to standard observations
+int	nstds		# Number of standard observations
+pointer	apertures	# Pointer to apertures array
+int	napertures	# Number of apertures
+
+int	i
+pointer	std
+
+begin
+	do i = 1, nstds {
+	    std = Memi[stds+i-1]
+	    if (STD_NWAVES(std) > 0) {
+		call mfree (STD_WAVES(std), TY_REAL)
+		call mfree (STD_FLUXES(std), TY_REAL)
+		call mfree (STD_DWAVES(std), TY_REAL)
+		call mfree (STD_COUNTS(std), TY_REAL)
+		call mfree (STD_SENS(std), TY_REAL)
+		call mfree (STD_FIT(std), TY_REAL)
+		call mfree (STD_WTS(std), TY_REAL)
+		call mfree (STD_IWTS(std), TY_REAL)
+		call mfree (STD_X(std), TY_REAL)
+		call mfree (STD_Y(std), TY_REAL)
+	    }
+	    call mfree (std, TY_STRUCT)
+	}
+	call mfree (stds, TY_INT)
+	call mfree (apertures, TY_INT)
+end
diff --git a/noao/onedspec/sensfunc/sftitle.x b/noao/onedspec/sensfunc/sftitle.x
new file mode 100644
index 00000000..50cece9a
--- /dev/null
+++ b/noao/onedspec/sensfunc/sftitle.x
@@ -0,0 +1,23 @@
+include	"sensfunc.h"
+
+
+# SF_TITLE -- Make title string for graphs.
+
+procedure sf_title (gp, aperture, function, order, npts, rms)
+
+pointer	gp
+int	aperture
+char	function[ARB]
+int	order
+int	npts
+real	rms
+
+begin
+	call sprintf (GP_TITLE(gp), GP_SZTITLE,
+	    "Aperture=%d  Function=%s  Order=%d  Points=%d  RMS=%6.4f")
+	    call pargi (aperture)
+	    call pargstr (function)
+	    call pargi (order)
+	    call pargi (npts)
+	    call pargr (rms)
+end
diff --git a/noao/onedspec/sensfunc/sfundelete.x b/noao/onedspec/sensfunc/sfundelete.x
new file mode 100644
index 00000000..25161cc3
--- /dev/null
+++ b/noao/onedspec/sensfunc/sfundelete.x
@@ -0,0 +1,141 @@
+include	<gset.h>
+include "sensfunc.h"
+
+
+# SF_UNDELETE -- Unelete point, star, or wavelength.
+
+procedure sf_undelete (gp, stds, nstds, key, istd, ipt)
+
+pointer	gp			# GIO pointer
+pointer	stds[nstds]		# Standard star data
+int	nstds			# Number of standard stars
+int	key			# Delete point, star, or wavelength
+int	istd			# Index of standard star
+int	ipt			# Index of point
+
+int	i, j, n, wcs, mark, mdel, color, stridx()
+real	wave, szmark, szmdel
+pointer	x, y, z, w, w1, gio
+
+begin
+	gio = GP_GIO(gp)
+	mdel = GP_MDEL(gp)
+	szmdel = GP_SZMDEL(gp)
+	szmark = GP_SZMARK(gp)
+
+	# Undelete points from each displayed graph.
+	for (wcs = 1; GP_GRAPHS(gp,wcs) != EOS; wcs = wcs + 1) {
+	    if (stridx (GP_GRAPHS(gp,wcs), "ars") == 0)
+	        next
+
+	    call gseti (gio, G_WCS, wcs)
+	    call gseti (gio, G_PMLTYPE, 0)
+	    call sf_data (stds, nstds, GP_GRAPHS(gp,wcs))
+	    switch (key) {
+	    case 'p':
+		if (istd != nstds-1) {
+		    mark = GP_MARK(gp)
+		    color = GP_CMARK(gp)
+		} else {
+		    mark = GP_MADD(gp)
+		    color = GP_CADD(gp)
+		}
+	        x = STD_X(stds[istd])+ipt-1
+	        y = STD_Y(stds[istd])+ipt-1
+	        w = STD_WTS(stds[istd])+ipt-1
+	        w1 = STD_IWTS(stds[istd])+ipt-1
+		call gseti (gio, G_PMLTYPE, 0)
+	        call gmark (gio, Memr[x], Memr[y], mdel, szmdel, szmdel)
+		call gseti (gio, G_PMLTYPE, 1)
+		call gseti (gio, G_PLCOLOR, color)
+	        call gmark (gio, Memr[x], Memr[y], mark, szmark , szmark)
+	    case 's':
+		if (istd != nstds-1) {
+		    mark = GP_MARK(gp)
+		    color = GP_CMARK(gp)
+		} else {
+		    mark = GP_MADD(gp)
+		    color = GP_CADD(gp)
+		}
+	        n = STD_NWAVES(stds[istd])
+	        x = STD_X(stds[istd])
+	        y = STD_Y(stds[istd])
+	        w = STD_WTS(stds[istd])
+		do j = 1, n {
+		    if (Memr[w] == 0.) {
+			call gseti (gio, G_PMLTYPE, 0)
+	        	call gmark (gio, Memr[x], Memr[y], mdel, szmdel, szmdel)
+			call gseti (gio, G_PMLTYPE, 1)
+			call gseti (gio, G_PLCOLOR, color)
+	        	call gmark (gio, Memr[x], Memr[y], mark, szmark, szmark)
+		    }
+		    x = x + 1
+		    y = y + 1
+		    w = w + 1
+	        }
+	    case 'w':
+	        wave = Memr[STD_WAVES(stds[istd])+ipt-1]
+	        do i = 1, nstds {
+	            if (STD_FLAG(stds[i]) != SF_INCLUDE)
+		        next
+		    if (i != nstds-1) {
+			mark = GP_MARK(gp)
+			color = GP_CMARK(gp)
+		    } else {
+			mark = GP_MADD(gp)
+			color = GP_CADD(gp)
+		    }
+	            n = STD_NWAVES(stds[i])
+	            x = STD_X(stds[i])
+	            y = STD_Y(stds[i])
+		    z = STD_WAVES(stds[i])
+	            w = STD_WTS(stds[i])
+		    do j = 1, n {
+		        if ((Memr[z] == wave) && (Memr[w] == 0.)) {
+			    call gseti (gio, G_PMLTYPE, 0)
+	    		    call gmark (gio, Memr[x], Memr[y], mdel, szmdel,
+				szmdel)
+			    call gseti (gio, G_PMLTYPE, 1)
+			    call gseti (gio, G_PLCOLOR, color)
+	    		    call gmark (gio, Memr[x], Memr[y], mark, szmark,
+				szmark)
+		        }
+			x = x + 1
+			y = y + 1
+			z = z + 1
+			w = w + 1
+		    }
+	        }
+	    }
+	}
+
+	# Now actually undelete the points by resetting the weights.
+	switch (key) {
+	case 'p':
+	    w = STD_WTS(stds[istd])+ipt-1
+	    w1 = STD_IWTS(stds[istd])+ipt-1
+	    Memr[w] = Memr[w1]
+	case 's':
+	    n = STD_NWAVES(stds[istd])
+	    w = STD_WTS(stds[istd])
+	    w1 = STD_IWTS(stds[istd])
+	    call amovr (Memr[w1], Memr[w], n)
+	case 'w':
+	    wave = Memr[STD_WAVES(stds[istd])+ipt-1]
+	    do i = 1, nstds {
+	        if (STD_FLAG(stds[i]) != SF_INCLUDE)
+		    next
+	        n = STD_NWAVES(stds[i])
+		z = STD_WAVES(stds[i])
+	        w = STD_WTS(stds[i])
+	        w1 = STD_IWTS(stds[i])
+		do j = 1, n {
+		    if (Memr[z] == wave)
+			Memr[w] = Memr[w1]
+		    z = z + 1
+		    w = w + 1
+		    w1 = w1 + 1
+		}
+	    }
+	}
+end
diff --git a/noao/onedspec/sensfunc/sfvstats.x b/noao/onedspec/sensfunc/sfvstats.x
new file mode 100644
index 00000000..add49da7
--- /dev/null
+++ b/noao/onedspec/sensfunc/sfvstats.x
@@ -0,0 +1,104 @@
+include	"sensfunc.h"
+
+# SF_VSTATS -- Verbose statistics output.
+
+procedure sf_vstats (fd, stds, nstds, cv, wextn, extn, nextn, ecv)
+
+int	fd		# Output file descriptor (may be STDOUT)
+pointer	stds[nstds]	# Standard star data
+int	nstds		# Number of standard stars
+pointer	cv		# Sensitivity function curve
+real	wextn[nextn]	# Extinction table wavelength
+real	extn[nextn]	# Extinction table values
+int	nextn		# Number of extinction table values
+pointer	ecv		# Residual extinction curve
+
+int	i, j, n, nwaves
+real	w, fit, ext, dext, cveval()
+double	sum, sum2, s
+pointer	sp, waves, sens, xp, yp, zp
+
+begin
+	nwaves = 0
+	do i = 1, nstds-1
+	    if (STD_FLAG(stds[i]) != SF_EXCLUDE)
+	        nwaves = nwaves + STD_NWAVES(stds[i])
+
+	call smark (sp)
+	call salloc (waves, nwaves, TY_REAL)
+	call salloc (sens, nwaves, TY_REAL)
+	    
+	nwaves = 0
+	do i = 1, nstds-1 {
+	    if (STD_FLAG(stds[i]) == SF_EXCLUDE)
+		next
+	    n = STD_NWAVES(stds[i])
+	    xp = STD_WAVES(stds[i])
+	    yp = STD_SENS(stds[i])
+	    zp = STD_WTS(stds[i])
+	    do j = 1, n {
+		if (Memr[zp] != 0.) {
+		    Memr[waves+nwaves] = Memr[xp]
+		    Memr[sens+nwaves] = Memr[yp]
+		    nwaves = nwaves + 1
+		}
+		xp = xp + 1
+		yp = yp + 1
+		zp = zp + 1
+	    }
+	}
+	call xt_sort2 (Memr[waves], Memr[sens], nwaves)
+
+	call fprintf (fd, "%8s %7s %7s %7s %7s %5s %7s %7s\n")
+	    call pargstr ("Lambda")
+	    call pargstr ("Fit")
+	    call pargstr ("Avg")
+	    call pargstr ("Resid")
+	    call pargstr ("SD Avg")
+	    call pargstr ("N")
+	    call pargstr ("Ext")
+	    call pargstr ("Dext")
+
+	dext = 0.
+	n = 0
+	sum = 0.
+	sum2 = 0.
+	do i = 0, nwaves-1 {
+	    w = Memr[waves+i]
+	    s = Memr[sens+i]
+	    n = n + 1
+	    sum = sum + s
+	    sum2 = sum2 + s * s
+
+	    if ((i < nwaves-1) && (w == Memr[waves+i+1]))
+		next
+
+	    sum = sum / n
+	    sum2 = sum2 / n - sum * sum
+	    if (sum2 > 0)
+	        sum2 = sqrt (sum2 / n)
+	    else
+		sum2 = 0.
+	    fit = cveval (cv, w)
+	    call intrp (1, wextn, extn, nextn, w, ext, j)
+	    if (ecv != NULL)
+	        dext = cveval (ecv, w)
+	    call fprintf (fd, "%8.2f %7.3f %7.3f %7.4f %7.4f %5d %7.4f %7.4f\n")
+		call pargr (w)
+		call pargr (fit)
+		call pargd (sum)
+		call pargd (sum - fit)
+		call pargd (sum2)
+		call pargi (n)
+		call pargr (ext)
+		call pargr (dext)
+	    n = 0
+	    sum = 0.
+	    sum2 = 0.
+	}
+
+	# Trailing spacer
+	call fprintf (fd, "\n")
+
+	call sfree (sp)
+end
diff --git a/noao/onedspec/sensfunc/sfweights.x b/noao/onedspec/sensfunc/sfweights.x
new file mode 100644
index 00000000..2ce24b1a
--- /dev/null
+++ b/noao/onedspec/sensfunc/sfweights.x
@@ -0,0 +1,51 @@
+include "sensfunc.h"
+
+
+# SF_WEIGHTS -- Change weights for point, star, or wavelength.
+# The original input weight is permanently lost.
+
+procedure sf_weights (stds, nstds, key, istd, ipt, weight)
+
+pointer	stds[nstds]		# Standard star data
+int	nstds			# Number of standard stars
+int	key			# Delete point, star, or wavelength
+int	istd			# Index of standard star
+int	ipt			# Index of point
+real	weight			# New weight
+
+int	i, j, n
+real	wave
+pointer	z, w, iw
+
+begin
+	switch (key) {
+	case 'p':
+	    Memr[STD_WTS(stds[istd])+ipt-1] = weight
+	    Memr[STD_IWTS(stds[istd])+ipt-1] = weight
+	case 's':
+	    n = STD_NWAVES(stds[istd])
+	    w = STD_WTS(stds[istd])
+	    iw = STD_IWTS(stds[istd])
+	    call amovkr (weight, Memr[w], n)
+	    call amovkr (weight, Memr[iw], n)
+	case 'w':
+	    wave = Memr[STD_WAVES(stds[istd])+ipt-1]
+	    do i = 1, nstds {
+	        if (STD_FLAG(stds[i]) != SF_INCLUDE)
+		    next
+	        n = STD_NWAVES(stds[i])
+		z = STD_WAVES(stds[i])
+	        w = STD_WTS(stds[i])
+	        iw = STD_IWTS(stds[i])
+		do j = 1, n {
+		    if (Memr[z] == wave) {
+	                Memr[w] = weight
+	                Memr[iw] = weight
+		    }
+		    w = w + 1
+		    iw = iw + 1
+		    z = z + 1
+		}
+	    }
+	}
+end
diff --git a/noao/onedspec/sensfunc/t_sensfunc.x b/noao/onedspec/sensfunc/t_sensfunc.x
new file mode 100644
index 00000000..82f18678
--- /dev/null
+++ b/noao/onedspec/sensfunc/t_sensfunc.x
@@ -0,0 +1,99 @@
+include	"sensfunc.h"
+
+
+# T_SENSFUNC -- Determine sensitivities and residual extinctions.
+# The input is a file of standard star produced by the task STANDARD.
+# The input data is read into an array of structures, one per standard
+# star.  The stars common to the aperture to be fit are flagged
+# and then the data is passed to the main routine SF_SENSFUNC.
+# This routine determines a sensitivity curve for the aperture which
+# is output by the procedure as well as some optional statistical
+# information.  It returns an optional residual extinction curve
+# for each aperture.  The residual extinctions curves are finally combined
+# and output as a revised extinction table.
+
+procedure t_sensfunc ()
+
+pointer	standards		# Input standard star data filename
+pointer	sensitivity		# Output root sensitivity function imagename
+pointer	aps			# Aperture list
+bool	ignoreaps		# Ignore apertures?	
+pointer	logfile			# Output log for statistics
+pointer	function		# Sensitivity function type
+int	order			# Order of sensitivity function
+int	interactive		# Interactive?
+
+int	i, j, aperture, nstds, napertures, nextn, clgeti()
+pointer	sp, str, stds, apertures, wextn, extn, ecvs, gp
+bool	clgetb()
+pointer	rng_open()
+errchk	sf_sensfunc
+
+begin
+	call smark (sp)
+	call salloc (standards, SZ_FNAME, TY_CHAR)
+	call salloc (sensitivity, SZ_FNAME, TY_CHAR) 
+	call salloc (str, SZ_LINE, TY_CHAR) 
+	call salloc (logfile, SZ_FNAME, TY_CHAR) 
+	call salloc (function, SZ_FNAME, TY_CHAR)
+
+	# CL parameter input.
+	call clgstr ("standards", Memc[standards], SZ_FNAME)
+	call clgstr ("sensitivity", Memc[sensitivity], SZ_FNAME) 
+	call clgstr ("apertures", Memc[str], SZ_LINE)
+	ignoreaps = clgetb ("ignoreaps")
+	call clgstr ("logfile", Memc[logfile], SZ_FNAME) 
+	call clgstr ("function", Memc[function], SZ_FNAME)
+	order = clgeti ("order")
+	if (clgetb ("interactive"))
+	    interactive = 2
+	else
+	    interactive = 3
+
+	# Decode aperture list.
+	iferr (aps = rng_open (Memc[str], INDEF, INDEF, INDEF))
+	    call error (0, "Bad aperture list")
+
+	# Get the standard star data, the aperture array, and the
+	# extinction table, and allocate and initialize an array of
+	# residual extinction curves for each aperture.
+
+	call sf_stds (Memc[standards], aps, ignoreaps, stds, nstds)
+	if (nstds == 0) {
+	    call sfree (sp)
+	    return
+	}
+	call sf_apertures (Memi[stds], nstds, apertures, napertures)
+	call ext_load (wextn, extn, nextn)
+	call salloc (ecvs, napertures, TY_INT)
+	call amovki (NULL, Memi[ecvs], napertures)
+
+	# For each aperture flag standard stars to be used and call sf_sensfunc.
+	gp = NULL
+	do j = 1, napertures {
+	    aperture = Memi[apertures+j-1]
+	    do i = 1, nstds - 2
+		if (STD_BEAM(Memi[stds+i-1]) == aperture)
+		    STD_FLAG(Memi[stds+i-1]) = SF_INCLUDE
+		else
+		    STD_FLAG(Memi[stds+i-1]) = SF_EXCLUDE
+	    call sf_sensfunc (gp, Memi[stds], nstds, Memr[wextn], Memr[extn],
+		nextn, Memc[sensitivity], Memc[logfile], Memi[ecvs+j-1],
+		Memc[function], order, ignoreaps, interactive)
+	}
+	call sf_gfree (gp)
+
+	# Output a revised extinction table by combining the residual
+	# extinction curves for the apertures.  The table name is obtained
+	# by this proceudre.
+
+	call sf_eout (Memr[wextn], Memr[extn], nextn, Memi[ecvs], napertures)
+
+	# Finish up.
+	call sf_free (stds, nstds, apertures, napertures)
+	call ext_free (wextn, extn)
+	do j = 1, napertures
+	    call cvfree (Memi[ecvs+j-1])
+	call rng_close (aps)
+	call sfree (sp)
+end