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