Initial commit

13 files changed, 5765 insertions, 0 deletions

diff --git a/noao/obsutil/src/starfocus/Revisions b/noao/obsutil/src/starfocus/Revisions
new file mode 100644
index 00000000..d3f93e90
--- /dev/null
+++ b/noao/obsutil/src/starfocus/Revisions
@@ -0,0 +1,162 @@
+.help revisions Nov01 obsutil
+.nf
+t_starfocus.x
+starfocus.par
+psfmeasure.par
+starfocus.hlp
+psfmeasure.hlp
+    Added a "wcs" parameter to allow specifying the coordinate system type
+    for cursor input.  (10/7/98, Valdes)
+
+stfprofile.x
+    The logic in STF_FIT for determining the points to fit and the point
+    to use for the initial width estimate was faulty allowing some bad
+    cases to get through.  (7/31/98, Valdes)
+
+t_starfocus.x
+    There was an error in the changes to allow focus sequences to go
+    up or down such that simple PSF measurement failed.
+    (6/11/97, Valdes)
+
+stfmeasure.x
+stfprofile.x
+    Changes were need to better deal with error conditions.
+    (6/11/97, Valdes)
+
+stfmeasure.x
+    Changes to deal with error conditions.
+    (6/11/97, Valdes)
+
+t_starfocus.x
+    Fixed a typo error in an earlier change which left out the pointer
+    variable in a macro.  (5/8/97, Valdes)
+
+t_starfocus.x
+    1.  Added a minimum radius of 3 to the input and to interative setting.
+    2.  The estimate for the next focus position is not based on the previous
+	center plus the step.
+    3.  Added commented out code to allow setting the selected star to always
+	be the top or bottom, etc.
+    (2/7/97, Valdes)
+
+stfprofile.x
+    1.  Added a minimum radius of 3 to stf_find.
+    2.  Added a error return if sum2 <= 0 in stf_find
+    3.  Added errchks in stf_widths for called procedures.
+    4.  Added error returns in stf_fit for bad profile or parameters.
+    (2/6/97, Valdes)
+
+t_starfocus.x
+    The parabola fitting routine would fail if the independent variable
+    (the focus values) became too large because a) use of real and b)
+    unscaled variables.  The routine was revised for both these problems.
+    (10/24/96, Valdes)
+
+stfprofile.x
+stfmeasure.x
+    Fixed bug in evaluation of enclosed flux profile in which the scaled
+    radius was used for the gaussian subtraction stage instead of pixels.
+    This does not currently affect IMEXAM because the scale is fixed
+    at 1.  (8/29/96, Valdes)
+
+stfmeasure.x
+stfprofile.x
+    Fixed incorrect datatype declaration "real np" -> "int np" in various
+    related places.  (4/9/96, Valdes)
+
+stfmeasure.x -
+    1. Restricted the peak pixel to be within the specified radius.
+    2. Renamed file and procedures to avoid library conflict.
+    (3/14/96, Valdes)
+
+stfprofile.x
+    Minor bug fix that does not affect anything.  (3/14/96, Valdes)
+
+kpnofocus.cl
+starfocus.par
+psfmeasure.par
+    Changed the defaults for radius=5, sbuffer=5, swidth=5 to try and make
+    IMEXAMINE and the PSF tasks measure the same thing by default.
+    (3/13/96, Valdes)
+
+stfprofile.x
+    Changed the centering routine to only use the data in the specified
+    radius rather than the extended region including the sky.  (11/7/95, Valdes)
+
+stfmeasure.x
+mkpkg
+    Added a routine which can be called separately to compute the
+    enclosed flux width or radius and the direct FWHM.  (11/6/95, Valdes)
+
+stfprofile.x
+starfocus.h
+    Added a computation of the direct FWHM.   (11/6/95, Valdes)
+
+t_starfocus.x
+    Make a change to avoid reloading the image display when the specified
+    image name and the image name in the display are same apart from
+    the image extension.  (6/26/95, Valdes)
+
+stfprofile.h
+t_starfocus.x
+stfprofile.x
+starfocus.par
+psfmeasure.par
+kpnofocus.cl
+starfocus.par
+psfmeasure.par
+kpnofocus.cl
+    Added saturation and ignore_sat parameters to allow flagging measurements
+    with saturated pixels and either ignoring them or including them but
+    showing them in the output log.  (10/28/94, Valdes)
+
+t_starfocus.x
+    Fixed bug in interpreting the focus parameter.  (9/15/94, Valdes)
+
+starfocus.h
+t_starfocus.x
+stfprofile.x
+stfgraph.x
+psfmeasure.par
+starfocus.par
+kpnofocus.cl
+psfmeasure.hlp
+starfocus.hlp
+kpnofocus.hlp
+    Added an "iteration" parameter that, if greater than 1, uses the
+    previous FWHM estimate to adjust the "radius" parameter.  (8/5/94, Valdes)
+
+psfmeasure.hlp
+starfocus.hlp
+    Clarified the description of the imagecur parameter.  (2/1/94, Valdes)
+
+t_starfocus.x
+starfocus.h
+starfocus.par
+psfmeasure.par
+kpnofocus.cl	+
+starfocus.hlp
+psfmeasure.hlp
+kpnofocus.hlp
+    1.  A new parameter, fstep, was added to STARFOCUS to allow specifying
+	the focus sequence as a starting value, by the focus parameters,
+	and a step.
+    2.  STARFOCUS was modified to allow specification of header keywords
+	for the focus (this was true previously), the focus step,
+	the number of exposures, and the multiple exposure shift.
+	This allows multiple exposure images to be completely header
+	drivien if the appropriate keywords are present.
+    3.  A new task KPNOFOCUS was added.  This is a script calling
+	STARFOCUS with parameters set specifically for Kitt Peak
+	headers containing the focus parameters in the header.
+	Many of the parameters are fixed in the script and the task
+	parameters are then simpler.
+    4.  STARFOCUS/PSFMEASURE were modified to search all display frames
+	for the requested image rather than just frame 1.
+    5.  A new parameters, frame, was added to STARFOCUS/PSFMEASURE
+	to specify the display frame to load if necessary.  Previously
+	it was always frame 1.
+    7.  In STARFOCUS the default number of steps was changed to 7 and
+	the default step size to 30.
+    (11/13/93, Valdes)
+.endhelp
diff --git a/noao/obsutil/src/starfocus/mkpkg b/noao/obsutil/src/starfocus/mkpkg
new file mode 100644
index 00000000..3feffc76
--- /dev/null
+++ b/noao/obsutil/src/starfocus/mkpkg
@@ -0,0 +1,22 @@
+# Make the STARFOCUS tasks.
+
+$checkout libpkg.a ../
+$update   libpkg.a
+$checkin  libpkg.a ../
+$exit
+
+standalone:
+	$update	libpkg.a
+	$omake	x_starfocus.x
+	$link	x_starfocus.o libpkg.a -lxtools -lnlfit -liminterp -lgsurfit\
+	    -lds -o xx_starfocus.e
+	;
+
+libpkg.a:
+	stfgraph.x	starfocus.h <error.h> <gset.h> <mach.h>
+	stfmeasure.x	starfocus.h <error.h> <imhdr.h> <imset.h>\
+			<math/nlfit.h>
+	stfprofile.x	starfocus.h <imhdr.h> <mach.h> <math.h>\
+			<math/iminterp.h> <math/nlfit.h>
+	t_starfocus.x	starfocus.h <error.h> <imhdr.h> <imset.h> <mach.h>
+	;
diff --git a/noao/obsutil/src/starfocus/psfhelp.key b/noao/obsutil/src/starfocus/psfhelp.key
new file mode 100644
index 00000000..9e617a62
--- /dev/null
+++ b/noao/obsutil/src/starfocus/psfhelp.key
@@ -0,0 +1,60 @@
+			PSFMEASURE COMMAND OPTIONS
+
+				SUMMARY
+
+? Help              m Magnitude         r Redraw            z Zoom
+a Spatial           n Normalize         s Mag symbols       <space> Next
+d Delete            o Offset            t Field radius      
+e Enclosed flux     p Radial profile    u Undelete          
+i Info              q Quit              x Delete            
+
+:level <val>        :radius <val>       :show <file>        :xcenter <val>
+:overplot <y|n>     :scale <val>        :size <type>        :ycenter <val>
+
+
+				CURSOR COMMANDS
+
+All plots may not be available depending on the number of stars.
+
+?	Page this help information
+a	Spatial plot
+d	Delete star nearest to cursor
+e	Enclosed flux for all stars
+i	Information about star nearest the cursor
+m	Size and ellipticity vs relative magnitude
+n	Normalize enclosed flux at x cursor position
+o	Offset enclosed flux to x,y cursor position by adjusting background
+p	Radial profiles for all stars
+          Profiles are determined from the derivatives of the enclosed flux
+q	Quit
+r	Redraw
+s	Toggle magnitude symbols in spatial plot
+t	Size and ellipticity vs radius from field center
+u	Undelete all deleted points
+x	Delete nearest point or star (selected by query)
+z	Zoom to a single measurement showing enclosed flux and radial profile
+<space> Step through different stars in some plots
+
+
+			COLON COMMANDS
+
+A command without a value generally shows the current value of the
+parameter while with a value it sets the value of the parameter.
+
+:level <val>	Level at which the size parameter is evaluated
+:overplot <y|n> Overplot the profiles from the narrowest profile?
+:radius <val>   Change profile radius(*)
+:show <file>	Page all information for the current set of objects
+:size <type>	Size type (Radius|FWHM|GFWHM|MFWHM) (**)
+:scale <val>	Pixel scale for size values
+:xcenter <val>	X field center for radius from field center plots
+:ycenter <val>	Y field center for radius from field center plots
+
+(*) The profile radius may not exceed the initial value set by the task
+    parameter.
+
+(**)
+Radius = radius enclosing the fraction of the flux specified by "level"
+FWHM   = Full-width at half-maximum based on the radially smoothed profile
+GFWHM  = Full-width at half-maximum of Gaussian function fit to enclosed flux
+MFWHM  = Full-width at half-maximum of Moffat function fit to enclosed flux
diff --git a/noao/obsutil/src/starfocus/psfmeasure.par b/noao/obsutil/src/starfocus/psfmeasure.par
new file mode 100644
index 00000000..16dea3e9
--- /dev/null
+++ b/noao/obsutil/src/starfocus/psfmeasure.par
@@ -0,0 +1,24 @@
+# PSFMEASURE -- Measure PSF.
+
+images,s,a,,,,"List of images"
+coords,s,h,"markall","center|mark1|markall",,"Object coordinates"
+wcs,s,h,"logical","logical|physical|world",,"Coordinate system"
+display,b,h,yes,,,"Display images?"
+frame,i,h,1,1,,"Display frame to use
+"
+level,r,h,0.5,,,"Measurement level (fraction or percent)"
+size,s,h,"FWHM","Radius|FWHM|GFWHM|MFWHM",,"Size to display"
+beta,r,h,INDEF,2.1,,Moffat beta parameter
+scale,r,h,1.,,,"Pixel scale"
+radius,r,h,5.,,,"Measurement radius (pixels)"
+sbuffer,r,h,5.,,,"Sky buffer (pixels)"
+swidth,r,h,5.,,,"Sky width (pixels)"
+saturation,r,h,INDEF,,,"Saturation level"
+ignore_sat,b,h,no,,,"Ignore objects with saturated pixels?"
+iterations,i,h,2,1,,"Number of radius adjustment iterations"
+xcenter,r,h,INDEF,,,X field center (pixels)
+ycenter,r,h,INDEF,,,X field center (pixels)
+logfile,s,h,"logfile",,,"Logfile
+"
+imagecur,*imcur,h,"",,,"Image cursor input"
+graphcur,*gcur,h,"",,,"Graphics cursor input"
diff --git a/noao/obsutil/src/starfocus/starfocus.h b/noao/obsutil/src/starfocus/starfocus.h
new file mode 100644
index 00000000..871e8c98
--- /dev/null
+++ b/noao/obsutil/src/starfocus/starfocus.h
@@ -0,0 +1,140 @@
+# STARFOCUS
+
+# Types of coordinates
+define	SF_TYPES	"|center|mark1|markall|"
+define	SF_CENTER	1		# Star at center of image
+define	SF_MARK1	2		# Mark stars in first image
+define	SF_MARKALL	3		# Mark stars in all images
+
+# Task type
+define	STARFOCUS	1
+define	PSFMEASURE	2
+
+# Radius types
+define	SF_WTYPES	"|Radius|FWHM|GFWHM|MFWHM|"
+
+define	SF_RMIN		16		# Minimum centering search radius
+define	MAX_FRAMES	8		# Maximum number of display frames
+
+# Data structures for STARFOCUS
+
+define	NBNDRYPIX	100		# Number of boundary pixels
+define	TYBNDRY		BT_REFLECT	# Type of boundary extension
+define	SAMPLE		.2		# Subpixel sampling size
+define	SF_SZFNAME	79		# Length of file names
+define	SF_SZWTYPE	7		# Length of width type string
+
+# Main data structure
+define	SF		40
+define	SF_TASK		Memi[$1]	# Task type
+define	SF_WTYPE	Memc[P2C($1+1)]	# Width type string
+define	SF_WCODE	Memi[$1+5]	# Width code
+define	SF_BETA		Memr[P2R($1+6)]	# Moffat beta
+define	SF_SCALE	Memr[P2R($1+7)]	# Pixel scale
+define	SF_LEVEL	Memr[P2R($1+8)]	# Profile measurement level
+define	SF_RADIUS	Memr[P2R($1+9)]	# Profile radius
+define	SF_SBUF		Memr[P2R($1+10)]# Sky region buffer
+define	SF_SWIDTH	Memr[P2R($1+11)]# Sky region width
+define	SF_SAT		Memr[P2R($1+12)]# Saturation
+define	SF_NIT		Memi[$1+13]	# Number of iterations for radius
+define	SF_OVRPLT	Memi[$1+14]	# Overplot the best profile?
+define	SF_NCOLS	Memi[$1+15]	# Number of image columns
+define	SF_NLINES	Memi[$1+16]	# Number of image lines
+define	SF_XF		Memr[P2R($1+17)]# X field center
+define	SF_YF		Memr[P2R($1+18)]# Y field center
+define	SF_GP		Memi[$1+19]	# GIO pointer
+define	SF_F		Memr[P2R($1+20)]# Best focus
+define	SF_W		Memr[P2R($1+21)]# Width at best focus
+define	SF_M		Memr[P2R($1+22)]# Brightest star magnitude
+define	SF_XP1		Memr[P2R($1+23)]# First derivative point to plot
+define	SF_XP2		Memr[P2R($1+24)]# Last derivative point to plot
+define	SF_YP1		Memr[P2R($1+25)]# Minimum of derivative profile
+define	SF_YP2		Memr[P2R($1+26)]# Maximum of derivative profile
+define	SF_N		Memi[$1+27]	# Number of points not deleted
+define	SF_NSFD		Memi[$1+28]	# Number of data points
+define	SF_SFDS		Memi[$1+29]	# Pointer to data structures
+define	SF_NS		Memi[$1+30]	# Number of stars not deleted
+define	SF_NSTARS	Memi[$1+31]	# Number of stars
+define	SF_STARS	Memi[$1+32]	# Pointer to star groups
+define	SF_NF		Memi[$1+33]	# Number of focuses not deleted
+define	SF_NFOCUS	Memi[$1+34]	# Number of different focus values
+define	SF_FOCUS	Memi[$1+35]	# Pointer to focus groups
+define	SF_NI		Memi[$1+36]	# Number of images not deleted
+define	SF_NIMAGES	Memi[$1+37]	# Number of images
+define	SF_IMAGES	Memi[$1+38]	# Pointer to image groups
+define	SF_BEST		Memi[$1+39]	# Pointer to best focus star
+
+define	SF_SFD		Memi[SF_SFDS($1)+$2-1]
+define	SF_SFS		Memi[SF_STARS($1)+$2-1]
+define	SF_SFF		Memi[SF_FOCUS($1)+$2-1]
+define	SF_SFI		Memi[SF_IMAGES($1)+$2-1]
+
+# Basic data structure.
+define	SFD		94
+define	SFD_IMAGE	Memc[P2C($1)]	# Image name
+define	SFD_DATA	Memi[$1+40]	# Pointer to real image raster
+define	SFD_RADIUS	Memr[P2R($1+41)]# Profile radius
+define	SFD_NP		Memi[$1+42]	# Number of profile points
+define	SFD_NPMAX	Memi[$1+43]	# Maximum number of profile points
+define	SFD_X1		Memi[$1+44]	# Image raster limits
+define	SFD_X2		Memi[$1+45]
+define	SFD_Y1		Memi[$1+46]
+define	SFD_Y2		Memi[$1+47]
+define	SFD_ID		Memi[$1+48]	# Star ID
+define	SFD_X		Memr[P2R($1+49)]# Star X position
+define	SFD_Y		Memr[P2R($1+50)]# Star Y position
+define	SFD_F		Memr[P2R($1+51)]# Focus
+define	SFD_W		Memr[P2R($1+52)]# Width to use
+define	SFD_M		Memr[P2R($1+53)]# Magnitude
+define	SFD_E		Memr[P2R($1+54)]# Ellipticity
+define	SFD_PA		Memr[P2R($1+55)]# Position angle
+define	SFD_R		Memr[P2R($1+56)]# Radius at given level
+define	SFD_DFWHM	Memr[P2R($1+57)]# Direct FWHM
+define	SFD_GFWHM	Memr[P2R($1+58)]# Gaussian FWHM
+define	SFD_MFWHM	Memr[P2R($1+59)]# Moffat FWHM
+define	SFD_ASI1	Memi[$1+60]	# Pointer to enclosed flux profile
+define	SFD_ASI2	Memi[$1+61]	# Pointer to derivative profile
+define	SFD_YP1		Memr[P2R($1+62)]# Minimum of derivative profile
+define	SFD_YP2		Memr[P2R($1+63)]# Maximum of derivative profile
+define	SFD_FWHM	Memr[P2R($1+$2+63)]	# FWHM vs level=0.5*i (i=1-19)
+define	SFD_BKGD	Memr[P2R($1+83)]# Background value
+define	SFD_BKGD1	Memr[P2R($1+84)]# Original background value
+define	SFD_MISO	Memr[P2R($1+85)]# Moment isophote
+define	SFD_SIGMA	Memr[P2R($1+86)]# Moffat alpha
+define	SFD_ALPHA	Memr[P2R($1+87)]# Moffat alpha
+define	SFD_BETA	Memr[P2R($1+88)]# Moffat beta
+define	SFD_STATUS	Memi[$1+89]	# Status
+define	SFD_NSAT	Memi[$1+90]	# Number of saturated pixels
+define	SFD_SFS		Memi[$1+91]	# Pointer to star group
+define	SFD_SFF		Memi[$1+92]	# Pointer to focus group
+define	SFD_SFI		Memi[$1+93]	# Pointer to image group
+
+
+# Structure grouping data by star.
+define	SFS		($1+7)
+define	SFS_ID		Memi[$1]	# Star ID
+define	SFS_F		Memr[P2R($1+1)]	# Best focus
+define	SFS_W		Memr[P2R($1+2)]	# Best width
+define	SFS_M		Memr[P2R($1+3)]	# Average magnitude
+define	SFS_N		Memi[$1+4]	# Number of points used
+define	SFS_NF		Memi[$1+5]	# Number of focuses
+define	SFS_NSFD	Memi[$1+6]	# Number of data points
+define	SFS_SFD		Memi[$1+$2+6]	# Array of data structures
+
+
+# Structure grouping stars by focus values.
+define	SFF		($1+5)
+define	SFF_F		Memr[P2R($1)]	# Focus
+define	SFF_W		Memr[P2R($1+1)]	# Average width
+define	SFF_N		Memi[$1+2]	# Number in average
+define	SFF_NI		Memi[$1+3]	# Number of images
+define	SFF_NSFD	Memi[$1+4]	# Number of data points
+define	SFF_SFD		Memi[$1+$2+4]	# Array of data structures
+
+
+# Structure grouping stars by image.
+define	SFI		($1+42)
+define	SFI_IMAGE	Memc[P2C($1)]	# Image
+define	SFI_N		Memi[$1+40]	# Number in imagE
+define	SFI_NSFD	Memi[$1+41]	# Number of data points
+define	SFI_SFD		Memi[$1+$2+41]	# Array of data structures
diff --git a/noao/obsutil/src/starfocus/starfocus.key b/noao/obsutil/src/starfocus/starfocus.key
new file mode 100644
index 00000000..119e74a6
--- /dev/null
+++ b/noao/obsutil/src/starfocus/starfocus.key
@@ -0,0 +1,15 @@
+			COMMAND OPTIONS
+
+?	Page this help information.
+g	Measure object and graph the results.
+m	Measure object.
+q	Quit object marking and go to next image.
+	At the end of all images go to analysis of all measurements.
+
+:show	Show the current results.
+
+When an object is measured the center and enclosed flux profile is determined.
+When using the "mark1" option typing 'q' will measure all remaining images
+at the same cursor positions while the "markall" option goes to the
+next image until the image list is finished.  If measuring only one
+object with the 'g' key then a 'q' will exit the program.
diff --git a/noao/obsutil/src/starfocus/starfocus.par b/noao/obsutil/src/starfocus/starfocus.par
new file mode 100644
index 00000000..b1b1d5b0
--- /dev/null
+++ b/noao/obsutil/src/starfocus/starfocus.par
@@ -0,0 +1,32 @@
+# STARFOCUS -- Determine best star focus.
+
+images,s,a,,,,"List of images"
+focus,s,h,"1x1",,,"Focus values"
+fstep,s,h,"",,,"Focus step
+"
+nexposures,s,h,"1",,,"Number of exposures"
+step,s,h,"30",,,"Step in pixels"
+direction,s,h,"-line","-line|+line|-column|+column",,"Step direction"
+gap,s,h,"end","none|beginning|end",,"Double step gap
+"
+coords,s,h,"mark1","center|mark1|markall",,"Object coordinates"
+wcs,s,h,"logical","logical|physical|world",,"Coordinate system"
+display,b,h,yes,,,"Display images?"
+frame,i,h,1,1,,"Display frame to use
+"
+level,r,h,0.5,,,"Measurement level (fraction or percent)"
+size,s,h,"FWHM","Radius|FWHM|GFWHM|MFWHM",,"Size to display"
+beta,r,h,INDEF,2.1,,Moffat beta parameter
+scale,r,h,1.,,,"Pixel scale"
+radius,r,h,5.,,,"Measurement radius (pixels)"
+sbuffer,r,h,5,,,"Sky buffer (pixels)"
+swidth,r,h,5.,,,"Sky width (pixels)"
+saturation,r,h,INDEF,,,"Saturation level"
+ignore_sat,b,h,no,,,"Ignore objects with saturated pixels?"
+iterations,i,h,2,1,,"Number of radius adjustment iterations"
+xcenter,r,h,INDEF,,,X field center (pixels)
+ycenter,r,h,INDEF,,,Y field center (pixels)
+logfile,s,h,"logfile",,,"Logfile
+"
+imagecur,*imcur,h,"",,,"Image cursor input"
+graphcur,*gcur,h,"",,,"Graphics cursor input"
diff --git a/noao/obsutil/src/starfocus/stfgraph.x b/noao/obsutil/src/starfocus/stfgraph.x
new file mode 100644
index 00000000..890dfc81
--- /dev/null
+++ b/noao/obsutil/src/starfocus/stfgraph.x
@@ -0,0 +1,2682 @@
+include	<error.h>
+include	<gset.h>
+include	<mach.h>
+include	"starfocus.h"
+
+# Interactive help files.  There is one for STARFOCUS and one for PSFMEASUE.
+define	STFHELP		"starfocus$stfhelp.key"
+define	PSFHELP		"starfocus$psfhelp.key"
+define	PROMPT		"Options"
+
+# View ports for all plots.
+define	VX1	.15			# Minimum X viewport for left graph
+define	VX2	.47			# Maximum X viewport for left graph
+define	VX3	.63			# Minimum X viewport for right graph
+define	VX4	.95			# Maximum X viewport for right graph
+define	VY1	.10			# Minimum Y viewport for bottom graph
+define	VY2	.44			# Minimum Y viewport for bottom graph
+define	VY3	.54			# Minimum Y viewport for top graph
+define	VY4	.88			# Maximum Y viewport for top graph
+
+# Miscellaneous graphics parameters.
+define	NMAX	5			# Maximum number of samples for labeling
+define	HLCOLOR	2			# Highlight color
+define	HLWIDTH	4.			# Highlight width
+define	GM_MARK	GM_CROSS		# Point marker
+define	GM_MAG	GM_PLUS+GM_CROSS	# Magnitude marker
+
+
+# STF_GRAPH -- Interactive graphing of results.
+
+procedure stf_graph (sf)
+
+pointer	sf			#I Starfocus structure
+
+real	wx, wy, x, y, r2, r2min, fa[8]
+int	i, j, ix, iy, nx, ny, wcs, key, pkey, skey, redraw, clgcur()
+pointer	sp, sysidstr, title, cmd, gp, gopen()
+pointer	sfd, sfs, sff, current, nearest
+
+data	fa/0.,1.,1.,0.,0.,0.,1.,1./
+
+begin
+	call smark (sp)
+	call salloc (sysidstr, SZ_LINE, TY_CHAR)
+	call salloc (title, SZ_LINE, TY_CHAR)
+	call salloc (cmd, SZ_LINE, TY_CHAR)
+
+	# Set system id label
+	call sysid (Memc[sysidstr], SZ_LINE)
+
+	# Open graphics and enter interactive graphics loop
+	SF_GP(sf) = gopen ("stdgraph", NEW_FILE, STDGRAPH)
+	gp = SF_GP(sf)
+	wcs = 0
+	if (SF_NF(sf) > 1)
+	    key = 'f'
+	else if (SF_NS(sf) > 1)
+	    key = 'a'
+	else
+	    key = 'z'
+	pkey = 0
+	skey = 1
+	current = SF_BEST(sf)
+	repeat {
+	    switch (key) {
+	    case 'q':	# Quit
+		break
+	    case '?':	# Help
+		if (SF_TASK(sf) == PSFMEASURE)
+		    call gpagefile (gp, PSFHELP, PROMPT)
+		else
+		    call gpagefile (gp, STFHELP, PROMPT)
+		next
+	    case ':':	# Colon commands
+		iferr (call stf_colon (sf, Memc[cmd], redraw))
+		    redraw = NO
+		if (redraw == NO)
+		    next
+	    case 'a', 'b', 'e', 'f', 'g', 'm', 'p', 't', 'z':	# Plots
+		# When there is not enough data for the requested plot
+		# map the key to another one.  This is done mostly to
+		# avoid redrawing the same graph when different keys
+		# map to the same pkey.  The 'e', 'g', and 'p' key may
+		# select a different object so the check for the same
+		# plot is deferred.
+
+		if (SF_NS(sf) > 1 && SF_NF(sf) > 1) {
+		    ;
+		} else if (SF_NS(sf) > 1) {
+		    if (key == 'b')
+			key = 'a'
+		    if (key == 'f')
+			key = 'm'
+		} else if (SF_NF(sf) > 1) {
+		    if (key == 'a' || key == 'b' || key == 'm' || key == 't')
+			key = 'f'
+		} else {
+		    key = 'z'
+		}
+
+		switch (key) {
+		case 'e', 'g', 'p':
+		    ;
+		default:
+		    if (key == pkey)
+			next
+		}
+	    case 's':	# Toggle plotting of magnitude symbols
+		if (pkey != 'a' && pkey != 'b')
+		    next
+		skey = mod (skey+1, 2)
+	    case 'u':	# Undelete all
+		j = 0
+		do i = 1, SF_NSFD(sf) {
+		    sfd = SF_SFD(sf,i)
+		    if (SFD_STATUS(sfd) != 0) {
+			SFD_STATUS(sfd) = 0
+			j = j + 1
+		    }
+		}
+		if (j == 0)
+		    next
+		call stf_fitfocus (sf)
+	    case 'd', 'n', 'o', 'r', 'x', 'i', ' ':	# Misc
+	        ;
+	    default:	# Unknown
+		call printf ("\007")
+		next
+	    }
+
+	    # Find the nearest or next object if needed.
+	    switch (key) {
+	    case 'r', 's', 'u', ':':	# Redraw last graph
+		pkey = pkey
+		nearest = current
+	    case 'n', 'o':		# Renormalize enclosed flux profile
+		if (wcs != 7 || pkey == 'p')
+		    next
+		pkey = pkey
+		nearest = current
+		if (key == 'n')
+		    call stf_norm (sf, nearest, wx, INDEF)
+		else
+		    call stf_norm (sf, nearest, wx, wy)
+		call stf_widths (sf, nearest)
+		call stf_fwhms (sf, nearest)
+		call stf_fitfocus (sf)
+	    case ' ':			# Select next focus or star
+		switch (pkey) {
+		case 'a', 'm', 't':
+		    sff = SFD_SFF(current)
+		    for (i=1; SF_SFF(sf,i)!=sff; i=i+1)
+			;
+		    j = SF_NFOCUS(sf)
+		    i = mod (i, j) + 1
+		    for (; SFF_N(SF_SFF(sf,i))==0; i=mod(i,j)+1)
+			;
+		    if (SF_SFF(sf,i) == sff)
+			next
+		    sff = SF_SFF(sf,i)
+		    do i = 1, SFF_NSFD(sff) {
+			nearest = SFF_SFD(sff,i)
+			if (SFD_STATUS(nearest) == 0)
+			    break
+		    }
+		case 'e', 'g', 'p', 'z':
+		    switch (wcs) {
+		    case 7, 8, 11:
+			for (i=1; SF_SFD(sf,i)!=current; i=i+1)
+			    ;
+			j = SF_NSFD(sf)
+			i = mod (i, j) + 1
+			for (; SFD_STATUS(SF_SFD(sf,i))!=0; i=mod(i,j)+1)
+			    ;
+			nearest = SF_SFD(sf,i)
+		    case 9:
+			sfs = SFD_SFS(current)
+			for (i=1; SFS_SFD(sfs,i)!=current; i=i+1)
+			    ;
+			j = SFS_NSFD(sfs)
+			i = mod (i, j) + 1
+			for (; SFD_STATUS(SFS_SFD(sfs,i))!=0; i=mod(i,j)+1)
+			    ;
+			nearest = SFS_SFD(sfs,i)
+			if (nearest == current)
+			    next
+		    case 10:
+			sff = SFD_SFF(current)
+			for (i=1; SFF_SFD(sff,i)!=current; i=i+1)
+			    ;
+			j = SFF_NSFD(sff)
+			i = mod (i, j) + 1
+			for (; SFD_STATUS(SFF_SFD(sff,i))!=0; i=mod(i,j)+1)
+			    ;
+			nearest = SFF_SFD(sff,i)
+			if (nearest == current)
+			    next
+		    }
+		default:
+		    next
+		}
+	    default:			# Select nearest to cursor
+		switch (pkey) {
+		case 'a':
+		    r2min = MAX_REAL
+		    call gctran (gp, wx, wy, wx, wy, wcs, 0)
+		    sff = SFD_SFF(current)
+		    do i = 1, SFF_NSFD(sff) {
+			sfd = SFF_SFD(sff,i)
+			if (SFD_STATUS(sfd) != 0)
+			    next
+			switch (wcs) {
+			case 1:
+			    x = SFD_X(sfd)
+			    y = SFD_Y(sfd)
+			case 2:
+			    x = SFD_X(sfd)
+			    y = SFD_W(sfd)
+			case 3:
+			    x = SFD_W(sfd)
+			    y = SFD_Y(sfd)
+			case 4:
+			    x = SFD_X(sfd)
+			    y = SFD_E(sfd)
+			case 5:
+			    x = SFD_E(sfd)
+			    y = SFD_Y(sfd)
+			}
+			call gctran (gp, x, y, x, y, wcs, 0)
+			r2 = (x-wx)**2 + (y-wy)**2
+			if (r2 < r2min) {
+			    r2min = r2
+			    nearest = sfd
+			}
+		    }
+		case 'b': 
+		    r2min = MAX_REAL
+		    call gctran (gp, wx, wy, wx, wy, wcs, 0)
+		    do i = 1, SF_NSTARS(sf) {
+			sfs = SF_SFS(sf,i)
+			if (SFS_N(sfs) == 0)
+			    next
+			switch (wcs) {
+			case 1:
+			    x = SFD_X(SFS_SFD(sfs,1))
+			    y = SFD_Y(SFS_SFD(sfs,1))
+			case 2:
+			    x = SFD_X(SFS_SFD(sfs,1))
+			    y = SFS_W(sfs)
+			case 3:
+			    x = SFS_W(sfs)
+			    y = SFD_Y(SFS_SFD(sfs,1))
+			case 4:
+			    x = SFD_X(SFS_SFD(sfs,1))
+			    y = SFS_F(sfs)
+			case 5:
+			    x = SFS_F(sfs)
+			    y = SFD_Y(SFS_SFD(sfs,1))
+			}
+			call gctran (gp, x, y, x, y, wcs, 0)
+			r2 = (x-wx)**2 + (y-wy)**2
+			if (r2 < r2min) {
+			    r2min = r2
+			    nearest = sfs
+			}
+		    }
+		    sfs = nearest
+		    r2min = MAX_REAL
+		    do i = 1, SFS_NSFD(sfs) {
+			sfd = SFS_SFD(sfs,i)
+			if (SFD_STATUS(sfd) != 0)
+			    next
+			r2 = SFD_W(sfd)
+			if (r2 < r2min) {
+			    r2min = r2
+			    nearest = sfd
+			}
+		    }
+		case 'e', 'g', 'p':
+		    switch (wcs) {
+		    case 9:
+			sfs = SFD_SFS(current)
+			i = SFS_N(sfs)
+			if (i < 4) {
+			    nx = i
+			    ny = 1
+			} else {
+			    nx = nint (sqrt (real (i)))
+			    if (mod (i-1, nx+1) >= mod (i-1, nx))
+				nx = nx + 1
+			    ny = (i - 1) / nx + 1
+			}
+			ix = max (1, min (nx, nint(wx)))
+			iy = max (1, min (ny, nint(wy)))
+
+			j = 0
+			do i = 1, SFS_NSFD(sfs) {
+			    sfd = SFS_SFD(sfs, i)
+			    if (SFD_STATUS(sfd) != 0)
+				next
+			    if (ix == 1 + mod (j, nx) && iy == 1 + j / nx) {
+				nearest = sfd
+				break
+			    }
+			    j = j + 1
+			}
+		    case 10:
+			sff = SFD_SFF(current)
+			i = SFF_N(sff)
+			if (i < 4) {
+			    nx = i
+			    ny = 1
+			} else {
+			    nx = nint (sqrt (real (i)))
+			    if (mod (i-1, nx+1) >= mod (i-1, nx))
+				nx = nx + 1
+			    ny = (i - 1) / nx + 1
+			}
+			ix = max (1, min (nx, nint(wx)))
+			iy = max (1, min (ny, nint(wy)))
+
+			j = 0
+			do i = 1, SFF_NSFD(sff) {
+			    sfd = SFF_SFD(sff, i)
+			    if (SFD_STATUS(sfd) != 0)
+				next
+			    if (ix == 1 + mod (j, nx) && iy == 1 + j / nx) {
+				nearest = sfd
+				break
+			    }
+			    j = j + 1
+			}
+		    }
+		    if (key == pkey && nearest == current)
+			next
+		default:
+		    switch (wcs) {
+		    case 1, 2:
+			r2min = MAX_REAL
+			call gctran (gp, wx, wy, wx, wy, wcs, 0)
+			do i = 1, SF_NSFD(sf) {
+			    sfd = SF_SFD(sf,i)
+			    if (SFD_STATUS(sfd) != 0)
+				next
+			    switch (wcs) {
+			    case 1:
+				x = SFD_F(sfd)
+				y = SFD_W(sfd)
+			    case 2:
+				x = SFD_F(sfd)
+				y = SFD_E(sfd)
+			    }
+			    call gctran (gp, x, y, x, y, wcs, 0)
+			    r2 = (x-wx)**2 + (y-wy)**2
+			    if (r2 < r2min) {
+				r2min = r2
+				nearest = sfd
+			    }
+			}
+		    case 3, 4, 5, 6:
+			r2min = MAX_REAL
+			call gctran (gp, wx, wy, wx, wy, wcs, 0)
+			sff = SFD_SFF(current)
+			do i = 1, SFF_NSFD(sff) {
+			    sfd = SFF_SFD(sff,i)
+			    if (SFD_STATUS(sfd) != 0)
+				next
+			    switch (wcs) {
+			    case 3:
+				x = -2.5 * log10 (SFS_M(SFD_SFS(sfd))/SF_M(sf))
+				y = SFD_W(sfd)
+			    case 4:
+				x = -2.5 * log10 (SFS_M(SFD_SFS(sfd))/SF_M(sf))
+				y = SFD_E(sfd)
+			    case 5:
+				x = sqrt ((SFD_X(sfd) - SF_XF(sf)) ** 2 +
+				    (SFD_Y(sfd) - SF_YF(sf)) ** 2)
+				y = SFD_W(sfd)
+			    case 6:
+				x = sqrt ((SFD_X(sfd) - SF_XF(sf)) ** 2 +
+				    (SFD_Y(sfd) - SF_YF(sf)) ** 2)
+				y = SFD_E(sfd)
+			    }
+			    call gctran (gp, x, y, x, y, wcs, 0)
+			    r2 = (x-wx)**2 + (y-wy)**2
+			    if (r2 < r2min) {
+				r2min = r2
+				nearest = sfd
+			    }
+			}
+		    default:
+			nearest = current
+		    }
+		}
+
+		# Act on selection for delete or info.
+		switch (key) {
+		case 'd':
+		    if (SF_NS(sf) > 1) {
+			sfs = SFD_SFS(nearest)
+			do i = 1, SFS_NSFD(sfs)
+			    SFD_STATUS(SFS_SFD(sfs,i)) = 1
+		    } else
+			SFD_STATUS(nearest) = 1
+		    call stf_fitfocus (sf)
+		case 'x':
+		    repeat {
+			switch (key) {
+			case 'f':
+			    sff = SFD_SFF(nearest)
+			    do i = 1, SFF_NSFD(sff)
+				SFD_STATUS(SFF_SFD(sff,i)) = 1
+			case 'i':
+			    sfd = SFD_SFI(nearest)
+			    do i = 1, SFI_NSFD(sfd)
+				SFD_STATUS(SFI_SFD(sfd,i)) = 1
+			case 'p':
+			    SFD_STATUS(nearest) = 1
+			case 's':
+			    sfs = SFD_SFS(nearest)
+			    do i = 1, SFS_NSFD(sfs)
+				SFD_STATUS(SFS_SFD(sfs,i)) = 1
+			default:
+			    call printf (
+			    "Delete image, star, focus, or point? (i|s|f|p)")
+			    next
+			}
+			call stf_fitfocus (sf)
+			break
+		    } until (clgcur ("graphcur",
+			wx, wy, wcs, key, Memc[cmd], SZ_LINE) == EOF)
+		case 'i':
+		    switch (pkey) {
+		    case 'b':
+			sfs = SFD_SFS(nearest)
+			call stf_title (sf, NULL, sfs, NULL,
+			    Memc[title], SZ_LINE)
+		    default:
+			call stf_title (sf, nearest, NULL, NULL,
+			    Memc[title], SZ_LINE)
+		    }
+		    call printf ("%s\n")
+			call pargstr (Memc[title])
+		    next
+		default:
+		    pkey = key
+		}
+	    }
+
+	    # If current object has been deleted select another.
+	    if (SFD_STATUS(nearest) == 0)
+		current = nearest
+	    else
+		current = SF_BEST(sf)
+
+	    # Make the graphs.  The graph depends on the number of stars
+	    # and number of focus values.  Note that the pkey has already
+	    # been mapped but all the keys are shown for clarity.
+
+	    call gclear (gp)
+	    call gseti (gp, G_FACOLOR, 0)
+
+	    if (SF_NS(sf) > 1 && SF_NF(sf) > 1) {
+		switch (pkey) {
+		case 'a':
+		    sff = SFD_SFF(current)
+		    call stf_title (sf, NULL, NULL, sff, Memc[title], SZ_LINE)
+		    call gseti (gp, G_WCS, 1)
+		    call gsview (gp, VX1, VX4, VY1, VY4)
+		    call stf_g11 (sf, current, skey, Memc[title])
+		case 'b':
+		    call sprintf (Memc[title], SZ_LINE,
+			"Best focus estimates for each star")
+		    call gseti (gp, G_WCS, 1)
+		    call gsview (gp, VX1, VX4, VY1, VY4)
+		    call stf_g12 (sf, current, skey, Memc[title])
+		case 'e':
+		    sfs = SFD_SFS(current)
+		    call sprintf (Memc[title], SZ_LINE,
+			"Star: x=%.2f, y=%.2f, m=%.2f")
+			call pargr (SFD_X(current))
+			call pargr (SFD_Y(current))
+			call pargr (-2.5 * log10 (SFS_M(sfs) / SF_M(sf)))
+		    call gseti (gp, G_WCS, 9)
+		    call gsview (gp, VX1, VX4, VY3, VY4)
+		    call stf_g2 (sf, current, Memc[title])
+		    sff = SFD_SFF(current)
+		    call stf_title (sf, NULL, NULL, sff, Memc[title], SZ_LINE)
+		    call gseti (gp, G_WCS, 10)
+		    call gsview (gp, VX1, VX4, VY1, VY2)
+		    call stf_g3 (sf, current, Memc[title])
+		case 'f':
+		    call gseti (gp, G_WCS, 1)
+		    call gsview (gp, VX1, VX4, VY3, VY4)
+		    call stf_g1 (sf, current, 'f', 'r', "", "", SF_WTYPE(sf))
+		    call gseti (gp, G_WCS, 2)
+		    call gsview (gp, VX1, VX4, VY1, VY2)
+		    call gfill (gp, fa, fa[5], 4, GF_SOLID)
+		    call stf_g1 (sf, current, 'f', 'e', "", "Focus",
+			"Ellipticity")
+		case 'g':
+		    sfs = SFD_SFS(current)
+		    call sprintf (Memc[title], SZ_LINE,
+			"Star: x=%.2f, y=%.2f, m=%.2f")
+			call pargr (SFD_X(current))
+			call pargr (SFD_Y(current))
+			call pargr (-2.5 * log10 (SFS_M(sfs) / SF_M(sf)))
+		    call gseti (gp, G_WCS, 9)
+		    call gsview (gp, VX1, VX4, VY3, VY4)
+		    call stf_g9 (sf, current, Memc[title])
+		    sff = SFD_SFF(current)
+		    call stf_title (sf, NULL, NULL, sff, Memc[title], SZ_LINE)
+		    call gseti (gp, G_WCS, 10)
+		    call gsview (gp, VX1, VX4, VY1, VY2)
+		    call stf_g10 (sf, current, Memc[title])
+		case 'm':
+		    sff = SFD_SFF(current)
+		    call stf_title (sf, NULL, NULL, sff, Memc[title], SZ_LINE)
+		    call gseti (gp, G_WCS, 3)
+		    call gsview (gp, VX1, VX4, VY3, VY4)
+		    call stf_g1 (sf, current, 'm', 'r', Memc[title],
+			"", SF_WTYPE(sf))
+		    call gseti (gp, G_WCS, 4)
+		    call gsview (gp, VX1, VX4, VY1, VY2)
+		    call gfill (gp, fa, fa[5], 4, GF_SOLID)
+		    call stf_g1 (sf, current, 'm', 'e', "", "Magnitude",
+			"Ellipticity")
+		case 'p':
+		    sfs = SFD_SFS(current)
+		    call sprintf (Memc[title], SZ_LINE,
+			"Star: x=%.2f, y=%.2f, m=%.2f")
+			call pargr (SFD_X(current))
+			call pargr (SFD_Y(current))
+			call pargr (-2.5 * log10 (SFS_M(sfs) / SF_M(sf)))
+		    call gseti (gp, G_WCS, 9)
+		    call gsview (gp, VX1, VX4, VY3, VY4)
+		    call stf_g4 (sf, current, Memc[title])
+		    sff = SFD_SFF(current)
+		    call stf_title (sf, NULL, NULL, sff, Memc[title], SZ_LINE)
+		    call gseti (gp, G_WCS, 10)
+		    call gsview (gp, VX1, VX4, VY1, VY2)
+		    call stf_g5 (sf, current, Memc[title])
+		case 't':
+		    sff = SFD_SFF(current)
+		    call stf_title (sf, NULL, NULL, sff, Memc[title], SZ_LINE)
+		    call gseti (gp, G_WCS, 5)
+		    call gsview (gp, VX1, VX4, VY3, VY4)
+		    call stf_g1 (sf, current, 't', 'r', Memc[title],
+			"", SF_WTYPE(sf))
+		    call gseti (gp, G_WCS, 6)
+		    call gsview (gp, VX1, VX4, VY1, VY2)
+		    call gfill (gp, fa, fa[5], 4, GF_SOLID)
+		    call stf_g1 (sf, current, 't', 'e', "", "Field radius",
+			"Ellipticity")
+		case 'z':
+		    call gseti (gp, G_WCS, 7)
+		    call gsview (gp, VX1, VX2, VY3, VY4)
+		    call stf_g6 (sf, current, "", "", "Enclosed flux")
+		    call gseti (gp, G_WCS, 8)
+		    call gsview (gp, VX1, VX2, VY1, VY2)
+		    call gfill (gp, fa, fa[5], 4, GF_SOLID)
+		    call stf_g7 (sf, current, "", "Radius", "Profile")
+		    call gseti (gp, G_WCS, 11)
+		    call gsview (gp, VX3, VX4, VY3, VY4)
+		    call gfill (gp, fa, fa[5], 4, GF_SOLID)
+		    call stf_g8 (sf, current, "", "Enclosed flux", "FWHM")
+
+		    call stf_title (sf, current, NULL, NULL, Memc[title],
+			SZ_LINE)
+		    call gseti (gp, G_WCS, 0)
+		    call gsetr (gp, G_PLWIDTH, 2.0)
+		    call gline (gp, 0., 0., 0., 0.)
+		    call gtext (gp, 0.5, 0.93, Memc[title], "h=c,v=t")
+		}
+	    } else if (SF_NS(sf) > 1) {
+		switch (pkey) {
+		case 'a', 'b':
+		    sff = SFD_SFF(current)
+		    call stf_title (sf, NULL, NULL, sff, Memc[title], SZ_LINE)
+		    call gseti (gp, G_WCS, 1)
+		    call gsview (gp, VX1, VX4, VY1, VY4)
+		    call stf_g11 (sf, current, skey, Memc[title])
+		case 'e':
+		    sff = SFD_SFF(current)
+		    call stf_title (sf, NULL, NULL, sff, Memc[title], SZ_LINE)
+		    call gseti (gp, G_WCS, 10)
+		    call gsview (gp, VX1, VX4, VY3, VY4)
+		    call stf_g3 (sf, current, Memc[title])
+		    call stf_title (sf, current, NULL, NULL, Memc[title],
+			SZ_LINE)
+		    call gseti (gp, G_WCS, 7)
+		    call gsview (gp, VX1, VX4, VY1, VY2)
+		    call gfill (gp, fa, fa[5], 4, GF_SOLID)
+		    call stf_g6 (sf, current, Memc[title], "Radius",
+			"Enclosed flux")
+		case 'f', 'm':
+		    sff = SFD_SFF(current)
+		    call stf_title (sf, NULL, NULL, sff, Memc[title], SZ_LINE)
+		    call gseti (gp, G_WCS, 3)
+		    call gsview (gp, VX1, VX4, VY3, VY4)
+		    call stf_g1 (sf, current, 'm', 'r', Memc[title], "",
+			SF_WTYPE(sf))
+		    call gseti (gp, G_WCS, 4)
+		    call gsview (gp, VX1, VX4, VY1, VY2)
+		    call gfill (gp, fa, fa[5], 4, GF_SOLID)
+		    call stf_g1 (sf, current, 'm', 'e', "", "Magnitude",
+			"Ellipticity")
+		case 'g':
+		    sff = SFD_SFF(current)
+		    call stf_title (sf, NULL, NULL, sff, Memc[title], SZ_LINE)
+		    call gseti (gp, G_WCS, 10)
+		    call gsview (gp, VX1, VX4, VY3, VY4)
+		    call stf_g10 (sf, current, Memc[title])
+		    call stf_title (sf, current, NULL, NULL, Memc[title],
+			SZ_LINE)
+		    call gseti (gp, G_WCS, 11)
+		    call gsview (gp, VX1, VX4, VY1, VY2)
+		    call gfill (gp, fa, fa[5], 4, GF_SOLID)
+		    call stf_g8 (sf, current, Memc[title], "Enclosed flux",
+			"FWHM")
+		case 'p':
+		    sff = SFD_SFF(current)
+		    call stf_title (sf, NULL, NULL, sff, Memc[title], SZ_LINE)
+		    call gseti (gp, G_WCS, 10)
+		    call gsview (gp, VX1, VX4, VY3, VY4)
+		    call stf_g5 (sf, current, Memc[title])
+		    call stf_title (sf, current, NULL, NULL, Memc[title],
+			SZ_LINE)
+		    call gseti (gp, G_WCS, 7)
+		    call gsview (gp, VX1, VX4, VY1, VY2)
+		    call gfill (gp, fa, fa[5], 4, GF_SOLID)
+		    call stf_g7 (sf, current, Memc[title], "Radius", "Profile")
+		case 't':
+		    sff = SFD_SFF(current)
+		    call stf_title (sf, NULL, NULL, sff, Memc[title], SZ_LINE)
+		    call gseti (gp, G_WCS, 5)
+		    call gsview (gp, VX1, VX4, VY3, VY4)
+		    call stf_g1 (sf, current, 't', 'r', Memc[title], "",
+			SF_WTYPE(sf))
+		    call gseti (gp, G_WCS, 6)
+		    call gsview (gp, VX1, VX4, VY1, VY2)
+		    call gfill (gp, fa, fa[5], 4, GF_SOLID)
+		    call stf_g1 (sf, current, 't', 'e', "", "Field radius",
+			"Ellipticity")
+		case 'z':
+		    call gseti (gp, G_WCS, 7)
+		    call gsview (gp, VX1, VX2, VY3, VY4)
+		    call stf_g6 (sf, current, "", "", "Enclosed flux")
+		    call gseti (gp, G_WCS, 8)
+		    call gsview (gp, VX1, VX2, VY1, VY2)
+		    call gfill (gp, fa, fa[5], 4, GF_SOLID)
+		    call stf_g7 (sf, current, "", "Radius", "Profile")
+		    call gseti (gp, G_WCS, 11)
+		    call gsview (gp, VX3, VX4, VY3, VY4)
+		    call gfill (gp, fa, fa[5], 4, GF_SOLID)
+		    call stf_g8 (sf, current, "", "Enclosed flux", "FWHM")
+
+		    call stf_title (sf, current, NULL, NULL, Memc[title],
+			SZ_LINE)
+		    call gseti (gp, G_WCS, 0)
+		    call gsetr (gp, G_PLWIDTH, 2.0)
+		    call gline (gp, 0., 0., 0., 0.)
+		    call gtext (gp, 0.5, 0.93, Memc[title], "h=c,v=t")
+		}
+	    } else if (SF_NF(sf) > 1) {
+		switch (pkey) {
+		case 'a', 'b', 'f', 'm', 't':
+		    call gseti (gp, G_WCS, 1)
+		    call gsview (gp, VX1, VX4, VY3, VY4)
+		    call stf_g1 (sf, current, 'f', 'r', "", "", SF_WTYPE(sf))
+		    call gseti (gp, G_WCS, 2)
+		    call gsview (gp, VX1, VX4, VY1, VY2)
+		    call gfill (gp, fa, fa[5], 4, GF_SOLID)
+		    call stf_g1 (sf, current, 'f', 'e', "", "Focus",
+			"Ellipticity")
+		case 'e':
+		    sfs = SFD_SFS(current)
+		    call sprintf (Memc[title], SZ_LINE,
+			"Star: x=%.2f, y=%.2f, m=%.2f")
+			call pargr (SFD_X(current))
+			call pargr (SFD_Y(current))
+			call pargr (-2.5 * log10 (SFS_M(sfs) / SF_M(sf)))
+		    call gseti (gp, G_WCS, 9)
+		    call gsview (gp, VX1, VX4, VY3, VY4)
+		    call stf_g2 (sf, current, Memc[title])
+		    call stf_title (sf, current, NULL, NULL, Memc[title],
+			SZ_LINE)
+		    call gseti (gp, G_WCS, 7)
+		    call gsview (gp, VX1, VX4, VY1, VY2)
+		    call gfill (gp, fa, fa[5], 4, GF_SOLID)
+		    call stf_g6 (sf, current, Memc[title], "Radius",
+			"Enclosed flux")
+		case 'g':
+		    sfs = SFD_SFS(current)
+		    call sprintf (Memc[title], SZ_LINE,
+			"Star: x=%.2f, y=%.2f, m=%.2f")
+			call pargr (SFD_X(current))
+			call pargr (SFD_Y(current))
+			call pargr (-2.5 * log10 (SFS_M(sfs) / SF_M(sf)))
+		    call gseti (gp, G_WCS, 9)
+		    call gsview (gp, VX1, VX4, VY3, VY4)
+		    call stf_g9 (sf, current, Memc[title])
+		    call stf_title (sf, current, NULL, NULL, Memc[title],
+			SZ_LINE)
+		    call gseti (gp, G_WCS, 11)
+		    call gsview (gp, VX1, VX4, VY1, VY2)
+		    call gfill (gp, fa, fa[5], 4, GF_SOLID)
+		    call stf_g8 (sf, current, Memc[title], "Enclosed flux",
+			"FWHM")
+		case 'p':
+		    sfs = SFD_SFS(current)
+		    call sprintf (Memc[title], SZ_LINE,
+			"Star: x=%.2f, y=%.2f, m=%.2f")
+			call pargr (SFD_X(current))
+			call pargr (SFD_Y(current))
+			call pargr (-2.5 * log10 (SFS_M(sfs) / SF_M(sf)))
+		    call gseti (gp, G_WCS, 9)
+		    call gsview (gp, VX1, VX4, VY3, VY4)
+		    call stf_g4 (sf, current, Memc[title])
+		    call stf_title (sf, current, NULL, NULL, Memc[title],
+			SZ_LINE)
+		    call gseti (gp, G_WCS, 7)
+		    call gsview (gp, VX1, VX4, VY1, VY2)
+		    call gfill (gp, fa, fa[5], 4, GF_SOLID)
+		    call stf_g7 (sf, current, Memc[title], "Radius",
+			"profile")
+		case 'z':
+		    call gseti (gp, G_WCS, 7)
+		    call gsview (gp, VX1, VX2, VY3, VY4)
+		    call stf_g6 (sf, current, "", "", "Enclosed flux")
+		    call gseti (gp, G_WCS, 8)
+		    call gsview (gp, VX1, VX2, VY1, VY2)
+		    call gfill (gp, fa, fa[5], 4, GF_SOLID)
+		    call stf_g7 (sf, current, "", "Radius", "Profile")
+		    call gseti (gp, G_WCS, 11)
+		    call gsview (gp, VX3, VX4, VY3, VY4)
+		    call gfill (gp, fa, fa[5], 4, GF_SOLID)
+		    call stf_g8 (sf, current, "", "Enclosed flux", "FWHM")
+
+		    call stf_title (sf, current, NULL, NULL, Memc[title],
+			SZ_LINE)
+		    call gseti (gp, G_WCS, 0)
+		    call gsetr (gp, G_PLWIDTH, 2.0)
+		    call gline (gp, 0., 0., 0., 0.)
+		    call gtext (gp, 0.5, 0.93, Memc[title], "h=c,v=t")
+		}
+	    } else {
+		switch (pkey) {
+		case 'a', 'b', 'f', 'm', 'p', 'z', 'e', 't':
+		    call gseti (gp, G_WCS, 7)
+		    call gsview (gp, VX1, VX2, VY3, VY4)
+		    call stf_g6 (sf, current, "", "", "Enclosed flux")
+		    call gseti (gp, G_WCS, 8)
+		    call gsview (gp, VX1, VX2, VY1, VY2)
+		    call gfill (gp, fa, fa[5], 4, GF_SOLID)
+		    call stf_g7 (sf, current, "", "Radius", "Profile")
+		    call gseti (gp, G_WCS, 11)
+		    call gsview (gp, VX3, VX4, VY3, VY4)
+		    call gfill (gp, fa, fa[5], 4, GF_SOLID)
+		    call stf_g8 (sf, current, "", "Enclosed flux", "FWHM")
+
+		    call stf_title (sf, current, NULL, NULL, Memc[title],
+			SZ_LINE)
+		    call gseti (gp, G_WCS, 0)
+		    call gsetr (gp, G_PLWIDTH, 2.0)
+		    call gline (gp, 0., 0., 0., 0.)
+		    call gtext (gp, 0.5, 0.93, Memc[title], "h=c,v=t")
+		}
+	    }
+
+	    # Add banner title.
+	    call stf_title (sf, NULL, NULL, NULL, Memc[title], SZ_LINE)
+	    call gseti (gp, G_WCS, 0)
+	    call gsetr (gp, G_PLWIDTH, 2.0)
+	    call gline (gp, 0., 0., 0., 0.)
+	    call gtext (gp, 0.5, 0.99, Memc[sysidstr], "h=c,v=t")
+	    call gtext (gp, 0.5, 0.96, Memc[title], "h=c,v=t")
+
+	    if (SF_NSFD(sf) == 1)
+		break
+
+	} until (clgcur ("graphcur", wx, wy, wcs, key, Memc[cmd], SZ_LINE)==EOF)
+
+	call gclose (gp)
+	call sfree (sp)
+end
+
+
+# List of colon commands.
+define	CMDS		"|show|level|size|scale|radius|xcenter|ycenter\
+			|overplot|beta|"
+define	SHOW		1	# Show current results
+define	LEVEL		2	# Measurement level
+define	SIZE		3	# Size type
+define	SCALE		4	# Pixel scale
+define	RADIUS		5	# Maximum radius
+define	XCENTER		6	# X field center
+define	YCENTER		7	# Y field center
+define	OVERPLOT	8	# Overplot best profile
+define	BETA		9	# Beta value for Moffat function
+
+# STF_COLON -- Respond to colon command.
+
+procedure stf_colon (sf, cmd, redraw)
+
+pointer	sf			#I Starfocus pointer
+char	cmd[ARB]		#I Colon command
+int	redraw			#O Redraw?
+
+bool	bval
+real	rval, stf_r2i()
+int	i, j, ncmd, nscan(), strdic(), open(), btoi()
+pointer	sp, str, sfd
+errchk	open, delete, stf_log, stf_norm, stf_radius, stf_fitfocus
+
+begin
+	call smark (sp)
+	call salloc (str, SZ_FNAME, TY_CHAR)
+
+	# Scan the command string and get the first word.
+	call sscan (cmd)
+	call gargwrd (Memc[str], SZ_FNAME)
+	ncmd = strdic (Memc[str], Memc[str], SZ_FNAME, CMDS)
+
+	switch (ncmd) {
+	case SHOW:
+	    call gargwrd (Memc[str], SZ_FNAME)
+	    iferr {
+		if (nscan() == 1) {
+		    call mktemp ("tmp$iraf", Memc[str], SZ_FNAME)
+		    i = open (Memc[str], APPEND, TEXT_FILE)
+		    call stf_log (sf, i)
+		    call close (i)
+		    call gpagefile (SF_GP(sf), Memc[str], "starfocus")
+		    call delete (Memc[str])
+		} else {
+		    i = open (Memc[str], APPEND, TEXT_FILE)
+		    call stf_log (sf, i)
+		    call close (i)
+		}
+	    } then
+		call erract (EA_WARN)
+	    redraw = NO
+	case LEVEL:
+	    call gargr (rval)
+	    if (nscan() == 2) {
+		if (rval > 1.)
+		    rval = rval / 100.
+		SF_LEVEL(sf) = max (0.05, min (0.95, rval))
+		do i = 1, SF_NSFD(sf) {
+		    sfd = SF_SFD(sf,i)
+		    call stf_radius (sf, sfd, SF_LEVEL(sf), SFD_R(sfd))
+		}
+		if (SF_WCODE(sf) == 1)
+		    call stf_fitfocus (sf) 
+		redraw = YES
+	    } else {
+		call printf ("level %g\n")
+		    call pargr (SF_LEVEL(sf))
+		redraw = NO
+	    }
+	case SIZE:
+	    call gargwrd (Memc[str], SZ_FNAME)
+	    if (nscan() == 2) {
+		ncmd = strdic (Memc[str], Memc[str], SZ_FNAME, SF_WTYPES)
+		if (ncmd == 0) {
+		    call eprintf ("Invalid size type\n")
+		    redraw = NO
+		} else {
+		    call strcpy (Memc[str], SF_WTYPE(sf), SF_SZWTYPE)
+		    SF_WCODE(sf) = ncmd
+		    do i = 1, SF_NSFD(sf) {
+			sfd = SF_SFD(sf,i)
+			switch (SF_WCODE(sf)) {
+			case 1:
+			    SFD_W(sfd) = SFD_R(sfd)
+			case 2:
+			    SFD_W(sfd) = SFD_DFWHM(sfd)
+			case 3:
+			    SFD_W(sfd) = SFD_GFWHM(sfd)
+			case 4:
+			    SFD_W(sfd) = SFD_MFWHM(sfd)
+			}
+			call stf_fwhms (sf, sfd)
+		    }
+		    call stf_fitfocus (sf) 
+		    redraw = YES
+		}
+	    } else {
+		call printf ("size %s\n")
+		    call pargstr (SF_WTYPE(sf))
+		redraw = NO
+	    }
+	case SCALE:
+	    call gargr (rval)
+	    if (nscan() == 2) {
+		rval = rval / SF_SCALE(sf)
+		SF_SCALE(sf) = SF_SCALE(sf) * rval
+		do i = 1, SF_NSFD(sf) {
+		    sfd = SF_SFD(sf,i)
+		    switch (SF_WCODE(sf)) {
+		    case 1:
+			SFD_R(sfd) = SFD_R(sfd) * rval
+			SFD_W(sfd) = SFD_R(sfd)
+		    case 2:
+			SFD_DFWHM(sfd) = SFD_DFWHM(sfd) * rval
+			SFD_W(sfd) = SFD_DFWHM(sfd)
+		    case 3:
+			SFD_SIGMA(sfd) = SFD_SIGMA(sfd) * rval
+			SFD_GFWHM(sfd) = SFD_GFWHM(sfd) * rval
+			SFD_W(sfd) = SFD_GFWHM(sfd)
+		    case 4:
+			SFD_ALPHA(sfd) = SFD_ALPHA(sfd) * rval
+			SFD_MFWHM(sfd) = SFD_MFWHM(sfd) * rval
+			SFD_W(sfd) = SFD_MFWHM(sfd)
+		    }
+		    do j = 1, 19
+			SFD_FWHM(sfd,j) = SFD_FWHM(sfd,j) * rval
+		}
+		do i = 1, SF_NSTARS(sf) {
+		    sfd = SF_SFS(sf,i)
+		    SFS_W(sfd) = SFS_W(sfd) * rval
+		}
+		do i = 1, SF_NFOCUS(sf) {
+		    sfd = SF_SFF(sf,i)
+		    SFF_W(sfd) = SFF_W(sfd) * rval
+		}
+		SF_W(sf) = SF_W(sf) * rval
+		redraw = YES
+	    } else {
+		call printf ("scale %g\n")
+		    call pargr (SF_SCALE(sf))
+		redraw = NO
+	    }
+	case RADIUS:
+	    call gargr (rval)
+	    if (nscan() == 2) {
+		j = stf_r2i (rval) + 1
+		SF_RADIUS(sf) = rval
+		do i = 1, SF_NSFD(sf) {
+		    sfd = SF_SFD(sf,i)
+		    if (j > SFD_NPMAX(sfd))
+			next
+		    SFD_NP(sfd) = j
+		    SFD_RADIUS(sf) = SF_RADIUS(sf)
+		    call stf_norm (sf, sfd, INDEF, INDEF)
+		    call stf_widths (sf, sfd)
+		    call stf_fwhms (sf, sfd)
+		}
+		call stf_fitfocus (sf)
+		redraw = YES
+	    } else {
+		call printf ("radius %g\n")
+		    call pargr (SF_RADIUS(sf))
+		redraw = NO
+	    }
+	case XCENTER:
+	    call gargr (rval)
+	    if (nscan() == 2) {
+		if (IS_INDEF(rval))
+		    SF_XF(sf) = (SF_NCOLS(sf) + 1) / 2.
+		else
+		    SF_XF(sf) = rval
+		redraw = NO
+	    } else {
+		call printf ("xcenter %g\n")
+		    call pargr (SF_XF(sf))
+		redraw = NO
+	    }
+	case YCENTER:
+	    call gargr (rval)
+	    if (nscan() == 2) {
+		if (IS_INDEF(rval))
+		    SF_YF(sf) = (SF_NLINES(sf) + 1) / 2.
+		else
+		    SF_YF(sf) = rval
+		redraw = NO
+	    } else {
+		call printf ("ycenter %g\n")
+		    call pargr (SF_YF(sf))
+		redraw = NO
+	    }
+	case OVERPLOT:
+	    call gargb (bval)
+	    if (nscan() == 2) {
+		SF_OVRPLT(sf) = btoi (bval)
+		redraw = YES
+	    } else {
+		call printf ("overplot %b\n")
+		    call pargi (SF_OVRPLT(sf))
+		redraw = NO
+	    }
+	case BETA:
+	    call gargr (rval)
+	    if (nscan() == 2) {
+		SF_BETA(sf) = rval
+		do i = 1, SF_NSFD(sf) {
+		    sfd = SF_SFD(sf,i)
+		    call stf_widths (sf, sfd)
+		    switch (SF_WCODE(sf)) {
+		    case 1:
+			SFD_W(sfd) = SFD_R(sfd)
+		    case 2:
+			SFD_W(sfd) = SFD_DFWHM(sfd)
+		    case 3:
+			SFD_W(sfd) = SFD_GFWHM(sfd)
+		    case 4:
+			SFD_W(sfd) = SFD_MFWHM(sfd)
+		    }
+		    call stf_fwhms (sf, sfd)
+		}
+		call stf_fitfocus (sf) 
+		redraw = YES
+	    } else {
+		call printf ("beta %g\n")
+		    call pargr (SF_BETA(sf))
+		redraw = NO
+	    }
+	default:
+	    call printf ("Unrecognized or ambiguous command\007")
+	    redraw = NO
+	}
+
+	call sfree (sp)
+end
+
+
+# STF_G1 -- Plot of size/ellip vs. focus/mag/radius.
+
+procedure stf_g1 (sf, current, xkey, ykey, title, xlabel, ylabel)
+
+pointer	sf			#I Starfocus pointer
+pointer	current			#I Current sfd pointer
+int	xkey			#I X axis key
+int	ykey			#I Y axis key
+char	title[ARB]		#I Title
+char	xlabel[ARB]		#I X label
+char	ylabel[ARB]		#I Y label
+
+int	i, j
+bool	hl
+real	x, x1, x2, dx, y, y1, y2, dy
+pointer	gp, sff, sfd
+
+begin
+	# Determine data range
+	x1 = MAX_REAL
+	x2 = -MAX_REAL
+	switch (ykey) {
+	case 'r':
+	    y1 = SF_W(sf)
+	    y2 = 1.5 * SF_W(sf)
+	case 'e':
+	    y1 = 0
+	    y2 = 1
+	}
+	do j = 1, SF_NFOCUS(sf) {
+	    sff = SF_SFF(sf,j)
+	    if (xkey != 'f' && sff != SFD_SFF(current))
+		next
+	    do i = 1, SFF_NSFD(sff) {
+		sfd = SFF_SFD(sff,i)
+		if (SFD_STATUS(sfd) == 0) {
+		    switch (xkey) {
+		    case 'f':
+			x = SFD_F(sfd)
+		    case 'm':
+			x = -2.5 * log10 (SFS_M(SFD_SFS(sfd)) / SF_M(sf))
+		    case 't':
+			x = sqrt ((SFD_X(sfd) - SF_XF(sf)) ** 2 +
+			    (SFD_Y(sfd) - SF_YF(sf)) ** 2)
+		    }
+		    switch (ykey) {
+		    case 'r':
+			y = SFD_W(sfd)
+		    case 'e':
+			y = SFD_E(sfd)
+		    }
+		    x1 = min (x1, x)
+		    x2 = max (x2, x)
+		    y1 = min (y1, y)
+		    y2 = max (y2, y)
+		}
+	    }
+	}
+
+	dx = (x2 - x1)
+	dy = (y2 - y1)
+	x1 = x1 - dx * 0.05
+	x2 = x2 + dx * 0.05
+	y1 = y1 - dy * 0.05
+	y2 = y2 + dy * 0.05
+	gp = SF_GP (sf)
+	call gswind (gp, x1, x2, y1, y2)
+	call glabax (gp, title, xlabel, ylabel)
+
+	do j = 1, SF_NFOCUS(sf) {
+	    sff = SF_SFF(sf,j)
+	    if (xkey != 'f' && sff != SFD_SFF(current))
+		next
+	    do i = 1, SFF_NSFD(sff) {
+		sfd = SFF_SFD(sff,i)
+		if (SFD_STATUS(sfd) == 0) {
+		    hl = false
+		    switch (xkey) {
+		    case 'f':
+			x = SFD_F(sfd)
+			#hl = (SFD_SFS(sfd) == SFD_SFS(current))
+		    case 'm':
+			x = -2.5 * log10 (SFS_M(SFD_SFS(sfd)) / SF_M(sf))
+			#hl = (SFD_SFF(sfd) != SFD_SFF(current))
+		    case 't':
+			x = sqrt ((SFD_X(sfd) - SF_XF(sf)) ** 2 +
+			    (SFD_Y(sfd) - SF_YF(sf)) ** 2)
+			#hl = (SFD_SFF(sfd) != SFD_SFF(current))
+		    }
+		    switch (ykey) {
+		    case 'r':
+			y = SFD_W(sfd)
+		    case 'e':
+			y = SFD_E(sfd)
+		    }
+		    if (hl) {
+			call gseti (gp, G_PLCOLOR, HLCOLOR)
+			if (sfd == current)
+			    call gmark (gp, x, y, GM_BOX, 3., 3.)
+			call gmark (gp, x, y, GM_PLUS, 3., 3.)
+			call gseti (gp, G_PLCOLOR, 1)
+		    } else
+			call gmark (gp, x, y, GM_MARK, 2., 2.)
+		}
+	    }
+	}
+
+	call gseti (gp, G_PLTYPE, 2)
+	if (xkey == 'f')
+	    call gline (gp, SF_F(sf), y1, SF_F(sf), y2)
+	if (ykey == 'r')
+	    call gline (gp, x1, SF_W(sf), x2, SF_W(sf))
+	call gseti (gp, G_PLTYPE, 1)
+end
+
+
+# STF_G2 -- Enclosed flux profiles for a given star.
+
+procedure stf_g2 (sf, current, title)
+
+pointer	sf			#I Starfocus pointer
+pointer	current			#I Current sfd pointer
+char	title[ARB]		#I Title
+
+int	i, j, np, np1, nx, ny, ix, iy
+real	vx, dvx, vy, dvy, x1, x2, y1, y2, z, z1, r, r1, r2, dr, fa[10]
+pointer	sp, str, gp, sfs, sfd, asi
+real	stf_i2r(), stf_r2i(), asieval()
+
+begin
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	gp = SF_GP(sf)
+	sfs = SFD_SFS(current)
+	np = SFD_NP(current)
+
+	# Set grid layout
+	i = SFS_N(sfs)
+	if (i < 4) {
+	    nx = i
+	    ny = 1
+	} else {
+	    nx = nint (sqrt (real (i)))
+	    if (mod (i-1, nx+1) >= mod (i-1, nx))
+		nx = nx + 1
+	    ny = (i - 1) / nx + 1
+	}
+
+        # Set subview port parameters
+        call ggview (gp, vx, dvx, vy, dvy)
+        dvx = (dvx - vx) / nx
+        dvy = (dvy - vy) / ny
+
+	# Set data window parameters
+	x1 = -0.05
+	x2 = 1.05
+	y1 = -0.15
+	y2 = 1.05
+	call gswind (gp, x1, x2, y1, y2)
+
+	# Set fill area
+        fa[1] = x1; fa[6] = y1
+        fa[2] = x2; fa[7] = y1
+        fa[3] = x2; fa[8] = y2
+        fa[4] = x1; fa[9] = y2
+        fa[5] = x1; fa[10] = y1
+
+	# Draw profiles.
+	j = 0
+	do i = 1, SFS_NSFD(sfs) {
+	    sfd = SFS_SFD(sfs, i)
+	    if (SFD_STATUS(sfd) != 0)
+		next
+	    np1 = SFD_NP(sfd)
+	    ix = 1 + mod (j, nx)
+	    iy = 1 + j / nx
+	    j = j + 1
+	    call gsview (gp, vx+dvx*(ix-1), vx+dvx*ix,
+		vy+dvy*(ny-iy), vy+(ny-iy+1)*dvy)
+	    call gfill (gp, fa, fa[6], 4, GF_SOLID)
+	    call gseti (gp, G_DRAWTICKS, NO)
+	    call glabax (gp, "", "", "")
+	    if (sfd == current) {
+		call gsview (gp, vx+dvx*(ix-1)+.005, vx+dvx*ix-.005,
+		    vy+dvy*(ny-iy)+.005, vy+(ny-iy+1)*dvy-.005)
+		call gsetr (gp, G_PLWIDTH, HLWIDTH)
+		call gseti (gp, G_PLCOLOR, HLCOLOR)
+		call gpline (gp, fa, fa[6], 5)
+		call gsetr (gp, G_PLWIDTH, 1.)
+		call gseti (gp, G_PLCOLOR, 1)
+		call gsview (gp, vx+dvx*(ix-1), vx+dvx*ix,
+		    vy+dvy*(ny-iy), vy+(ny-iy+1)*dvy)
+	    }
+
+	    asi = SFD_ASI1(sfd)
+	    r2 = stf_i2r (real(np))
+	    call gamove (gp, 0., 0.)
+	    for (z = 1.; z <= np1; z = z + 0.1)
+		call gadraw (gp, stf_i2r(z)/r2, asieval(asi,z))
+	    if (SF_OVRPLT(sf) == YES && sfd != SF_BEST(sf)) {
+		call gseti (gp, G_PLCOLOR, HLCOLOR)
+		np1 = SFD_NP(SF_BEST(sf))
+		asi = SFD_ASI1(SF_BEST(sf))
+		r1 = stf_i2r (1.)
+		r2 = stf_i2r (real(np))
+		dr = 0.05 * (r2 - r1)
+		for (r = r1; r <= r2; r = r + dr) {
+		    z = stf_r2i (r)
+		    z1 = stf_r2i (r+0.7*dr)
+		    if (z > 1. && z1 <= np1)
+			call gline (gp, r/r2, asieval(asi,z),
+			    (r+0.7*dr)/r2, asieval(asi,z1))
+		}
+		call gseti (gp, G_PLCOLOR, 1)
+	    }
+
+	    call sprintf (Memc[str], SZ_LINE, "%.3g")
+		call pargr (SFD_W(sfd))
+	    call gtext (gp, 0.95, -0.1, Memc[str], "h=r;v=b")
+	    if (nx < NMAX && ny < NMAX) {
+		call sprintf (Memc[str], SZ_LINE, "%.4g")
+		    call pargr (SFD_F(sfd))
+		call gtext (gp, 0.05, -0.1, Memc[str], "h=l;v=b")
+	    }
+	}
+
+	call gsview (gp, vx, vx+nx*dvx, vy, vy+ny*dvy)
+	call gswind (gp, 0.5, 0.5+nx, 0.5+ny, 0.5)
+	call gamove (gp, 1., 1.)
+
+        # Draw label
+        call gseti (gp, G_DRAWAXES, 0)
+        call glabax (gp, title, "", "")
+        call gseti (gp, G_DRAWAXES, 3)
+
+	call sfree (sp)
+end
+
+
+# STF_G3 -- Enclosed flux profiles for a given focus.
+
+procedure stf_g3 (sf, current, title)
+
+pointer	sf			#I Starfocus pointer
+pointer	current			#I Current sfd pointer
+char	title[ARB]		#I Title
+
+int	i, j, np, np1, nx, ny, ix, iy
+real	vx, dvx, vy, dvy, x1, x2, y1, y2, z, z1, r, r1, r2, dr, fa[10]
+pointer	sp, str, gp, sff, sfd, asi
+real	stf_i2r(), stf_r2i(), asieval()
+
+begin
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	gp = SF_GP(sf)
+	sff = SFD_SFF(current)
+	np = SFD_NP(current)
+
+	# Set grid layout
+	i = SFF_N(sff)
+	if (i < 4) {
+	    nx = i
+	    ny = 1
+	} else {
+	    nx = nint (sqrt (real (i)))
+	    if (mod (i-1, nx+1) >= mod (i-1, nx))
+		nx = nx + 1
+	    ny = (i - 1) / nx + 1
+	}
+
+        # Set subview port parameters
+        call ggview (gp, vx, dvx, vy, dvy)
+        dvx = (dvx - vx) / nx
+        dvy = (dvy - vy) / ny
+
+	# Set data window parameters
+	x1 = -0.05
+	x2 = 1.05
+	y1 = -0.2
+	y2 = 1.05
+	call gswind (gp, x1, x2, y1, y2)
+
+	# Set fill area
+        fa[1] = x1; fa[6] = y1
+        fa[2] = x2; fa[7] = y1
+        fa[3] = x2; fa[8] = y2
+        fa[4] = x1; fa[9] = y2
+        fa[5] = x1; fa[10] = y1
+
+	# Draw profiles.
+	j = 0
+	do i = 1, SFF_NSFD(sff) {
+	    sfd = SFF_SFD(sff, i)
+	    if (SFD_STATUS(sfd) != 0)
+		next
+	    np1 = SFD_NP(sfd)
+	    ix = 1 + mod (j, nx)
+	    iy = 1 + j / nx
+	    j = j + 1
+	    call gsview (gp, vx+dvx*(ix-1), vx+dvx*ix,
+		vy+dvy*(ny-iy), vy+(ny-iy+1)*dvy)
+	    call gfill (gp, fa, fa[6], 4, GF_SOLID)
+	    call gseti (gp, G_DRAWTICKS, NO)
+	    call glabax (gp, "", "", "")
+	    if (sfd == current) {
+		call gsview (gp, vx+dvx*(ix-1)+.005, vx+dvx*ix-.005,
+		    vy+dvy*(ny-iy)+.005, vy+(ny-iy+1)*dvy-.005)
+		call gsetr (gp, G_PLWIDTH, HLWIDTH)
+		call gseti (gp, G_PLCOLOR, HLCOLOR)
+		call gpline (gp, fa, fa[6], 5)
+		call gsetr (gp, G_PLWIDTH, 1.)
+		call gseti (gp, G_PLCOLOR, 1)
+		call gsview (gp, vx+dvx*(ix-1), vx+dvx*ix,
+		    vy+dvy*(ny-iy), vy+(ny-iy+1)*dvy)
+	    }
+
+	    asi = SFD_ASI1(sfd)
+	    r2 = stf_i2r (real(np))
+	    call gamove (gp, 0., 0.)
+	    for (z = 1.; z <= np1; z = z + 0.1)
+		call gadraw (gp, stf_i2r(z)/r2, asieval(asi,z))
+	    if (SF_OVRPLT(sf) == YES && sfd != SF_BEST(sf)) {
+		call gseti (gp, G_PLCOLOR, HLCOLOR)
+		np1 = SFD_NP(SF_BEST(sf))
+		asi = SFD_ASI1(SF_BEST(sf))
+		r1 = stf_i2r (1.)
+		r2 = stf_i2r (real(np))
+		dr = 0.05 * (r2 - r1)
+		for (r = r1; r <= r2; r = r + dr) {
+		    z = stf_r2i (r)
+		    z1 = stf_r2i (r+0.7*dr)
+		    if (z > 1. && z1 <= np1)
+			call gline (gp, r/r2, asieval(asi,z),
+			    (r+0.7*dr)/r2, asieval(asi,z1))
+		}
+		call gseti (gp, G_PLCOLOR, 1)
+	    }
+
+	    call sprintf (Memc[str], SZ_LINE, "%.3g")
+		call pargr (SFD_W(sfd))
+	    call gtext (gp, 0.95, -.1, Memc[str], "h=r;v=b")
+	    if (nx < NMAX && ny < NMAX) {
+		call sprintf (Memc[str], SZ_LINE, "%d %d")
+		    call pargr (SFD_X(sfd))
+		    call pargr (SFD_Y(sfd))
+		call gtext (gp, 0.05, -.1, Memc[str], "h=l;v=b")
+	    }
+	}
+
+	call gsview (gp, vx, vx+nx*dvx, vy, vy+ny*dvy)
+	call gswind (gp, 0.5, 0.5+nx, 0.5+ny, 0.5)
+	call gamove (gp, 1., 1.)
+
+        # Draw label
+        call gseti (gp, G_DRAWAXES, 0)
+        call glabax (gp, title, "", "")
+        call gseti (gp, G_DRAWAXES, 3)
+
+	call sfree (sp)
+end
+
+
+# STF_G4 -- Radial profiles (derivative of enclosed flux) for a given star.
+
+procedure stf_g4 (sf, current, title)
+
+pointer	sf			#I Starfocus pointer
+pointer	current			#I Current sfd pointer
+char	title[ARB]		#I Title
+
+int	i, j, np, np1, nx, ny, ix, iy
+real	vx, dvx, vy, dvy, x1, x2, y1, y2, z, z1, r, r1, r2, dr, rmax, fa[10]
+pointer	sp, str, gp, sfs, sfd, asi
+real	stf_i2r(), stf_r2i(), asieval()
+
+begin
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	gp = SF_GP(sf)
+	sfs = SFD_SFS(current)
+	np = SFD_NP(current)
+
+	# Set grid layout
+	i = SFS_N(sfs)
+	if (i < 4) {
+	    nx = i
+	    ny = 1
+	} else {
+	    nx = nint (sqrt (real (i)))
+	    if (mod (i-1, nx+1) >= mod (i-1, nx))
+		nx = nx + 1
+	    ny = (i - 1) / nx + 1
+	}
+
+        # Set subview port parameters
+        call ggview (gp, vx, dvx, vy, dvy)
+        dvx = (dvx - vx) / nx
+        dvy = (dvy - vy) / ny
+
+	# Set data window parameters
+	x1 = -0.05
+	x2 = 1.05
+	z = SF_YP2(sf) - SF_YP1(sf)
+	y1 = SF_YP1(sf) - 0.05 * z
+	y2 = SF_YP2(sf) + 0.15 * z
+
+	# Set fill area
+        fa[1] = x1; fa[6] = y1
+        fa[2] = x2; fa[7] = y1
+        fa[3] = x2; fa[8] = y2
+        fa[4] = x1; fa[9] = y2
+        fa[5] = x1; fa[10] = y1
+
+	# Draw profiles.
+	j = 0
+	do i = 1, SFS_NSFD(sfs) {
+	    sfd = SFS_SFD(sfs, i)
+	    if (SFD_STATUS(sfd) != 0)
+		next
+	    np1 = SFD_NP(sfd)
+	    ix = 1 + mod (j, nx)
+	    iy = 1 + j / nx
+	    j = j + 1
+	    call gsview (gp, vx+dvx*(ix-1), vx+dvx*ix,
+		vy+dvy*(ny-iy), vy+(ny-iy+1)*dvy)
+	    call gswind (gp, x1, x2, y1, y2)
+	    call gfill (gp, fa, fa[6], 4, GF_SOLID)
+	    call gseti (gp, G_DRAWTICKS, NO)
+	    call glabax (gp, "", "", "")
+	    if (sfd == current) {
+		call gsview (gp, vx+dvx*(ix-1)+.005, vx+dvx*ix-.005,
+		    vy+dvy*(ny-iy)+.005, vy+(ny-iy+1)*dvy-.005)
+		call gsetr (gp, G_PLWIDTH, HLWIDTH)
+		call gseti (gp, G_PLCOLOR, HLCOLOR)
+		call gpline (gp, fa, fa[6], 5)
+		call gsetr (gp, G_PLWIDTH, 1.)
+		call gseti (gp, G_PLCOLOR, 1)
+		call gsview (gp, vx+dvx*(ix-1), vx+dvx*ix,
+		    vy+dvy*(ny-iy), vy+(ny-iy+1)*dvy)
+	    }
+
+	    asi = SFD_ASI2(sfd)
+	    rmax = stf_i2r (real(np))
+	    z = SF_XP1(sf)
+	    call gamove (gp, stf_i2r(z)/rmax, asieval(asi,z))
+	    for (; z <= SF_XP2(sf); z = z + 0.1)
+		call gadraw (gp, stf_i2r(z)/rmax, asieval(asi,z))
+	    if (SF_OVRPLT(sf) == YES && sfd != SF_BEST(sf)) {
+		call gseti (gp, G_PLCOLOR, HLCOLOR)
+		np1 = SFD_NP(SF_BEST(sf))
+		asi = SFD_ASI2(SF_BEST(sf))
+		rmax = stf_i2r (real(np))
+		r1 = stf_i2r (SF_XP1(sf))
+		r2 = stf_i2r (SF_XP2(sf))
+		dr = 0.05 * (rmax - stf_i2r(1.))
+		for (r = r1; r <= r2; r = r + dr) {
+		    z = stf_r2i (r)
+		    z1 = stf_r2i (r+0.7*dr)
+		    if (z > 1. && z1 <= np1)
+			call gline (gp, r/rmax, asieval(asi,z),
+			    (r+0.7*dr)/rmax, asieval(asi,z1))
+		}
+		call gseti (gp, G_PLCOLOR, 1)
+	    }
+
+	    call gswind (gp, 0., 1., 0., 1.)
+	    call sprintf (Memc[str], SZ_LINE, "%.3g")
+		call pargr (SFD_W(sfd))
+	    call gtext (gp, 0.95, 0.98, Memc[str], "h=r;v=t")
+	    if (nx < NMAX && ny < NMAX) {
+		call sprintf (Memc[str], SZ_LINE, "%.4g")
+		    call pargr (SFD_F(sfd))
+		call gtext (gp, 0.05, 0.98, Memc[str], "h=l;v=t")
+	    }
+	}
+
+	call gsview (gp, vx, vx+nx*dvx, vy, vy+ny*dvy)
+	call gswind (gp, 0.5, 0.5+nx, 0.5+ny, 0.5)
+	call gamove (gp, 1., 1.)
+
+        # Draw label
+        call gseti (gp, G_DRAWAXES, 0)
+        call glabax (gp, title, "", "")
+        call gseti (gp, G_DRAWAXES, 3)
+
+	call sfree (sp)
+end
+
+
+# STF_G5 -- Radial profiles (derivative of enclosed flux) for a given focus.
+
+procedure stf_g5 (sf, current, title)
+
+pointer	sf			#I Starfocus pointer
+pointer	current			#I Current sfd pointer
+char	title[ARB]		#I Title
+
+int	i, j, np, np1, nx, ny, ix, iy
+real	vx, dvx, vy, dvy, x1, x2, y1, y2, z, z1, r, r1, r2, dr, rmax, fa[10]
+pointer	sp, str, gp, sff, sfd, asi
+real	stf_i2r(), stf_r2i(), asieval()
+
+begin
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	gp = SF_GP(sf)
+	sff = SFD_SFF(current)
+	np = SFD_NP(current)
+
+	# Set grid layout
+	i = SFF_N(sff)
+	if (i < 4) {
+	    nx = i
+	    ny = 1
+	} else {
+	    nx = nint (sqrt (real (i)))
+	    if (mod (i-1, nx+1) >= mod (i-1, nx))
+		nx = nx + 1
+	    ny = (i - 1) / nx + 1
+	}
+
+        # Set subview port parameters
+        call ggview (gp, vx, dvx, vy, dvy)
+        dvx = (dvx - vx) / nx
+        dvy = (dvy - vy) / ny
+
+	# Set data window parameters
+	x1 = -0.05
+	x2 = 1.05
+	z = SF_YP2(sf) - SF_YP1(sf)
+	y1 = SF_YP1(sf) - 0.05 * z
+	y2 = SF_YP2(sf) + 0.15 * z
+
+	# Set fill area
+        fa[1] = x1; fa[6] = y1
+        fa[2] = x2; fa[7] = y1
+        fa[3] = x2; fa[8] = y2
+        fa[4] = x1; fa[9] = y2
+        fa[5] = x1; fa[10] = y1
+
+	# Draw profiles.
+	j = 0
+	do i = 1, SFF_NSFD(sff) {
+	    sfd = SFF_SFD(sff, i)
+	    if (SFD_STATUS(sfd) != 0)
+		next
+	    np1 = SFD_NP(sfd)
+	    ix = 1 + mod (j, nx)
+	    iy = 1 + j / nx
+	    j = j + 1
+	    call gsview (gp, vx+dvx*(ix-1), vx+dvx*ix,
+		vy+dvy*(ny-iy), vy+(ny-iy+1)*dvy)
+	    call gswind (gp, x1, x2, y1, y2)
+	    call gfill (gp, fa, fa[6], 4, GF_SOLID)
+	    call gseti (gp, G_DRAWTICKS, NO)
+	    call glabax (gp, "", "", "")
+	    if (sfd == current) {
+		call gsview (gp, vx+dvx*(ix-1)+.005, vx+dvx*ix-.005,
+		    vy+dvy*(ny-iy)+.005, vy+(ny-iy+1)*dvy-.005)
+		call gsetr (gp, G_PLWIDTH, HLWIDTH)
+		call gseti (gp, G_PLCOLOR, HLCOLOR)
+		call gpline (gp, fa, fa[6], 5)
+		call gsetr (gp, G_PLWIDTH, 1.)
+		call gseti (gp, G_PLCOLOR, 1)
+		call gsview (gp, vx+dvx*(ix-1), vx+dvx*ix,
+		    vy+dvy*(ny-iy), vy+(ny-iy+1)*dvy)
+	    }
+
+	    asi = SFD_ASI2(sfd)
+	    rmax = stf_i2r (real(np))
+	    z = SF_XP1(sf)
+	    call gamove (gp, stf_i2r(z)/rmax, asieval(asi,z))
+	    for (; z <= SF_XP2(sf); z = z + 0.1)
+		call gadraw (gp, stf_i2r(z)/rmax, asieval(asi,z))
+	    if (SF_OVRPLT(sf) == YES && sfd != SF_BEST(sf)) {
+		call gseti (gp, G_PLCOLOR, HLCOLOR)
+		np1 = SFD_NP(SF_BEST(sf))
+		asi = SFD_ASI2(SF_BEST(sf))
+		rmax = stf_i2r (real(np))
+		r1 = stf_i2r (SF_XP1(sf))
+		r2 = stf_i2r (SF_XP2(sf))
+		dr = 0.05 * (rmax - stf_i2r (1.))
+		for (r = r1; r <= r2; r = r + dr) {
+		    z = stf_r2i (r)
+		    z1 = stf_r2i (r+0.7*dr)
+		    if (z > 1. && z1 <= np1)
+			call gline (gp, r/rmax, asieval(asi,z),
+			    (r+0.7*dr)/rmax, asieval(asi,z1))
+		}
+		call gseti (gp, G_PLCOLOR, 1)
+	    }
+
+	    call gswind (gp, 0., 1., 0., 1.)
+	    call sprintf (Memc[str], SZ_LINE, "%.3g")
+		call pargr (SFD_W(sfd))
+	    call gtext (gp, 0.95, 0.98, Memc[str], "h=r;v=t")
+	    if (nx < NMAX && ny < NMAX) {
+		call sprintf (Memc[str], SZ_LINE, "%d %d")
+		    call pargr (SFD_X(sfd))
+		    call pargr (SFD_Y(sfd))
+		call gtext (gp, 0.05, 0.98, Memc[str], "h=l;v=t")
+	    }
+	}
+
+	call gsview (gp, vx, vx+nx*dvx, vy, vy+ny*dvy)
+	call gswind (gp, 0.5, 0.5+nx, 0.5+ny, 0.5)
+	call gamove (gp, 1., 1.)
+
+        # Draw label
+        call gseti (gp, G_DRAWAXES, 0)
+        call glabax (gp, title, "", "")
+        call gseti (gp, G_DRAWAXES, 3)
+
+	call sfree (sp)
+end
+
+
+# STF_G6 -- Enclosed flux profile of a star.
+
+procedure stf_g6 (sf, current, title, xlabel, ylabel)
+
+pointer	sf			#I Starfocus pointer
+pointer	current			#I Star pointer
+char	title[ARB]		#I Title
+char	xlabel[ARB]		#I X label
+char	ylabel[ARB]		#I Y label
+
+int	np, np1
+real	scale, level, radius, flux, profile
+pointer	gp, asi
+
+real	x1, x2, y1, y2, z, z1, r, r1, r2, dr
+real	stf_i2r(), stf_r2i(), asieval()
+
+begin
+	gp = SF_GP(sf)
+	level = SF_LEVEL(sf)
+	scale = SF_SCALE(sf)
+	np = SFD_NP(current)
+	asi = SFD_ASI1(current)
+
+	x1 = -0.5 * scale
+	x2 = (stf_i2r (real(np)) + 0.5) * scale
+	y1 = -0.05
+	y2 = 1.05
+	call gswind (gp, x1, x2, y1, y2)
+
+        call gseti (gp, G_DRAWTICKS, YES)
+	call gseti (gp, G_XNMAJOR, 6)
+	call gseti (gp, G_XNMINOR, 4)
+	call gseti (gp, G_YNMAJOR, 6)
+	call gseti (gp, G_YNMINOR, 4)
+	call glabax (gp, title, xlabel, ylabel)
+
+	# Draw profiles.
+	if (SFD_STATUS(current) == 0) {
+	    call gseti (gp, G_PLCOLOR, 1)
+	    for (z = 1.; z <= np; z = z + 1)
+		call gmark (gp, stf_i2r(z)*scale, asieval(asi,z),
+		    GM_PLUS, 2., 2.)
+	    call gamove (gp, 0., 0.)
+	    for (z = 1.; z <= np; z = z + 0.1)
+		call gadraw (gp, stf_i2r(z)*scale, asieval(asi,z))
+	    switch (SF_WCODE(sf)) {
+	    case 1:
+		radius = SFD_W(current)
+		call gseti (gp, G_PLTYPE, 2)
+		call gline (gp, x1, level, radius, level)
+		call gline (gp, radius, level, radius, y1)
+		call gseti (gp, G_PLTYPE, 1)
+	    default:
+		radius = SFD_W(current) / 2.
+		call gseti (gp, G_PLTYPE, 2)
+		call gline (gp, radius, y1, radius, y2)
+		call gseti (gp, G_PLTYPE, 1)
+	    }
+
+	    call gseti (gp, G_PLCOLOR, HLCOLOR)
+	    call stf_model (sf, current, 0., profile, flux)
+	    call gamove (gp, 0., flux)
+	    for (z = 1.; z <= np; z = z + 0.1) {
+		r = stf_i2r(z) * scale
+		call stf_model (sf, current, r, profile, flux)
+		call gadraw (gp, r, flux)
+	    }
+	    call gseti (gp, G_PLCOLOR, 1)
+	    if (SF_OVRPLT(sf) == YES && current != SF_BEST(sf)) {
+		call gseti (gp, G_PLCOLOR, HLCOLOR)
+		np1 = SFD_NP(SF_BEST(sf))
+		asi = SFD_ASI1(SF_BEST(sf))
+		r1 = stf_i2r(1.)
+		r2 = stf_i2r (real(np))
+		dr = 0.05 * (r2 - r1)
+		for (r = r1; r <= r2; r = r + dr) {
+		    z = stf_r2i (r)
+		    z1 = stf_r2i (r+0.7*dr)
+		    if (z > 1. && z1 <= np1)
+			call gline (gp, r*scale, asieval(asi,z),
+			    (r+0.7*dr)*scale, asieval(asi,z1))
+		}
+		call gseti (gp, G_PLCOLOR, 1)
+	    }
+	}
+end
+
+
+# STF_G7 -- Radial profile (derivative of enclosed flux) for a star.
+
+procedure stf_g7 (sf, current, title, xlabel, ylabel)
+
+pointer	sf			#I Starfocus pointer
+pointer	current			#I Star pointer
+char	title[ARB]		#I Title
+char	xlabel[ARB]		#I X label
+char	ylabel[ARB]		#I Y label
+
+int	np, np1
+real	scale, level, radius, profile, flux
+pointer	gp, asi
+
+real	x1, x2, y1, y2, z, z1, r, r1, r2, dr
+real	stf_i2r(), stf_r2i(), asieval()
+
+begin
+	gp = SF_GP(sf)
+	level = SF_LEVEL(sf)
+	scale = SF_SCALE(sf)
+	np = SFD_NP(current)
+	asi = SFD_ASI2(current)
+
+	x1 = -0.5 * scale
+	x2 = (stf_i2r (real(np)) + 0.5) * scale
+	z = SFD_YP2(current) - SFD_YP1(current)
+	y1 = SFD_YP1(current) - 0.05 * z
+	y2 = SFD_YP2(current) + 0.05 * z
+	call gswind (gp, x1, x2, y1, y2)
+
+        call gseti (gp, G_XDRAWTICKS, YES)
+        call gseti (gp, G_YDRAWTICKS, NO)
+	call gseti (gp, G_XNMAJOR, 6)
+	call gseti (gp, G_XNMINOR, 4)
+	call gseti (gp, G_YNMAJOR, 6)
+	call gseti (gp, G_YNMINOR, 4)
+	call glabax (gp, title, xlabel, ylabel)
+
+	# Draw profile
+	call gseti (gp, G_PLCOLOR, 1)
+	for (z = SF_XP1(sf); z <= SF_XP2(sf); z = z + 1)
+	    call gmark (gp, stf_i2r(z)*scale, asieval(asi,z),
+		GM_PLUS, 2., 2.)
+	z = SF_XP1(sf)
+	call gamove (gp, stf_i2r(z)*scale, asieval(asi,z))
+	for (; z <= SF_XP2(sf); z = z + 0.1)
+	    call gadraw (gp, stf_i2r (z)*scale, asieval(asi,z))
+
+	switch (SF_WCODE(sf)) {
+	case 1:
+	    radius = SFD_W(current)
+	default:
+	    radius = SFD_W(current) / 2.
+	}
+	call gseti (gp, G_PLTYPE, 2)
+	call gline (gp, radius, y1, radius, y2)
+	call gseti (gp, G_PLTYPE, 1)
+
+	call gseti (gp, G_PLCOLOR, HLCOLOR)
+	z = SF_XP1(sf)
+	r = stf_i2r(z) * scale
+	call stf_model (sf, current, r, profile, flux)
+	call gamove (gp, r, profile)
+	for (; z <= np; z = z + 0.1) {
+	    r = stf_i2r(z) * scale
+	    call stf_model (sf, current, r, profile, flux)
+	    call gadraw (gp, r, profile)
+	}
+	call gseti (gp, G_PLCOLOR, 1)
+	if (SF_OVRPLT(sf) == YES && current != SF_BEST(sf)) {
+	    call gseti (gp, G_PLCOLOR, HLCOLOR)
+	    np1 = SFD_NP(SF_BEST(sf))
+	    asi = SFD_ASI2(SF_BEST(sf))
+	    r1 = stf_i2r (SF_XP1(sf))
+	    r2 = stf_i2r (SF_XP2(sf))
+	    dr = 0.05 * (r2 - r1)
+	    for (r = r1; r <= r2; r = r + dr) {
+		z = stf_r2i (r)
+		z1 = stf_r2i (r+0.7*dr)
+		if (z > 1. && z1 <= np1)
+		    call gline (gp, r*scale, asieval(asi,z),
+			(r+0.7*dr)*scale, asieval(asi,z1))
+	    }
+	    call gseti (gp, G_PLCOLOR, 1)
+	}
+end
+
+
+# STF_G8 -- FWHM vs level.
+
+procedure stf_g8 (sf, current, title, xlabel, ylabel)
+
+pointer	sf			#I Starfocus pointer
+pointer	current			#I Star pointer
+char	title[ARB]		#I Title
+char	xlabel[ARB]		#I X label
+char	ylabel[ARB]		#I Y label
+
+real	y1, y2, level, fwhm
+pointer	gp
+
+begin
+	level = SF_LEVEL(sf)
+	if (SF_WCODE(sf) == 1)
+	    fwhm = SFD_MFWHM(current)
+	else
+	    fwhm = SFD_W(current)
+
+	call alimr (SFD_FWHM(current,2), 17, y1, y2)
+	y2 = y2 - y1
+	y1 = y1 - 0.05 * y2
+	y2 = y1 + 1.10 * y2
+	y1 = min (y1, 0.9 * fwhm)
+	y2 = max (y2, 1.1 * fwhm)
+
+	gp = SF_GP(sf)
+        call gseti (gp, G_DRAWTICKS, YES)
+	call gseti (gp, G_XNMAJOR, 6)
+	call gseti (gp, G_XNMINOR, 4)
+	call gseti (gp, G_YNMAJOR, 6)
+	call gseti (gp, G_YNMINOR, 4)
+	call gswind (gp, 0., 1., y1, y2)
+	call glabax (gp, title, xlabel, ylabel)
+	call gvline (gp, SFD_FWHM(current,2), 17, 0.1, 0.9)
+	call gvmark (gp, SFD_FWHM(current,2), 17, 0.1, 0.9, GM_PLUS, 2., 2.)
+
+	switch (SF_WCODE(sf)) {
+	case 1:
+	    call gseti (gp, G_PLTYPE, 2)
+	    call gline (gp, 0., fwhm, level, fwhm)
+	    call gline (gp, level, y1, level, fwhm)
+	    call gseti (gp, G_PLTYPE, 1)
+	default:
+	    call gseti (gp, G_PLTYPE, 2)
+	    call gline (gp, 0., fwhm, 1., fwhm)
+	    call gseti (gp, G_PLTYPE, 1)
+	}
+end
+
+
+# STF_G9 -- FWHM vs level for a given star.
+
+procedure stf_g9 (sf, current, title)
+
+pointer	sf			#I Starfocus pointer
+pointer	current			#I Current sfd pointer
+char	title[ARB]		#I Title
+
+int	i, j, nx, ny, ix, iy
+real	level, fwhm, vx, dvx, vy, dvy, x1, x2, y1, y2, fa[10]
+pointer	sp, str, gp, sfs, sfd
+
+begin
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	gp = SF_GP(sf)
+	sfs = SFD_SFS(current)
+	level = SF_LEVEL(sf)
+	if (SF_WCODE(sf) == 1)
+	    fwhm = SFD_MFWHM(current)
+	else
+	    fwhm = SFD_W(current)
+
+	# Set grid layout
+	i = SFS_N(sfs)
+	if (i < 4) {
+	    nx = i
+	    ny = 1
+	} else {
+	    nx = nint (sqrt (real (i)))
+	    if (mod (i-1, nx+1) >= mod (i-1, nx))
+		nx = nx + 1
+	    ny = (i - 1) / nx + 1
+	}
+
+        # Set subview port parameters
+        call ggview (gp, vx, dvx, vy, dvy)
+        dvx = (dvx - vx) / nx
+        dvy = (dvy - vy) / ny
+
+	# Set data window parameters
+	y1 = 0.9 * fwhm
+	y2 = 1.1 * fwhm
+	do i = 1, SFS_NSFD(sfs) {
+	    sfd = SFS_SFD(sfs,i)
+	    if (SFD_STATUS(sfd) != 0)
+		next
+	    call alimr (SFD_FWHM(sfd,2), 17, x1, x2)
+	    x2 = x2 - x1
+	    x1 = x1 - 0.05 * x2
+	    x2 = x1 + 1.10 * x2
+	    y1 = min (x1, y1)
+	    y2 = max (x2, y2)
+	}
+	x2 = y2 - y1
+	y1 = min (y1, fwhm - 0.2 * x2)
+	y2 = max (y2, fwhm + 0.2 * x2)
+	x1 = 0.
+	x2 = 1.
+	call gswind (gp, x1, x2, y1, y2)
+
+	# Set fill area
+        fa[1] = x1; fa[6] = y1
+        fa[2] = x2; fa[7] = y1
+        fa[3] = x2; fa[8] = y2
+        fa[4] = x1; fa[9] = y2
+        fa[5] = x1; fa[10] = y1
+
+	# Draw profiles.
+	j = 0
+	do i = 1, SFS_NSFD(sfs) {
+	    sfd = SFS_SFD(sfs, i)
+	    if (SFD_STATUS(sfd) != 0)
+		next
+	    ix = 1 + mod (j, nx)
+	    iy = 1 + j / nx
+	    j = j + 1
+	    call gsview (gp, vx+dvx*(ix-1), vx+dvx*ix,
+		vy+dvy*(ny-iy), vy+(ny-iy+1)*dvy)
+	    call gfill (gp, fa, fa[6], 4, GF_SOLID)
+	    call gseti (gp, G_DRAWTICKS, NO)
+	    call glabax (gp, "", "", "")
+	    if (sfd == current) {
+		call gsview (gp, vx+dvx*(ix-1)+.005, vx+dvx*ix-.005,
+		    vy+dvy*(ny-iy)+.005, vy+(ny-iy+1)*dvy-.005)
+		call gsetr (gp, G_PLWIDTH, HLWIDTH)
+		call gseti (gp, G_PLCOLOR, HLCOLOR)
+		call gpline (gp, fa, fa[6], 5)
+		call gsetr (gp, G_PLWIDTH, 1.)
+		call gseti (gp, G_PLCOLOR, 1)
+		call gsview (gp, vx+dvx*(ix-1), vx+dvx*ix,
+		    vy+dvy*(ny-iy), vy+(ny-iy+1)*dvy)
+	    }
+
+	    call gvline (gp, SFD_FWHM(sfd,2), 17, 0.1, 0.9)
+	    #call gseti (gp, G_PLTYPE, 2)
+	    #call gline (gp, x1, fwhm, x2, fwhm)
+	    #call gseti (gp, G_PLTYPE, 1)
+
+	    call sprintf (Memc[str], SZ_LINE, "%.3g")
+		call pargr (SFD_W(sfd))
+	    call gtext (gp, 0.95, 0.95*y2+0.05*y1, Memc[str], "h=r;v=t")
+	    if (nx < NMAX && ny < NMAX) {
+		call sprintf (Memc[str], SZ_LINE, "%.4g")
+		    call pargr (SFD_F(sfd))
+		call gtext (gp, 0.05, 0.95*y2+0.05*y1, Memc[str], "h=l;v=t")
+	    }
+	}
+
+	call gsview (gp, vx, vx+nx*dvx, vy, vy+ny*dvy)
+	call gswind (gp, 0.5, 0.5+nx, 0.5+ny, 0.5)
+	call gamove (gp, 1., 1.)
+
+        # Draw label
+        call gseti (gp, G_DRAWAXES, 0)
+        call glabax (gp, title, "", "")
+        call gseti (gp, G_DRAWAXES, 3)
+
+	call sfree (sp)
+end
+
+
+# STF_G10 -- FWHM vs level for a given focus.
+
+procedure stf_g10 (sf, current, title)
+
+pointer	sf			#I Starfocus pointer
+pointer	current			#I Current sfd pointer
+char	title[ARB]		#I Title
+
+int	i, j, nx, ny, ix, iy
+real	level, fwhm, vx, dvx, vy, dvy, x1, x2, y1, y2, fa[10]
+pointer	sp, str, gp, sff, sfd
+
+begin
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	gp = SF_GP(sf)
+	sff = SFD_SFF(current)
+	level = SF_LEVEL(sf)
+	if (SF_WCODE(sf) == 1)
+	    fwhm = SFD_MFWHM(current)
+	else
+	    fwhm = SFD_W(current)
+
+	# Set grid layout
+	i = SFF_N(sff)
+	if (i < 4) {
+	    nx = i
+	    ny = 1
+	} else {
+	    nx = nint (sqrt (real (i)))
+	    if (mod (i-1, nx+1) >= mod (i-1, nx))
+		nx = nx + 1
+	    ny = (i - 1) / nx + 1
+	}
+
+        # Set subview port parameters
+        call ggview (gp, vx, dvx, vy, dvy)
+        dvx = (dvx - vx) / nx
+        dvy = (dvy - vy) / ny
+
+	# Set data window parameters
+	y1 = 0.9 * fwhm
+	y2 = 1.1 * fwhm
+	do i = 1, SFF_NSFD(sff) {
+	    sfd = SFF_SFD(sff,i)
+	    if (SFD_STATUS(sfd) != 0)
+		next
+	    call alimr (SFD_FWHM(sfd,2), 17, x1, x2)
+	    x2 = x2 - x1
+	    x1 = x1 - 0.05 * x2
+	    x2 = x1 + 1.10 * x2
+	    y1 = min (x1, y1)
+	    y2 = max (x2, y2)
+	}
+	x2 = y2 - y1
+	y1 = min (y1, fwhm - 0.2 * x2)
+	y2 = max (y2, fwhm + 0.2 * x2)
+	x1 = 0.
+	x2 = 1.
+	call gswind (gp, x1, x2, y1, y2)
+
+	# Set fill area
+        fa[1] = x1; fa[6] = y1
+        fa[2] = x2; fa[7] = y1
+        fa[3] = x2; fa[8] = y2
+        fa[4] = x1; fa[9] = y2
+        fa[5] = x1; fa[10] = y1
+
+	# Draw plots.
+	j = 0
+	do i = 1, SFF_NSFD(sff) {
+	    sfd = SFF_SFD(sff, i)
+	    if (SFD_STATUS(sfd) != 0)
+		next
+	    ix = 1 + mod (j, nx)
+	    iy = 1 + j / nx
+	    j = j + 1
+	    call gsview (gp, vx+dvx*(ix-1), vx+dvx*ix,
+		vy+dvy*(ny-iy), vy+(ny-iy+1)*dvy)
+	    call gfill (gp, fa, fa[6], 4, GF_SOLID)
+	    call gseti (gp, G_DRAWTICKS, NO)
+	    call glabax (gp, "", "", "")
+	    if (sfd == current) {
+		call gsview (gp, vx+dvx*(ix-1)+.005, vx+dvx*ix-.005,
+		    vy+dvy*(ny-iy)+.005, vy+(ny-iy+1)*dvy-.005)
+		call gsetr (gp, G_PLWIDTH, HLWIDTH)
+		call gseti (gp, G_PLCOLOR, HLCOLOR)
+		call gpline (gp, fa, fa[6], 5)
+		call gsetr (gp, G_PLWIDTH, 1.)
+		call gseti (gp, G_PLCOLOR, 1)
+		call gsview (gp, vx+dvx*(ix-1), vx+dvx*ix,
+		    vy+dvy*(ny-iy), vy+(ny-iy+1)*dvy)
+	    }
+
+	    call gvline (gp, SFD_FWHM(sfd,2), 17, 0.1, 0.9)
+	    #call gseti (gp, G_PLTYPE, 2)
+	    #call gline (gp, x1, fwhm, x2, fwhm)
+	    #call gseti (gp, G_PLTYPE, 1)
+
+	    call sprintf (Memc[str], SZ_LINE, "%.3g")
+		call pargr (SFD_W(sfd))
+	    call gtext (gp, 0.95, 0.95*y2+0.05*y1, Memc[str], "h=r;v=t")
+	    if (nx < NMAX && ny < NMAX) {
+		call sprintf (Memc[str], SZ_LINE, "%d %d")
+		    call pargr (SFD_X(sfd))
+		    call pargr (SFD_Y(sfd))
+		call gtext (gp, 0.05, 0.95*y2+0.05*y1, Memc[str], "h=l;v=t")
+	    }
+	}
+
+	call gsview (gp, vx, vx+nx*dvx, vy, vy+ny*dvy)
+	call gswind (gp, 0.5, 0.5+nx, 0.5+ny, 0.5)
+	call gamove (gp, 1., 1.)
+
+        # Draw label
+        call gseti (gp, G_DRAWAXES, 0)
+        call glabax (gp, title, "", "")
+        call gseti (gp, G_DRAWAXES, 3)
+
+	call sfree (sp)
+end
+
+
+# STF_G11 -- Spatial plot at one focus.
+
+procedure stf_g11 (sf, current, key, title)
+
+pointer	sf			#I Starfocus pointer
+pointer	current			#I Current sfd pointer
+int	key			#I Plot magnitude symbol?
+char	title[ARB]		#I Title
+
+int	i
+real	x, y, z, x1, x2, y1, y2, rbest, rmin, rmax, emin, emax
+real	vx[3,2], vy[3,2], dvx, dvy, fa[8]
+pointer	gp, sfd, sff
+
+data	fa/0.,1.,1.,0.,0.,0.,1.,1./
+
+begin
+	gp = SF_GP(sf)
+	sff = SFD_SFF(current)
+
+	# Range of X, Y, R, E.
+	x1 = 1.
+	y1 = 1.
+	x2 = SF_NCOLS(sf)
+	y2 = SF_NLINES(sf)
+
+	rbest = SFD_W(SF_BEST(sf))
+	rmin = SF_W(sf)
+	rmax = 1.5 * SF_W(sf)
+	emin = 0
+	emax = 1
+	do i = 1, SFF_NSFD(sff) {
+	    sfd = SFF_SFD(sff,i)
+	    if (SFD_STATUS(sfd) != 0)
+		next
+	    rmin = min (rmin, SFD_W(sfd))
+	    rmax = max (rmax, SFD_W(sfd))
+	    emin = min (emin, SFD_E(sfd))
+	    emax = max (emax, SFD_E(sfd))
+	}
+	z = rmax - rmin
+	rmin = rmin - 0.1 * z
+	rmax = rmax + 0.1 * z
+
+	# Set view ports
+	call ggview (gp, vx[1,1], vx[3,2], vy[1,1], vy[3,2])
+	dvx = vx[3,2] - vx[1,1]
+	dvy = vy[3,2] - vy[1,1]
+	vx[1,1] = vx[1,1] + 0.00 * dvx
+	vx[1,2] = vx[1,1] + 0.20 * dvx
+	vx[2,1] = vx[1,1] + 0.25 * dvx
+	vx[2,2] = vx[1,1] + 0.75 * dvx
+	vx[3,1] = vx[1,1] + 0.80 * dvx
+	vx[3,2] = vx[1,1] + 1.00 * dvx
+	vy[1,1] = vy[1,1] + 0.00 * dvy
+	vy[1,2] = vy[1,1] + 0.20 * dvy
+	vy[2,1] = vy[1,1] + 0.25 * dvy
+	vy[2,2] = vy[1,1] + 0.75 * dvy
+	vy[3,1] = vy[1,1] + 0.80 * dvy
+	vy[3,2] = vy[1,1] + 1.00 * dvy
+
+	# (X,R)
+	call gseti (gp, G_WCS, 2)
+	call gseti (gp, G_DRAWAXES, 3)
+	call gseti (gp, G_XLABELTICKS, YES)
+	call gseti (gp, G_YLABELTICKS, YES)
+	call gseti (gp, G_XNMAJOR, 6)
+	call gseti (gp, G_XNMINOR, 4)
+	call gseti (gp, G_YNMAJOR, 4)
+	call gseti (gp, G_YNMINOR, 0)
+	call gsview (gp, vx[2,1], vx[2,2], vy[1,1], vy[1,2])
+	call gswind (gp, 0., 1., 0., 1.)
+	call gfill (gp, fa, fa[5], 4, GF_SOLID)
+
+	call gswind (gp, x1, x2, rmin, rmax)
+	call glabax (gp, "", "Column", "")
+
+	do i = 1, SFF_NSFD(sff) {
+	    sfd = SFF_SFD(sff,i)
+	    if (SFD_STATUS(sfd) != 0)
+		next
+	    x = SFD_X(sfd)
+	    y = SFD_W(sfd)
+	    if (key == 1) {
+		z = sqrt (SFS_M(SFD_SFS(sfd)) / SF_M(sf))
+		z = max (0.005, 0.03 * z)
+		call gmark (gp, x, y, GM_MAG, z, z)
+	    }
+	    if (SFD_W(sfd) < SF_W(sf))
+		call gseti (gp, G_PLCOLOR, 2)
+	    else
+		call gseti (gp, G_PLCOLOR, 3)
+	    z = min (2., SFD_W(sfd) / rbest)
+	    z = 0.010 * (1 + (z - 1) * 5)
+	    call gmark (gp, x, y, GM_CIRCLE, z, z)
+	    call gseti (gp, G_PLCOLOR, 1)
+	}
+
+	call gseti (gp, G_PLTYPE, 2)
+	call gline (gp, x1, SF_W(sf), x2, SF_W(sf))
+	call gseti (gp, G_PLTYPE, 1)
+
+	# (R,Y)
+	call gseti (gp, G_WCS, 3)
+	call gseti (gp, G_XLABELTICKS, YES)
+	call gseti (gp, G_YLABELTICKS, YES)
+	call gseti (gp, G_XNMAJOR, 4)
+	call gseti (gp, G_XNMINOR, 0)
+	call gseti (gp, G_YNMAJOR, 6)
+	call gseti (gp, G_YNMINOR, 4)
+	call gsview (gp, vx[1,1], vx[1,2], vy[2,1], vy[2,2])
+	call gswind (gp, 0., 1., 0., 1.)
+	call gfill (gp, fa, fa[5], 4, GF_SOLID)
+
+	call gswind (gp, rmin, rmax, y1, y2)
+	call glabax (gp, "", SF_WTYPE(sf), "Line")
+
+	do i = 1, SFF_NSFD(sff) {
+	    sfd = SFF_SFD(sff,i)
+	    if (SFD_STATUS(sfd) != 0)
+		next
+	    x = SFD_W(sfd)
+	    y = SFD_Y(sfd)
+	    if (key == 1) {
+		z = sqrt (SFS_M(SFD_SFS(sfd)) / SF_M(sf))
+		z = max (0.005, 0.03 * z)
+		call gmark (gp, x, y, GM_MAG, z, z)
+	    }
+	    if (SFD_W(sfd) < SF_W(sf))
+		call gseti (gp, G_PLCOLOR, 2)
+	    else
+		call gseti (gp, G_PLCOLOR, 3)
+	    z = min (2., SFD_W(sfd) / rbest)
+	    z = 0.010 * (1 + (z - 1) * 5)
+	    call gmark (gp, x, y, GM_CIRCLE, z, z)
+	    call gseti (gp, G_PLCOLOR, 1)
+	}
+
+	call gseti (gp, G_PLTYPE, 2)
+	call gline (gp, SF_W(sf), y1, SF_W(sf), y2)
+	call gseti (gp, G_PLTYPE, 1)
+
+	# (E,R)
+	call gseti (gp, G_WCS, 4)
+	call gseti (gp, G_DRAWAXES, 3)
+	call gseti (gp, G_XLABELTICKS, NO)
+	call gseti (gp, G_YLABELTICKS, YES)
+	call gseti (gp, G_XNMAJOR, 6)
+	call gseti (gp, G_XNMINOR, 4)
+	call gseti (gp, G_YNMAJOR, 4)
+	call gseti (gp, G_YNMINOR, 0)
+	call gsview (gp, vx[2,1], vx[2,2], vy[3,1], vy[3,2])
+	call gswind (gp, 0., 1., 0., 1.)
+	call gfill (gp, fa, fa[5], 4, GF_SOLID)
+
+	call gswind (gp, x1, x2, emin, emax)
+	call glabax (gp, "", "", "Ellip")
+
+	do i = 1, SFF_NSFD(sff) {
+	    sfd = SFF_SFD(sff,i)
+	    if (SFD_STATUS(sfd) != 0)
+		next
+	    x = SFD_X(sfd)
+	    y = SFD_E(sfd)
+	    if (key == 1) {
+		z = sqrt (SFS_M(SFD_SFS(sfd)) / SF_M(sf))
+		z = max (0.005, 0.03 * z)
+		call gmark (gp, x, y, GM_MAG, z, z)
+	    }
+	    if (SFD_W(sfd) < SF_W(sf))
+		call gseti (gp, G_PLCOLOR, 2)
+	    else
+		call gseti (gp, G_PLCOLOR, 3)
+	    z = min (2., SFD_W(sfd) / rbest)
+	    z = 0.010 * (1 + (z - 1) * 5)
+	    call gmark (gp, x, y, GM_CIRCLE, z, z)
+	    call gseti (gp, G_PLCOLOR, 1)
+	}
+
+	# (E,Y)
+	call gseti (gp, G_WCS, 5)
+	call gseti (gp, G_XLABELTICKS, YES)
+	call gseti (gp, G_YLABELTICKS, NO)
+	call gseti (gp, G_XNMAJOR, 4)
+	call gseti (gp, G_XNMINOR, 0)
+	call gseti (gp, G_YNMAJOR, 6)
+	call gseti (gp, G_YNMINOR, 4)
+	call gsview (gp, vx[3,1], vx[3,2], vy[2,1], vy[2,2])
+	call gswind (gp, 0., 1., 0., 1.)
+	call gfill (gp, fa, fa[5], 4, GF_SOLID)
+
+	call gswind (gp, emin, emax, y1, y2)
+	call glabax (gp, "", "Ellip", "")
+
+	do i = 1, SFF_NSFD(sff) {
+	    sfd = SFF_SFD(sff,i)
+	    if (SFD_STATUS(sfd) != 0)
+		next
+	    x = SFD_E(sfd)
+	    y = SFD_Y(sfd)
+	    if (key == 1) {
+		z = sqrt (SFS_M(SFD_SFS(sfd)) / SF_M(sf))
+		z = max (0.005, 0.03 * z)
+		call gmark (gp, x, y, GM_MAG, z, z)
+	    }
+	    if (SFD_W(sfd) < SF_W(sf))
+		call gseti (gp, G_PLCOLOR, 2)
+	    else
+		call gseti (gp, G_PLCOLOR, 3)
+	    z = min (2., SFD_W(sfd) / rbest)
+	    z = 0.010 * (1 + (z - 1) * 5)
+	    call gmark (gp, x, y, GM_CIRCLE, z, z)
+	    call gseti (gp, G_PLCOLOR, 1)
+	}
+
+	# Label window.
+	call gseti (gp, G_WCS, 1)
+	call gseti (gp, G_DRAWAXES, 0)
+	call gsview (gp, vx[1,1], vx[3,2], vy[1,1], vy[3,2])
+	call glabax (gp, title, "", "")
+
+	# (X,Y)
+	call gseti (gp, G_DRAWAXES, 3)
+	call gseti (gp, G_LABELTICKS, NO)
+	call gseti (gp, G_XNMAJOR, 6)
+	call gseti (gp, G_XNMINOR, 4)
+	call gseti (gp, G_YNMAJOR, 6)
+	call gseti (gp, G_YNMINOR, 4)
+	call gsview (gp, vx[2,1], vx[2,2], vy[2,1], vy[2,2])
+	call gswind (gp, 0., 1., 0., 1.)
+	call gfill (gp, fa, fa[5], 4, GF_SOLID)
+
+	call gswind (gp, x1, x2, y1, y2)
+	call glabax (gp, "", "", "")
+
+	do i = 1, SFF_NSFD(sff) {
+	    sfd = SFF_SFD(sff,i)
+	    if (SFD_STATUS(sfd) != 0)
+		next
+	    x = SFD_X(sfd)
+	    y = SFD_Y(sfd)
+	    if (key == 1) {
+		z = sqrt (SFS_M(SFD_SFS(sfd)) / SF_M(sf))
+		z = max (0.005, 0.03 * z)
+		call gmark (gp, x, y, GM_MAG, z, z)
+	    }
+	    if (SFD_W(sfd) < SF_W(sf))
+		call gseti (gp, G_PLCOLOR, 2)
+	    else
+		call gseti (gp, G_PLCOLOR, 3)
+	    z = min (2., SFD_W(sfd) / rbest)
+	    z = 0.010 * (1 + (z - 1) * 5)
+	    call gmark (gp, x, y, GM_CIRCLE, z, z)
+	    call gseti (gp, G_PLCOLOR, 1)
+	}
+end
+
+
+# STF_G9 -- Spatial plots at best focus.
+
+procedure stf_g12 (sf, current, key, title)
+
+pointer	sf			#I Starfocus pointer
+pointer	current			#I Current sfd pointer
+int	key			#I Plot magnitude symbol?
+char	title[ARB]		#I Title
+
+int	i
+real	x, y, z, x1, x2, y1, y2, fmin, fmax, rbest, rmin, rmax
+real	vx[3,2], vy[3,2], dvx, dvy, fa[8]
+pointer	gp, sfs, sfd
+
+data	fa/0.,1.,1.,0.,0.,0.,1.,1./
+
+begin
+	gp = SF_GP(sf)
+
+	# Range of X, Y, R, F.
+	x1 = 1.
+	y1 = 1.
+	x2 = SF_NCOLS(sf)
+	y2 = SF_NLINES(sf)
+
+	rbest = SFD_W(SF_BEST(sf))
+	fmin = MAX_REAL
+	fmax = -MAX_REAL
+	rmin = SF_W(sf)
+	rmax = 1.5 * SF_W(sf)
+	do i = 1, SF_NSTARS(sf) {
+	    sfs = SF_SFS(sf,i)
+	    if (SFS_N(sfs) == 0)
+		next
+	    fmin = min (fmin, SFS_F(sfs))
+	    fmax = max (fmax, SFS_F(sfs))
+	    rmin = min (rmin, SFS_W(sfs))
+	    rmax = max (rmax, SFS_W(sfs))
+	}
+	z = fmax - fmin
+	fmin = fmin - 0.1 * z
+	fmax = fmax + 0.1 * z
+	z = rmax - rmin
+	rmin = rmin - 0.1 * z
+	rmax = rmax + 0.1 * z
+
+	# Set view ports
+	call ggview (gp, vx[1,1], vx[3,2], vy[1,1], vy[3,2])
+	dvx = vx[3,2] - vx[1,1]
+	dvy = vy[3,2] - vy[1,1]
+	vx[1,1] = vx[1,1] + 0.00 * dvx
+	vx[1,2] = vx[1,1] + 0.20 * dvx
+	vx[2,1] = vx[1,1] + 0.25 * dvx
+	if (SF_NF(sf) > 1) {
+	    vx[2,2] = vx[1,1] + 0.75 * dvx
+	    vx[3,1] = vx[1,1] + 0.80 * dvx
+	    vx[3,2] = vx[1,1] + 1.00 * dvx
+	} else {
+	    vx[2,2] = vx[1,1] + 1.00 * dvx
+	    vx[3,1] = vx[1,1] + 1.00 * dvx
+	    vx[3,2] = vx[1,1] + 1.00 * dvx
+	}
+	vy[1,1] = vy[1,1] + 0.00 * dvy
+	vy[1,2] = vy[1,1] + 0.20 * dvy
+	vy[2,1] = vy[1,1] + 0.25 * dvy
+	if (SF_NF(sf) > 1) {
+	    vy[2,2] = vy[1,1] + 0.75 * dvy
+	    vy[3,1] = vy[1,1] + 0.80 * dvy
+	    vy[3,2] = vy[1,1] + 1.00 * dvy
+	} else {
+	    vy[2,2] = vy[1,1] + 1.00 * dvy
+	    vy[3,1] = vy[1,1] + 1.00 * dvy
+	    vy[3,2] = vy[1,1] + 1.00 * dvy
+	}
+
+	dvx = vx[2,1] - vx[2,2]
+	dvy = vy[1,1] - vy[1,2]
+	if (abs (dvx) > 0.01 && abs (dvy) > 0.01) {
+	    # (X,R)
+	    call gseti (gp, G_WCS, 2)
+	    call gseti (gp, G_DRAWAXES, 3)
+	    call gseti (gp, G_XLABELTICKS, YES)
+	    call gseti (gp, G_YLABELTICKS, YES)
+	    call gseti (gp, G_XNMAJOR, 6)
+	    call gseti (gp, G_XNMINOR, 4)
+	    call gseti (gp, G_YNMAJOR, 4)
+	    call gseti (gp, G_YNMINOR, 0)
+	    call gsview (gp, vx[2,1], vx[2,2], vy[1,1], vy[1,2])
+	    call gswind (gp, 0., 1., 0., 1.)
+	    call gfill (gp, fa, fa[5], 4, GF_SOLID)
+
+	    call gswind (gp, x1, x2, rmin, rmax)
+	    call glabax (gp, "", "Column", "")
+
+	    do i = 1, SF_NSTARS(sf) {
+		sfs = SF_SFS(sf,i)
+		if (SFS_N(sfs) == 0)
+		    next
+		x = SFD_X(SFS_SFD(sfs,1))
+		y = SFS_W(sfs)
+		if (key == 1) {
+		    z = sqrt (SFS_M(sfs) / SF_M(sf))
+		    z = max (0.005, 0.03 * z)
+		    call gmark (gp, x, y, GM_MAG, z, z)
+		}
+		if (SFS_F(sfs) < SF_F(sf))
+		    call gseti (gp, G_PLCOLOR, 2)
+		else
+		    call gseti (gp, G_PLCOLOR, 3)
+		z = min (2., SFS_W(sfs) / rbest)
+		z = 0.010 * (1 + (z - 1) * 5)
+		call gmark (gp, x, y, GM_CIRCLE, z, z)
+		call gseti (gp, G_PLCOLOR, 1)
+	    }
+
+	    call gseti (gp, G_PLTYPE, 2)
+	    call gline (gp, x1, SF_W(sf), x2, SF_W(sf))
+	    call gseti (gp, G_PLTYPE, 1)
+	}
+
+	dvx = vx[1,1] - vx[1,2]
+	dvy = vy[2,1] - vy[2,2]
+	if (abs (dvx) > 0.01 && abs (dvy) > 0.01) {
+	    # (R,Y)
+	    call gseti (gp, G_WCS, 3)
+	    call gseti (gp, G_XLABELTICKS, YES)
+	    call gseti (gp, G_YLABELTICKS, YES)
+	    call gseti (gp, G_XNMAJOR, 4)
+	    call gseti (gp, G_XNMINOR, 0)
+	    call gseti (gp, G_YNMAJOR, 6)
+	    call gseti (gp, G_YNMINOR, 4)
+	    call gsview (gp, vx[1,1], vx[1,2], vy[2,1], vy[2,2])
+	    call gswind (gp, 0., 1., 0., 1.)
+	    call gfill (gp, fa, fa[5], 4, GF_SOLID)
+
+	    call gswind (gp, rmin, rmax, y1, y2)
+	    call glabax (gp, "", SF_WTYPE(sf), "Line")
+
+	    do i = 1, SF_NSTARS(sf) {
+		sfs = SF_SFS(sf,i)
+		if (SFS_N(sfs) == 0)
+		    next
+		x = SFS_W(sfs)
+		y = SFD_Y(SFS_SFD(sfs,1))
+		if (key == 1) {
+		    z = sqrt (SFS_M(sfs) / SF_M(sf))
+		    z = max (0.005, 0.03 * z)
+		    call gmark (gp, x, y, GM_MAG, z, z)
+		}
+		if (SFS_F(sfs) < SF_F(sf))
+		    call gseti (gp, G_PLCOLOR, 2)
+		else
+		    call gseti (gp, G_PLCOLOR, 3)
+		z = min (2., SFS_W(sfs) / rbest)
+		z = 0.010 * (1 + (z - 1) * 5)
+		call gmark (gp, x, y, GM_CIRCLE, z, z)
+		call gseti (gp, G_PLCOLOR, 1)
+	    }
+
+	    call gseti (gp, G_PLTYPE, 2)
+	    call gline (gp, SF_W(sf), y1, SF_W(sf), y2)
+	    call gseti (gp, G_PLTYPE, 1)
+	}
+
+	dvx = vx[2,1] - vx[2,2]
+	dvy = vy[3,1] - vy[3,2]
+	if (abs (dvx) > 0.01 && abs (dvy) > 0.01) {
+	    # (X,F)
+	    call gseti (gp, G_WCS, 4)
+	    call gseti (gp, G_XLABELTICKS, NO)
+	    call gseti (gp, G_YLABELTICKS, YES)
+	    call gseti (gp, G_XNMAJOR, 6)
+	    call gseti (gp, G_XNMINOR, 4)
+	    call gseti (gp, G_YNMAJOR, 4)
+	    call gseti (gp, G_YNMINOR, 0)
+	    call gsview (gp, vx[2,1], vx[2,2], vy[3,1], vy[3,2])
+	    call gswind (gp, 0., 1., 0., 1.)
+	    call gfill (gp, fa, fa[5], 4, GF_SOLID)
+
+	    call gswind (gp, x1, x2, fmin, fmax)
+	    call glabax (gp, "", "", "Focus")
+
+	    do i = 1, SF_NSTARS(sf) {
+		sfs = SF_SFS(sf,i)
+		if (SFS_N(sfs) == 0)
+		    next
+		x = SFD_X(SFS_SFD(sfs,1))
+		y = SFS_F(sfs)
+		if (key == 1) {
+		    z = sqrt (SFS_M(sfs) / SF_M(sf))
+		    z = max (0.005, 0.03 * z)
+		    call gmark (gp, x, y, GM_MAG, z, z)
+		}
+		if (SFS_F(sfs) < SF_F(sf))
+		    call gseti (gp, G_PLCOLOR, 2)
+		else
+		    call gseti (gp, G_PLCOLOR, 3)
+		z = min (2., SFS_W(sfs) / rbest)
+		z = 0.010 * (1 + (z - 1) * 5)
+		call gmark (gp, x, y, GM_CIRCLE, z, z)
+		call gseti (gp, G_PLCOLOR, 1)
+	    }
+
+	    call gseti (gp, G_PLTYPE, 2)
+	    call gline (gp, x1, SF_F(sf), x2, SF_F(sf))
+	    call gseti (gp, G_PLTYPE, 1)
+	}
+
+	dvx = vx[3,1] - vx[3,2]
+	dvy = vy[2,1] - vy[2,2]
+	if (abs (dvx) > 0.01 && abs (dvy) > 0.01) {
+	    # (F,Y)
+	    call gseti (gp, G_WCS, 5)
+	    call gseti (gp, G_XLABELTICKS, YES)
+	    call gseti (gp, G_YLABELTICKS, NO)
+	    call gseti (gp, G_XNMAJOR, 4)
+	    call gseti (gp, G_XNMINOR, 0)
+	    call gseti (gp, G_YNMAJOR, 6)
+	    call gseti (gp, G_YNMINOR, 4)
+	    call gsview (gp, vx[3,1], vx[3,2], vy[2,1], vy[2,2])
+	    call gswind (gp, 0., 1., 0., 1.)
+	    call gfill (gp, fa, fa[5], 4, GF_SOLID)
+
+	    call gswind (gp, fmin, fmax, y1, y2)
+	    call glabax (gp, "", "Focus", "")
+
+	    do i = 1, SF_NSTARS(sf) {
+		sfs = SF_SFS(sf,i)
+		if (SFS_N(sfs) == 0)
+		    next
+		x = SFS_F(sfs)
+		y = SFD_Y(SFS_SFD(sfs,1))
+		if (key == 1) {
+		    z = sqrt (SFS_M(sfs) / SF_M(sf))
+		    z = max (0.005, 0.03 * z)
+		    call gmark (gp, x, y, GM_MAG, z, z)
+		}
+		if (SFS_F(sfs) < SF_F(sf))
+		    call gseti (gp, G_PLCOLOR, 2)
+		else
+		    call gseti (gp, G_PLCOLOR, 3)
+		z = min (2., SFS_W(sfs) / rbest)
+		z = 0.010 * (1 + (z - 1) * 5)
+		call gmark (gp, x, y, GM_CIRCLE, z, z)
+		call gseti (gp, G_PLCOLOR, 1)
+	    }
+
+	    call gseti (gp, G_PLTYPE, 2)
+	    call gline (gp, SF_F(sf), y1, SF_F(sf), y2)
+	    call gseti (gp, G_PLTYPE, 1)
+	}
+
+	# Label window.
+	call gseti (gp, G_WCS, 1)
+	call gseti (gp, G_DRAWAXES, 0)
+	call gsview (gp, vx[1,1], vx[3,2], vy[1,1], vy[3,2])
+	call glabax (gp, title, "", "")
+
+	dvx = vx[2,1] - vx[2,2]
+	dvy = vy[2,1] - vy[2,2]
+	if (abs (dvx) > 0.01 && abs (dvy) > 0.01) {
+	    # (X,Y)
+	    call gseti (gp, G_DRAWAXES, 3)
+	    call gseti (gp, G_LABELTICKS, NO)
+	    call gseti (gp, G_XNMAJOR, 6)
+	    call gseti (gp, G_XNMINOR, 4)
+	    call gseti (gp, G_YNMAJOR, 6)
+	    call gseti (gp, G_YNMINOR, 4)
+	    call gsview (gp, vx[2,1], vx[2,2], vy[2,1], vy[2,2])
+	    call gswind (gp, 0., 1., 0., 1.)
+	    call gfill (gp, fa, fa[5], 4, GF_SOLID)
+
+	    call gswind (gp, x1, x2, y1, y2)
+	    call glabax (gp, "", "", "")
+
+	    do i = 1, SF_NSTARS(sf) {
+		sfs = SF_SFS(sf,i)
+		if (SFS_N(sfs) == 0)
+		    next
+		sfd = SFS_SFD(sfs,1)
+		x = SFD_X(sfd)
+		y = SFD_Y(sfd)
+		if (key == 1) {
+		    z = sqrt (SFS_M(sfs) / SF_M(sf))
+		    z = max (0.005, 0.03 * z)
+		    call gmark (gp, x, y, GM_MAG, z, z)
+		}
+		if (SFS_F(sfs) < SF_F(sf))
+		    call gseti (gp, G_PLCOLOR, 2)
+		else
+		    call gseti (gp, G_PLCOLOR, 3)
+		z = min (2., SFS_W(sfs) / rbest)
+		z = 0.010 * (1 + (z - 1) * 5)
+		call gmark (gp, x, y, GM_CIRCLE, z, z)
+		call gseti (gp, G_PLCOLOR, 1)
+	    }
+	}
+end
diff --git a/noao/obsutil/src/starfocus/stfhelp.key b/noao/obsutil/src/starfocus/stfhelp.key
new file mode 100644
index 00000000..948bb9b7
--- /dev/null
+++ b/noao/obsutil/src/starfocus/stfhelp.key
@@ -0,0 +1,63 @@
+			STARFOCUS COMMAND OPTIONS
+
+				SUMMARY
+
+? Help              f Focus             p Radial profile    u Undelete
+a Spatial           i Info              q Quit              x Delete        
+b Best              m Magnitude         r Redraw            z Zoom
+d Delete            n Normalize         s Mag symbols       <space> Next
+e Enclosed flux     o Offset            t Field radius      
+
+:level <val>        :radius <val>       :show <file>        :xcenter <val>
+:overplot <y|n>     :scale <val>        :size <type>        :ycenter <val>
+
+
+				CURSOR COMMANDS
+
+All plots may not be available depending on the number of focus values and
+the number of stars.
+
+?	Page this help information
+a	Spatial plot at a single focus
+b	Spatial plot of best focus values
+d	Delete star nearest to cursor
+e	Enclosed flux for all stars at one focus and all focus for one star
+f	Size and ellipticity vs focus for all data
+i	Information about point nearest the cursor
+m	Size and ellipticity vs relative magnitude at one focus
+n	Normalize enclosed flux at x cursor position
+o	Offset enclosed flux to x,y cursor position by adjusting background
+p	Radial profiles for all stars at one focus and all focus for one star
+          The profiles are determined from the derivatives of the enclosed flux
+q	Quit
+r	Redraw
+s	Toggle magnitude symbols in spatial plots
+t	Size and ellipticity vs radius from field center at one focus
+u	Undelete all deleted points
+x	Delete nearest point, star, or focus (selected by query)
+z	Zoom to a single measurement showing enclosed flux and radial profile
+<space> Step through different focus or stars in current plot type
+
+
+			COLON COMMANDS
+
+A command without a value generally shows the current value of the
+parameter while with a value it sets the value of the parameter.
+
+:level <val>	Level at which the size parameter is evaluated
+:overplot <y|n> Overplot the profiles from the narrowest profile?
+:radius <val>   Change profile radius(*)
+:show <file>	Page all information for the current set of objects
+:size <type>	Size type (Radius|FWHM|GFWHM|MFWHM) (**)
+:scale <val>	Pixel scale for size values
+:xcenter <val>	X field center for radius from field center plots
+:ycenter <val>	Y field center for radius from field center plots
+
+(*) The profile radius may not exceed the initial value set by the task
+    parameter.
+
+(**)
+Radius = radius enclosing the fraction of the flux specified by "level"
+FWHM   = Full-width at half-maximum based on the radially smoothed profile
+GFWHM  = Full-width at half-maximum of Gaussian function fit to enclosed flux
+MFWHM  = Full-width at half-maximum of Moffat function fit to enclosed flux
diff --git a/noao/obsutil/src/starfocus/stfmeasure.x b/noao/obsutil/src/starfocus/stfmeasure.x
new file mode 100644
index 00000000..808b6058
--- /dev/null
+++ b/noao/obsutil/src/starfocus/stfmeasure.x
@@ -0,0 +1,134 @@
+include	<error.h>
+include	<imhdr.h>
+include	<imset.h>
+include	<math/nlfit.h>
+include	"starfocus.h"
+
+
+# STF_MEASURE -- PSF measuring routine.
+# This is a stand-alone routine that can be called to return the FWHM.
+# It is a greatly abbreviated version of starfocus.
+
+procedure stf_measure (im, xc, yc, beta, level, radius, nit,
+	sbuffer, swidth, saturation, gp, logfd,
+	bkg, renclosed, dfwhm, gfwhm, mfwhm)
+
+pointer	im		#I Image pointer
+real	xc		#I Initial X center
+real	yc		#I Initial Y center
+real	beta		#I Moffat beta
+real	level		#I Measurement level
+real	radius		#U Profile radius
+int	nit		#I Number of iterations on radius
+real	sbuffer		#I Sky buffer (pixels)
+real	swidth		#I Sky width (pixels)
+real	saturation	#I Saturation
+pointer	gp		#I Graphics output if not NULL
+int	logfd		#I Log output if not NULL
+real	bkg		#O Background used
+real	renclosed	#O Enclosed flux radius
+real	dfwhm		#O Direct FWHM
+real	gfwhm		#O Gaussian FWHM
+real	mfwhm		#O Moffat FWHM
+
+int	i
+bool	ignore_sat
+pointer	sp, str, sf, sfd, sfds
+
+int	strdic()
+real	stf_r2i()
+errchk	stf_find, stf_bkgd, stf_profile, stf_widths, stf_fwhms, stf_organize
+
+begin
+	call smark (sp)
+	call salloc (str, SZ_FNAME, TY_CHAR)
+	call salloc (sf, SF, TY_STRUCT)
+	call salloc (sfd, SFD, TY_STRUCT)
+	call salloc (sfds, 1, TY_POINTER)
+	call aclri (Memi[sf], SF)
+	call aclri (Memi[sfd], SFD)
+	Memi[sfds] = sfd
+
+	# Initialize parameters.
+	SF_TASK(sf) = PSFMEASURE
+	SF_WCODE(sf) = strdic ("FWHM", SF_WTYPE(sf), SF_SZWTYPE, SF_WTYPES)
+	SF_SCALE(sf) = 1.
+	SF_LEVEL(sf) = level
+	SF_BETA(sf) = beta
+	SF_RADIUS(sf) = radius
+	SF_SBUF(sf) = sbuffer
+	SF_SWIDTH(sf) = swidth
+	SF_SAT(sf) = saturation
+	SF_NIT(sf) = nit
+	SF_OVRPLT(sf) = NO
+	SF_NCOLS(sf) = IM_LEN(im,1)
+	SF_NLINES(sf) = IM_LEN(im,2)
+	SF_XF(sf) = (IM_LEN(im,1) + 1) / 2.
+	SF_YF(sf) = (IM_LEN(im,2) + 1) / 2.
+	ignore_sat = false
+
+	call imstats (im, IM_IMAGENAME, SFD_IMAGE(sfd), SF_SZFNAME)
+	SFD_ID(sfd) = 1
+	SFD_X(sfd) = xc
+	SFD_Y(sfd) = yc
+	SFD_F(sfd) = INDEF
+	SFD_STATUS(sfd) = 0
+	SFD_SFS(sfd) = NULL
+	SFD_SFF(sfd) = NULL
+	SFD_SFI(sfd) = NULL
+
+	if (SF_LEVEL(sf) > 1.)
+	    SF_LEVEL(sf) = SF_LEVEL(sf) / 100.
+	SF_LEVEL(sf) = max (0.05, min (0.95, SF_LEVEL(sf)))
+
+	# Evaluate PSF data.
+	iferr {
+	    do i = 1, SF_NIT(sf) {
+		if (i == 1)
+		    SFD_RADIUS(sfd) = SF_RADIUS(sf)
+		else
+		    SFD_RADIUS(sfd) = 3. * SFD_DFWHM(sfd)
+		SFD_NPMAX(sfd) = stf_r2i (SFD_RADIUS(sfd)) + 1
+		SFD_NP(sfd) = SFD_NPMAX(sfd)
+		call stf_find (sf, sfd, im)
+		call stf_bkgd (sf, sfd)
+		if (SFD_NSAT(sfd) > 0 && i == 1) {
+		    if (ignore_sat)
+			call error (0,
+			"Saturated pixels found - ignoring object")
+		    else
+			call eprintf (
+			    "WARNING: Saturated pixels found.\n")
+		}
+		call stf_profile (sf, sfd)
+		call stf_widths (sf, sfd)
+		call stf_fwhms (sf, sfd)
+	    }
+
+	    # Set output results.
+	    radius = SFD_RADIUS(sfd)
+	    bkg = SFD_BKGD(sfd)
+	    renclosed = SFD_R(sfd)
+	    dfwhm = SFD_DFWHM(sfd)
+	    mfwhm = SFD_MFWHM(sfd)
+	    gfwhm = SFD_GFWHM(sfd)
+
+	    # Optional graph and log output.  Note that the gp pointer is only
+	    # used to indicate whether to make a graph.  The stf_graph
+	    # procedure opens its own graphics stream.
+
+	    call stf_organize (sf, sfds, 1)
+	    if (gp != NULL)
+		call stf_graph (sf)
+	    if (logfd != NULL)
+		call stf_log (sf, logfd)
+
+	    call asifree (SFD_ASI1(sfd))
+	    call asifree (SFD_ASI2(sfd))
+	} then
+	    call erract (EA_WARN)
+
+	# Finish up
+	call stf_free (sf)
+	call sfree (sp)
+end
diff --git a/noao/obsutil/src/starfocus/stfprofile.x b/noao/obsutil/src/starfocus/stfprofile.x
new file mode 100644
index 00000000..d26c085d
--- /dev/null
+++ b/noao/obsutil/src/starfocus/stfprofile.x
@@ -0,0 +1,1189 @@
+include	<imhdr.h>
+include	<mach.h>
+include	<math.h>
+include	<math/iminterp.h>
+include	<math/nlfit.h>
+include	"starfocus.h"
+
+
+# STF_FIND      -- Find the object and return the data raster and object center.
+# STF_BKGD      -- Compute the background.
+# STF_PROFILE   -- Compute enclosed flux profile, derivative, and moments.
+# STF_NORM      -- Renormalized enclosed flux profile
+# STF_WIDTHS    -- Set widths.
+# STF_I2R       -- Radius from sample index.
+# STF_R2I       -- Sample index from radius.
+# STF_R2N       -- Number of subsamples from radius.
+# STF_MODEL     -- Return model values.
+# STF_DFWHM     -- Direct FWHM from profile.
+# STF_FWHMS     -- Measure FWHM vs level.
+# STF_RADIUS    -- Measure the radius at the specified level.
+# STF_FIT       -- Fit model.
+# STF_GAUSS1    -- Gaussian function used in NLFIT.
+# STF_GAUSS2    -- Gaussian function and derivatives used in NLFIT.
+# STF_MOFFAT1   -- Moffat function used in NLFIT.
+# STF_MOFFAT2   -- Moffat function and derivatives used in NLFIT.
+
+
+# STF_FIND -- Find the object and return the data raster and object center.
+# Centering uses centroid of marginal distributions of data above the mean.
+
+procedure stf_find (sf, sfd, im)
+
+pointer	sf			#I Starfocus pointer
+pointer	sfd			#I Object pointer
+pointer	im			#I Image pointer
+
+long	lseed
+int	i, j, k, x1, x2, y1, y2, nx, ny, npts
+real	radius, buffer, width, xc, yc, xlast, ylast, r1, r2
+real	mean, sum, sum1, sum2, sum3, asumr(), urand()
+pointer	data, ptr, imgs2r()
+errchk	imgs2r
+
+begin
+	radius = max (3., SFD_RADIUS(sfd))
+	buffer = SF_SBUF(sf)
+	width = SF_SWIDTH(sf)
+
+	xc = SFD_X(sfd)
+	yc = SFD_Y(sfd)
+	r1 = radius + buffer + width
+	r2 = radius
+
+	# Iterate on the center finding.
+	do k = 1, 3 {
+
+	    # Extract region around current center.
+	    xlast = xc
+	    ylast = yc
+
+	    x1 = max (1-NBNDRYPIX, nint (xc - r2))
+	    x2 = min (IM_LEN(im,1)+NBNDRYPIX, nint (xc + r2))
+	    nx = x2 - x1 + 1
+	    y1 = max (1-NBNDRYPIX, nint (yc - r2))
+	    y2 = min (IM_LEN(im,2)+NBNDRYPIX, nint (yc + r2))
+	    ny = y2 - y1 + 1
+	    npts = nx * ny
+	    data = imgs2r (im, x1, x2, y1, y2)
+
+	    # Find center of gravity of marginal distributions above mean.
+	    npts = nx * ny
+	    sum = asumr (Memr[data], npts)
+	    mean = sum / nx
+	    sum1 = 0.
+	    sum2 = 0.
+
+	    do i = x1, x2 {
+		ptr = data + i - x1
+		sum3 = 0.
+		do j = y1, y2 {
+		    sum3 = sum3 + Memr[ptr]
+		    ptr = ptr + nx
+		}
+		sum3 = sum3 - mean
+		if (sum3 > 0.) {
+		    sum1 = sum1 + i * sum3
+		    sum2 = sum2 + sum3
+		}
+	    }
+	    if (sum2 <= 0)
+		call error (1, "Centering failed to converge")
+	    xc = sum1 / sum2
+	    if (xlast - xc > 0.2 * nx)
+		xc = xlast - 0.2 * nx
+	    if (xc - xlast > 0.2 * nx)
+		xc = xlast + 0.2 * nx
+
+	    ptr = data
+	    mean = sum / ny
+	    sum1 = 0.
+	    sum2 = 0.
+	    do j = y1, y2 {
+		sum3 = 0.
+		do i = x1, x2 {
+		    sum3 = sum3 + Memr[ptr]
+		    ptr = ptr + 1
+		}
+		sum3 = sum3 - mean
+		if (sum3 > 0.) {
+		    sum1 = sum1 + j * sum3
+		    sum2 = sum2 + sum3
+		}
+	    }
+	    if (sum2 <= 0)
+		call error (1, "Centering failed to converge")
+	    yc = sum1 / sum2
+	    if (ylast - yc > 0.2 * ny)
+		yc = ylast - 0.2 * ny
+	    if (yc - ylast > 0.2 * ny)
+		yc = ylast + 0.2 * ny
+
+	    if (nint(xc) == nint(xlast) && nint(yc) == nint(ylast))
+		break
+	}
+
+	# Get a new centered raster if necessary.
+	if (nint(xc) != nint(xlast) || nint(yc) != nint(ylast) || r2 < r1) {
+	    x1 = max (1-NBNDRYPIX, nint (xc - r1))
+	    x2 = min (IM_LEN(im,1)+NBNDRYPIX, nint (xc + r1))
+	    nx = x2 - x1 + 1
+	    y1 = max (1-NBNDRYPIX, nint (yc - r1))
+	    y2 = min (IM_LEN(im,2)+NBNDRYPIX, nint (yc + r1))
+	    ny = y2 - y1 + 1
+	    npts = nx * ny
+	    data = imgs2r (im, x1, x2, y1, y2)
+	}
+
+	# Add a dither for integer data.  The random numbers are always
+	# the same to provide reproducibility.
+
+	i = IM_PIXTYPE(im)
+	if (i == TY_SHORT || i == TY_INT || i == TY_LONG) {
+	    lseed = 1
+	    do i = 0, npts-1
+		Memr[data+i] = Memr[data+i] + urand(lseed) - 0.5
+	}
+
+	SFD_DATA(sfd) = data
+	SFD_X1(sfd) = x1
+	SFD_X2(sfd) = x2
+	SFD_Y1(sfd) = y1
+	SFD_Y2(sfd) = y2
+	SFD_X(sfd) = xc
+	SFD_Y(sfd) = yc
+end
+
+
+# STF_BKGD -- Compute the background.
+# A mode is estimated from the minimum slope in the sorted background pixels
+# with a bin width of 5%.
+
+procedure stf_bkgd (sf, sfd)
+
+pointer	sf			#I Parameter structure
+pointer	sfd			#I Star structure
+
+int	i, j, x1, x2, y1, y2, xc, yc, nx, ny, npts, ns, nsat
+real	sat, bkgd, miso
+real	r, r1, r2, r3, dx, dy, dz
+pointer	sp, data, bdata, ptr
+
+begin
+	data = SFD_DATA(sfd)
+	x1 = SFD_X1(sfd)
+	x2 = SFD_X2(sfd)
+	y1 = SFD_Y1(sfd)
+	y2 = SFD_Y2(sfd)
+	xc = SFD_X(sfd)
+	yc = SFD_Y(sfd)
+
+	nx = x2 - x1 + 1
+	ny = y2 - y1 + 1
+	npts = nx * ny
+
+	ns = 0
+	nsat = 0
+	r1 = SFD_RADIUS(sfd) ** 2
+	r2 = (SFD_RADIUS(sfd) + SF_SBUF(sf)) ** 2
+	r3 = (SFD_RADIUS(sfd) + SF_SBUF(sf) + SF_SWIDTH(sf)) ** 2
+	sat = SF_SAT(sf)
+	if (IS_INDEF(sat))
+	    sat = MAX_REAL
+
+	call smark (sp)
+	call salloc (bdata, npts, TY_REAL)
+
+	ptr = data
+	do j = y1, y2 {
+	    dy = (yc - j) ** 2
+	    do i = x1, x2 {
+		dx = (xc - i) ** 2
+		r = dx + dy
+		if (r <= r1) {
+		    if (Memr[ptr] >= sat)
+			nsat = nsat + 1
+		} else if (r >= r2 && r <= r3) {
+		    Memr[bdata+ns] = Memr[ptr]
+		    ns = ns + 1
+		}
+		ptr = ptr + 1
+	    }
+	}
+
+	if (ns > 9) {
+	    call asrtr (Memr[bdata], Memr[bdata], ns)
+	    r = Memr[bdata+ns-1] - Memr[bdata]
+	    bkgd = Memr[bdata] + r / 2
+	    miso = r / 2
+
+	    j = 1 + 0.50 * ns
+	    do i = 0, ns - j {
+		dz = Memr[bdata+i+j-1] - Memr[bdata+i]
+		if (dz < r) {
+		    r = dz
+		    bkgd = Memr[bdata+i] + dz / 2
+		    miso = dz / 2
+		}
+	    }
+	} else {
+	    bkgd = 0.
+	    miso = 0.
+	}
+
+	SFD_BKGD1(sfd) = bkgd
+	SFD_BKGD(sfd) = bkgd
+	SFD_MISO(sfd) = miso
+	SFD_NSAT(sfd) = nsat
+
+	call sfree (sp)
+end
+
+
+# STF_PROFILE -- Compute enclosed flux profile, derivative, direct FWHM, and
+# profile moments..
+# 1.  The flux profile is normalized at the maximum value.
+# 2.  The radial profile is computed from the numerical derivative of the
+#     enclose flux profile.
+
+procedure stf_profile (sf, sfd)
+
+pointer	sf			#I Parameter structure
+pointer	sfd			#I Star structure
+
+int	np
+real	radius, xc, yc
+
+int	i, j, k, l, m, ns, nx, ny, x1, x2, y1, y2
+real	bkgd, miso, sigma, peak
+real	r, r1, r2, r3, dx, dy, dx1, dx2, dy1, dy2, dz, xx, yy, xy, ds, da
+pointer	sp, data, profile, ptr, asi, msi, gs
+int	stf_r2n()
+real	asieval(), msieval(), gseval(), stf_i2r(), stf_r2i()
+errchk	asiinit, asifit, msiinit, msifit, gsrestore
+
+real	gsdata[24]
+data	gsdata/ 1., 4., 4., 1., 0., 0.6726812, 1., 2., 1.630641, 0.088787,
+		0.00389378, -0.001457133, 0.3932125, -0.1267456, -0.004864541,
+		0.00249941, 0.03078612, 0.02731274, -4.875850E-4, 2.307464E-4,
+		-0.002134843, 0.007603908, -0.002552385, -8.010564E-4/
+
+begin
+	data = SFD_DATA(sfd)
+	x1 = SFD_X1(sfd)
+	x2 = SFD_X2(sfd)
+	y1 = SFD_Y1(sfd)
+	y2 = SFD_Y2(sfd)
+	xc = SFD_X(sfd)
+	yc = SFD_Y(sfd)
+	bkgd = SFD_BKGD(sfd)
+	miso = SFD_MISO(sfd)
+	radius = SFD_RADIUS(sfd)
+	np = SFD_NP(sfd)
+
+	nx = x2 - x1 + 1
+	ny = y2 - y1 + 1
+
+	# Use an image interpolator fit to the data.
+	call msiinit (msi, II_BISPLINE3)
+	call msifit (msi, Memr[data], nx, ny, nx)
+
+	# To avoid trying to interpolate outside the center of the
+	# edge pixels, a requirement of the interpolator functions,
+	# we reset the data limits.
+	x1 = x1 + 1
+	x2 = x2 - 1
+	y1 = y1 + 1
+	y2 = y2 - 1
+
+	# Compute the enclosed flux profile, its derivative, and moments.
+	call smark (sp)
+	call salloc (profile, np, TY_REAL)
+	call aclrr (Memr[profile], np)
+
+	xx = 0.
+	yy = 0.
+	xy = 0.
+	do j = y1, y2 {
+	    ptr = data + (j-y1+1)*nx + 1
+	    dy = j - yc
+	    do i = x1, x2 {
+		dx = i - xc
+
+		# Set the subpixel sampling which may be a function of radius.
+		r = sqrt (dx * dx + dy * dy)
+		ns = stf_r2n (r)
+		ds = 1. / ns
+		da = ds * ds
+		dz = 0.5 + 0.5 * ds
+
+		# Sum the interpolator values over the subpixels and compute
+		# an offset to give the correct total for the pixel.
+
+		r2 = 0.
+		dy1 = dy - dz
+		do l = 1, ns {
+		    dy1 = dy1 + ds
+		    dy2 = dy1 * dy1
+		    dx1 = dx - dz
+		    do k = 1, ns {
+			dx1 = dx1 + ds
+			dx2 = dx1 * dx1
+			r1 = msieval (msi, dx1+xc-x1+2, dy1+yc-y1+2)
+			r2 = r2 + r1
+		    }
+		}
+
+		r1 = Memr[ptr] - bkgd
+		ptr = ptr + 1
+		r2 = r1 - r2 * da
+
+		# Accumulate the enclosed flux over the sub pixels.
+		dy1 = dy - dz
+		do l = 1, ns {
+		    dy1 = dy1 + ds
+		    dy2 = dy1 * dy1
+		    dx1 = dx - dz
+		    do k = 1, ns {
+			dx1 = dx1 + ds
+			dx2 = dx1 * dx1
+			r = max (0., sqrt (dx2 + dy2) - ds / 2)
+			if (r < radius) {
+			    r1 = da * (msieval (msi, dx1+xc-x1+2, dy1+yc-y1+2) +
+				r2)
+
+			    # Use approximation for fractions of a subpixel.
+			    for (m=stf_r2i(r)+1; m<=np; m=m+1) {
+				r3 = (stf_i2r (real(m)) - r) / ds
+				if (r3 >= 1.)
+				    break
+				Memr[profile+m-1] = Memr[profile+m-1] + r3 * r1
+			    }
+
+			    # The subpixel is completely within these radii.
+			    for (; m<=np; m=m+1)
+				Memr[profile+m-1] = Memr[profile+m-1] + r1
+
+			    # Accumulate the moments above an isophote.
+			    if (r1 > miso) {
+				xx = xx + dx2 * r1
+				yy = yy + dy2 * r1
+				xy = xy + dx1 * dy1 * r1
+			    }
+			}
+		    }
+		}
+	    }
+	}
+
+	call msifree (msi)
+
+	# Compute the ellipticity and position angle from the moments.
+	r = (xx + yy)
+	if (r > 0.) {
+	    r1 = (xx - yy) / r
+	    r2 = 2 * xy / r
+	    SFD_E(sfd) = sqrt (r1**2 + r2**2)
+	    SFD_PA(sfd) = RADTODEG (atan2 (r2, r1) / 2.)
+	} else {
+	    SFD_E(sfd) = 0.
+	    SFD_PA(sfd) = 0.
+	}
+
+	# The magnitude and profile normalization is from the max enclosed flux.
+	call alimr (Memr[profile], np, r, SFD_M(sfd))
+	if (SFD_M(sfd) <= 0.)
+	    call error (1, "Invalid flux profile")
+	call adivkr (Memr[profile], SFD_M(sfd), Memr[profile], np)
+
+	# Fit interpolator to the enclosed flux profile.
+	call asiinit (asi, II_SPLINE3)
+	call asifit (asi, Memr[profile], np)
+	SFD_ASI1(sfd) = asi
+
+	# Estimate a gaussian sigma (actually sqrt(2)*sigma) and if it is
+	# it is small subtract the gaussian so that the image interpolator
+	# can more accurately estimate subpixel values.
+
+	#call stf_radius (sf, sfd, SF_LEVEL(sf), r)
+	#sigma = r / sqrt (log (1/(1-SF_LEVEL(sf))))
+	call stf_radius (sf, sfd, 0.8, r)
+	r = r / SF_SCALE(sf)
+	sigma = 2 * r * sqrt (log(2.) / log (1/(1-0.8)))
+	if (sigma < 5.) {
+	    if (sigma <= 2.) {
+		call gsrestore (gs, gsdata)
+		dx = xc - nint (xc)
+		dy = yc - nint (yc)
+		r = sqrt (dx * dx + dy * dy)
+		dx = 1.
+		ds = abs (sigma - gseval (gs, r, dx))
+		for (da = 1.; da <= 2.; da = da + .01) {
+		    dz = abs (sigma - gseval (gs, r, da))
+		    if (dz < ds) {
+			ds = dz
+			dx = da
+		    }
+		}
+		sigma = dx
+		call gsfree (gs)
+	    }
+
+	    sigma = sigma / (2 * sqrt (log(2.)))
+	    sigma = sigma * sigma
+
+	    # Compute the peak that gives the correct central pixel value.
+	    i = nint (xc)
+	    j = nint (yc)
+	    dx = i - xc
+	    dy = j - yc
+	    r = sqrt (dx * dx + dy * dy)
+	    ns = stf_r2n (r)
+	    ds = 1. / ns
+	    da = ds * ds
+	    dz = 0.5 + 0.5 * ds
+
+	    r1 = 0.
+	    dy1 = dy - dz
+	    do l = 1, ns {
+		dy1 = dy1 + ds
+		dy2 = dy1 * dy1
+		dx1 = dx - dz
+		do k = 1, ns {
+		    dx1 = dx1 + ds
+		    dx2 = dx1 * dx1
+		    r2 = (dx2 + dy2) / sigma
+		    if (r2 < 25.)
+			r1 = r1 + exp (-r2)
+		}
+	    }
+	    ptr = data + (j - y1 + 1) * nx + (i - x1 + 1)
+	    peak = (Memr[ptr] - bkgd) / (r1 * da)
+
+	    # Subtract the gaussian from the data.
+	    do j = y1, y2 {
+		ptr = data + (j - y1 + 1) * nx + 1
+		dy = j - yc
+		do i = x1, x2 {
+		    dx = i - xc
+		    r = sqrt (dx * dx + dy * dy)
+		    ns = stf_r2n (r)
+		    ds = 1. / ns
+		    da = ds * ds
+		    dz = 0.5 + 0.5 * ds
+
+		    r1 = 0.
+		    dy1 = dy - dz
+		    do l = 1, ns {
+			dy1 = dy1 + ds
+			dy2 = dy1 * dy1
+			dx1 = dx - dz
+			do k = 1, ns {
+			    dx1 = dx1 + ds
+			    dx2 = dx1 * dx1
+			    r2 = (dx2 + dy2) / sigma
+			    if (r2 < 25.)
+				r1 = r1 + peak * exp (-r2)
+			}
+		    }
+		    Memr[ptr] = Memr[ptr] - r1 * da
+		    ptr = ptr + 1
+		}
+	    }
+
+	    # Fit the image interpolator to the residual data.
+	    call msiinit (msi, II_BISPLINE3)
+	    call msifit (msi, Memr[data], nx, ny, nx)
+
+	    # Recompute the enclosed flux profile and moments
+	    # using the gaussian plus image interpolator fit to the residuals.
+
+	    call aclrr (Memr[profile], np)
+
+	    xx = 0.
+	    yy = 0.
+	    xy = 0.
+	    do j = y1, y2 {
+		ptr = data + (j - y1 + 1) * nx + 1
+		dy = j - yc
+		do i = x1, x2 {
+		    dx = i - xc
+		    r = sqrt (dx * dx + dy * dy)
+		    ns = stf_r2n (r)
+		    ds = 1. / ns
+		    da = ds * ds
+		    dz = 0.5 + 0.5 * ds
+
+		    # Compute interpolator correction.
+		    r2 = 0.
+		    dy1 = dy - dz
+		    do l = 1, ns {
+			dy1 = dy1 + ds
+			dx1 = dx - dz
+			do k = 1, ns {
+			    dx1 = dx1 + ds
+			    r1 = msieval (msi, dx1+xc-x1+2, dy1+yc-y1+2)
+			    r2 = r2 + r1
+			}
+		    }
+
+		    r1 = Memr[ptr] - bkgd
+		    ptr = ptr + 1
+		    r2 = r1 - r2 * da
+
+		    # Accumulate the enclosed flux and moments.
+		    dy1 = dy - dz
+		    do l = 1, ns {
+			dy1 = dy1 + ds
+			dy2 = dy1 * dy1
+			dx1 = dx - dz
+			do k = 1, ns {
+			    dx1 = dx1 + ds
+			    dx2 = dx1 * dx1
+			    r3 = (dx2 + dy2) / sigma
+			    if (r3 < 25.)
+				r3 = peak * exp (-r3)
+			    else
+				r3 = 0.
+			    r = max (0., sqrt (dx2 + dy2) - ds / 2)
+			    if (r < radius) {
+				r1 = msieval (msi, dx1+xc-x1+2, dy1+yc-y1+2)
+				r1 = da * (r1 + r2 + r3)
+
+				for (m=stf_r2i(r)+1; m<=np; m=m+1) {
+				    r3 = (stf_i2r (real(m)) - r) / ds
+				    if (r3 >= 1.)
+					break
+				    Memr[profile+m-1] = Memr[profile+m-1] +
+					r3 * r1
+				}
+				for (; m<=np; m=m+1)
+				    Memr[profile+m-1] = Memr[profile+m-1] + r1
+
+				if (r1 > miso) {
+				    xx = xx + dx2 * r1
+				    yy = yy + dy2 * r1
+				    xy = xy + dx1 * dy1 * r1
+				}
+			    }
+			}
+		    }
+		}
+	    }
+
+	    call msifree (msi)
+
+	    # Recompute the moments, magnitude, normalized flux, and interp.
+	    r = (xx + yy)
+	    if (r > 0.) {
+		r1 = (xx - yy) / r
+		r2 = 2 * xy / r
+		SFD_E(sfd) = sqrt (r1**2 + r2**2)
+		SFD_PA(sfd) = RADTODEG (atan2 (r2, r1) / 2.)
+	    } else {
+		SFD_E(sfd) = 0.
+		SFD_PA(sfd) = 0.
+	    }
+
+	    call alimr (Memr[profile], np, r, SFD_M(sfd))
+	    if (SFD_M(sfd) <= 0.)
+		call error (1, "Invalid flux profile")
+	    call adivkr (Memr[profile], SFD_M(sfd), Memr[profile], np)
+
+	    call asifit (asi, Memr[profile], np)
+	    SFD_ASI1(sfd) = asi
+	}
+
+	# Compute derivative of enclosed flux profile and fit an image
+	# interpolator.
+
+	dx = 0.25
+	Memr[profile] = 0.
+	ns = 0
+	do i = 1, np {
+	    r = stf_i2r (real(i))
+	    r2 = stf_r2i (r + dx)
+	    if (r2 > np) {
+		k = i
+		break
+	    }
+	    r1 = stf_r2i (r - dx)
+	    if (r1 < 1) {
+		if (i > 1) {
+		    dy = asieval (asi, real(i)) / r**2
+		    Memr[profile] = (ns * Memr[profile] + dy) / (ns + 1)
+		    ns = ns + 1
+		}
+		j = i
+	    } else {
+		dy = (asieval (asi, r2) - asieval (asi, r1)) /
+		    (4 * r * dx)
+		Memr[profile+i-1] = dy
+	    }
+	}
+	do i = 2, j
+	    Memr[profile+i-1] = (Memr[profile+j] - Memr[profile]) / j *
+		(i - 1) + Memr[profile]
+	do i = k, np
+	    Memr[profile+i-1] = Memr[profile+k-2]
+
+	call adivkr (Memr[profile], SF_SCALE(sf)**2, Memr[profile], np)
+	call alimr (Memr[profile], np, SFD_YP1(sfd), SFD_YP2(sfd))
+	call asiinit (asi, II_SPLINE3)
+	call asifit (asi, Memr[profile], np)
+	SFD_ASI2(sfd) = asi
+	#SF_XP1(sf) = j+1
+	SF_XP1(sf) = 1
+	SF_XP2(sf) = k-1
+
+	call sfree (sp)
+end
+
+
+# STF_NORM -- Renormalize the enclosed flux profile.
+
+procedure stf_norm (sf, sfd, x, y)
+
+pointer	sf			#I Parameter structure
+pointer	sfd			#I Star structure
+real	x			#I Radius
+real	y			#I Flux
+
+int	npmax, np
+pointer	asi
+
+int	i, j, k
+real	r, r1, r2, dx, dy
+pointer	sp, profile
+real	asieval(), stf_i2r(), stf_r2i()
+errchk	asifit
+
+begin
+	npmax = SFD_NPMAX(sfd)
+	np = SFD_NP(sfd)
+	asi = SFD_ASI1(sfd)
+
+	call smark (sp)
+	call salloc (profile, npmax, TY_REAL)
+
+	# Renormalize the enclosed flux profile.
+	if (IS_INDEF(x) || x <= 0.) {
+	    dy = SFD_BKGD(sfd) - SFD_BKGD1(sfd)
+	    SFD_BKGD(sfd) = SFD_BKGD(sfd) - dy
+	    do i = 1, npmax
+		Memr[profile+i-1] = asieval (asi, real(i)) +
+		    dy * stf_i2r(real(i)) ** 2
+	    call alimr (Memr[profile], np, r1, r2)
+	    call adivkr (Memr[profile], r2, Memr[profile], npmax)
+	} else if (IS_INDEF(y)) {
+	    r = max (1., min (real(np), stf_r2i (x)))
+	    r2 = asieval (asi, r)
+	    if (r2 <= 0.)
+		return
+	    do i = 1, npmax
+		Memr[profile+i-1] = asieval (asi, real(i))
+	    call adivkr (Memr[profile], r2, Memr[profile], npmax)
+	} else {
+	    r = max (1., min (real(np), stf_r2i (x)))
+	    r1 = asieval (asi, r)
+	    dy = (y - r1) / x ** 2
+	    SFD_BKGD(sfd) = SFD_BKGD(sfd) - dy
+	    do i = 1, npmax
+		Memr[profile+i-1] = asieval (asi, real(i)) +
+		    dy * stf_i2r(real(i)) ** 2
+	}
+
+	call asifit (asi, Memr[profile], npmax)
+	SFD_ASI1(sfd) = asi
+
+	# Compute derivative of enclosed flux profile and fit an image
+	# interpolator.
+
+	dx = 0.25
+	do i = 1, npmax {
+	    r = stf_i2r (real(i))
+	    r2 = stf_r2i (r + dx)
+	    if (r2 > np) {
+		k = i
+		break
+	    }
+	    r1 = stf_r2i (r - dx)
+	    if (r1 < 1) {
+		if (i > 1) {
+		    dy = asieval (asi, real(i)) / r**2
+		    Memr[profile] = dy
+		}
+		j = i
+	    } else {
+		dy = (asieval (asi, r2) - asieval (asi, r1)) /
+		    (4 * r * dx)
+		Memr[profile+i-1] = dy
+	    }
+	}
+	do i = 2, j
+	    Memr[profile+i-1] = (Memr[profile+j] - Memr[profile]) / j *
+		(i - 1) + Memr[profile]
+	do i = k, npmax
+	    Memr[profile+i-1] = Memr[profile+k-2]
+
+	call adivkr (Memr[profile], SF_SCALE(sf)**2, Memr[profile], np)
+	call alimr (Memr[profile], np, SFD_YP1(sfd), SFD_YP2(sfd))
+	asi = SFD_ASI2(sfd)
+	call asifit (asi, Memr[profile], np)
+	SFD_ASI2(sfd) = asi
+	#SF_XP1(sf) = min (j+1, np)
+	SF_XP1(sf) = 1
+	SF_XP2(sf) = min (k-1, np)
+
+	call sfree (sp)
+end
+
+
+# STF_WIDTHS -- Set the widhts.
+
+procedure stf_widths (sf, sfd)
+
+pointer	sf			#I Main data structure
+pointer	sfd			#I Star data structure
+
+errchk	stf_radius, stf_dfwhm, stf_fit
+
+begin
+	call stf_radius (sf, sfd, SF_LEVEL(sf), SFD_R(sfd))
+	call stf_dfwhm (sf, sfd)
+	call stf_fit (sf, sfd)
+
+	switch (SF_WCODE(sf)) {
+	case 1:
+	    SFD_W(sfd) = SFD_R(sfd)
+	case 2:
+	    SFD_W(sfd) = SFD_DFWHM(sfd)
+	case 3:
+	    SFD_W(sfd) = SFD_GFWHM(sfd)
+	case 4:
+	    SFD_W(sfd) = SFD_MFWHM(sfd)
+	}
+end
+
+
+# STF_I2R -- Compute radius from sample index.
+
+real procedure stf_i2r (i)
+
+real	i			#I Index
+real	r			#O Radius
+
+begin
+	if (i < 20)
+	    r = 0.05 * i
+	else if (i < 30)
+	    r = 0.1 * i - 1
+	else if (i < 40)
+	    r = 0.2 * i - 4
+	else if (i < 50)
+	    r = 0.5 * i - 16
+	else
+	    r = i - 41
+	return (r)
+end
+
+
+# STF_R2I -- Compute sample index from radius.
+
+real procedure stf_r2i (r)
+
+real	r			#I Radius
+real	i			#O Index
+
+begin
+	if (r < 1)
+	    i = 20 * r
+	else if (r < 2)
+	    i = 10 * (r + 1)
+	else if (r < 4)
+	    i = 5 * (r + 4)
+	else if (r < 9)
+	    i = 2 * (r + 16)
+	else
+	    i = r + 41
+	return (i)
+end
+
+
+# STF_R2N -- Compute number of subsamples from radius.
+
+int procedure stf_r2n (r)
+
+real	r			#I Radius
+int	n			#O Number of subsamples
+
+begin
+	if (r < 1)
+	    n = 20
+	else if (r < 2)
+	    n = 10
+	else if (r < 4)
+	    n = 5
+	else if (r < 9)
+	    n = 2
+	else
+	    n = 1
+	return (n)
+end
+
+
+# STF_MODEL -- Return model value.
+
+procedure stf_model (sf, sfd, r, profile, flux)
+
+pointer	sf			#I Main data structure
+pointer	sfd			#I Star data structure
+real	r			#I Radius at level
+real	profile			#I Profile value
+real	flux			#I Enclosed flux value
+
+real	x, x1, x2, r1, r2, dr
+
+begin
+	dr = 0.25 * SF_SCALE(sf)
+	r1 = r - dr
+	r2 = r + dr
+	if (r1 < 0.) {
+	    r1 = dr
+	    r2 = r1 + dr
+	}
+
+	switch (SF_WCODE(sf)) {
+	case 3:
+	    x = r**2 / (2. * SFD_SIGMA(sfd)**2)
+	    if (x < 20.)
+		flux = 1 - exp (-x)
+	    else
+		flux = 0.
+
+	    x1 = r1**2 / (2. * SFD_SIGMA(sfd)**2)
+	    x2 = r2**2 / (2. * SFD_SIGMA(sfd)**2)
+	    if (x2 < 20.) {
+		x1 = 1 - exp (-x1)
+		x2 = 1 - exp (-x2)
+	    } else {
+		x1 = 1.
+		x2 = 1.
+	    }
+	    if (r <= dr) {
+		x1 = x1 / dr ** 2
+		x2 = x2 / (4 * dr ** 2)
+		profile = (x2 - x1) / dr * r + x1
+	    } else {
+		profile = (x2 - x1) / (4 * r * dr)
+	    }
+	default:
+	    x = 1 + (r / SFD_ALPHA(sfd)) ** 2
+	    flux = 1 - x ** (1 - SFD_BETA(sfd))
+
+	    x1 = 1 + (r1 / SFD_ALPHA(sfd)) ** 2
+	    x2 = 1 + (r2 / SFD_ALPHA(sfd)) ** 2
+	    x1 = 1 - x1 ** (1 - SFD_BETA(sfd))
+	    x2 = 1 - x2 ** (1 - SFD_BETA(sfd))
+	    if (r <= dr) {
+		x1 = x1 / dr ** 2
+		x2 = x2 / (4 * dr ** 2)
+		profile = (x2 - x1) / dr * r + x1
+	    } else {
+		profile = (x2 - x1) / (4 * r * dr)
+	    }
+	}
+end
+
+
+# STF_DFWHM -- Direct FWHM from profile.
+
+procedure stf_dfwhm (sf, sfd)
+
+pointer	sf			#I Main data structure
+pointer	sfd			#I Star data structure
+
+int	np
+real	r, rpeak, profile, peak, asieval(), stf_i2r()
+pointer	asi
+
+begin
+	asi = SFD_ASI2(sfd)
+	np = SFD_NP(sfd)
+
+	rpeak = 1.
+	peak = 0.
+	for (r=1.; r <= np; r = r + 0.01) {
+	    profile = asieval (asi, r)
+	    if (profile > peak) {
+		rpeak = r
+		peak = profile
+	    }
+	}
+	
+	peak = peak / 2.
+	for (r=rpeak; r <= np && asieval (asi, r) > peak; r = r + 0.01)
+	    ;
+
+	SFD_DFWHM(sfd) = 2 * stf_i2r (r) * SF_SCALE(sf)
+end
+
+
+# STF_FWHMS -- Measure FWHM vs level.
+
+procedure stf_fwhms (sf, sfd)
+
+pointer	sf			#I Main data structure
+pointer	sfd			#I Star data structure
+
+int	i
+real	level, r
+
+begin
+	do i = 1, 19 {
+	    level = i * 0.05
+	    call stf_radius (sf, sfd, level, r)
+	    switch (SF_WCODE(sf)) {
+	    case 3:
+		SFD_FWHM(sfd,i) = 2 * r * sqrt (log (2.) / log (1/(1-level)))
+	    default:
+		r = r / sqrt ((1.-level)**(1./(1.-SFD_BETA(sfd))) - 1.)
+		SFD_FWHM(sfd,i) = 2 * r * sqrt (2.**(1./SFD_BETA(sfd))-1.)
+	    }
+	}
+end
+
+
+# STF_RADIUS -- Measure the radius at the specified level.
+
+procedure stf_radius (sf, sfd, level, r)
+
+pointer	sf			#I Main data structure
+pointer	sfd			#I Star data structure
+real	level			#I Level to measure
+real	r			#O Radius
+
+int	np
+pointer	asi
+real	f, fmax, rmax, asieval(), stf_i2r()
+
+begin
+	np = SFD_NP(sfd)
+	asi = SFD_ASI1(sfd)
+
+	for (r=1; r <= np && asieval (asi, r) < level; r = r + 0.01)
+	    ;
+	if (r > np) {
+	    fmax = 0.
+	    rmax = 0.
+	    for (r=1; r <= np; r = r + 0.01) {
+		f = asieval (asi, r)
+		if (f > fmax) {
+		    fmax = f
+		    rmax = r
+		}
+	    }
+	    r = rmax
+	}
+	r = stf_i2r (r) * SF_SCALE(sf)
+end
+
+
+# STF_FIT -- Fit models to enclosed flux.
+
+procedure stf_fit (sf, sfd)
+
+pointer	sf			#I Main data structure
+pointer	sfd			#I Star data structure
+
+int	i, j, n, np, pfit[2]
+real	beta, z, params[3]
+pointer	asi, nl
+pointer	sp, x, y, w
+
+int	locpr()
+real	asieval(), stf_i2r()
+extern	stf_gauss1(), stf_gauss2(), stf_moffat1(), stf_moffat2()
+errchk	nlinitr, nlfitr
+
+data	pfit/2,3/
+
+begin
+	np = SFD_NP(sfd)
+	asi = SFD_ASI1(sfd)
+
+	call smark (sp)
+	call salloc (x, np, TY_REAL)
+	call salloc (y, np, TY_REAL)
+	call salloc (w, np, TY_REAL)
+
+	n = 0
+	j = 0
+	do i = 1, np {
+	    z = 1. - max (0., asieval (asi, real(i)))
+	    if (n > np/3 && z < 0.5)
+		break
+	    if ((n < np/3 && z > 0.01) || z > 0.5)
+		j = n
+
+	    Memr[x+n] = stf_i2r (real(i)) * SF_SCALE(sf)
+	    Memr[y+n] = z
+	    Memr[w+n] = 1.
+	    n = n + 1
+	}
+
+	# Gaussian.
+	np = 1
+	params[2] = Memr[x+j] / sqrt (2. * log (1./min(0.99,Memr[y+j])))
+	params[1] = 1
+	call nlinitr (nl, locpr (stf_gauss1), locpr (stf_gauss2),
+	    params, params, 2, pfit, np, .001, 100)
+	call nlfitr (nl, Memr[x], Memr[y], Memr[w], n, 1, WTS_USER, i)
+	if (i != SINGULAR && i != NO_DEG_FREEDOM) {
+	    call nlpgetr (nl, params, i)
+	    if (params[2] < 0.)
+		params[2] = Memr[x+j] / sqrt (2. * log (1./min(0.99,Memr[y+j])))
+	}
+	SFD_SIGMA(sfd) = params[2]
+	SFD_GFWHM(sfd) = 2 * SFD_SIGMA(sfd) * sqrt (2. * log (2.))
+
+	# Moffat.
+	if (SF_BETA(sf) < 1.1) {
+	    call nlfreer (nl)
+	    call sfree (sp)
+	    call error (1, "Cannot measure FWHM - Moffat beta too small")
+	}
+
+	beta = SF_BETA(sf)
+	if (IS_INDEFR(beta)) {
+	    beta = 2.5
+	    np = 2
+	} else {
+	    np = 1
+	}
+	params[3] = 1 - beta
+	params[2] = Memr[x+j] / sqrt (min(0.99,Memr[y+j])**(1./params[3]) - 1.)
+	params[1] = 1
+	call nlinitr (nl, locpr (stf_moffat1), locpr (stf_moffat2),
+	    params, params, 3, pfit, np, .001, 100)
+	call nlfitr (nl, Memr[x], Memr[y], Memr[w], n, 1, WTS_USER, i)
+	if (i != SINGULAR && i != NO_DEG_FREEDOM) {
+	    call nlpgetr (nl, params, i)
+	    if (params[2] < 0.) {
+		params[3] = 1. - beta
+		params[2] = Memr[x+j] /
+		    sqrt (min(0.99,Memr[y+j])**(1./params[3]) - 1.)
+	    }
+	}
+	SFD_ALPHA(sfd) = params[2]
+	SFD_BETA(sfd) = 1 - params[3]
+	SFD_MFWHM(sfd) = 2 * SFD_ALPHA(sfd) * sqrt (2.**(1./SFD_BETA(sfd))-1.)
+
+	call nlfreer (nl)
+	call sfree (sp)
+end
+
+
+# STF_GAUSS1 -- Gaussian function used in NLFIT.  The parameters are the
+# amplitude and sigma and the input variable is the radius.
+
+procedure stf_gauss1 (x, nvars, p, np, z)
+
+real	x[nvars]		#I Input variables
+int	nvars			#I Number of variables
+real	p[np]			#I Parameter vector
+int	np			#I Number of parameters
+real	z			#O Function return
+
+real	r2
+
+begin
+	r2 = x[1]**2 / (2 * p[2]**2)
+	if (abs (r2) > 20.)
+	    z = 0.
+	else
+	    z = p[1] * exp (-r2)
+end
+
+
+# STF_GAUSS2 -- Gaussian function and derivatives used in NLFIT.  The parameters
+# are the amplitude and sigma and the input variable is the radius.
+
+procedure stf_gauss2 (x, nvars, p, dp, np, z, der)
+
+real	x[nvars]		#I Input variables
+int	nvars			#I Number of variables
+real	p[np]			#I Parameter vector
+real	dp[np]			#I Dummy array of parameters increments
+int	np			#I Number of parameters
+real	z			#O Function return
+real	der[np]			#O Derivatives
+
+real	r2
+
+begin
+	r2 = x[1]**2 / (2 * p[2]**2)
+	if (abs (r2) > 20.) {
+	    z = 0.
+	    der[1] = 0.
+	    der[2] = 0.
+	} else {
+	    der[1] = exp (-r2)
+	    z = p[1] * der[1]
+	    der[2] = z * 2 * r2 / p[2]
+	}
+end
+
+
+# STF_MOFFAT1 -- Moffat function used in NLFIT.  The parameters are the
+# amplitude, alpha squared, and beta and the input variable is the radius.
+
+procedure stf_moffat1 (x, nvars, p, np, z)
+
+real	x[nvars]		#I Input variables
+int	nvars			#I Number of variables
+real	p[np]			#I Parameter vector
+int	np			#I Number of parameters
+real	z			#O Function return
+
+real	y
+
+begin
+	y = 1 + (x[1] / p[2]) ** 2
+	if (abs (y) > 20.)
+	    z = 0.
+	else
+	    z = p[1] * y ** p[3]
+end
+
+
+# STF_MOFFAT2 -- Moffat function and derivatives used in NLFIT.  The
+# parameters are the amplitude, alpha squared, and beta and the input
+# variable is the radius.
+
+procedure stf_moffat2 (x, nvars, p, dp, np, z, der)
+
+real	x[nvars]		#I Input variables
+int	nvars			#I Number of variables
+real	p[np]			#I Parameter vector
+real	dp[np]			#I Dummy array of parameters increments
+int	np			#I Number of parameters
+real	z			#O Function return
+real	der[np]			#O Derivatives
+
+real	y
+
+begin
+	y = 1 + (x[1] / p[2]) ** 2
+	if (abs (y) > 20.) {
+	    z = 0.
+	    der[1] = 0.
+	    der[2] = 0.
+	    der[3] = 0.
+	} else {
+	    der[1] = y ** p[3]
+	    z = p[1] * der[1]
+	    der[2] = -2 * z / y * p[3] / p[2] * (x[1] / p[2]) ** 2
+	    der[3] = z * log (y)
+	}
+end
diff --git a/noao/obsutil/src/starfocus/t_starfocus.x b/noao/obsutil/src/starfocus/t_starfocus.x
new file mode 100644
index 00000000..2c3213e9
--- /dev/null
+++ b/noao/obsutil/src/starfocus/t_starfocus.x
@@ -0,0 +1,1240 @@
+include	<error.h>
+include	<imhdr.h>
+include	<imset.h>
+include	<mach.h>
+include	"starfocus.h"
+
+define	HELP		"nmisc$src/starfocus.key"
+define	PROMPT		"Options"
+
+
+# T_STARFOCUS -- Stellar focusing task.
+
+procedure t_starfocus ()
+
+begin
+	call starfocus (STARFOCUS)
+end
+
+
+# T_PSFMEASURE -- PSF measuring task.
+
+procedure t_psfmeasure ()
+
+begin
+	call starfocus (PSFMEASURE)
+end
+
+
+# STARFOCUS -- Stellar focusing and PSF measuring main routine.
+
+procedure starfocus (type)
+
+int	type		#I Task type
+
+int	list		# List of images
+pointer	fvals		# Focus values
+pointer	fstep		# Focus step
+pointer	nexposures	# Number of exposures
+pointer	step		# step in pixels
+int	direction	# Step direction
+int	gap		# Double step gap
+int	coords		# Type of image data
+bool	display		# Display images?
+int	frame		# Display frame
+int	logfd		# Log file descriptor
+bool	ignore_sat	# Ignore saturation?
+
+real	wx, wy, f, df, xstep, ystep
+int	i, i1, i2, i3, j, k, l, ip, wcs, key, id, ncols, nlines
+int	nexp, nsfd, nimages, nstars, ngraph, nmark
+pointer	sp, sf, image, system, cmd, rg, mark, im, mw, ct
+pointer	sfds, sfd
+
+bool	clgetb(), streq()
+real	clgetr(), imgetr(), stf_r2i()
+int	clgeti(), clgwrd(), clgcur(), imtopenp(), imtgetim(), imgeti()
+int	nowhite(), open(), rng_index(), strdic(), ctoi(), ctor()
+pointer	rng_open(), immap(), mw_openim(), mw_sctran()
+errchk	immap, open, imgetr, imgeti, mw_openim, mw_sctran
+errchk	stf_find, stf_bkgd, stf_profile, stf_widths, stf_fwhms, stf_radius
+errchk	stf_organize, stf_graph, stf_display
+
+begin
+	call smark (sp)
+	call salloc (sf, SF, TY_STRUCT)
+	call salloc (image, SZ_FNAME, TY_CHAR)
+	call salloc (fvals, SZ_LINE, TY_CHAR)
+	call salloc (fstep, SZ_LINE, TY_CHAR)
+	call salloc (nexposures, SZ_LINE, TY_CHAR)
+	call salloc (step, SZ_LINE, TY_CHAR)
+	call salloc (system, SZ_LINE, TY_CHAR)
+	call salloc (cmd, SZ_LINE, TY_CHAR)
+
+	call aclri (Memi[sf], SF)
+	SF_TASK(sf) = type
+
+	# Set task parameters.
+	switch (type) {
+	case STARFOCUS:
+	    call clgstr ("focus", Memc[fvals], SZ_LINE)
+	    call clgstr ("fstep", Memc[fstep], SZ_LINE)
+	    call clgstr ("nexposures", Memc[nexposures], SZ_LINE)
+	    call clgstr ("step", Memc[step], SZ_LINE)
+
+	    direction = clgwrd ("direction", Memc[cmd], SZ_LINE,
+		"|-line|+line+|-column|+column|")
+	    gap = clgwrd ("gap", Memc[cmd], SZ_LINE, "|none|beginning|end|") 
+
+	    if (nowhite (Memc[fvals], Memc[fvals], SZ_LINE) != 0) {
+		iferr (rg = rng_open (Memc[fvals], -MAX_REAL, MAX_REAL, 1.))
+		    rg = NULL
+	    } else
+		rg = NULL
+	case PSFMEASURE:
+	    Memc[fvals] = EOS
+	    rg = NULL
+	    nexp = 1
+	}
+
+	list = imtopenp ("images")
+	display = clgetb ("display")
+	frame = clgeti ("frame")
+	coords = clgwrd ("coords", Memc[cmd], SZ_LINE, SF_TYPES)
+	call clgstr ("wcs", Memc[system], SZ_LINE)
+
+	SF_XF(sf) = clgetr ("xcenter")
+	SF_YF(sf) = clgetr ("ycenter")
+	SF_LEVEL(sf) = clgetr ("level")
+	SF_WCODE(sf) = clgwrd ("size", SF_WTYPE(sf), SF_SZWTYPE, SF_WTYPES)
+	SF_BETA(sf) = clgetr ("beta")
+	SF_SCALE(sf) = clgetr ("scale")
+	SF_RADIUS(sf) = max (3., clgetr ("radius"))
+	SF_NIT(sf) = clgeti ("iterations")
+	SF_SBUF(sf) = clgetr ("sbuffer")
+	SF_SWIDTH(sf) = clgetr ("swidth")
+	SF_SAT(sf) = clgetr ("saturation")
+	ignore_sat = clgetb ("ignore_sat")
+	SF_OVRPLT(sf) = NO
+
+	if (SF_LEVEL(sf) > 1.)
+	    SF_LEVEL(sf) = SF_LEVEL(sf) / 100.
+	SF_LEVEL(sf) = max (0.05, min (0.95, SF_LEVEL(sf)))
+
+	# Accumulate the psf/focus data.
+	key = 'm'
+	mark = NULL
+	nstars = 0
+	nmark = 0
+	ngraph = 0
+	nimages = 0
+	nsfd = 0
+	while (imtgetim (list, Memc[image], SZ_FNAME) != EOF) {
+	    im = immap (Memc[image], READ_ONLY, 0)
+	    call imseti (im, IM_TYBNDRY, TYBNDRY)
+	    call imseti (im, IM_NBNDRYPIX, NBNDRYPIX)
+	    if (streq (Memc[system], "logical")) {
+		mw = NULL
+		ct = NULL
+	    } else {
+		mw = mw_openim (im)
+		ct = mw_sctran (mw, Memc[system], "logical", 03)
+	    }
+	    ncols = IM_LEN(im,1)
+	    nlines = IM_LEN(im,2)
+	    nimages = nimages + 1
+	    if (nimages == 1) {
+		SF_NCOLS(sf) = ncols
+		SF_NLINES(sf) = nlines
+		if (IS_INDEF(SF_XF(sf)))
+		    SF_XF(sf) = (SF_NCOLS(sf) + 1) / 2.
+		if (IS_INDEF(SF_YF(sf)))
+		    SF_YF(sf) = (SF_NLINES(sf) + 1) / 2.
+	    } else if (ncols!=SF_NCOLS(sf)||nlines!=SF_NLINES(sf))
+		call eprintf ("WARNING: Images have different sizes\n")
+
+	    # Display the image if needed.
+	    if (display) {
+		switch (coords) {
+		case SF_MARK1:
+		    if (nimages == 1)
+			call stf_display (Memc[image], frame)
+		case SF_MARKALL:
+		    call stf_display (Memc[image], frame)
+		}
+		if (nimages == 1) {
+		    call printf (
+		 "** Select stars to measure with 'm' and finish with 'q'.\n")
+		    call printf (
+		 "** Additional options are '?', 'g', and :show.\n")
+		call flush (STDOUT)
+		}
+	    }
+
+	    # Accumulate objects.
+	    repeat {
+		switch (coords) {
+		case SF_CENTER:
+		    if (nstars == nimages)
+			break
+		    if (rg == NULL && Memc[fvals] == EOS)
+			id = nstars
+		    else
+			id = 0
+		    wx = 1 + (ncols - 1) / 2.
+		    wy = 1 + (nlines - 1) / 2.
+		    key = 0
+		case SF_MARK1:
+		    if (nimages == 1) {
+			if (clgcur ("imagecur", wx, wy, wcs, key,
+			    Memc[cmd], SZ_LINE) == EOF)
+			    break
+			switch (key) {
+			case '?':
+			    call pagefile (HELP, PROMPT)
+			    next
+			case ':':
+			    if (strdic (Memc[cmd], Memc[cmd], SZ_LINE,
+				"|show|") == 1) {
+				if (nsfd > 0) {
+				    call stf_organize (sf, sfds, nsfd)
+				    call mktemp ("tmp$iraf", Memc[cmd], SZ_LINE)
+				    logfd = open (Memc[cmd], APPEND, TEXT_FILE)
+				    call stf_log (sf, logfd)
+				    call close (logfd)
+				    call pagefile (Memc[cmd], "starfocus")
+				    call delete (Memc[cmd])
+				}
+			    }
+			    next
+			case 'q':
+			    break
+			}
+			id = nstars
+			if (mark == NULL)
+			    call malloc (mark, 3*10, TY_REAL)
+			else if (mod (nmark, 10) == 0)
+			    call realloc (mark, 3*(nmark+10), TY_REAL)
+			Memr[mark+3*nmark] = id
+			Memr[mark+3*nmark+1] = wx
+			Memr[mark+3*nmark+2] = wy
+			nmark = nmark+1
+		    } else {
+			if (nmark == 0)
+			    break
+			if (nstars / nmark == nimages)
+			    break
+			i = mod (nstars, nmark)
+			id = Memr[mark+3*i]
+			wx = Memr[mark+3*i+1]
+			wy = Memr[mark+3*i+2]
+			key = 0
+		    }
+		    if (ct != NULL)
+			call mw_c2tranr (ct, wx, wy, wx, wy)
+		case SF_MARKALL:
+		    if (clgcur ("imagecur", wx, wy, wcs, key,
+			Memc[cmd], SZ_LINE) == EOF)
+			break
+		    switch (key) {
+		    case '?':
+			call pagefile (HELP, PROMPT)
+			next
+		    case ':':
+			if (strdic(Memc[cmd],Memc[cmd],SZ_LINE,"|show|")==1) {
+			    if (nsfd > 0) {
+				call stf_organize (sf, sfds, nsfd)
+				call mktemp ("tmp$iraf", Memc[cmd], SZ_LINE)
+				logfd = open (Memc[cmd], APPEND, TEXT_FILE)
+				call stf_log (sf, logfd)
+				call close (logfd)
+				call pagefile (Memc[cmd], "starfocus")
+				call delete (Memc[cmd])
+			    }
+			}
+			next
+		    case 'q':
+			break
+		    }
+		    id = nstars
+		    if (ct != NULL)
+			call mw_c2tranr (ct, wx, wy, wx, wy)
+		}
+
+		if (type == STARFOCUS) {
+		    ip = 1
+		    if (ctoi (Memc[nexposures], ip, nexp) == 0)
+			nexp = imgeti (im, Memc[nexposures])
+		    ip = 1
+		    if (ctor (Memc[step], ip, f) == 0)
+			f = imgetr (im, Memc[step])
+
+		    xstep = 0.
+		    ystep = 0.
+		    switch (direction) {
+		    case 1:
+			ystep = -f
+		    case 2:
+			ystep = f
+		    case 3:
+			xstep = -f
+		    case 4:
+			xstep = f
+		    }
+
+		    # Set the steps and order of evaluation.
+		    # This depends on which star in the sequence is marked.
+		    # Below we assume the minimum x or maximum y is marked.
+
+		    i1 = 1; i2 = nexp; i3 = 1
+#		    if (xstep < 0.) {
+#			i1 = nexp; i2 = 1; i3 = -1; xstep = -xstep
+#		    }
+#		    if (ystep > 0.) {
+#			i1 = nexp; i2 = 1; i3 = -1; ystep = -ystep
+#		    }
+		} else {
+		    i1 = 1; i2 = 1; i3 = 1
+		}
+
+		k = nsfd
+		do i = i1, i2, i3 {
+		    if (i > 1) {
+			wx = wx + xstep
+			wy = wy + ystep
+			switch (gap) {
+			case 2:
+			    if ((i==2 && i3==1) || (i==1 && i3==-1)) {
+				wx = wx + xstep
+				wy = wy + ystep
+			    }
+			case 3:
+			    if ((i==nexp && i3==1) || (i==nexp-1 && i3==-1)) {
+				wx = wx + xstep
+				wy = wy + ystep
+			    }
+			}
+		    }
+				
+		    if (wx < SF_RADIUS(sf)-NBNDRYPIX ||
+			wx > IM_LEN(im,1)-SF_RADIUS(sf)+NBNDRYPIX ||
+		        wy < SF_RADIUS(sf)-NBNDRYPIX ||
+			wy > IM_LEN(im,2)-SF_RADIUS(sf)+NBNDRYPIX)
+			next
+		    if (nexp == 1)
+			j = nimages
+		    else
+			j = i
+		    if (nsfd == 0)
+			call malloc (sfds, 10, TY_POINTER)
+		    else if (mod (nsfd, 10) == 0)
+			call realloc (sfds, nsfd+10, TY_POINTER)
+		    call malloc (sfd, SFD, TY_STRUCT)
+		    call strcpy (Memc[image], SFD_IMAGE(sfd), SF_SZFNAME)
+		    SFD_ID(sfd) = id
+		    SFD_X(sfd) = wx
+		    SFD_Y(sfd) = wy
+		    if (Memc[fvals] == EOS)
+			f = INDEF
+		    else if (Memc[fstep] != EOS) {
+			ip = 1
+			if (ctor (Memc[fvals], ip, f) == 0)
+			    f = imgetr (im, Memc[fvals])
+			ip = 1
+			if (ctor (Memc[fstep], ip, df) == 0)
+			    df = imgetr (im, Memc[fstep])
+			f = f + (i - 1) * df
+		    } else if (rg != NULL) {
+			if (rng_index (rg, j, f) == EOF)
+			    call error (1, "Focus list ended prematurely")
+		    } else
+			f = imgetr (im, Memc[fvals])
+		    SFD_F(sfd) = f
+		    SFD_STATUS(sfd) = 0
+		    SFD_SFS(sfd) = NULL
+		    SFD_SFF(sfd) = NULL
+		    SFD_SFI(sfd) = NULL
+
+		    iferr {
+			do l = 1, SF_NIT(sf) {
+			    if (l == 1)
+				SFD_RADIUS(sfd) = max (3., SF_RADIUS(sf))
+			    else
+				SFD_RADIUS(sfd) = max (3., 3. * SFD_DFWHM(sfd))
+			    SFD_NPMAX(sfd) = stf_r2i (SFD_RADIUS(sfd)) + 1
+			    SFD_NP(sfd) = SFD_NPMAX(sfd)
+			    call stf_find (sf, sfd, im)
+			    call stf_bkgd (sf, sfd)
+			    if (SFD_NSAT(sfd) > 0 && l == 1) {
+				if (ignore_sat)
+				    call error (0,
+				    "Saturated pixels found - ignoring object")
+				else
+				    call eprintf (
+					"WARNING: Saturated pixels found.\n")
+			    }
+			    call stf_profile (sf, sfd)
+			    call stf_widths (sf, sfd)
+			    call stf_fwhms (sf, sfd)
+			}
+			Memi[sfds+nsfd] = sfd
+			nsfd = nsfd + 1
+			wx = SFD_X(sfd)
+			wy = SFD_Y(sfd)
+		    } then {
+			call erract (EA_WARN)
+			call mfree (sfd, TY_STRUCT)
+		    }
+		}
+		if (nsfd > k) {
+		    nstars = nstars + 1
+		    if (key == 'g') {
+			if (nsfd - k > 0) {
+			    call stf_organize (sf, sfds+k, nsfd-k)
+			    call stf_graph (sf)
+			    ngraph = ngraph + 1
+			}
+		    }
+		}
+	    }
+	    if (mw != NULL)
+		call mw_close (mw)
+	    call imunmap (im)
+	}
+
+	if (nsfd == 0)
+	    call error (1, "No input data")
+
+	# Organize the objects, graph the data, and log the results.
+	if (nstars > 1 || ngraph != nstars) {
+	    call stf_organize (sf, sfds, nsfd)
+	    call stf_graph (sf)
+	}
+	call stf_log (sf, STDOUT)
+	call clgstr ("logfile", Memc[image], SZ_FNAME)
+	ifnoerr (logfd = open (Memc[image], APPEND, TEXT_FILE)) {
+	    call stf_log (sf, logfd)
+	    call close (logfd)
+	}
+
+	# Finish up
+	call rng_close (rg)
+	call imtclose (list)
+	call stf_free (sf)
+	do i = 1, SF_NSFD(sf) {
+	    sfd = SF_SFD(sf,i)
+	    call asifree (SFD_ASI1(sfd))
+	    call asifree (SFD_ASI2(sfd))
+	    call mfree (sfd, TY_STRUCT)
+	}
+	call mfree (SF_SFDS(sf), TY_POINTER)
+	call mfree (mark, TY_REAL)
+	call sfree (sp)
+end
+
+
+# STF_FREE -- Free the starfocus data structures.
+
+procedure stf_free (sf)
+
+pointer	sf			#I Starfocus structure
+int	i
+
+begin
+	do i = 1, SF_NSTARS(sf)
+	    call mfree (SF_SFS(sf,i), TY_STRUCT)
+	do i = 1, SF_NFOCUS(sf)
+	    call mfree (SF_SFF(sf,i), TY_STRUCT)
+	do i = 1, SF_NIMAGES(sf)
+	    call mfree (SF_SFI(sf,i), TY_STRUCT)
+	call mfree (SF_STARS(sf), TY_POINTER)
+	call mfree (SF_FOCUS(sf), TY_POINTER)
+	call mfree (SF_IMAGES(sf), TY_POINTER)
+	SF_NSTARS(sf) = 0
+	SF_NFOCUS(sf) = 0
+	SF_NIMAGES(sf) = 0
+end
+
+
+# STF_ORGANIZE -- Organize the individual object structures by star, focus,
+# and image.  Compute focus, radius, and magnitude by group and over all
+# data.
+
+procedure stf_organize (sf, sfds, nsfd)
+
+pointer	sf			#I Starfocus structure
+pointer	sfds			#I Pointer to array of object structures
+int	nsfd			#I Number of object structures
+
+int	i, j, k, nstars, nfocus, nimages, key
+real	f
+pointer	stars, focus, images, sfd, sfs, sff, sfi
+pointer	sp, image
+bool	streq()
+errchk	malloc
+
+int	stf_focsort(), stf_magsort()
+extern	stf_focsort, stf_magsort
+
+begin
+	call smark (sp)
+	call salloc (image, SZ_FNAME, TY_CHAR)
+
+	# Free previous structures.
+	call stf_free (sf)
+
+	# Organize sfds by star.
+	nstars = 0
+	for (i = 0; i < nsfd; i = i + 1) {
+	    key = SFD_ID(Memi[sfds+i])
+	    for (j = 0; SFD_ID(Memi[sfds+j]) != key; j = j + 1)
+		;
+	    if (j == i)
+		nstars = nstars + 1
+	}
+	call malloc (stars, nstars, TY_POINTER)
+
+	nstars = 0
+	for (i = 0; i < nsfd; i = i + 1) {
+	    key = SFD_ID(Memi[sfds+i])
+	    for (j = 0; j < nstars; j = j + 1)
+		if (SFS_ID(Memi[stars+j]) == key)
+		    break
+	    if (j == nstars) {
+		k = 0
+		for (j = i; j < nsfd; j = j + 1)
+		    if (SFD_ID(Memi[sfds+j]) == key)
+			k = k + 1
+		call malloc (sfs, SFS(k), TY_STRUCT)
+		SFS_ID(sfs) = key
+		SFS_NSFD(sfs) = k
+		k = 0
+		for (j = i; j < nsfd; j = j + 1) {
+		    sfd = Memi[sfds+j]
+		    if (SFD_ID(sfd) == key) {
+			k = k + 1
+			SFD_SFS(sfd) = sfs
+			SFS_SFD(sfs,k) = sfd
+		    }
+		}
+		Memi[stars+nstars] = sfs
+		nstars = nstars + 1
+	    }
+	}
+
+	# Organize sfds by focus values.  Sort by magnitude.
+	nfocus = 0
+	for (i = 0; i < nsfd; i = i + 1) {
+	    f = SFD_F(Memi[sfds+i])
+	    for (j = 0; SFD_F(Memi[sfds+j]) != f; j = j + 1)
+		;
+	    if (j == i)
+		nfocus = nfocus + 1
+	}
+	call malloc (focus, nfocus, TY_POINTER)
+
+	nfocus = 0
+	for (i = 0; i < nsfd; i = i + 1) {
+	    f = SFD_F(Memi[sfds+i])
+	    for (j = 0; j < nfocus; j = j + 1)
+		if (SFF_F(Memi[focus+j]) == f)
+		    break
+	    if (j == nfocus) {
+		k = 0
+		for (j = i; j < nsfd; j = j + 1)
+		    if (SFD_F(Memi[sfds+j]) == f)
+			k = k + 1
+		call malloc (sff, SFF(k), TY_STRUCT)
+		SFF_F(sff) = f
+		SFF_NSFD(sff) = k
+		k = 0
+		for (j = i; j < nsfd; j = j + 1) {
+		    sfd = Memi[sfds+j]
+		    if (SFD_F(sfd) == f) {
+			k = k + 1
+			SFD_SFF(sfd) = sff
+			SFF_SFD(sff,k) = sfd
+		    }
+		}
+		Memi[focus+nfocus] = sff
+		nfocus = nfocus + 1
+	    }
+	}
+
+	# Organize sfds by image.
+	nimages = 0
+	for (i = 0; i < nsfd; i = i + 1) {
+	    call strcpy (SFD_IMAGE(Memi[sfds+i]), Memc[image], SZ_FNAME)
+	    for (j = 0; !streq (SFD_IMAGE(Memi[sfds+j]), Memc[image]); j = j+1)
+		;
+	    if (j == i)
+		nimages = nimages + 1
+	}
+	call malloc (images, nimages, TY_POINTER)
+
+	nimages = 0
+	for (i = 0; i < nsfd; i = i + 1) {
+	    call strcpy (SFD_IMAGE(Memi[sfds+i]), Memc[image], SZ_FNAME)
+	    for (j = 0; j < nimages; j = j + 1)
+		if (streq (SFI_IMAGE(Memi[images+j]), Memc[image]))
+		    break
+	    if (j == nimages) {
+		k = 0
+		for (j = i; j < nsfd; j = j + 1)
+		    if (streq (SFD_IMAGE(Memi[sfds+j]), Memc[image]))
+			k = k + 1
+		call malloc (sfi, SFI(k), TY_STRUCT)
+		call strcpy (Memc[image], SFI_IMAGE(sfi), SF_SZFNAME)
+		SFI_NSFD(sfi) = k
+		k = 0
+		for (j = i; j < nsfd; j = j + 1) {
+		    sfd = Memi[sfds+j]
+		    if (streq (SFD_IMAGE(sfd), Memc[image])) {
+			k = k + 1
+			SFD_SFI(sfd) = sfi
+			SFI_SFD(sfi,k) = sfd
+		    }
+		}
+		Memi[images+nimages] = sfi
+		nimages = nimages + 1
+	    }
+	}
+		    
+	SF_NSFD(sf) = nsfd
+	SF_SFDS(sf) = sfds
+	SF_NSTARS(sf) = nstars
+	SF_STARS(sf) = stars
+	SF_NFOCUS(sf) = nfocus
+	SF_FOCUS(sf) = focus
+	SF_NIMAGES(sf) = nimages
+	SF_IMAGES(sf) = images
+
+	# Find the average and best focus values.  Sort the focus groups
+	# by magnitude and the star groups by focus.
+
+	call stf_fitfocus (sf)
+	do i = 1, SF_NFOCUS(sf) {
+	    sff = SF_SFF(sf,i)
+	    call qsort (SFF_SFD(sff,1), SFF_NSFD(sff), stf_magsort)
+	}
+	do i = 1, SF_NSTARS(sf) {
+	    sfs = SF_SFS(sf,i)
+	    call qsort (SFS_SFD(sfs,1), SFS_NSFD(sfs), stf_focsort)
+	}
+
+	call sfree (sp)
+end
+
+
+# STF_LOG -- Print log of results
+
+procedure stf_log (sf, fd)
+
+pointer	sf			#I Main data structure
+int	fd			#I File descriptor
+
+int	i, j, n
+pointer	sp, str, sfd, sfs, sff, sfi
+
+begin
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	# Print banner and overall result.
+	call sysid (Memc[str], SZ_LINE)
+	call fprintf (fd, "%s\n\n")
+	    call pargstr (Memc[str])
+
+	# Print each individual object organized by image.
+	call fprintf (fd, "%15.15s %7s %7s %7s")
+	    call pargstr ("Image")
+	    call pargstr ("Column")
+	    call pargstr ("Line")
+	    call pargstr ("Mag")
+	if (IS_INDEF(SF_F(sf))) {
+	    call fprintf (fd, " %7s")
+		call pargstr (SF_WTYPE(sf))
+	} else {
+	    call fprintf (fd, " %7s %7s")
+		call pargstr ("Focus")
+		call pargstr (SF_WTYPE(sf))
+	}
+	if (SF_WCODE(sf) == 4) {
+	    call fprintf (fd, " %4s")
+		call pargstr ("Beta")
+	}
+	call fprintf (fd, " %7s %7s %3s\n")
+	    call pargstr ("Ellip")
+	    call pargstr ("PA")
+	    call pargstr ("SAT")
+
+	do i = 1, SF_NIMAGES(sf) {
+	    sfi = SF_SFI(sf,i)
+	    n = 0
+	    do j = 1, SFI_NSFD(sfi) {
+		sfd = SFI_SFD(sfi,j)
+		if (SFD_STATUS(sfd) != 0)
+		    next
+		if (n == 0) {
+		    call fprintf (fd, "%15.15s")
+			call pargstr (SFD_IMAGE(sfd))
+		} else
+		    call fprintf (fd, "%15w")
+		call fprintf (fd, " %7.2f %7.2f %7.2f")
+		    call pargr (SFD_X(sfd))
+		    call pargr (SFD_Y(sfd))
+		    call pargr (-2.5*log10 (SFD_M(sfd) / SF_M(sf)))
+		if (IS_INDEF(SFD_F(sfd))) {
+		    call fprintf (fd,
+			" %7.3f")
+			call pargr (SFD_W(sfd))
+		} else {
+		    call fprintf (fd, " %7.6g %7.3f")
+			call pargr (SFD_F(sfd))
+			call pargr (SFD_W(sfd))
+		}
+		if (SF_WCODE(sf) == 4) {
+		    call fprintf (fd, " %4.1f")
+			call pargr (SFD_BETA(sfd))
+		}
+		call fprintf (fd, " %7.2f %7d")
+		    call pargr (SFD_E(sfd))
+		    call pargr (SFD_PA(sfd))
+		if (SFD_NSAT(sfd) == 0)
+		    call fprintf (fd, "\n")
+		else
+		    call fprintf (fd, "  *\n")
+		n = n + 1
+	    }
+	}
+	if (n > 0)
+	    call fprintf (fd, "\n")
+
+	# Print best values by star.
+	if (SF_NS(sf) > 1) {
+	    n = 0
+	    do i = 1, SF_NSTARS(sf) {
+		sfs = SF_SFS(sf,i)
+		if (SFS_NF(sfs) > 1 || SFS_N(sfs) > 1) {
+		    call stf_title (sf, NULL, sfs, NULL, Memc[str], SZ_LINE)
+		    call fprintf (fd, "  %s\n")
+			call pargstr (Memc[str])
+		    n = n + 1
+		}
+	    }
+	    if (n > 0)
+		call fprintf (fd, "\n")
+	}
+
+	# Print averages at each focus.
+	if (SF_NF(sf) > 1) {
+	    n = 0
+	    do i = 1, SF_NFOCUS(sf) {
+		sff = SF_SFF(sf,i)
+		if (SFF_N(sff) > 1) {
+		    call stf_title (sf, NULL, NULL, sff, Memc[str], SZ_LINE)
+		    call fprintf (fd, "  %s\n")
+			call pargstr (Memc[str])
+		    n = n + 1
+		}
+	    }
+	    if (n > 0)
+		call fprintf (fd, "\n")
+	}
+
+	call stf_title (sf, NULL, NULL, NULL, Memc[str], SZ_LINE)
+	call fprintf (fd, "%s\n")
+	    call pargstr (Memc[str])
+end
+
+
+# STF_TITLE -- Return result title string.
+# The title is dependent on whether an overall title, a title for a star
+# group, for a focus group, or for an indivdual object is desired.
+# The title also is adjusted for the select size type and the number
+# of objects in a group.
+
+procedure stf_title (sf, sfd, sfs, sff, title, sz_title)
+
+pointer	sf			#I Starfocus pointer
+pointer	sfd			#I Data pointer
+pointer	sfs			#I Star pointer
+pointer	sff			#I Focus pointer
+char	title[sz_title]		#O Title string 
+int	sz_title		#I Size of title string 
+
+pointer	ptr
+int	i, fd, stropen()
+errchk	stropen
+
+begin
+	fd = stropen (title, sz_title, WRITE_ONLY)
+
+	if (sfd != NULL) {
+	    call fprintf (fd, "%s @ (%.2f, %.2f):")
+		call pargstr (SFD_IMAGE(sfd))
+		call pargr (SFD_X(sfd))
+		call pargr (SFD_Y(sfd))
+	    switch (SF_WCODE(sf)) {
+	    case 4:
+		call fprintf (fd, " %s=%.2f (%3.1f), e=%.2f, pa=%d")
+		    call pargstr (SF_WTYPE(sf))
+		    call pargr (SFD_W(sfd))
+		    call pargr (SFD_BETA(sfd))
+		    call pargr (SFD_E(sfd))
+		    call pargr (SFD_PA(sfd))
+	    default:
+		call fprintf (fd, " %s=%.2f, e=%.2f, pa=%d")
+		    call pargstr (SF_WTYPE(sf))
+		    call pargr (SFD_W(sfd))
+		    call pargr (SFD_E(sfd))
+		    call pargr (SFD_PA(sfd))
+	    }
+	    if (SFD_SFS(sfd) != NULL) {
+		if (SFS_M(SFD_SFS(sfd)) != SF_M(sf)) {
+		    call fprintf (fd, " , m=%.2f")
+			call pargr (-2.5*log10 (SFS_M(SFD_SFS(sfd)) / SF_M(sf)))
+		}
+	    }
+	    if (!IS_INDEF(SFD_F(sfd))) {
+		call fprintf (fd, ", f=%.4g")
+		    call pargr (SFD_F(sfd))
+	    }
+	} else if (sfs != NULL) {
+	    ptr = SFS_SFD(sfs,1)
+	    call fprintf (fd, "%s")
+	    if (SFS_NF(sfs) > 1)
+		call pargstr ("Best focus estimate")
+	    else if (SFS_N(sfs) > 1)
+		call pargstr ("Average star")
+	    else {
+		for (i=1; SFD_STATUS(SFS_SFD(sfs,i))!=0; i=i+1)
+		    ;
+		call pargstr (SFD_IMAGE(SFS_SFD(sfs,i)))
+	    }
+	    call fprintf (fd, " @ (%.2f, %.2f): %s=%.2f")
+		    call pargr (SFD_X(ptr))
+		    call pargr (SFD_Y(ptr))
+		    call pargstr (SF_WTYPE(sf))
+		    call pargr (SFS_W(sfs))
+	    #if (SFS_M(sfs) != SF_M(sf)) {
+		call fprintf (fd, ", m=%.2f")
+		    call pargr (-2.5 * log10 (SFS_M(sfs) / SF_M(sf)))
+	    #}
+	    if (!IS_INDEF(SFS_F(sfs))) {
+		call fprintf (fd, ", f=%.4g")
+		    call pargr (SFS_F(sfs))
+	    }
+	} else if (sff != NULL) {
+	    if (SFF_NI(sff) == 1) {
+		for (i=1; SFD_STATUS(SFF_SFD(sff,i))!=0; i=i+1)
+		    ;
+		call fprintf (fd, "%s")
+		    call pargstr (SFD_IMAGE(SFF_SFD(sff,i)))
+		if (!IS_INDEF(SFF_F(sff))) {
+		    call fprintf (fd, " at focus %.4g")
+			call pargr (SFF_F(sff))
+		}
+		call fprintf (fd, " with average")
+	    } else {
+		if (IS_INDEF(SFF_F(sff)))
+		    call fprintf (fd, "Average")
+		else {
+		    call fprintf (fd, "Focus %.4g with average")
+			call pargr (SFF_F(sff))
+		}
+	    }
+	    call fprintf (fd, " %s of %.2f")
+		call pargstr (SF_WTYPE(sf))
+		call pargr (SFF_W(sff))
+	} else {
+	    if (IS_INDEF(SF_F(sf))) {
+		if (SF_WCODE(sf) == 1) {
+		    call fprintf (fd, "  %s%d%% enclosed flux radius of ")
+			if (SF_N(sf) > 1)
+			    call pargstr ("Average ")
+			else
+			    call pargstr ("")
+			call pargr (100 * SF_LEVEL(sf))
+		} else {
+		    if (SF_N(sf) > 1)
+			call fprintf (fd,
+			    "  Average full width at half maximum (%s) of ")
+		    else
+			call fprintf (fd,
+			    "  Full width at half maximum (%s) of ")
+		    call pargstr (SF_WTYPE(sf))
+		}
+		call fprintf (fd, "%.4f")
+		    call pargr (SF_W(sf))
+	    } else {
+		call fprintf (fd, "  %s of %.6g with ")
+		    if (SF_NS(sf) > 1) {
+			if (SF_NF(sf) > 1)
+			    call pargstr ("Average best focus")
+			else
+			    call pargstr ("Average focus")
+		    } else {
+			if (SF_NF(sf) > 1)
+			    call pargstr ("Best focus")
+			else
+			    call pargstr ("Focus")
+		    }
+		    call pargr (SF_F(sf))
+		if (SF_WCODE(sf) == 1) {
+		    call fprintf (fd, "%d%% enclosed flux radius of ")
+			call pargr (100 * SF_LEVEL(sf))
+		} else {
+		    call fprintf (fd, "%s of ")
+			call pargstr (SF_WTYPE(sf))
+		}
+		call fprintf (fd, "%.2f")
+		    call pargr (SF_W(sf))
+	    }
+	}
+
+	call strclose (fd)
+end
+
+
+# STF_FITFOCUS -- Find the best focus.
+
+procedure stf_fitfocus (sf)
+
+pointer	sf			#I Starfocus pointer
+
+int	i, j, k, n, jmin
+pointer	x, y, sfd, sfs, sff, sfi
+real	f, r, m, wr, wf
+bool	fp_equalr()
+
+begin
+	# Set number of valid points, stars, focuses, images.
+	SF_N(sf) = 0
+	SF_YP1(sf) = 0
+	SF_YP2(sf) = 0
+	do i = 1, SF_NSFD(sf) {
+	    sfd = SF_SFD(sf,i)
+	    if (SFD_STATUS(sfd) == 0) {
+		SF_N(sf) = SF_N(sf) + 1
+		SF_YP1(sf) = min (SF_YP1(sf), SFD_YP1(sfd))
+		SF_YP2(sf) = max (SF_YP2(sf), SFD_YP2(sfd))
+	    }
+	}
+	SF_NS(sf) = 0
+	do i = 1, SF_NSTARS(sf) {
+	    sfs = SF_SFS(sf,i)
+	    SFS_N(sfs) = 0
+	    SFS_M(sfs) = 0.
+	    SFS_NF(sfs) = 0
+	    do j = 1, SFS_NSFD(sfs) {
+		sfd = SFS_SFD(sfs,j)
+		if (SFD_STATUS(SFS_SFD(sfs,j)) != 0)
+		    next
+		SFS_N(sfs) = SFS_N(sfs) + 1
+		SFS_M(sfs) = SFS_M(sfs) + SFD_M(sfd)
+		sff = SFD_SFF(sfd)
+		for (k = 1; SFD_SFF(SFS_SFD(sfs,k)) != sff; k = k + 1)
+		    ;
+		if (k == j)
+		    SFS_NF(sfs) = SFS_NF(sfs) + 1
+	    }
+	    if (SFS_N(sfs) > 0) {
+		SFS_M(sfs) = SFS_M(sfs) / SFS_N(sfs)
+		SF_NS(sf) = SF_NS(sf) + 1
+	    }
+	}
+	SF_NF(sf) = 0
+	do i = 1, SF_NFOCUS(sf) {
+	    sff = SF_SFF(sf,i)
+	    SFF_W(sff) = 0.
+	    SFF_N(sff) = 0
+	    SFF_NI(sff) = 0
+	    wr = 0
+	    do j = 1, SFF_NSFD(sff) {
+		sfd = SFF_SFD(sff,j)
+		if (SFD_STATUS(sfd) != 0)
+		    next
+		m = SFS_M(SFD_SFS(sfd))
+		wr = wr + m
+		SFF_W(sff) = SFF_W(sff) + m * SFD_W(sfd)
+		SFF_N(sff) = SFF_N(sff) + 1
+		sfi = SFD_SFI(sfd)
+		for (k = 1; SFD_SFI(SFF_SFD(sff,k)) != sfi; k = k + 1)
+		    ;
+		if (k == j)
+		    SFF_NI(sff) = SFF_NI(sff) + 1
+	    }
+	    if (SFF_N(sff) > 0) {
+		SFF_W(sff) = SFF_W(sff) / wr
+		SF_NF(sf) = SF_NF(sf) + 1
+	    }
+	}
+	SF_NI(sf) = 0
+	do i = 1, SF_NIMAGES(sf) {
+	    sfi = SF_SFI(sf,i)
+	    SFI_N(sfi) = 0
+	    do j = 1, SFI_NSFD(sfi)
+		if (SFD_STATUS(SFI_SFD(sfi,j)) == 0)
+		    SFI_N(sfi) = SFI_N(sfi) + 1
+	    if (SFI_N(sfi) > 0)
+		SF_NI(sf) = SF_NI(sf) + 1
+	}
+
+	# Find the average magnitude, best focus, and radius for each star.
+	# Find the brightest magnitude and average best focus and radius
+	# over all stars.
+
+	SF_BEST(sf) = SF_SFD(sf,1)
+	SF_F(sf) = 0.
+	SF_W(sf) = 0.
+	SF_M(sf) = 0.
+	SF_NS(sf) = 0
+	wr = 0.
+	wf = 0.
+	do i = 1, SF_NSTARS(sf) {
+	    sfs = SF_SFS(sf,i)
+	    call malloc (x, SFS_NSFD(sfs), TY_REAL)
+	    call malloc (y, SFS_NSFD(sfs), TY_REAL)
+	    k = 0
+	    n = 0
+	    do j = 1, SFS_NSFD(sfs) {
+		sfd = SFS_SFD(sfs,j)
+		if (SFD_STATUS(sfd) != 0)
+		    next
+		r = SFD_W(sfd)
+		f = SFD_F(sfd)
+		if (!IS_INDEF(f))
+		    k = k + 1
+		Memr[x+n] = f
+		Memr[y+n] = r
+		n = n + 1
+		if (r < SFD_W(SF_BEST(sf)))
+		    SF_BEST(sf) = sfd
+	    }
+
+	    # Find the best focus and radius.
+	    if (n == 0) {
+		SFS_F(sfs) = INDEF
+		SFS_W(sfs) = INDEF
+		SFS_M(sfs) = INDEF
+		SFS_N(sfs) = 0
+	    } else if (k == 0) {
+		call alimr (Memr[y], n, f, r)
+		f = INDEF
+		m = SFS_M(sfs)
+		wr = wr + m
+		SFS_F(sfs) = f
+		SFS_W(sfs) = r
+		SFS_M(sfs) = m
+		SFS_N(sfs) = n
+		SF_W(sf) = SF_W(sf) + m * r
+		SF_M(sf) = max (SF_M(sf), m)
+		SF_NS(sf) = SF_NS(sf) + 1
+	    } else {
+		SFS_N(sfs) = n
+		if (k < n) {
+		    k = 0
+		    do j = 0, n-1 {
+			if (!IS_INDEF(Memr[x+j])) {
+			    Memr[x+k] = Memr[x+j]
+			    Memr[y+k] = Memr[y+j]
+			    k = k + 1
+			}
+		    }
+		}
+		call xt_sort2 (Memr[x], Memr[y], k)
+		n = 0
+		do j = 1, k-1 {
+		    if (fp_equalr (Memr[x+j], Memr[x+n])) {
+			if (Memr[y+j] < Memr[y+n])
+			    Memr[y+n] = Memr[y+j]
+		    } else {
+			n = n + 1
+			Memr[x+n] = Memr[x+j]
+			Memr[y+n] = Memr[y+j]
+		    }
+		}
+		n = n + 1
+
+		# Find the minimum radius
+		jmin = 0
+		do j = 0, n-1
+		    if (Memr[y+j] < Memr[y+jmin])
+			jmin = j
+
+		# Use parabolic interpolation to find the best focus
+		if (jmin == 0 || jmin == n-1) {
+		    f = Memr[x+jmin]
+		    r = Memr[y+jmin]
+		} else
+		    call stf_parab (Memr[x+jmin-1], Memr[y+jmin-1], f, r)
+
+		m = SFS_M(sfs)
+		wr = wr + m
+		wf = wf + m
+		SFS_F(sfs) = f
+		SFS_W(sfs) = r
+		SFS_M(sfs) = m
+		SF_F(sf) = SF_F(sf) + m * f
+		SF_W(sf) = SF_W(sf) + m * r
+		SF_M(sf) = max (SF_M(sf), m)
+		SF_NS(sf) = SF_NS(sf) + 1
+	    }
+	    call mfree (x, TY_REAL)
+	    call mfree (y, TY_REAL)
+	}
+
+	if (wr > 0.)
+	    SF_W(sf) = SF_W(sf) / wr
+	else {
+	    SF_W(sf) = INDEF
+	    SF_M(sf) = INDEF
+	}
+	if (wf > 0.)
+	    SF_F(sf) = SF_F(sf) / wf
+	else
+	    SF_F(sf) = INDEF
+end
+
+
+# STF_PARAB -- Find the minimum of a parabolic fit to three points.
+
+procedure stf_parab (x, y, xmin, ymin)
+
+real	x[3]
+real	y[3]
+real	xmin
+real	ymin
+
+double	x12, x13, x23, x213, x223, y13, y23, a, b, c
+
+begin
+	x12 = x[1] - x[2]
+	x13 = x[1] - x[3]
+	x23 = x[2] - x[3]
+	x213 = x13 * x13
+	x223 = x23 * x23
+	y13 = y[1] - y[3]
+	y23 = y[2] - y[3]
+	c = (y13 - y23 * x13 / x23) / (x213 - x223 * x13 / x23)
+	b = (y23 - c * x223) / x23
+	a = y[3]
+	xmin = -b / (2 * c)
+	ymin = a + b * xmin + c * xmin * xmin
+	xmin = xmin + x[3]
+end
+
+
+# STF_MAGSORT -- Compare two star structures by average magnitude.
+
+int procedure stf_magsort (sfd1, sfd2)
+
+pointer	sfd1, sfd2		# Structures to compare
+pointer	sfs1, sfs2		# Star structures for magnitudes
+
+begin
+	sfs1 = SFD_SFS(sfd1)
+	sfs2 = SFD_SFS(sfd2)
+	if (SFS_M(sfs1) > SFS_M(sfs2))
+	    return (-1)
+	else if (SFS_M(sfs1) < SFS_M(sfs2))
+	    return (1)
+	else
+	    return (0)
+end
+
+
+# STF_FOCSORT -- Compare two star structures by focus.
+
+int procedure stf_focsort (sfd1, sfd2)
+
+pointer	sfd1, sfd2		# Structures to compare
+
+begin
+	if (SFD_F(sfd1) < SFD_F(sfd2))
+	    return (-1)
+	else if (SFD_F(sfd1) > SFD_F(sfd2))
+	    return (1)
+	else
+	    return (0)
+end
+
+
+# STF_DISPLAY -- Display image if necessary.
+# The user is required to display the first image separately.
+
+procedure stf_display (image, frame)
+
+char	image[ARB]		#I Image to display
+int	frame			#I Display frame to use
+
+int	i, status
+pointer	sp, dname, ds, iw, imd_mapframe(), iw_open()
+bool	xt_imnameeq()
+errchk	clcmdw
+
+begin
+	call smark (sp)
+	call salloc (dname, SZ_LINE, TY_CHAR)
+
+	ds = imd_mapframe (1, READ_WRITE, NO)
+	do i = 1, MAX_FRAMES {
+	    iferr (iw = iw_open (ds, i, Memc[dname], SZ_LINE, status))
+		next
+	    call iw_close (iw)
+	    if (xt_imnameeq (image, Memc[dname]))
+		break
+	}
+	call imunmap (ds)
+
+	if (!xt_imnameeq (image, Memc[dname])) {
+	    call sprintf (Memc[dname], SZ_LINE, "display %s frame=%d fill=yes")
+		call pargstr (image)
+		call pargi (frame)
+	    call clcmdw (Memc[dname])
+	}
+
+	call sfree (sp)
+end
+
+
+# STFCUR -- Debugging routine.
+# Replace calls to clgcur with stfcur so that an text file containing the
+# cursor coordinates may be specified when running standalone (such as
+# under a debugger).
+
+int procedure stfcur (cur, wx, wy, wcs, key, cmd, sz_cmd)
+
+char	cur[ARB]		# Cursor name
+real	wx, wy			# Cursor coordinate
+int	wcs			# WCS
+int	key			# Key
+char	cmd[sz_cmd]		# Command
+int	sz_cmd			# Size of command
+
+int	fd, stat, open(), fscan()
+pointer	fname
+errchk	open
+
+begin
+	if (fd == NULL) {
+	    call malloc (fname, SZ_FNAME, TY_CHAR) 
+	    call clgstr (cur, Memc[fname], SZ_FNAME)
+	    fd = open (Memc[fname], READ_ONLY, TEXT_FILE)
+	    call mfree (fname, TY_CHAR)
+	}
+
+	stat = fscan (fd)
+	if (stat == EOF) {
+	    call close (fd)
+	    return (stat)
+	}
+
+	call gargr (wx)
+	call gargr (wy)
+	call gargi (wcs)
+	call gargi (key)
+
+	return (stat)
+end
diff --git a/noao/obsutil/src/starfocus/x_starfocus.x b/noao/obsutil/src/starfocus/x_starfocus.x
new file mode 100644
index 00000000..a3642bae
--- /dev/null
+++ b/noao/obsutil/src/starfocus/x_starfocus.x
@@ -0,0 +1,2 @@
+task    psfmeasure	= t_psfmeasure,
+	starfocus	= t_starfocus