Initial commit

8 files changed, 2620 insertions, 0 deletions

diff --git a/noao/obsutil/src/specfocus/Revisions b/noao/obsutil/src/specfocus/Revisions
new file mode 100644
index 00000000..3104b9f7
--- /dev/null
+++ b/noao/obsutil/src/specfocus/Revisions
@@ -0,0 +1,9 @@
+.help revisions Nov01 obsutil
+.nf
+spfgraph.x
+    Fixed case where if the average of the minimum and maximum focus values
+    is zero then a floating divide by zero would occur.  (2/14/97, Valdes)
+
+t_specfocus.x
+    Fixed bug in interpreting the focus parameter.  (6/26/95, Valdes)
+.endhelp
diff --git a/noao/obsutil/src/specfocus/mkpkg b/noao/obsutil/src/specfocus/mkpkg
new file mode 100644
index 00000000..c41d9a34
--- /dev/null
+++ b/noao/obsutil/src/specfocus/mkpkg
@@ -0,0 +1,19 @@
+# Make the SPECFOCUS task.
+
+$checkout libpkg.a ../
+$update   libpkg.a
+$checkin  libpkg.a ../
+$exit
+
+standalone:
+	$update	libpkg.a
+	$omake	x_specfocus.x
+	$link	x_specfocus.o libpkg.a -lxtools -lcurfit -liminterp\
+	    -o xx_specfocus.e
+	;
+
+libpkg.a:
+	spfgraph.x	specfocus.h <gset.h> <mach.h>
+	t_specfocus.x	specfocus.h <error.h> <imhdr.h> <mach.h> <math.h>\
+			<math/curfit.h> <math/iminterp.h>
+	;
diff --git a/noao/obsutil/src/specfocus/specfocus.h b/noao/obsutil/src/specfocus/specfocus.h
new file mode 100644
index 00000000..d79f7c6e
--- /dev/null
+++ b/noao/obsutil/src/specfocus/specfocus.h
@@ -0,0 +1,33 @@
+# Data structures for SPECFOCUS
+
+define	SZ_SFFNAME	79		  # Length of file names
+define	LEN_SF		54		  # Length of image data structure
+
+define	SF_IMAGE	Memc[P2C($1)]	  # Image name
+define	SF_FOCUS	Memr[P2R($1+40)]  # Focus
+define	SF_WIDTH	Memr[P2R($1+41)]  # Width
+define	SF_LEVEL	Memr[P2R($1+42)]  # Level of width
+define	SF_AXIS		Memi[$1+43]	  # Dispersion axis
+define	SF_X1		Memi[$1+44]	  # Start of dispersion sampling
+define	SF_DX		Memi[$1+45]	  # Dispersion sampling step
+define	SF_NX		Memi[$1+46]	  # Number of dispersion samples
+define	SF_Y1		Memi[$1+47]	  # Start of cross-dispersion sampling
+define	SF_DY		Memi[$1+48]	  # Cross-dispersion sampling step
+define	SF_NY		Memi[$1+49]	  # Number of cross-dispersion samples
+define	SF_SFD		Memi[$1+50]	  # Pointer to data structures
+define	SF_NSFD		Memi[$1+51]	  # Number of data structures
+define	SF_DATA		Memi[$1+52]	  # Pointer to spectrum data 
+define	SF_NPIX		Memi[$1+53]	  # Number of pixels per spectrum
+
+define	LEN_SFD		8		  # Length of spectrum data structure
+
+define	SF_X		Memr[P2R($1)]	  # Dispersion axis coordinate
+define	SF_Y		Memr[P2R($1+1)]	  # Spatial axis coordinate
+define	SF_SPEC		Memi[$1+2]	  # Pointer to spectrum
+define	SF_ASI		Memi[$1+3]	  # Pointer to correlation profile
+define	SF_FOC		Memr[P2R($1+4)]	  # Focus
+define	SF_WID		Memr[P2R($1+5)]	  # Width
+define	SF_POS		Memr[P2R($1+6)]	  # Position
+define	SF_DEL		Memi[$1+7]	  # Deleted?
+
+define	SFD		Memi[SF_SFD($1)+($3-1)*SF_NX($1)+$2-1]
diff --git a/noao/obsutil/src/specfocus/specfocus.par b/noao/obsutil/src/specfocus/specfocus.par
new file mode 100644
index 00000000..bf21416c
--- /dev/null
+++ b/noao/obsutil/src/specfocus/specfocus.par
@@ -0,0 +1,13 @@
+images,s,a,,,,List of images
+focus,s,h,"",,,Focus values
+corwidth,i,h,20,,,Correlation width
+level,r,h,0.5,,,Percent or fraction of peak for width measurement
+shifts,b,h,yes,,,"Compute shifts across the dispersion?
+"
+dispaxis,i,h,2,1,2,Dispersion axis (long slit only)
+nspectra,i,h,1,1,,Number of spectral samples (long slit only)
+ndisp,i,h,1,1,,Number of dispersion samples
+slit1,i,h,INDEF,,,Lower slit edge
+slit2,i,h,INDEF,,,"Upper slit edge
+"
+logfile,s,h,"logfile",,,"Logfile"
diff --git a/noao/obsutil/src/specfocus/spfgraph.x b/noao/obsutil/src/specfocus/spfgraph.x
new file mode 100644
index 00000000..0430d494
--- /dev/null
+++ b/noao/obsutil/src/specfocus/spfgraph.x
@@ -0,0 +1,1637 @@
+include	<gset.h>
+include	<mach.h>
+include	"specfocus.h"
+
+define	HELP		"specfocus$specfocus.key"
+define	PROMPT		"specfocus options"
+
+define	VX1	.15		# Minimum X viewport
+define	VX2	.95		# Maximum X viewport
+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
+
+define	NMAX	5		# Maximum number of samples
+define	HLCOLOR	2		# Highlight color
+define	HLWIDTH	4.		# Highlight width
+
+
+# SPF_GRAPH -- Interactive graphing of results
+
+procedure spf_graph (sfavg, sfbest, sfs, nimages, lag)
+
+pointer	sfavg			# Average image
+pointer	sfbest			# Best image
+pointer	sfs[nimages]		# Images
+int	nimages			# Number of images
+int	lag			# Maximum lag
+
+char	cmd[10]
+real	wx, wy, f, w, r2, r2min, fa[8]
+int	i, j, k, l, i1, j1, nx, ny, nxgrid, nygrid, nsfd
+int	wcs, key, pkey, del, clgcur()
+pointer	sp, sysidstr, title, gp, gopen()
+pointer	sf, sfd, sfcur
+
+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)
+
+	# Set system id label
+	call sysid (Memc[sysidstr], SZ_LINE)
+
+	# Set current image and sample
+	nsfd = SF_NSFD(sfavg)
+	nx = SF_NX(sfavg)
+	ny = SF_NY(sfavg)
+	i = (nx + 1) / 2
+	j = (ny + 1) / 2
+	for (k=1; k<nimages && sfs[k]!=sfbest; k=k+1)
+	    ;
+
+	# Set grid layout
+	if (nimages < 3) {
+	    nxgrid = nimages
+	    nygrid = 1
+	} else {
+	    nxgrid = nint (sqrt (real (nimages)))
+	    if (mod (nimages, nxgrid+1) >= mod (nimages, nxgrid))
+		nxgrid = nxgrid + 1
+	    nygrid = (nimages-1) / nxgrid + 1
+	}
+
+	# Open graphics and enter interactive graphics loop
+	gp = gopen ("stdgraph", NEW_FILE, STDGRAPH)
+	key = 'b'
+	wcs = 0
+	repeat {
+	    # Verify keys, check for '?', and 'q', set redraw.
+	    switch (key) {
+	    case '?':
+		call gpagefile (gp, HELP, PROMPT)
+		next
+	    case 'b', 'p', 's', 'w', 'z':
+		pkey = key
+		del = NO
+	    case 'q':
+		break
+	    case 'r':
+		key = pkey
+		del = NO
+	    case ' ', 'd':
+		del = NO
+	    case 'u':
+		del = YES
+	    default:
+		call printf ("\007")
+		next
+	    }
+
+	    # Map the cursor position to an image and sample.
+	    switch (wcs) {
+	    case 1:
+		k = 1
+		call spf_sample (sfavg, 1, del, wx, wy, i, j, k)
+		f = SF_FOC(SFD(sfavg,i,j))
+		r2min = MAX_REAL
+		do l = 1, nimages {
+		    sf = sfs[l]
+		    sfd = SFD(sf,i,j)
+		    if (SF_DEL(sfd) == del) {
+			r2 = abs (f - SF_FOC(sfd))
+			if (r2 < r2min) {
+			    r2min = r2
+			    k = l
+			}
+		    }
+		}
+		call spf_sample (sfs, nimages, del, wx, wy, i, j, k)
+	    case 2, 6:
+		call spf_sample (sfs, nimages, del, wx, wy, i, j, k)
+	    case 3:
+		r2min = MAX_REAL
+		call gctran (gp, wx, wy, wx, wy, wcs, 0)
+		do l = 1, nimages {
+		    sf = sfs[l]
+		    do j1 = 1, ny {
+			do i1 = 1, nx {
+			    sfd = SFD(sf,i1,j1)
+			    if (SF_DEL(sfd) == del) {
+				f = SF_FOC(sfd)
+				w = SF_WID(sfd)
+				call gctran (gp, f, w, f, w, wcs, 0)
+				r2 = (f-wx)**2 + (w-wy)**2
+				if (r2 < r2min) {
+				    r2min = r2
+				    i = i1
+				    j = j1
+				    k = l
+				}
+			    }
+			}
+		    }
+		}
+	    case 4, 5, 8:
+		i1 = max (1, min (nxgrid, nint(wx)))
+		j1 = max (1, min (nygrid, nint(wy)))
+		k = max (1, min (nimages, (j1-1) * nxgrid + i1))
+		call spf_sample (sfs, nimages, del, real(i), real(j), i, j, k)
+		if (wcs == 8) {
+		    wx = nx * (wx - nint (wx) + 0.5) + 0.5
+		    wy = -(ny+1) * (wy - nint (wy) + 0.5) + (ny+1.5)
+		    call spf_sample (sfs, nimages, del, wx, wy, i, j, k)
+		}
+	    }
+
+	    # Switch on action key
+	    switch (key) {
+	    case ' ':
+		if (wcs == 1)
+		    sf = sfavg
+		else
+		    sf = sfs[k]
+		sfd = SFD(sf,i,j)
+		call printf (
+		    "Image %s at (%d, %d), Focus = %.3g, Width = %.2f")
+		    call pargstr (SF_IMAGE(sf))
+		    call pargr (SF_X(sfd))
+		    call pargr (SF_Y(sfd))
+		    call pargr (SF_FOC(sfd))
+		    call pargr (SF_WID(sfd))
+		if (abs (SF_POS(sfd)) > .01) {
+		    call printf (", Shift = %.2f")
+			call pargr (SF_POS(sfd))
+		}
+		call printf ("\n")
+		next
+	    case 'd':
+		repeat {
+		    switch (key) {
+		    case 'i':
+			do j1 = 1, ny
+			   do i1 = 1, nx
+				SF_DEL(SFD(sfs[k],i1,j1)) = YES
+			call spf_sample (sfs, nimages, del, real(i), real(j),
+			    i, j, k)
+		    case 's':
+			do l = 1, nimages {
+			    sfd = SFD(sfs[l],i,j)
+			    SF_DEL(sfd) = YES
+			}
+			call spf_sample (sfs, nimages, del, real(i), real(j),
+			    i, j, k)
+		    case 'p':
+			sfd = SFD(sfs[k],i,j)
+			SF_DEL(sfd) = YES
+			call spf_sample (sfs, nimages, del, real(i), real(j),
+			    i, j, k)
+		    default:
+			call printf ("Delete image, sample, or point?")
+			next
+		    }
+		    call spf_fitfocus (sfs, nimages, sfavg, sfbest)
+		    key = pkey
+		    break
+		} until (clgcur ("gcur", wx, wy, wcs, key, cmd, 10) == EOF)
+	    case 'u':
+		repeat {
+		    switch (key) {
+		    case 'i':
+			call spf_sample (sfs, nimages, del, real(i), real(j),
+			    i, j, k)
+			do j1 = 1, ny
+			   do i1 = 1, nx
+				SF_DEL(SFD(sfs[k],i1,j1)) = NO
+		    case 's':
+			call spf_sample (sfs, nimages, del, real(i), real(j),
+			    i, j, k)
+			do l = 1, nimages {
+			    sfd = SFD(sfs[l],i,j)
+			    SF_DEL(sfd) = NO
+			}
+		    case 'p':
+			call spf_sample (sfs, nimages, del, real(i), real(j),
+			    i, j, k)
+			sfd = SFD(sfs[k],i,j)
+			SF_DEL(sfd) = NO
+		    default:
+			call printf ("Undelete image, sample or point?")
+			next
+		    }
+		    call spf_fitfocus (sfs, nimages, sfavg, sfbest)
+		    key = pkey
+		    break
+		} until (clgcur ("gcur", wx, wy, wcs, key, cmd, 10) == EOF)
+	    }
+	    sfcur = sfs[k]
+	    sfd = SFD(sfcur,i,j)
+
+	    # Make the graphs.
+	    call gclear (gp)
+	    call gseti (gp, G_FACOLOR, 0)
+
+	    if (nimages > 1 && nsfd > 1) {
+		switch (key) {
+		case 'p':
+		    call gseti (gp, G_WCS, 4)
+		    call gsview (gp, VX1, VX2, VY3, VY4)
+		    call spf_g4 (gp, sfavg, sfbest, sfcur, sfs, nimages, i, j,
+			nxgrid, nygrid, lag)
+		    if (nx > NMAX || ny > NMAX) {
+			call gseti (gp, G_WCS, 1)
+			call gsview (gp, VX1, VX2, VY1, VY2)
+			call gfill (gp, fa, fa[5], 4, GF_SOLID)
+			call spf_g1 (gp, sfcur, i, j, NO, NO)
+		    } else {
+			call gseti (gp, G_WCS, 2)
+			call gsview (gp, VX1, VX2, VY1, VY2)
+			call gfill (gp, fa, fa[5], 4, GF_SOLID)
+			call spf_g2 (gp, sfavg, sfbest, sfcur, i, j, lag)
+		    }
+		case 's':
+		    call gseti (gp, G_WCS, 5)
+		    call gsview (gp, VX1, VX2, VY3, VY4)
+		    call spf_g5 (gp, sfcur, sfs, nimages, i, j, nxgrid, nygrid)
+		    if (nx > NMAX || ny > NMAX) {
+			call gseti (gp, G_WCS, 1)
+			call gsview (gp, VX1, VX2, VY1, VY2)
+			call gfill (gp, fa, fa[5], 4, GF_SOLID)
+			call spf_g1 (gp, sfcur, i, j, NO, NO)
+		    } else {
+			call gseti (gp, G_WCS, 6)
+			call gsview (gp, VX1, VX2, VY1, VY2)
+			call gfill (gp, fa, fa[5], 4, GF_SOLID)
+			call spf_g6 (gp, sfcur, i, j)
+		    }
+		case 'w':
+		    call gseti (gp, G_WCS, 3)
+		    call gsview (gp, VX1, VX2, VY3, VY4)
+		    call spf_g3 (gp, sfavg, sfcur, sfs, nimages, i, j)
+		    if (nx > NMAX || ny > NMAX) {
+			call gseti (gp, G_WCS, 1)
+			call gsview (gp, VX1, VX2, VY1, VY2)
+			call gfill (gp, fa, fa[5], 4, GF_SOLID)
+			call spf_g1 (gp, sfcur, i, j, NO, NO)
+		    } else {
+			call gseti (gp, G_WCS, 8)
+			call gsview (gp, VX1, VX2, VY1, VY2)
+			call gfill (gp, fa, fa[5], 4, GF_SOLID)
+			call spf_g8 (gp, sfavg, sfcur, sfd, sfs, nimages,
+			    nxgrid, nygrid)
+		    }
+		case 'z':
+		    call gseti (gp, G_WCS, 7)
+		    call gsview (gp, VX1, VX2, VY3, VY4)
+		    call spf_g7 (gp, sfbest, sfcur, i, j, lag)
+		    call gseti (gp, G_WCS, 9)
+		    call gsview (gp, VX1, VX2, VY1, VY2)
+		    call gfill (gp, fa, fa[5], 4, GF_SOLID)
+		    call spf_g9 (gp, sfcur, i, j)
+		default:
+		    call gseti (gp, G_WCS, 1)
+		    call gsview (gp, VX1, VX2, VY1, VY4)
+		    call spf_g1 (gp, sfavg, i, j, YES, YES)
+		}
+	    } else if (nimages > 1) {
+		switch (key) {
+		case 'z':
+		    call gseti (gp, G_WCS, 7)
+		    call gsview (gp, VX1, VX2, VY3, VY4)
+		    call spf_g7 (gp, sfbest, sfcur, i, j, lag)
+		    call gseti (gp, G_WCS, 9)
+		    call gsview (gp, VX1, VX2, VY1, VY2)
+		    call gfill (gp, fa, fa[5], 4, GF_SOLID)
+		    call spf_g9 (gp, sfcur, i, j)
+		case 's':
+		    call gseti (gp, G_WCS, 3)
+		    call gsview (gp, VX1, VX2, VY3, VY4)
+		    call spf_g3 (gp, sfavg, sfcur, sfs, nimages, i, j)
+		    call gseti (gp, G_WCS, 5)
+		    call gsview (gp, VX1, VX2, VY1, VY2)
+		    call gfill (gp, fa, fa[5], 4, GF_SOLID)
+		    call spf_g5 (gp, sfcur, sfs, nimages, i, j, nxgrid, nygrid)
+		default:
+		    call gseti (gp, G_WCS, 3)
+		    call gsview (gp, VX1, VX2, VY3, VY4)
+		    call spf_g3 (gp, sfavg, sfcur, sfs, nimages, i, j)
+		    call gseti (gp, G_WCS, 4)
+		    call gsview (gp, VX1, VX2, VY1, VY2)
+		    call spf_g4 (gp, sfavg, sfbest, sfcur, sfs, nimages, i, j,
+			nxgrid, nygrid, lag)
+		}
+	    } else if (nsfd > 1) {
+		switch (key) {
+		case 'z':
+		    call gseti (gp, G_WCS, 7)
+		    call gsview (gp, VX1, VX2, VY3, VY4)
+		    call spf_g7 (gp, sfbest, sfcur, i, j, lag)
+		    call gseti (gp, G_WCS, 9)
+		    call gsview (gp, VX1, VX2, VY1, VY2)
+		    call gfill (gp, fa, fa[5], 4, GF_SOLID)
+		    call spf_g9 (gp, sfcur, i, j)
+		default:
+		    call gseti (gp, G_WCS, 2)
+		    call gsview (gp, VX1, VX2, VY3, VY4)
+		    call spf_g2 (gp, sfavg, sfbest, sfcur, i, j, lag)
+		    call gseti (gp, G_WCS, 6)
+		    call gsview (gp, VX1, VX2, VY1, VY2)
+		    call gfill (gp, fa, fa[5], 4, GF_SOLID)
+		    call spf_g6 (gp, sfcur, i, j)
+		}
+	    } else {
+		call gseti (gp, G_WCS, 7)
+		call gsview (gp, VX1, VX2, VY3, VY4)
+		call spf_g7 (gp, sfbest, sfcur, i, j, lag)
+		call gseti (gp, G_WCS, 9)
+		call gsview (gp, VX1, VX2, VY1, VY2)
+		call gfill (gp, fa, fa[5], 4, GF_SOLID)
+		call spf_g9 (gp, sfcur, i, j)
+	    }
+
+	    call sprintf (Memc[title], SZ_LINE,
+		"Best Average Focus at %.3g with Width of %.3g at %d%% of Peak")
+		call pargr (SF_FOCUS(sfavg))
+		call pargr (SF_WIDTH(sfavg))
+		call pargr (100 * SF_LEVEL(sfavg))
+	    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")
+
+	    pkey = key
+	} until (clgcur ("gcur", wx, wy, wcs, key, cmd, 10) == EOF)
+
+	call gclose (gp)
+	call sfree (sp)
+end
+
+
+# SPF_G1 -- Best Focus at each sample
+
+procedure spf_g1 (gp, sf, ix, iy, focplot, posplot)
+
+pointer	gp		# GIO pointer
+pointer	sf		# SF pointer
+int	ix, iy		# Sample
+int	focplot		# Focus plot?
+int	posplot		# Position plot?
+
+int	i, j, n, nx, ny
+real	wx1, wx2, wy1, wy2, ww1, ww2, wf1, wf2, wp1, wp2
+real	vx[3,2], vy[3,2], dvx, dvy, fa[8]
+real	x, y, z, last
+pointer	sp, str, sfd
+
+data	fa/0.,1.,1.,0.,0.,0.,1.,1./
+
+begin
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	nx = SF_NX(sf)
+	ny = SF_NY(sf)
+	wx1 = SF_X1(sf)
+	wy1 = SF_Y1(sf)
+	wx2 = wx1 + (nx + .3) * SF_DX(sf)
+	wy2 = wy1 + ny * SF_DY(sf)
+
+	# Determine the range of WID, FOC, and POS.
+	ww1 = MAX_REAL
+	wf1 = MAX_REAL
+	wp1 = MAX_REAL
+	ww2 = -MAX_REAL
+	wf2 = -MAX_REAL
+	wp2 = -MAX_REAL
+	do j = 1, ny {
+	    do i = 1, nx {
+		sfd = SFD(sf,i,j)
+		if (SF_DEL(sfd) == NO) {
+		    ww1 = min (ww1, SF_WID(sfd))
+		    wf1 = min (wf1, SF_FOC(sfd))
+		    wp1 = min (wp1, SF_POS(sfd))
+		    ww2 = max (ww2, SF_WID(sfd))
+		    wf2 = max (wf2, SF_FOC(sfd))
+		    wp2 = max (wp2, SF_POS(sfd))
+		}
+	    }
+	}
+	ww2 = max (2., ww2 - ww1)
+
+	# 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,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
+	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
+
+	z = abs ((wf1 + wf2) / 2.)
+	if (z == 0.)
+	    z = max (abs(wf1), abs(wf2))
+	if (z == 0.)
+	    z = 1.
+	if (focplot == NO || abs (wf2 - wf1) / z <= 0.01) {
+	    vx[2,1] = vx[1,1]
+	    vy[2,1] = vy[1,1]
+	}
+	if (posplot == NO || abs (wp2 - wp1) <= .05) {
+	    vx[2,2] = vx[3,2]
+	    vy[2,2] = vy[3,2]
+	}
+	if (nx == 1) {
+	    vy[2,1] = vy[1,1]
+	    vy[2,2] = vy[3,2]
+	}
+	if (ny == 1) {
+	    vx[2,1] = vx[1,1]
+	    vx[2,2] = vx[3,2]
+	}
+
+	call gseti (gp, G_DRAWAXES, 3)
+	call gseti (gp, G_DRAWTICKS, YES)
+
+	# FOC plot
+	if (focplot == YES && abs (wf2 - wf1) / z > 0.01) {
+	    z = wf2 - wf1
+	    wf1 = wf1 - 0.05 * z
+	    wf2 = wf2 + 0.05 * z
+
+	    if (nx > 1) {
+		call gseti (gp, G_LABELTICKS, YES)
+		call gseti (gp, G_NMAJOR, 6)
+		call gseti (gp, G_NMINOR, 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, wx1, wx2, wf1, wf2)
+		if (SF_AXIS(sf) == 1)
+		    call glabax (gp, "", "Column (Dispersion)", "Focus")
+		else
+		    call glabax (gp, "", "Column (Slit)", "Focus")
+
+		call gswind (gp, 0.5, 0.8+nx, wf1, wf2)
+		last = INDEF
+		do i = 1, nx {
+		    x = i
+		    do j = 1, ny {
+			sfd = SFD(sf,i,j)
+			if (SF_DEL(sfd) == NO) {
+			    z = SF_WID(sfd)
+			    z = 0.008 * (1 + (z - ww1) / ww2 * 3)
+			    call gmark (gp, x, SF_FOC(sfd), GM_CIRCLE, z, z)
+			}
+		    }
+
+		    n = 0
+		    z = 0.
+		    do j = 1, ny {
+			sfd = SFD(sf,i,j)
+			if (SF_DEL(sfd) == NO) {
+			    z = z + SF_FOC(SFD(sf,i,j))
+			    n = n + 1
+			}
+		    }
+		    if (n > 0) {
+			if (!IS_INDEF(last))
+			    call gline (gp, last, y, x, z/n)
+			y = z / n
+			last = x
+		    }
+		}
+
+		call gseti (gp, G_PLTYPE, 2)
+		call gline (gp, 0.5, SF_FOCUS(sf), 0.8+nx, SF_FOCUS(sf))
+		call gseti (gp, G_PLTYPE, 1)
+	    }
+
+	    if (ny > 1) {
+		call gseti (gp, G_LABELTICKS, YES)
+		call gseti (gp, G_NMAJOR, 6)
+		call gseti (gp, G_NMINOR, 4)
+		call gseti (gp, G_XNMAJOR, 4)
+		call gseti (gp, G_XNMINOR, 0)
+		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, wf1, wf2, wy1, wy2)
+		if (SF_AXIS(sf) == 1)
+		    call glabax (gp, "", "Focus", "Line (Slit)")
+		else
+		    call glabax (gp, "", "Focus", "Line (Dispersion)")
+
+		call gswind (gp, wf1, wf2, 0.5, 0.5+ny)
+		last = INDEF
+		do j = 1, ny {
+		    y = j
+		    do i = 1, nx {
+			sfd = SFD(sf,i,j)
+			if (SF_DEL(sfd) == NO) {
+			    z = SF_WID(sfd)
+			    z = 0.008 * (1 + (z - ww1) / ww2 * 3)
+			    call gmark (gp, SF_FOC(sfd), y, GM_CIRCLE, z, z)
+			}
+		    }
+
+		    n = 0
+		    z = 0.
+		    do i = 1, nx {
+			sfd = SFD(sf,i,j)
+			if (SF_DEL(sfd) == NO) {
+			    z = z + SF_FOC(SFD(sf,i,j))
+			    n = n + 1
+			}
+		    }
+		    if (n > 0) {
+			if (!IS_INDEF(last))
+			    call gline (gp, x, last, z/n, y)
+			x = z / n
+			last = y
+		    }
+		}
+
+		call gseti (gp, G_PLTYPE, 2)
+		call gline (gp, SF_FOCUS(sf), 0.5, SF_FOCUS(sf), 0.5+ny)
+		call gseti (gp, G_PLTYPE, 1)
+	    }
+	}
+
+	# POS plot
+	if (posplot == YES && abs (wp2 - wp1) > .05) {
+	    z = wp2 - wp1
+	    wp1 = wp1 - 0.05 * z
+	    wp2 = wp2 + 0.05 * z
+
+	    if (nx > 1) {
+		call gseti (gp, G_XLABELTICKS, NO)
+		call gseti (gp, G_YLABELTICKS, YES)
+		call gseti (gp, G_NMAJOR, 6)
+		call gseti (gp, G_NMINOR, 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, wx1, wx2, wp1, wp2)
+		call glabax (gp, "", "", "Shift")
+
+		call gswind (gp, 0.5, 0.8+nx, wp1, wp2)
+		last = INDEF
+		do i = 1, nx {
+		    x = i
+		    do j = 1, ny {
+			sfd = SFD(sf,i,j)
+			if (SF_DEL(sfd) == NO) {
+			    z = SF_WID(sfd)
+			    z = 0.008 * (1 + (z - ww1) / ww2 * 3)
+			    call gmark (gp, x, SF_POS(sfd), GM_CIRCLE, z, z)
+			}
+		    }
+
+		    n = 0
+		    z = 0.
+		    do j = 1, ny {
+			sfd = SFD(sf,i,j)
+			if (SF_DEL(sfd) == NO) {
+			    z = z + SF_POS(SFD(sf,i,j))
+			    n = n + 1
+			}
+		    }
+		    if (n > 0) {
+			if (!IS_INDEF(last))
+			    call gline (gp, last, y, x, z/n)
+			y = z / n
+			last = x
+		    }
+		}
+
+		call gseti (gp, G_PLTYPE, 2)
+		call gline (gp, 0.5, 0., 0.8+nx, 0.)
+		call gseti (gp, G_PLTYPE, 1)
+	    }
+
+	    if (ny > 1) {
+		call gseti (gp, G_XLABELTICKS, YES)
+		call gseti (gp, G_YLABELTICKS, NO)
+		call gseti (gp, G_NMAJOR, 6)
+		call gseti (gp, G_NMINOR, 4)
+		call gseti (gp, G_XNMAJOR, 4)
+		call gseti (gp, G_XNMINOR, 0)
+		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, wp1, wp2, wy1, wy2)
+		call glabax (gp, "", "Shift", "")
+
+		call gswind (gp, wp1, wp2, 0.5, 0.5+ny)
+		last = INDEF
+		do j = 1, ny {
+		    y = j
+		    do i = 1, nx {
+			sfd = SFD(sf,i,j)
+			if (SF_DEL(sfd) == NO) {
+			    z = SF_WID(sfd)
+			    z = 0.008 * (1 + (z - ww1) / ww2 * 3)
+			    call gmark (gp, SF_POS(sfd), y, GM_CIRCLE, z, z)
+			}
+		    }
+
+		    n = 0
+		    z = 0.
+		    do i = 1, nx {
+			sfd = SFD(sf,i,j)
+			if (SF_DEL(sfd) == NO) {
+			    z = z + SF_POS(sfd)
+			    n = n + 1
+			}
+		    }
+		    if (n > 0) {
+			if (!IS_INDEF(last))
+			    call gline (gp, x, last, z/n, y)
+			x = z / n
+			last = y
+		    }
+		}
+
+		call gseti (gp, G_PLTYPE, 2)
+		call gline (gp, 0., 0.5, 0., 0.5+ny)
+		call gseti (gp, G_PLTYPE, 1)
+	    }
+	}
+
+	# Spatial plot
+	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, wx1, wx2, wy1, wy2)
+	call gseti (gp, G_NMAJOR, 6)
+	call gseti (gp, G_NMINOR, 4)
+	if (vx[2,1] == vx[1,1] && vy[2,1] == vy[1,1]) {
+	    call gseti (gp, G_LABELTICKS, YES)
+	    if (SF_AXIS(sf) == 1)
+		call glabax (gp, "", "Column (Dispersion)", "Line (Slit)")
+	    else
+		call glabax (gp, "", "Column (Slit)", "Line (Dispersion)")
+	} else if (vx[2,1] == vx[1,1]) {
+	    call gseti (gp, G_XLABELTICKS, NO)
+	    call gseti (gp, G_YLABELTICKS, YES)
+	    if (SF_AXIS(sf) == 1)
+		call glabax (gp, "", "", "Line (Slit)")
+	    else
+		call glabax (gp, "", "", "Line (Dispersion)")
+	} else if (vy[2,1] == vy[1,1]) {
+	    call gseti (gp, G_XLABELTICKS, YES)
+	    call gseti (gp, G_YLABELTICKS, NO)
+	    if (SF_AXIS(sf) == 1)
+		call glabax (gp, "", "Column (Dispersion)", "")
+	    else
+		call glabax (gp, "", "Column (Slit)", "")
+	} else {
+	    call gseti (gp, G_LABELTICKS, NO)
+	    call glabax (gp, "", "", "")
+	}
+
+	call gswind (gp, 0.5, 0.8+nx, 0.5, 0.5+ny)
+
+	do j = 1, ny {
+	    do i = 1, nx {
+		sfd = SFD(sf,i,j)
+		if (SF_DEL(sfd) == NO) {
+		    x = i
+		    y = j
+		    z = SF_WID(sfd)
+		    z = 0.008 * (1 + (z - ww1) / ww2 * 3)
+		    call gmark (gp, x, y, GM_CIRCLE, z, z)
+		    if (i == ix && j == iy) {
+			call gseti (gp, G_PLCOLOR, HLCOLOR)
+			call gmark (gp, x, y, GM_BOX, -.8, -.8)
+			call gseti (gp, G_PLCOLOR, 1)
+		    }
+		    if (abs (wp2 - wp1) > .05) {
+			if (SF_AXIS(sf) == 1) {
+			    z =  0.25 * SF_POS(sfd) / max (abs(wp1), abs(wp2))
+			    call gadraw (gp, x+z, y)
+			} else {
+			    z = 0.25 * SF_POS(sfd) / max (abs(wp1), abs(wp2))
+			    call gadraw (gp, x, y+z)
+			}
+		    }
+		}
+	    }
+	}
+
+	if (nx <= 3 && ny <= 3) {
+	    call gsetr (gp, G_PLWIDTH, 2.)
+	    call gline (gp, wx1, wy1, wx1, wy1)
+	    do j = 1, ny {
+		do i = 1, nx {
+		    sfd = SFD(sf,i,j)
+		    if (SF_DEL(sfd) == NO) {
+			x = i
+			y = j
+			
+			call sprintf (Memc[str], SZ_LINE, "%.3g")
+			    call pargr (SF_FOC(sfd))
+			call gtext (gp, x+0.2, y+0.2, Memc[str], "h=l;v=c")
+			call sprintf (Memc[str], SZ_LINE, "%.2f")
+			    call pargr (SF_WID(sfd))
+			call gtext (gp, x+0.2, y-0.2, Memc[str], "h=l;v=c")
+			if (abs (SF_POS(sfd)) >= .01) {
+			    call sprintf (Memc[str], SZ_LINE, "%.2f")
+				call pargr (SF_POS(sfd))
+			    call gtext (gp, x+0.2, y, Memc[str], "h=l;v=c")
+			}
+		    }
+		}
+	    }
+	    call gsetr (gp, G_PLWIDTH, 1.)
+	}
+
+	if (SF_DATA(sf) == NULL) {
+	    call strcpy ("Best Focus at Each Sample", Memc[str], SZ_LINE)
+	} else {
+	    call sprintf (Memc[str], SZ_LINE, "Image %s with Focus %.3g")
+		call pargstr (SF_IMAGE(sf))
+		call pargr (SF_FOCUS(sf))
+	}
+	call gseti (gp, G_DRAWAXES, 0)
+	call gsview (gp, vx[1,1], vx[3,2], vy[1,1], vy[3,2])
+	call glabax (gp, Memc[str], "", "")
+
+	call sfree (sp)
+end
+
+
+# SPF_G2 -- Profiles at each sample for a given image
+
+procedure spf_g2 (gp, sfavg, sfbest, sfcur, ix, iy, lag)
+
+pointer	gp		# GIO pointer
+pointer	sfavg		# Average image
+pointer	sfbest		# Best image
+pointer	sfcur		# Current image
+int	ix, iy		# Sample
+int	lag		# Maximum lag
+
+int	i, j, nx, ny
+real	x1, x2, y1, y2, z1, z2, dz, vx, dvx, vy, dvy, p, x, fa[10], asieval()
+pointer	sp, str, sfd, asi
+
+begin
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	# Set range and draw axes
+	i = min (lag, nint (3 * SF_WIDTH(sfavg)))
+	x1 = -i
+	x2 = i
+	y1 = -0.05
+	y2 = 1.55
+	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
+	call gswind (gp, x1, x2, y1, y2)
+
+	call gseti (gp, G_DRAWTICKS, NO)
+	call sprintf (Memc[str], SZ_LINE, "Image %s with Focus %.3g")
+	    call pargstr (SF_IMAGE(sfcur))
+	    call pargr (SF_FOCUS(sfcur))
+	call glabax (gp, Memc[str], "", "")
+
+	# Set subviewport
+	nx = SF_NX(sfcur)
+	ny = SF_NY(sfcur)
+	call ggview (gp, vx, dvx, vy, dvy)
+	dvx = (dvx - vx) / nx
+	dvy = (dvy - vy) / ny
+
+	# Draw correlation profiles
+	do j = 1, ny {
+	    do i = 1, nx {
+		call gsview (gp, vx+(i-1)*dvx, vx+i*dvx, vy+(j-1)*dvy, vy+j*dvy)
+		call glabax (gp, "", "", "")
+
+		if (i == ix && j == iy) {
+		    call gsview (gp, vx+(i-1)*dvx+.005, vx+i*dvx-.005,
+			vy+(j-1)*dvy+.005, vy+j*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)
+		}
+
+		sfd = SFD(sfcur,i,j)
+		if (SF_DEL(sfd) == NO) {
+		    asi = SF_ASI(sfd)
+		    p = sqrt (2.) * SF_POS(sfd)
+		    z1 = max (x1, x1-p)
+		    z2 = min (x2, x2-p)
+		    dz = 1
+		    for (x=nint(z1); x<=nint(z2); x=x+dz)
+			call gmark (gp, x+p, asieval (asi, x+lag+1),
+			    GM_PLUS, 2., 2.)
+		    call gamove (gp, z1+p, asieval (asi, z1+lag+1))
+		    dz = .1
+		    for (x=z1+dz; x<=z2; x=x+dz)
+			call gadraw (gp, x+p, asieval (asi, x+lag+1))
+		    if (sfcur != sfbest) {
+			asi = SF_ASI(SFD(sfbest,i,j))
+			call gamove (gp, z1+p, asieval (asi, z1+lag+1))
+			for (x=z1+dz; x<=z2; x=x+dz)
+			    call gadraw (gp, x+p, asieval (asi, x+lag+1))
+		    }
+
+		    call gseti (gp, G_PLTYPE, 3)
+		    call gline (gp, 0., 0., 0., 1.)
+		    call gseti (gp, G_PLTYPE, 1)
+
+		    if (nx <= NMAX && ny <= NMAX)
+			call spf_label (gp, sfd, 21, 2)
+		}
+	    }
+	}
+
+	call gswind (gp, 0.5, 0.5+nx, 0.5, 0.5+ny)
+	call gsview (gp, vx, vx+nx*dvx, vy, vy+ny*dvy)
+	call gamove (gp, 1., 1.)
+
+	call sfree (sp)
+end
+
+
+# SPF_G3 -- Width vs. Focus
+
+procedure spf_g3 (gp, sfavg, sfcur, sfs, nimages, ix, iy)
+
+pointer	gp			# GIO pointer
+pointer	sfavg			# Average image
+pointer	sfcur			# Current image
+pointer	sfs[nimages]		# Images
+int	nimages			# Number of images
+int	ix, iy			# Current sample
+
+int	i, j, k, mark
+real	x, x1, x2, dx, y, y1, y2, dy, size
+pointer	sf, sfd
+
+begin
+	# Determine data range
+	x1 = MAX_REAL
+	y1 = MAX_REAL
+	x2 = -MAX_REAL
+	y2 = -MAX_REAL
+	do i = 1, nimages {
+	    sf = sfs[i]
+	    do j = 1, SF_NSFD(sf) {
+		sfd = SFD(sf,j,1)
+		if (SF_DEL(sfd) == NO) {
+		    x = SF_FOC(sfd)
+		    y = SF_WID(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
+	call gswind (gp, x1, x2, y1, y2)
+	call glabax (gp, "Profile Width vs. Focus", "", "")
+
+	do k = 1, nimages {
+	    sf = sfs[k]
+	    do j = 1, SF_NY(sf) {
+		do i = 1, SF_NX(sf) {
+		    call gseti (gp, G_PLCOLOR, 1)
+		    if (sf == sfcur)
+			mark = GM_PLUS
+		    else
+			mark = GM_CROSS
+		    size = 2.
+		    if (sf == sfcur && i == ix && j == iy) {
+			call gseti (gp, G_PLCOLOR, HLCOLOR)
+			mark = mark + GM_BOX
+			size = 3.
+		    }
+		    sfd = SFD(sf,i,j)
+		    if (SF_DEL(sfd) == NO) {
+			x = SF_FOC(sfd)
+			y = SF_WID(sfd)
+			call gmark (gp, x, y, mark, size, size)
+		    }
+		}
+	    }
+	}
+
+	call gseti (gp, G_PLTYPE, 3)
+	x = INDEF
+	do k = 1, nimages {
+	    sf = sfs[k]
+	    sfd = SFD(sf,ix,iy)
+	    if (SF_DEL(sfd) == NO) {
+		if (!IS_INDEF(x))
+		    call gline (gp, x, y, SF_FOC(sfd), SF_WID(sfd))
+		x = SF_FOC(sfd)
+		y = SF_WID(sfd)
+	    }
+	}
+
+	call gseti (gp, G_PLTYPE, 2)
+	call gline (gp, SF_FOCUS(sfavg), y1, SF_FOCUS(sfavg), y2)
+	call gline (gp, x1, SF_WIDTH(sfavg), x2, SF_WIDTH(sfavg))
+	call gseti (gp, G_PLTYPE, 1)
+end
+
+
+# SPF_G4 -- Profiles at a given sample
+
+procedure spf_g4 (gp, sfavg, sfbest, sfcur, sfs, nimages, ix, iy, nx, ny, lag)
+
+pointer	gp		# GIO pointer
+pointer	sfavg		# Average image
+pointer	sfbest		# Best image
+pointer	sfcur		# Current image
+pointer	sfs[nimages]	# Images
+int	nimages		# Number of images
+int	ix, iy		# Sample
+int	nx, ny		# Grid layout
+int	lag		# Maximum lag
+
+int	i, j, k
+real	x1, x2, y1, y2, z1, z2, dz, p, x, fa[10], asieval()
+real	vx, dvx, vy, dvy
+pointer	sp, str, sf, sfd, asi
+
+begin
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	if (SF_NSFD(sfcur) > 1) {
+	    call sprintf (Memc[str], SZ_LINE,
+		"All Images at Sample [%d:%d,%d:%d]")
+		call pargi (SF_X1(sfcur)+(ix-1)*SF_DX(sfcur))
+		call pargi (SF_X1(sfcur)+ix*SF_DX(sfcur)-1)
+		call pargi (SF_Y1(sfcur)+(iy-1)*SF_DY(sfcur))
+		call pargi (SF_Y1(sfcur)+iy*SF_DY(sfcur)-1)
+	} else
+	     Memc[str] = EOS
+
+	# Set windows
+	i = min (lag, nint (3 * SF_WIDTH(sfavg)))
+	x1 = -i
+	x2 = i
+	y1 = -0.05
+	y2 = 1.55
+	call gswind (gp, x1, x2, y1, y2)
+	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
+
+	# Set subview port
+	call ggview (gp, vx, dvx, vy, dvy)
+	dvx = (dvx - vx) / nx
+	dvy = (dvy - vy) / ny
+
+	# Draw correlation profiles
+	k = 0
+	do i =  1, ny {
+	    do j = 1, nx {
+		k = k + 1
+		if (k > nimages)
+		    break
+		sf = sfs[k]
+		sfd = SFD(sf,ix,iy)
+
+		call gsview (gp, vx+(j-1)*dvx, vx+j*dvx,
+		    vy+(ny-i)*dvy, vy+(ny-i+1)*dvy)
+		call gfill (gp, fa, fa[6], 4, GF_SOLID)
+		if (sf == sfcur) {
+		    call gsview (gp, vx+(j-1)*dvx+.005, vx+j*dvx-.005,
+			vy+(ny-i)*dvy+.005, vy+(ny-i+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+(j-1)*dvx, vx+j*dvx,
+			vy+(ny-i)*dvy, vy+(ny-i+1)*dvy)
+		}
+		call gseti (gp, G_DRAWAXES, 3)
+		call gseti (gp, G_DRAWTICKS, NO)
+		call glabax (gp, "", "", "")
+
+
+		if (SF_DEL(sfd) == NO) {
+		    asi = SF_ASI(sfd)
+		    p = sqrt (2.) * SF_POS(sfd)
+		    z1 = max (x1, x1-p)
+		    z2 = min (x2, x2-p)
+		    dz = 1
+		    for (x=nint(z1); x<=nint(z2); x=x+dz)
+			call gmark (gp, x+p, asieval (asi, x+lag+1),
+			    GM_PLUS, 2., 2.)
+		    call gamove (gp, z1+p, asieval (asi, z1+lag+1))
+		    dz = .1
+		    for (x=z1+dz; x<=z2; x=x+dz)
+			call gadraw (gp, x+p, asieval (asi, x+lag+1))
+		    if (sf != sfbest) {
+			asi = SF_ASI(SFD(sfbest,ix,iy))
+			call gamove (gp, z1+p, asieval (asi, z1+lag+1))
+			for (x=z1+dz; x<=z2; x=x+dz)
+			    call gadraw (gp, x+p, asieval (asi, x+lag+1))
+		    }
+
+		    call gseti (gp, G_PLTYPE, 3)
+		    call gline (gp, 0., 0., 0., 1.)
+		    call gseti (gp, G_PLTYPE, 1)
+
+		    if (nx <= NMAX && ny <= NMAX)
+			call spf_label (gp, sfd, 31, 2)
+		}
+	    }
+	}
+
+	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, Memc[str], "", "")
+
+	call sfree (sp)
+end
+
+
+# SPF_G5 -- Spectra at a given sample
+
+procedure spf_g5 (gp, sfcur, sfs, nimages, ix, iy, nx, ny)
+
+pointer	gp		# GIO pointer
+pointer	sfcur		# Current image
+pointer	sfs[nimages]	# Images
+int	nimages		# Number of images
+int	ix, iy		# Sample
+int	nx, ny		# Grid layout
+
+int	i, j, k, npts
+real	x1, x2, y1, y2, vx, dvx, vy, dvy, fa[10]
+pointer	sp, str, sf, sfd, spec
+
+begin
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	# Set subviewport parameters
+	call ggview (gp, vx, dvx, vy, dvy)
+	dvx = (dvx - vx) / nx
+	dvy = (dvy - vy) / ny
+
+	# Draw bounding box and label
+	if (SF_NSFD(sfcur) > 1) {
+	    call sprintf (Memc[str], SZ_LINE,
+		"All Images at Sample [%d:%d,%d:%d]")
+		call pargi (SF_X1(sfcur)+(ix-1)*SF_DX(sfcur))
+		call pargi (SF_X1(sfcur)+ix*SF_DX(sfcur)-1)
+		call pargi (SF_Y1(sfcur)+(iy-1)*SF_DY(sfcur))
+		call pargi (SF_Y1(sfcur)+iy*SF_DY(sfcur)-1)
+	} else
+		 Memc[str] = EOS
+
+	call gseti (gp, G_DRAWAXES, 0)
+	call gseti (gp, G_DRAWTICKS, NO)
+	call glabax (gp, Memc[str], "", "")
+	call gseti (gp, G_DRAWAXES, 3)
+
+	# Draw spectra
+	k = 0
+	do j = 1, ny {
+	    do i = 1, nx {
+		k = k + 1
+		if (k > nimages)
+		    break
+
+		sf = sfs[k]
+		sfd = SFD(sf, ix, iy)
+		spec = SF_SPEC(sfd)
+		if (SF_AXIS(sf) == 1)
+		    npts = SF_DX(sf)
+		else
+		    npts = SF_DY(sf)
+
+		x1 = 0.; x2 = 1.
+		call alimr (Memr[spec], npts, y1, y2)
+		y2 = y1 + (y2 - y1) * 1.5
+		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
+		call gswind (gp, x1, x2, y1, y2)
+
+		if (sf == sfcur) {
+		    call gsview (gp, vx+(i-1)*dvx+.005, vx+i*dvx-.005,
+			vy+(ny-j)*dvy+.005, vy+(ny-j+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+(i-1)*dvx, vx+i*dvx,
+		    vy+(ny-j)*dvy, vy+(ny-j+1)*dvy)
+		call glabax (gp, "", "", "")
+
+		if (SF_DEL(sfd) == NO) {
+		    call gvline (gp, Memr[spec], npts, x1, x2)
+
+		    if (nx <= NMAX && ny <= NMAX)
+			call spf_label (gp, sfd, 31, 2)
+		}
+	    }
+	}
+
+	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.)
+
+	call sfree (sp)
+end
+
+
+# SPF_G6 -- Spectra at a given image
+
+procedure spf_g6 (gp, sfcur, ix, iy)
+
+pointer	gp		# GIO pointer
+pointer	sfcur		# Current image
+int	ix, iy		# Sample
+
+int	i, j, nx, ny, npts
+real	x1, x2, y1, y2, vx, dvx, vy, dvy, fa[10]
+pointer	sp, str, sfd, spec
+
+begin
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	nx = SF_NX(sfcur)
+	ny = SF_NY(sfcur)
+	if (SF_AXIS(sfcur) == 1)
+	    npts = SF_DX(sfcur)
+	else
+	    npts = SF_DY(sfcur)
+
+	# Set subviewport parameters
+	call ggview (gp, vx, dvx, vy, dvy)
+	dvx = (dvx - vx) / nx
+	dvy = (dvy - vy) / ny
+
+	# Draw bounding box and label
+	call sprintf (Memc[str], SZ_LINE, "Image %s with Focus %.3g")
+	    call pargstr (SF_IMAGE(sfcur))
+	    call pargr (SF_FOCUS(sfcur))
+
+	call gseti (gp, G_DRAWTICKS, NO)
+	call glabax (gp, Memc[str], "", "")
+
+	# Draw spectra
+	do j = 1, ny {
+	    do i = 1, nx {
+		sfd = SFD(sfcur,i,j)
+		spec = SF_SPEC(sfd)
+
+		x1 = 0.; x2 = 1.
+		call alimr (Memr[spec], npts, y1, y2)
+		y2 = y1 + (y2 - y1) * 1.5
+		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
+		call gswind (gp, x1, x2, y1, y2)
+
+		if (i == ix && j == iy) {
+		    call gsview (gp, vx+(i-1)*dvx+.005, vx+i*dvx-.005,
+			vy+(j-1)*dvy+.005, vy+j*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+(i-1)*dvx, vx+i*dvx, vy+(j-1)*dvy, vy+j*dvy)
+		call glabax (gp, "", "", "")
+
+		if (SF_DEL(sfd) == NO) {
+		    call gvline (gp, Memr[spec], npts, x1, x2)
+
+		    if (nx <= NMAX && ny <= NMAX)
+			call spf_label (gp, sfd, 21, 2)
+		}
+	    }
+	}
+
+	call gsview (gp, vx, vx+nx*dvx, vy, vy+ny*dvy)
+	call gswind (gp, 0.5, 0.5+nx, 0.5, 0.5+ny)
+	call gamove (gp, 1., 1.)
+
+	call sfree (sp)
+end
+
+
+# SPF_G7 -- Profile at a given image and sample
+
+procedure spf_g7 (gp, sfbest, sfcur, ix, iy, lag)
+
+pointer	gp		# GIO pointer
+pointer	sfbest		# Best image
+pointer	sfcur		# Current image
+int	ix, iy		# Sample
+int	lag		# Maximum lag
+
+real	x1, x2, y1, y2, z1, z2, dz, x, p, s, asieval()
+pointer	sp, str, sf, sfd, asi
+
+begin
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	sf = sfcur
+	sfd = SFD(sf,ix,iy)
+	asi = SF_ASI(sfd)
+
+	s = 1 / sqrt (2.)
+	x2 = s * min (lag, nint (3 * SF_WIDTH(sf)))
+	x1 = -x2
+	y1 = -0.05
+	y2 = 1.05
+	call gswind (gp, x1, x2, y1, y2)
+
+	if (abs (SF_POS(sfd)) > .01) {
+	    call sprintf (Memc[str], SZ_LINE,
+	    "%s at (%d, %d), Focus %.3g, Width %.2f, Shift %.2f")
+	    call pargstr (SF_IMAGE(sf))
+	    call pargr (SF_X(sfd))
+	    call pargr (SF_Y(sfd))
+	    call pargr (SF_FOC(sfd))
+	    call pargr (SF_WID(sfd))
+	    call pargr (SF_POS(sfd))
+	} else {
+	    call sprintf (Memc[str], SZ_LINE,
+		"%s at (%d, %d), Focus %.3g, Width %.2f")
+	    call pargstr (SF_IMAGE(sf))
+	    call pargr (SF_X(sfd))
+	    call pargr (SF_Y(sfd))
+	    call pargr (SF_FOC(sfd))
+	    call pargr (SF_WID(sfd))
+	}
+	call glabax (gp, Memc[str], "Pixel", "Correlation")
+
+	# Draw correlation profiles
+	if (SF_DEL(sfd) == NO) {
+	    p = SF_POS(sfd)
+	    z1 = max (x1, x1-p)
+	    z2 = min (x2, x2-p)
+	    dz = s
+	    for (x=s*nint(z1/s); x<=s*nint(z2/s); x=x+dz)
+		call gmark (gp, x+p, asieval (asi, x/s+lag+1), GM_PLUS, 2., 2.)
+	    call gamove (gp, z1+p, asieval (asi, z1/s+lag+1))
+	    dz = .1 * s
+	    for (x=z1+dz; x<=z2; x=x+dz)
+		call gadraw (gp, x+p, asieval (asi, x/s+lag+1))
+	    if (sf != sfbest) {
+		asi = SF_ASI(SFD(sfbest,ix,iy))
+		call gamove (gp, z1+p, asieval (asi, z1/s+lag+1))
+		for (x=z1+dz; x<=z2; x=x+dz)
+		    call gadraw (gp, x+p, asieval (asi, x/s+lag+1))
+	    }
+
+	    call gseti (gp, G_PLTYPE, 3)
+	    call gline (gp, 0., y1, 0., y2)
+	    call gseti (gp, G_PLTYPE, 1)
+	}
+
+	call sfree (sp)
+end
+
+
+# SPF_G8 -- Spatial distribution of widths
+
+procedure spf_g8 (gp, sfavg, sfcur, sfdcur, sfs, nimages, nxgrid, nygrid)
+
+pointer	gp		# GIO pointer
+pointer	sfavg		# Average image
+pointer	sfcur		# Current image
+pointer	sfdcur		# Current sample
+pointer	sfs[nimages]	# Images
+int	nimages		# Number of images
+int	nxgrid, nygrid	# Grid layout
+
+int	nx, ny
+int	i, j, k, l, m
+real	x, y, z, x1, x2, y1, y2, z1, z2, fa[10]
+pointer	sp, str, sf, sfd, kmin, kptr
+
+begin
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	# Set subviewport parameters
+	nx = SF_NX(sfavg)
+	ny = SF_NY(sfavg)
+	call ggview (gp, x1, x2, y1, y2)
+	x2 = (x2 - x1) / nxgrid
+	y2 = (y2 - y1) / nygrid
+	fa[1] = 0.; fa[6] = 0.
+	fa[2] = 1.; fa[7] = 0.
+	fa[3] = 1.; fa[8] = 1.
+	fa[4] = 0.; fa[9] = 1.
+	fa[5] = 0.; fa[10] = 0.
+
+	call gseti (gp, G_DRAWTICKS, NO)
+
+	# Find best focus at each sample and range of WID
+	call salloc (kmin, nx*ny, TY_INT)
+	kptr = kmin
+	z1 = MAX_REAL
+	z2 = -MAX_REAL
+	do j = 1, ny {
+	    do i = 1, nx {
+		l = 0
+		sfd = SFD(sfavg,i,j)
+		if (SF_DEL(sfd) == NO) {
+		    x = SF_FOC(sfd)
+		    y = MAX_REAL
+		    do k = 1, nimages {
+			sfd = SFD(sfs[k],i,j)
+			if (SF_DEL(sfd) == NO) {
+			    z1 = min (z1, SF_WID(sfd))
+			    z2 = max (z2, SF_WID(sfd))
+			    z = abs (SF_FOC(sfd) - x)
+			    if (z < y) {
+				l = k
+				y = z
+			    }
+			}
+		    }
+		}
+		Memi[kptr] = l
+		kptr = kptr + 1
+	    }
+	}
+	z2 = max (2., z2 - z1)
+
+	# Make graphs
+	k = 0
+	do j = 1, nygrid {
+	    do i = 1, nxgrid {
+		k = k + 1
+		if (k > nimages)
+		    break
+
+		sf = sfs[k]
+		nx = SF_NX(sfs[1])
+		ny = SF_NY(sfs[1])
+
+		if (sf == sfcur) {
+		    call gsview (gp, x1+(i-1)*x2+.005, x1+i*x2-.005,
+			y1+(nygrid-j)*y2+.005, y1+(nygrid-j+1)*y2-.005)
+		    call gswind (gp, 0., 1., 0., 1.)
+		    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, x1+(i-1)*x2, x1+i*x2,
+		    y1+(nygrid-j)*y2, y1+(nygrid-j+1)*y2)
+		call gswind (gp, 0.5, 0.5+nx, 0.5, 1.5+ny)
+		call glabax (gp, "", "", "")
+
+		kptr = kmin
+		do m = 1, ny {
+		    do l = 1, nx {
+			sfd = SFD(sf,l,m)
+			if (SF_DEL(sfd) == NO) {
+			    x = l
+			    y = m
+			    z = SF_WID(sfd)
+			    z = 0.008 * (1 + (z - z1) / z2 * 3)
+			    call gmark (gp, x, y, GM_CIRCLE, z, z)
+			    if (Memi[kptr] == k)
+				call gmark (gp, x, y, GM_PLUS, z, z)
+			    if (sfd == sfdcur) {
+				call gseti (gp, G_PLCOLOR, HLCOLOR)
+				call gmark (gp, x, y, GM_BOX, -.8, -.8)
+				call gseti (gp, G_PLCOLOR, 1)
+			    }
+			}
+			kptr = kptr + 1
+		    }
+		}
+			
+		if (nxgrid <= NMAX && nygrid <= NMAX)
+		    call spf_label (gp, sfd, 11, 2)
+	    }
+	}
+
+	call gsview (gp, x1, x1+nxgrid*x2, y1, y1+nygrid*y2)
+	call gswind (gp, 0.5, 0.5+nxgrid, 0.5+nygrid, 0.5)
+	call gamove (gp, 1., 1.)
+
+	call sfree (sp)
+end
+
+
+# SPF_G9 -- Spectrum at a given image and sample
+
+procedure spf_g9 (gp, sfcur, ix, iy)
+
+pointer	gp		# GIO pointer
+pointer	sfcur		# Current image
+int	ix, iy		# Sample
+
+int	npts
+real	x1, x2
+pointer	sf, sfd, spec
+
+begin
+	sf = sfcur
+	sfd = SFD(sf,ix,iy)
+	spec = SF_SPEC(sfd)
+	if (SF_AXIS(sf) == 1) {
+	    npts = SF_DX(sf)
+	    x1 = SF_X1(sf) + (ix - 1) * npts
+	} else {
+	    npts = SF_DY(sf)
+	    x1 = SF_Y1(sf) + (iy - 1) * npts
+	}
+	x2 = x1 + npts - 1
+	call gswind (gp, x1, x2, INDEF, INDEF)
+	call gascale (gp, Memr[spec], npts, 2)
+
+	if (SF_AXIS(sf) == 1)
+	    call glabax (gp, "", "Column", "")
+	else
+	    call glabax (gp, "", "Line", "")
+	if (SF_DEL(sfd) == NO)
+	    call gvline (gp, Memr[spec], npts, x1, x2)
+end
+
+
+# SPF_LABEL -- Label
+
+procedure spf_label (gp, sfd, type, width)
+
+pointer	gp		# GIO pointer
+pointer	sfd		# Sample pointer
+int	type		# Label type
+int	width		# Line width
+
+int	bkup, gstati()
+real	x1, x2, y1, y2, xs, ys, x, y
+pointer	sp, str
+
+begin
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	bkup = gstati (gp, G_PLWIDTH)
+	call gseti (gp, G_PLWIDTH, width)
+	call ggwind (gp, x1, x2, y1, y2)
+	call ggscale (gp, x1, y1, xs, ys)
+	call gline (gp, x1, y1, x1, y1)
+
+	switch (type) {
+	case 11: # 1 across the top
+	    x = x2 - 0.01 * xs; y = y2 - 0.01 * ys
+	    call sprintf (Memc[str], SZ_LINE, "%.3g")
+		call pargr (SF_FOC(sfd))
+	    call gtext (gp, x, y, Memc[str], "h=r;v=t")
+	case 21: # 2 across the top
+	    x = x1 + 0.01 * xs; y = y2 - 0.01 * ys
+	    call sprintf (Memc[str], SZ_LINE, "%.2f")
+		call pargr (SF_WID(sfd))
+	    call gtext (gp, x, y, Memc[str], "h=l;v=t")
+	    if (abs (SF_POS(sfd)) >= .01) {
+		x = x2 - 0.01 * xs; y = y2 - 0.01 * ys
+		call sprintf (Memc[str], SZ_LINE, "%.2f")
+		    call pargr (SF_POS(sfd))
+		call gtext (gp, x, y, Memc[str], "h=r;v=t")
+	    }
+	case 31: # 3 across the top
+	    x = x1 + 0.01 * xs; y = y2 - 0.01 * ys
+	    call sprintf (Memc[str], SZ_LINE, "%.2f")
+		call pargr (SF_WID(sfd))
+	    call gtext (gp, x, y, Memc[str], "h=l;v=t")
+	    x = x2 - 0.01 * xs; y = y2 - 0.01 * ys
+	    call sprintf (Memc[str], SZ_LINE, "%.3g")
+		call pargr (SF_FOC(sfd))
+	    call gtext (gp, x, y, Memc[str], "h=r;v=t")
+	    if (abs (SF_POS(sfd)) >= .01) {
+		x = (x1 + x2) / 2; y = y2 - 0.01 * ys
+		call sprintf (Memc[str], SZ_LINE, "%.2f")
+		    call pargr (SF_POS(sfd))
+		call gtext (gp, x, y, Memc[str], "h=c;v=t")
+	    }
+	case 12: # 2 along the left
+	    x = x1 + 0.02 * xs; y = y2 - 0.02 * ys
+	    call sprintf (Memc[str], SZ_LINE, "%.2f")
+		call pargr (SF_WID(sfd))
+	    call gtext (gp, x, y, Memc[str], "h=l;v=t")
+	    if (abs (SF_POS(sfd)) >= .01) {
+		x = x1 + 0.02 * xs; y = y2 - 0.06 * ys
+		call sprintf (Memc[str], SZ_LINE, "%.2f")
+		    call pargr (SF_POS(sfd))
+		call gtext (gp, x, y, Memc[str], "h=l;v=t")
+	    }
+	case 13: # 3 along the left
+	    x = x1 + 0.02 * xs; y = y2 - 0.02 * ys
+	    call sprintf (Memc[str], SZ_LINE, "%.2f")
+		call pargr (SF_WID(sfd))
+	    call gtext (gp, x, y, Memc[str], "h=l;v=t")
+	    x = x1 + 0.02 * xs; y = y2 - 0.10 * ys
+	    call sprintf (Memc[str], SZ_LINE, "%.3g")
+		call pargr (SF_FOC(sfd))
+	    call gtext (gp, x, y, Memc[str], "h=l;v=t")
+	    if (abs (SF_POS(sfd)) >= .01) {
+		x = x1 + 0.02 * xs; y = y2 - 0.06 * ys
+		call sprintf (Memc[str], SZ_LINE, "%.2f")
+		    call pargr (SF_POS(sfd))
+		call gtext (gp, x, y, Memc[str], "h=l;v=t")
+	    }
+	}
+
+	call gseti (gp, G_PLWIDTH, bkup)
+
+	call sfree (sp)
+end
+
+
+# SPF_SAMPLE -- Find the nearest sample to the cursor position
+
+procedure spf_sample (sfs, nimages, del, wx, wy, i, j, k)
+
+pointer	sfs[nimages]	#I Images
+int	nimages		#I Number of images
+int	del		#I Deletion flag
+real	wx, wy		#I Cursor coordinate
+int	i, j, k		#O Nearest sample and image
+
+int	i1, j1, k1, k2
+real	r, rmin
+
+begin
+	rmin = MAX_REAL
+	k1 = k
+	do k2 = 0, 2 * nimages {
+	    if (mod (k2, 2) == 0)
+		k1 = k1 + k2
+	    else
+		k1 = k1 - k2
+	    if (k1 < 1 || k1 > nimages)
+		next
+	    do j1 = 1, SF_NY(sfs[k1]) {
+		do i1 = 1, SF_NX(sfs[k1]) {
+		    if (SF_DEL(SFD(sfs[k1],i1,j1)) == del) {
+			r = (i1 - wx) ** 2 + (j1 - wy) ** 2
+			if (r < rmin) {
+			    i = i1
+			    j = j1
+			    k = k1
+			    rmin = r
+			}
+		    }
+		}
+	    }
+
+	    if (rmin < MAX_REAL)
+		break
+	}
+end
diff --git a/noao/obsutil/src/specfocus/t_specfocus.x b/noao/obsutil/src/specfocus/t_specfocus.x
new file mode 100644
index 00000000..7f176de5
--- /dev/null
+++ b/noao/obsutil/src/specfocus/t_specfocus.x
@@ -0,0 +1,762 @@
+include	<error.h>
+include	<imhdr.h>
+include	<mach.h>
+include	<math.h>
+include	<math/curfit.h>
+include	<math/iminterp.h>
+include	"specfocus.h"
+
+
+# T_SPECFOCUS -- Spectral focusing task
+
+procedure t_specfocus ()
+
+int	list		# List of images
+pointer	fvals		# List of focus values
+int	dispaxis	# Default dispersion axis
+int	amin		# Lower edge of data along slit
+int	amax		# Upper edge of data along slit
+int	nspec		# Number of spectra to subdivide width
+int	ndisp		# Number of dispersion samples
+int	lag		# Maximum lag
+real	level		# Level for width
+bool	shifts		# Measure shifts?
+int	log		# Log file descriptor
+
+int	i, j, k, l, nimages, npix, nprep
+int	aaxis, a1, da, na
+int	baxis, b1, db, nb
+int	c1, dc, nc
+int	l1, dl, nl
+pointer	sp, image, sys
+pointer	rg, sfs, sf, sfd, sfavg, sfbest, im, mw, data, buf1, buf2
+pointer	rng_open(), immap(), imgl2r(), mw_openim()
+int	clgeti(), imgeti(), imtopenp(), imtlen(), imtgetim(), nowhite()
+int	rng_index(), open()
+real	rval, clgetr(), imgetr(), asumr()
+bool	ms, clgetb(), streq()
+errchk	immap, spf_width
+
+int	spf_compare()
+extern	spf_compare
+
+begin
+	call smark (sp)
+	call salloc (fvals, SZ_LINE, TY_CHAR)
+	call salloc (image, SZ_FNAME, TY_CHAR)
+	call salloc (sys, SZ_FNAME, TY_CHAR)
+
+	# Get task parameters (except log file)
+	list = imtopenp ("images")
+	call clgstr ("focus", Memc[fvals], SZ_LINE)
+	dispaxis = clgeti ("dispaxis")
+	amin = clgeti ("slit1")
+	amax = clgeti ("slit2")
+	nspec = clgeti ("nspectra")
+	ndisp = clgeti ("ndisp")
+	lag = (clgeti ("corwidth") + 1) / 2
+	level = clgetr ("level")
+	shifts = clgetb ("shifts")
+
+	if (level > 1.)
+	    level = level / 100.
+	level = max (0.05, min (0.95, level))
+
+	# Initialize focus values
+	if (nowhite (Memc[fvals], Memc[fvals], SZ_LINE) == 0)
+	    call strcpy ("1x1", Memc[fvals], SZ_LINE)
+	iferr (rg = rng_open (Memc[fvals], -MAX_REAL, MAX_REAL, 1.))
+	    rg = NULL
+
+	# Allocate array for the image focus data structure pointers
+	nimages = imtlen (list)
+	call malloc (sfs, nimages, TY_POINTER)
+
+	# Accumulate the focus data
+	nimages = 0
+	while (imtgetim (list, Memc[image], SZ_FNAME) != EOF) {
+	    im = immap (Memc[image], READ_ONLY, 0)
+	    mw = mw_openim (im)
+	    call mw_gsystem (mw, Memc[sys], SZ_FNAME)
+	    ms = streq (Memc[sys], "multispec")
+	    call mw_close (mw)
+
+	    # Set the focus value
+	    if (rg != NULL) {
+		if (rng_index (rg, nimages+1, rval) == EOF)
+		    call error (1, "Focus list ended prematurely")
+	    } else
+		rval = imgetr (im, Memc[fvals])
+
+	    # Set dispersion and cross dispersion axes
+	    if (ms) {
+		baxis = 1
+		aaxis = 2
+
+		# Set sampling across the dispersion axis
+		if (IS_INDEFI (amin))
+		    i = 1
+		else
+		    i = amin
+		if (IS_INDEFI (amax))
+		    j = IM_LEN(im,aaxis)
+		else
+		    j = amax
+		a1 = max (1, min (i, j))
+		da = min (IM_LEN(im,aaxis), max (i, j)) - a1 + 1
+		if (da < 1)
+		    call error (1, "Error in slit limits")
+		na = da
+		da = da / na
+
+		# Set sampling along the dispersion axis
+		npix = IM_LEN(im,baxis)
+		nb = min (ndisp, npix / 100)
+		db = npix / nb
+		b1 = 1 + (npix - nb * db) / 2
+
+		# Set sampling along the columns and lines
+		c1 = b1; dc = db; nc = nb
+		l1 = a1; dl = da; nl = na
+	    } else {
+		iferr (baxis = imgeti (im, "dispaxis"))
+		    baxis = dispaxis
+		aaxis = 3 - baxis
+
+		# Set sampling across the dispersion axis
+		if (IS_INDEFI (amin))
+		    i = 1
+		else
+		    i = amin
+		if (IS_INDEFI (amax))
+		    j = IM_LEN(im,aaxis)
+		else
+		    j = amax
+		a1 = max (1, min (i, j))
+		da = min (IM_LEN(im,aaxis), max (i, j)) - a1 + 1
+		if (da < 1)
+		    call error (1, "Error in slit limits")
+		na = min (nspec, da)
+		da = da / na
+
+		# Set sampling along the dispersion axis
+		npix = IM_LEN(im,baxis)
+		nb = min (ndisp, npix / 100)
+		db = npix / nb
+		b1 = 1 + (npix - nb * db) / 2
+
+		# Set sampling along the columns and lines
+		if (baxis == 1) {
+		    c1 = b1; dc = db; nc = nb
+		    l1 = a1; dl = da; nl = na
+		} else {
+		    c1 = a1; dc = da; nc = na
+		    l1 = b1; dl = db; nl = nb
+		}
+	    }
+
+	    # Check for consistency
+	    if (nimages > 0) {
+		if (baxis!=SF_AXIS(sf) ||
+		    c1!=SF_X1(sf) || dc!=SF_DX(sf) || nc!=SF_NX(sf) ||
+		    l1!=SF_Y1(sf) || dl!=SF_DY(sf) || nl!=SF_NY(sf))
+		    call error (1, "Input images have different formats")
+	    }
+
+	    # Allocate the focus data structure for the image
+	    call spf_alloc (sf, Memc[image], rval, level, baxis, npix, na,
+		c1, dc, nc, l1, dl, nl)
+
+	    # Get the spectrum samples
+	    if (baxis == 1) {
+		do i = 1, na {
+		    k = a1 + (i - 1) * da
+		    l = k + da - 1
+		    data = SF_DATA(sf) + (i - 1) * npix
+		    do j = k, l
+			call aaddr (Memr[imgl2r(im,j)], Memr[data], Memr[data],
+			    npix)
+		}
+	    } else {
+		do j = 1, npix {
+		    buf1 = imgl2r (im, j)
+		    data = SF_DATA(sf) + j - 1
+		    do i = 1, na {
+			k = a1 + (i - 1) * da
+			Memr[data] = asumr (Memr[buf1+k-1], da)
+			data = data + npix
+		    }
+		}
+	    }
+	    Memi[sfs+nimages] = sf
+	    nimages = nimages + 1
+
+	    call imunmap (im)
+	}
+
+	if (nimages == 0)
+	    call error (1, "No input data")
+
+	# Sort the structures
+	call qsort (Memi[sfs], nimages, spf_compare)
+
+	# Allocate structure for the best focus
+	call spf_alloc (sfavg, "Best", INDEF, level, baxis, 0, 0, c1, dc, nc,
+	    l1, dl, nl)
+
+	# Compute the correlations and profile width and position
+	nprep = db + 2 * lag
+	call malloc (buf1, nprep, TY_REAL)
+	call malloc (buf2, nprep, TY_REAL)
+	l = (na + 1) / 2
+	do k = 1, nimages {
+	    sf = Memi[sfs+k-1]
+	    do i = 1, nb {
+		if (baxis == 1)
+		    sfd = SFD(sf,i,l)
+		else
+		    sfd = SFD(sf,l,i)
+		call spf_prep (Memr[SF_SPEC(sfd)], db, Memr[buf1], nprep)
+		call spf_corr (Memr[buf1], Memr[buf1], nprep, lag,
+		    SF_ASI(sfd), SF_POS(sfd), SF_WID(sfd), SF_LEVEL(sf))
+		do j = 1, na {
+		    if (j != l) {
+			if (baxis == 1)
+			    sfd = SFD(sf,i,j)
+			else
+			    sfd = SFD(sf,j,i)
+			call spf_prep (Memr[SF_SPEC(sfd)], db, Memr[buf2],
+			    nprep)
+			call spf_corr (Memr[buf2], Memr[buf2], nprep, lag,
+			    SF_ASI(sfd), SF_POS(sfd), SF_WID(sfd), SF_LEVEL(sf))
+			if (shifts)
+			    call spf_corr (Memr[buf2], Memr[buf1], nprep,
+				lag, SF_ASI(sfd), SF_POS(sfd), rval,
+				SF_LEVEL(sf))
+		    }
+		}
+	    }
+	}
+	call mfree (buf1, TY_REAL)
+	call mfree (buf2, TY_REAL)
+
+	# Set the averages
+	call spf_fitfocus (Memi[sfs], nimages, sfavg, sfbest)
+
+	# Graph the results
+	call spf_graph (sfavg, sfbest, Memi[sfs], nimages, lag)
+
+	# Log the results
+	call spf_log (sfavg, sfbest, Memi[sfs], nimages, shifts, STDOUT)
+	call clgstr ("logfile", Memc[image], SZ_FNAME)
+	ifnoerr (log = open (Memc[image], APPEND, TEXT_FILE)) {
+	    call spf_log (sfavg, sfbest, Memi[sfs], nimages, shifts, log)
+	    call close (log)
+	}
+
+	# Finish up
+	do i = 1, nimages
+	    call spf_free (Memi[sfs+i-1])
+	call spf_free (sfavg)
+	call mfree (sfs, TY_POINTER)
+	call rng_close (rg)
+	call imtclose (list)
+	call sfree (sp)
+end
+
+
+# SPF_ALLOC -- Allocate a focus data structure for an image
+
+procedure spf_alloc (sf, image, focus, level, axis, ndisp, nspec,
+    x1, dx, nx, y1, dy, ny)
+
+pointer	sf		# Image focus data structure
+char	image[ARB]	# Image name
+real	focus		# Focus value
+real	level		# Level for width
+int	axis		# Dispersion axis
+int	ndisp		# Number of pixels (along dispersion)
+int	nspec		# Number of spectra (across dispersion)
+int	x1, dx, nx	# X sampling
+int	y1, dy, ny	# Y sampling
+
+int	i, j
+pointer	data, sfd
+
+begin
+	call calloc (sf, LEN_SF, TY_STRUCT)
+
+	call strcpy (image, SF_IMAGE(sf), SZ_SFFNAME)
+	SF_FOCUS(sf) = focus
+	SF_WIDTH(sf) = INDEF
+	SF_LEVEL(sf) = level
+	SF_AXIS(sf) = axis
+	SF_X1(sf) = x1
+	SF_DX(sf) = dx
+	SF_NX(sf) = nx
+	SF_Y1(sf) = y1
+	SF_DY(sf) = dy
+	SF_NY(sf) = ny
+	call malloc (SF_SFD(sf), nx*ny, TY_POINTER)
+	SF_NSFD(sf) = nx*ny
+	SF_NPIX(sf) = ndisp
+	if (ndisp > 0)
+	    call calloc (SF_DATA(sf), ndisp * nspec, TY_REAL)
+
+	data = SF_DATA(sf)
+	do j = 1, ny {
+	    do i = 1, nx {
+		call calloc (sfd, LEN_SFD, TY_STRUCT)
+		SFD(sf,i,j) = sfd
+		SF_X(sfd) = x1 + (i - 0.5) * dx
+		SF_Y(sfd) = y1 + (j - 0.5) * dy
+		if (ndisp > 0) {
+		    if (axis == 1)
+			SF_SPEC(sfd) = data + (j-1)*ndisp + (i-1)*dx + x1-1
+		    else
+			SF_SPEC(sfd) = data + (i-1)*ndisp + (j-1)*dy + y1-1
+		}
+		call asiinit (SF_ASI(sfd), II_SPLINE3)
+		SF_FOC(sfd) = focus
+		SF_WID(sfd) = INDEF
+		SF_POS(sfd) = INDEF
+		SF_DEL(sfd) = NO
+	    }
+	}
+end
+
+
+# SPF_FREE -- Free a focus image data structure
+
+procedure spf_free (sf)
+
+pointer	sf		# Image focus data structure
+
+int	i
+pointer	sfd
+
+begin
+	do i = 1, SF_NSFD(sf) {
+	    sfd = SFD(sf,i,1)
+	    call asifree (SF_ASI(sfd))
+	    call mfree (sfd, TY_STRUCT)
+	}
+	call mfree (SF_DATA(sf), TY_REAL)
+	call mfree (SF_SFD(sf), TY_POINTER)
+	call mfree (sf, TY_STRUCT)
+end
+
+
+# SPF_PREP -- Prepare spectra for correlation: fit continuum, subtract, taper
+
+procedure spf_prep (in, nin, out, nout)
+
+real	in[nin]			# Input spectrum
+int	nin			# Number of pixels in input spectrum
+real	out[nout]		# Output spectrum
+int	nout			# Number of pixels output spectrum (nin+2*lag)
+
+int	i, lag
+real	cveval()
+pointer	sp, x, w, ic, cv
+
+begin
+	call smark (sp)
+	call salloc (x, nin, TY_REAL)
+	call salloc (w, nin, TY_REAL)
+
+	call ic_open (ic)
+	call ic_pstr (ic, "function", "chebyshev")
+	call ic_puti (ic, "order", 3)
+	call ic_putr (ic, "low", 3.)
+	call ic_putr (ic, "high", 1.)
+	call ic_puti (ic, "niterate", 5)
+	call ic_putr (ic, "grow", 1.)
+	call ic_putr (ic, "xmin", 1.)
+	call ic_putr (ic, "xmax", real(nin))
+
+	do i = 1, nin {
+	    Memr[x+i-1] = i
+	    Memr[w+i-1] = 1
+	}
+	call ic_fit (ic, cv, Memr[x], in, Memr[w], nin, YES, YES, YES, YES)
+
+	lag = (nout - nin) / 2
+	do i = 1-lag, 0
+	    out[i+lag] = 0.
+	do i = 1, lag-1
+	    out[i+lag] = (1-cos (PI*i/lag))/2 * (in[i] - cveval (cv, real(i))) 
+	do i = lag, nin-lag+1
+	    out[i+lag] = (in[i] - cveval (cv, real(i))) 
+	do i = nin-lag+2, nin
+	    out[i+lag] = (1-cos (PI*(nin+1-i)/lag))/2 *
+		(in[i] - cveval (cv, real(i))) 
+	do i = nin+1, nin+lag
+	    out[i+lag] = 0.
+
+	call cvfree (cv)
+	call ic_closer (ic)
+	call sfree (sp)
+end
+
+
+# SPF_CORR -- Correlate spectra, fit profile, and measure center/width
+
+procedure spf_corr (spec1, spec2, npix, lag, asi, center, width, level)
+
+real	spec1[npix]		# First spectrum
+real	spec2[npix]		# Second spectrum
+int	npix			# Number of pixels in spectra
+int	lag			# Maximum correlation lag
+pointer	asi			# Pointer to correlation profile interpolator
+real	center			# Center of profile
+real	width			# Width of profile
+real	level			# Level at which width is determined
+
+int	i, j, nprof
+real	x, p, pmin, pmax, asieval()
+pointer	sp, prof
+
+begin
+	nprof = 2 * lag + 1
+
+	call smark (sp)
+	call salloc (prof, nprof, TY_REAL)
+
+	do j = -lag, lag {
+	    p = 0.
+	    do i = 1+lag, npix-lag
+		p = p + spec1[i] * spec2[i-j]
+	    Memr[prof+j+lag] = p
+	}
+
+	# Fit interpolator
+	call asifit (asi, Memr[prof], nprof)
+
+	# Find the minimum and maximum
+	center = 1.
+	pmin = asieval (asi, 1.)
+	pmax = pmin
+	for (x=1; x<=nprof; x=x+.01) {
+	    p = asieval (asi, x)
+	    if (p < pmin)
+		pmin = p
+	    if (p > pmax) {
+		pmax = p
+		center = x
+	    }
+	}
+
+	# Normalize
+	pmax = pmax - pmin
+	do i = 0, nprof-1
+	    Memr[prof+i] = (Memr[prof+i] - pmin) / pmax
+
+	call asifit (asi, Memr[prof], nprof)
+
+	# Find the equal flux points
+	for (x=center; x>=1 && asieval (asi,x)>level; x=x-0.01)
+	    ;
+	width = x
+	for (x=center; x<=nprof && asieval (asi,x)>level; x=x+0.01)
+	    ;
+	width = (x - width - 0.01) / sqrt (2.)
+	center = center - lag - 1
+
+	call sfree (sp)
+end
+
+
+# SPF_FITFOCUS -- Find the best focus at each sample and the average over all
+# samples.
+
+procedure spf_fitfocus (sfs, nimages, sfavg, sfbest)
+
+pointer	sfs[nimages]		#I Images
+int	nimages			#I Number of images
+pointer	sfavg			#U Average image
+pointer	sfbest			#U Best focus image
+
+int	i, j, n, jmin, nims
+pointer	sp, x, y, z, sfd
+real	focus, fwhm, pos, foc
+bool	fp_equalr()
+
+define	avg_	10
+
+begin
+	call smark (sp)
+	call salloc (x, nimages, TY_REAL)
+	call salloc (y, nimages, TY_REAL)
+	call salloc (z, nimages, TY_REAL)
+
+	do i = 1, SF_NSFD(sfavg) {
+	    # Collect the focus values
+	    nims = 0
+	    do j = 1, nimages {
+		sfd = SFD(sfs[j],i,1)
+		if (SF_DEL(sfd) == NO) {
+		    Memr[x+nims] = SF_FOC(sfd)
+		    Memr[y+nims] = SF_WID(sfd)
+		    Memr[z+nims] = SF_POS(sfd)
+		    nims = nims + 1
+		}
+	    }
+	    sfd = SFD(sfavg,i,1)
+
+	    # Take the smallest width at each unique focus.
+	    if (nims > 0) {
+		call xt_sort3 (Memr[x], Memr[y], Memr[z], nims)
+		n = 0
+		do j = 1, nims-1 {
+		    if (fp_equalr (Memr[x+n], Memr[x+j])) {
+			if (Memr[y+n] > Memr[y+j]) {
+			    Memr[y+n] = Memr[y+j]
+			    Memr[z+n] = Memr[z+j]
+			}
+		    } else {
+			n = n + 1
+			Memr[x+n] = Memr[x+j]
+			Memr[y+n] = Memr[y+j]
+			Memr[z+n] = Memr[z+j]
+		    }
+		}
+
+		# Find the minimum width
+		jmin = 0
+		do j = 1, n
+		    if (Memr[y+j] < Memr[y+jmin])
+			jmin = j
+
+		# Use parabolic interpolation to find the best focus
+		if (jmin == 0 || jmin == n) {
+		    focus = Memr[x+jmin]
+		    fwhm = Memr[y+jmin]
+		    pos = Memr[z+jmin]
+		} else
+		    call spf_parab (Memr[x+jmin-1], Memr[y+jmin-1],
+			Memr[z+jmin-1], focus, fwhm, pos)
+
+		SF_FOC(sfd) = focus
+		SF_WID(sfd) = fwhm
+		SF_POS(sfd) = pos
+		SF_DEL(sfd) = NO
+	    } else {
+		SF_FOC(sfd) = INDEF
+		SF_WID(sfd) = INDEF
+		SF_POS(sfd) = INDEF
+		SF_DEL(sfd) = YES
+	    }
+	}
+
+	call sfree (sp)
+
+avg_
+	# Set the averages over all samples
+	n = 0
+	focus = 0.
+	fwhm = 0.
+	do i = 1, SF_NSFD(sfavg) {
+	    sfd = SFD(sfavg,i,1)
+	    if (SF_DEL(sfd) == NO) {
+		focus = focus + SF_FOC(sfd)
+		fwhm = fwhm + SF_WID(sfd)
+		n = n + 1
+	    }
+	}
+
+	if (n > 0) {
+	    SF_FOCUS(sfavg) = focus / n
+	    SF_WIDTH(sfavg) = fwhm / n
+	} else {
+	    SF_FOCUS(sfavg) = INDEF
+	    SF_WIDTH(sfavg) = INDEF
+	}
+
+	do j = 1, nimages {
+	    n = 0
+	    focus = 0.
+	    fwhm = 0.
+	    do i = 1, SF_NSFD(sfs[j]) {
+		sfd = SFD(sfs[j],i,1)
+		if (SF_DEL(sfd) == NO) {
+		    fwhm = fwhm + SF_WID(sfd)
+		    n = n + 1
+		}
+	    }
+
+	    if (n > 0)
+		SF_WIDTH(sfs[j]) = fwhm / n
+	    else
+		SF_WIDTH(sfs[j]) = INDEF
+	}
+
+	# Set the best focus image
+	sfbest = NULL
+	focus = SF_FOCUS(sfavg)
+	if (!IS_INDEF(focus)) {
+	    pos = MAX_REAL
+	    do j = 1, nimages {
+		foc = SF_FOCUS(sfs[j])
+		if (!IS_INDEF(foc)) {
+		    fwhm = abs (focus - foc)
+		    if (fwhm < pos) {
+			pos = fwhm
+			sfbest = sfs[j]
+		    }
+		}
+	    }
+	}
+end
+
+
+# SPF_PARAB -- Find the minimum of a parabolic fit to three points.
+
+procedure spf_parab (x, y, z, xmin, ymin, zmin)
+
+real	x[3]
+real	y[3]
+real	z[3]
+real	xmin
+real	ymin
+real	zmin
+
+real	x12, x13, x23, x212, x213, x223, y12, y13, y23, a, b, c
+
+begin
+	x12 = x[1] - x[2]
+	x13 = x[1] - x[3]
+	x23 = x[2] - x[3]
+	x212 = x[1] * x[1] - x[2] * x[2]
+	x213 = x[1] * x[1] - x[3] * x[3]
+	x223 = x[2] * x[2] - x[3] * x[3]
+	y12 = y[1] - y[2]
+	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] - b * x[3] - c * x[3] * x[3]
+	xmin = -b / (2 * c)
+	ymin = a + b * xmin + c * xmin * xmin
+
+	if (xmin < x[2])
+	    zmin = z[2] + (z[1] - z[2]) / (x[1] - x[2]) * (xmin - x[2])
+	else
+	    zmin = z[2] + (z[3] - z[2]) / (x[3] - x[2]) * (xmin - x[2])
+end
+
+
+# SPF_COMPARE -- Compare two structures by focus values
+
+int procedure spf_compare (sf1, sf2)
+
+pointer	sf1, sf2	# Structures to be compared.
+
+begin
+	if (SF_FOCUS[sf1] < SF_FOCUS[sf2])
+	    return (-1)
+	else if (SF_FOCUS[sf1] > SF_FOCUS[sf2])
+	    return (1)
+	else
+	    return (0)
+end
+
+
+# SPF_LOG -- Print log of results
+
+procedure spf_log (sfavg, sfbest, sfs, nimages, shifts, log)
+
+pointer	sfavg			# Average image
+pointer	sfbest			# Best focus image
+pointer	sfs[nimages]		# Images
+int	nimages			# Number of images
+bool	shifts			# Measure shifts?
+int	log			# Log file descriptor
+
+int	i, j
+pointer	sp, str, sfd
+
+begin
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	call sysid (Memc[str], SZ_LINE)
+	call fprintf (log, "SPECFOCUS: %s\n")
+	    call pargstr (Memc[str])
+
+	call fprintf (log,
+    "  Best average focus at %g with average width of %.2f at %d%% of peak\n\n")
+	    call pargr (SF_FOCUS(sfavg))
+	    call pargr (SF_WIDTH(sfavg))
+	    call pargr (100 * SF_LEVEL(sfavg))
+
+	call fprintf (log, "  -- Average Over All Samples\n\n")
+	call fprintf (log, "\t%25wImage  Focus  Width\n")
+	do i = 1, nimages {
+	    call fprintf (log, "\t%30s  %5.3g  %5.2f\n")
+		call pargstr (SF_IMAGE(sfs[i]))
+		call pargr (SF_FOCUS(sfs[i]))
+		call pargr (SF_WIDTH(sfs[i]))
+	}
+	call fprintf (log, "\n")
+
+	call fprintf (log, "  -- Image %s at Focus %g --\n")
+	    call pargstr (SF_IMAGE(sfbest))
+	    call pargr (SF_FOCUS(sfbest))
+
+	if (SF_NSFD(sfbest) > 1) {
+	    call fprintf (log, "\n\n\tWidth at %d%% of Peak:\n")
+		call pargr (100 * SF_LEVEL(sfavg))
+	    call fprintf (log, "\n\t%9w    Columns\n\t%9w  ")
+	    do i = 1, SF_NX(sfbest) {
+		call fprintf (log, "  %4d-%-4d")
+		    call pargi (SF_X1(sfbest)+(i-1)*SF_DX(sfbest))
+		    call pargi (SF_X1(sfbest)+i*SF_DX(sfbest)-1)
+	    }
+	    call fprintf (log, "\n\t   Lines  +")
+	    do i = 1, SF_NX(sfbest)
+		call fprintf (log, "-----------")
+	    do j = 1, SF_NY(sfbest) {
+		if (SF_DY(sfbest) > 1.) {
+		    call fprintf (log, "\n\t%4d-%-4d |")
+			call pargi (SF_Y1(sfbest)+(j-1)*SF_DY(sfbest))
+			call pargi (SF_Y1(sfbest)+j*SF_DY(sfbest)-1)
+		} else {
+		    call fprintf (log, "\n\t%9d |")
+			call pargi (SF_Y1(sfbest)+(j-1)*SF_DY(sfbest))
+			call pargi (SF_Y1(sfbest)+j*SF_DY(sfbest)-1)
+		}
+		do i = 1, SF_NX(sfbest) {
+		    sfd = SFD(sfbest,i,j)
+		    call fprintf (log, "   %5.2f   ")
+			call pargr (SF_WID(sfd))
+		}
+	    }
+	    if (shifts) {
+		call fprintf (log, 
+		    "\n\n\tPosition Shifts Relative To Central Sample:\n")
+		call fprintf (log, "\n\t%9w    Columns\n\t%9w  ")
+		do i = 1, SF_NX(sfbest) {
+		    call fprintf (log, "  %4d-%-4d")
+			call pargi (SF_X1(sfbest)+(i-1)*SF_DX(sfbest))
+			call pargi (SF_X1(sfbest)+i*SF_DX(sfbest)-1)
+		}
+		call fprintf (log, "\n\t   Lines  +")
+		do i = 1, SF_NX(sfbest)
+		    call fprintf (log, "-----------")
+		do j = 1, SF_NY(sfbest) {
+		    call fprintf (log, "\n\t%4d-%-4d |")
+			call pargi (SF_Y1(sfbest)+(j-1)*SF_DY(sfbest))
+			call pargi (SF_Y1(sfbest)+j*SF_DY(sfbest)-1)
+		    do i = 1, SF_NX(sfbest) {
+			sfd = SFD(sfbest,i,j)
+			call fprintf (log, "   %5.2f   ")
+			    call pargr (SF_POS(sfd))
+		    }
+		}
+	    }
+	}
+	call fprintf (log, "\n\n")
+	    
+	call sfree (sp)
+end
diff --git a/noao/obsutil/src/specfocus/x_specfocus.f b/noao/obsutil/src/specfocus/x_specfocus.f
new file mode 100644
index 00000000..12500041
--- /dev/null
+++ b/noao/obsutil/src/specfocus/x_specfocus.f
@@ -0,0 +1,146 @@
+      integer function sysruk (task, cmd, rukarf, rukint)
+      integer rukarf
+      integer rukint
+      integer*2 task(*)
+      integer*2 cmd(*)
+      integer i
+      integer ntasks
+      integer lmarg
+      integer rmarg
+      integer maxch
+      integer ncol
+      integer rukean
+      integer envgei
+      integer envscn
+      logical streq
+      logical xerpop
+      logical xerflg
+      common /xercom/ xerflg
+      integer iyy
+      integer dp(2)
+      integer*2 dict(10)
+      integer*2 st0001(9)
+      integer*2 st0002(6)
+      integer*2 st0003(3)
+      integer*2 st0004(6)
+      integer*2 st0005(6)
+      integer*2 st0006(4)
+      integer*2 st0007(6)
+      integer*2 st0008(2)
+      integer*2 st0009(29)
+      integer*2 st0010(25)
+      save
+      data (dict(iyy),iyy= 1, 8) /115,112,101, 99,102,111, 99,117/
+      data (dict(iyy),iyy= 9,10) /115, 0/
+      data (st0001(iyy),iyy= 1, 8) /116,116,121,110, 99,111,108,115/
+      data (st0001(iyy),iyy= 9, 9) / 0/
+      data st0002 / 99,104,100,105,114, 0/
+      data st0003 / 99,100, 0/
+      data st0004 /104,111,109,101, 36, 0/
+      data st0005 / 72, 79, 77, 69, 36, 0/
+      data st0006 /115,101,116, 0/
+      data st0007 /114,101,115,101,116, 0/
+      data st0008 / 9, 0/
+      data (st0009(iyy),iyy= 1, 8) /105,110,118, 97,108,105,100, 32/
+      data (st0009(iyy),iyy= 9,16) /115,101,116, 32,115,116, 97,116/
+      data (st0009(iyy),iyy=17,24) /101,109,101,110,116, 58, 32, 39/
+      data (st0009(iyy),iyy=25,29) / 37,115, 39, 10, 0/
+      data (st0010(iyy),iyy= 1, 8) /105,110,118, 97,108,105,100, 32/
+      data (st0010(iyy),iyy= 9,16) / 83, 69, 84, 32,105,110, 32, 73/
+      data (st0010(iyy),iyy=17,24) / 82, 65, 70, 32, 77, 97,105,110/
+      data (st0010(iyy),iyy=25,25) / 0/
+      data (dp(iyy),iyy= 1, 2) / 1, 0/
+      data lmarg /5/, maxch /0/, ncol /0/, rukean /3/
+      data ntasks /0/
+         if (.not.(ntasks .eq. 0)) goto 110
+            i=1
+120         if (.not.(dp(i) .ne. 0)) goto 122
+121            i=i+1
+               goto 120
+122         continue
+            ntasks = i - 1
+110      continue
+         if (.not.(task(1) .eq. 63)) goto 130
+            call xerpsh
+            rmarg = envgei (st0001)
+            if (.not.xerpop()) goto 140
+               rmarg = 80
+140         continue
+            call strtbl (4, dict, dp, ntasks, lmarg, rmarg, maxch, ncol)
+            sysruk = (0)
+            goto 100
+130      continue
+         if (.not.(streq(task,st0002) .or. streq(task,st0003))) goto 150
+            call xerpsh
+            if (.not.(cmd(rukarf) .eq. 0)) goto 170
+               call xerpsh
+               call xfchdr(st0004)
+               if (.not.xerpop()) goto 180
+                  call xfchdr(st0005)
+180            continue
+               goto 171
+170         continue
+               call xfchdr(cmd(rukarf))
+171         continue
+162         if (.not.xerpop()) goto 160
+               if (.not.(rukint .eq. 1)) goto 190
+                  call erract (rukean)
+                  if (xerflg) goto 100
+                  goto 191
+190            continue
+191            continue
+160         continue
+            sysruk = (0)
+            goto 100
+150      continue
+         if (.not.(streq(task,st0006) .or. streq(task,st0007))) goto 200
+            call xerpsh
+            if (.not.(cmd(rukarf) .eq. 0)) goto 220
+               call envlit (4, st0008, 1)
+               call xffluh(4)
+               goto 221
+220         continue
+            if (.not.(envscn (cmd) .le. 0)) goto 230
+               if (.not.(rukint .eq. 1)) goto 240
+                  call eprinf (st0009)
+                  call pargsr (cmd)
+                  goto 241
+240            continue
+                  goto 91
+241            continue
+230         continue
+221         continue
+212         if (.not.xerpop()) goto 210
+               if (.not.(rukint .eq. 1)) goto 250
+                  call erract (rukean)
+                  if (xerflg) goto 100
+                  goto 251
+250            continue
+91                call syspac (0, st0010)
+251            continue
+210         continue
+            sysruk = (0)
+            goto 100
+200      continue
+151      continue
+131      continue
+         if (.not.(streq (task, dict(dp(1))))) goto 260
+            call tspecs
+            sysruk = (0)
+            goto 100
+260      continue
+         sysruk = (-1)
+         goto 100
+100      return
+      end
+c     rukint  ruk_interact
+c     sysruk  sys_runtask
+c     envscn  envscan
+c     tspecs  t_specfocus
+c     envgei  envgeti
+c     syspac  sys_panic
+c     eprinf  eprintf
+c     rukarf  ruk_argoff
+c     rukean  ruk_eawarn
+c     pargsr  pargstr
+c     envlit  envlist
diff --git a/noao/obsutil/src/specfocus/x_specfocus.x b/noao/obsutil/src/specfocus/x_specfocus.x
new file mode 100644
index 00000000..063f2ce6
--- /dev/null
+++ b/noao/obsutil/src/specfocus/x_specfocus.x
@@ -0,0 +1 @@
+task    specfocus	= t_specfocus