Initial commit

28 files changed, 3383 insertions, 0 deletions

diff --git a/noao/digiphot/apphot/center/apbcenter.x b/noao/digiphot/apphot/center/apbcenter.x
new file mode 100644
index 00000000..9a607ee3
--- /dev/null
+++ b/noao/digiphot/apphot/center/apbcenter.x
@@ -0,0 +1,96 @@
+include <fset.h>
+include "../lib/apphot.h"
+include "../lib/display.h"
+
+# APBCENTER -- Procedure to determine accurate centers for a list of objects
+# in batch mode using a simple text file list.
+
+procedure apbcenter (ap, im, cl, out, id, ld, mgd, gid, interactive)
+
+pointer ap			# pointer to the apphot structure
+pointer	im			# pointer to the IRAF image
+int	cl			# the coordinate file descriptor
+int	out			# output file descriptor
+int	id			# output file sequence number
+int	ld			# the coordinate file list number
+pointer	mgd			# pointer to the metacode file stream
+pointer	gid			# pointer to the image display stream
+int	interactive		# mode of use
+
+int	stdin, ier, ild
+pointer	sp, str
+real	wx, wy
+int	apfitcenter(), fscan(), nscan(), strncmp(), apstati()
+
+begin
+	call smark (sp)
+	call salloc (str, SZ_FNAME, TY_CHAR)
+	call fstats (cl, F_FILENAME, Memc[str], SZ_FNAME)
+
+	# Initialize.
+	ild = ld
+
+        # Print query.
+	if (strncmp ("STDIN", Memc[str], 5) == 0) {
+	    stdin = YES
+	    call printf ("Type object x and y coordinates (^D or ^Z to end): ")
+	    call flush (STDOUT)
+	} else
+	    stdin = NO
+
+	# Loop over the coordinate file.
+	while (fscan (cl) != EOF) {
+
+	    # Get and store the coordinates.
+	    call gargr (wx)
+	    call gargr (wy)
+	    if (nscan () != 2) {
+	        if (stdin == YES) {
+		    call printf (
+		        "Type object x and y coordinates (^D or ^Z to end): ")
+		    call flush (STDOUT)
+	        } 
+		next
+	    }
+
+	    # Transform the input coordinates.
+	    switch (apstati(ap,WCSIN)) {
+	    case WCS_WORLD, WCS_PHYSICAL:
+		call ap_itol (ap, wx, wy, wx, wy, 1)
+	    case WCS_TV:
+		call ap_vtol (im, wx, wy, wx, wy, 1)
+	    default:
+		;
+	    }
+
+	    call apsetr (ap, CWX, wx)
+	    call apsetr (ap, CWY, wy)
+
+	    # Fit the center.
+	    ier = apfitcenter (ap, im, wx, wy)
+
+	    # Print and/or plot the results.
+	    if (interactive == YES) {
+		call ap_qcenter (ap, ier)
+		if (gid != NULL)
+		    call apmark (ap, gid, apstati (ap, MKCENTER), NO, NO)
+	    }
+	    call ap_cplot (ap, id, mgd, YES)
+	    if (id == 1)
+	        call ap_param (ap, out, "center")
+	    call ap_pcenter (ap, out, id, ild, ier)
+
+	    # Setup for the next object.
+	    id = id + 1
+	    ild = ild + 1
+	    call apsetr (ap, WX, wx)
+	    call apsetr (ap, WY, wy)
+	    if (stdin == YES) {
+		call printf (
+		    "Type object x and y coordinates (^D or ^Z to end): ")
+		call flush (STDOUT)
+	    } 
+	}
+
+	call sfree (sp)
+end
diff --git a/noao/digiphot/apphot/center/apcconfirm.x b/noao/digiphot/apphot/center/apcconfirm.x
new file mode 100644
index 00000000..6c570bd4
--- /dev/null
+++ b/noao/digiphot/apphot/center/apcconfirm.x
@@ -0,0 +1,75 @@
+include "../lib/apphot.h"
+include "../lib/noise.h"
+include "../lib/center.h"
+
+# AP_CCONFIRM -- Procedure to confirm the critical center parameters.
+
+procedure ap_cconfirm (ap, out, stid)
+
+pointer	ap		# pointer to the apphot structure
+int	out		# output file descriptor
+int	stid		# output file sequence number
+
+pointer	sp, str
+real	fwhmpsf, capert, skysigma, datamin, datamax
+int	apstati()
+real	apstatr(), ap_vfwhmpsf(), ap_vcapert(), ap_vsigma()
+real	ap_vdatamin(), ap_vdatamax()
+
+begin
+	call smark (sp)
+	call salloc (str, SZ_FNAME, TY_CHAR)
+
+	call printf ("\n")
+
+	# Confirm the centering algorithm.
+	call ap_vcstring (ap, Memc[str], SZ_FNAME)
+
+	if (apstati (ap, CENTERFUNCTION) != AP_NONE) {
+
+	    # Confirm the fwhmpsf.
+	    if (apstati (ap, CENTERFUNCTION) != AP_CENTROID1D) {
+		fwhmpsf = ap_vfwhmpsf (ap)
+	    } else
+	    	fwhmpsf = apstatr (ap, FWHMPSF)
+
+	    # Confirm the centering box.
+	    capert = 2.0 * ap_vcapert (ap)
+
+	    # Confirm the sky sigma parameter.
+	    skysigma = ap_vsigma (ap)
+
+	} else {
+
+	    fwhmpsf = apstatr (ap, FWHMPSF)
+	    capert = 2.0 * apstatr (ap, CAPERT)
+	    skysigma = apstatr (ap, SKYSIGMA)
+
+	}
+
+	# Confirm the sky sigma parameter.
+
+	# Confirm the minimum and maximum good data values.
+	datamin = ap_vdatamin (ap)
+	datamax = ap_vdatamax (ap)
+
+	call printf ("\n")
+
+	# Update the database file.
+	if (out != NULL && stid > 1) {
+	    call ap_sparam (out, KY_CSTRING, Memc[str], UN_CALGORITHM,
+		"centering algorithm")
+	    call ap_rparam (out, KY_FWHMPSF, fwhmpsf, UN_ASCALEUNIT,
+	        "full width half maximum of the psf")
+	    call ap_rparam (out, KY_CAPERT, capert, UN_CSCALEUNIT,
+	        "centering box width")
+	    call ap_rparam (out, KY_SKYSIGMA, skysigma, UN_NCOUNTS,
+	        "standard deviation of 1 sky pixel")
+	    call ap_rparam (out, KY_DATAMIN, datamin, UN_ACOUNTS,
+	        "minimum good data value")
+	    call ap_rparam (out, KY_DATAMAX, datamax, UN_ACOUNTS,
+	        "maximum good data value")
+	}
+
+	call sfree (sp)
+end
diff --git a/noao/digiphot/apphot/center/apcencolon.x b/noao/digiphot/apphot/center/apcencolon.x
new file mode 100644
index 00000000..9f2f40f8
--- /dev/null
+++ b/noao/digiphot/apphot/center/apcencolon.x
@@ -0,0 +1,329 @@
+include "../lib/apphot.h"
+include "../lib/display.h"
+include "../lib/noise.h"
+include "../lib/center.h"
+
+# APCENTERCOLON -- Process colon commands from the centering task.
+
+procedure apcentercolon (ap, im, cl, out, stid, ltid, cmdstr, newimage,
+	newbuf, newfit)
+
+pointer	ap		# pointer to the apphot structure
+pointer	im		# pointer to the iraf image
+int	cl		# coordinate file descriptor
+int	out		# output file descriptor
+int	stid		# output file sequence number
+int	ltid		# input list sequence number
+char	cmdstr[ARB]	# command string
+int	newimage	# new image ?
+int	newbuf		# new center buffer ?
+int	newfit		# new center fit ?
+
+int	junk
+pointer	sp, incmd, outcmd
+int	strdic()
+
+begin
+	# Allocate working space.
+	call smark (sp)
+	call salloc (incmd, SZ_LINE, TY_CHAR)
+	call salloc (outcmd, SZ_LINE, TY_CHAR)
+
+	# Get the command.
+	call sscan (cmdstr)
+	call gargwrd (Memc[incmd], SZ_LINE)
+	if (Memc[incmd] == EOS) {
+	    call sfree (sp)
+	    return
+	}
+
+	# Process the command.
+	if (strdic (Memc[incmd], Memc[outcmd], SZ_LINE, CCMDS) != 0)
+	    call ap_ccolon (ap, out, stid, cmdstr, newbuf, newfit)
+	else if (strdic (Memc[incmd], Memc[outcmd], SZ_LINE, APCMDS) != 0)
+	    call ap_apcolon (ap, im, cl, out, stid, ltid, cmdstr, newimage,
+	        newbuf, newfit, junk, junk, junk, junk)
+	else if (strdic (Memc[incmd], Memc[outcmd], SZ_LINE, NCMDS) != 0)
+	    call ap_nscolon (ap, im, out, stid, cmdstr, newbuf, newfit,
+		junk, junk, junk, junk)
+	else
+	    call ap_cimcolon (ap, cmdstr) 
+
+	call sfree (sp)
+end
+
+
+# AP_CCOLON -- Process colon commands affecting the centering parameters.
+
+procedure ap_ccolon (ap, out, stid, cmdstr, newbuf, newfit)
+
+pointer	ap		# pointer to the apphot structure
+int	out		# output file descriptor
+int	stid		# output list id number
+char	cmdstr		# command string
+int	newbuf		# new sky buffer
+int	newfit		# new sky fit
+
+bool	bval
+int	ival, ncmd, stat
+pointer	sp, cmd, str
+real	rval
+
+bool	itob()
+int	nscan(), strdic(), btoi(), apstati()
+real	apstatr()
+
+begin
+	call smark (sp)
+	call salloc (cmd, SZ_LINE, TY_CHAR)
+	call salloc (str, SZ_FNAME, TY_CHAR)
+
+	# Get the command.
+	call sscan (cmdstr)
+	call gargwrd (Memc[cmd], SZ_LINE)
+	if (Memc[cmd] == EOS) {
+	    call sfree (sp)
+	    return
+	}
+
+	# Process the command.
+	ncmd = strdic (Memc[cmd], Memc[cmd], SZ_LINE, CCMDS)
+	switch (ncmd) {
+
+	case CCMD_CBOXWIDTH:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("%s = %g %s\n")
+		    call pargstr (KY_CAPERT)
+		    call pargr (2.0 * apstatr (ap, CAPERT))
+		    call pargstr (UN_CSCALEUNIT)
+	    } else {
+		call apsetr (ap, CAPERT, (rval / 2.0))
+		if (stid > 1)
+		    call ap_rparam (out, KY_CAPERT, rval, UN_CSCALEUNIT,
+			"centering subraster width")
+		newbuf = YES
+		newfit = YES
+	    }
+
+	case CCMD_CALGORITHM:
+	    call gargwrd (Memc[cmd], SZ_LINE)
+	    if (Memc[cmd] == EOS) {
+		call apstats (ap, CSTRING, Memc[str], SZ_FNAME)
+	        call printf ("%s = %s %s\n")
+		    call pargstr (KY_CSTRING)
+		    call pargstr (Memc[str])
+		    call pargstr (UN_CALGORITHM)
+	    } else {
+		stat = strdic (Memc[cmd], Memc[cmd], SZ_LINE, CFUNCS)
+		if (stat > 0) {
+		    call apseti (ap, CENTERFUNCTION, stat)
+		    call apsets (ap, CSTRING, Memc[cmd])
+		    if (stid > 1)
+		        call ap_sparam (out, KY_CSTRING, Memc[cmd],
+			    UN_CALGORITHM, "centering algorihtm")
+		    newbuf = YES
+		    newfit = YES
+		}
+	    }
+
+	case CCMD_CTHRESHOLD:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("%s = %g %s\n")
+		    call pargstr (KY_CTHRESHOLD)
+		    call pargr (apstatr (ap, CTHRESHOLD))
+		    call pargstr (UN_CSIGMA)
+	    } else {
+		call apsetr (ap, CTHRESHOLD, rval)
+		if (stid > 1)
+		    call ap_rparam (out, KY_CTHRESHOLD, rval, UN_CSIGMA,
+			"threshold for centering")
+		newbuf = YES
+		newfit = YES
+	    }
+
+	case CCMD_MAXSHIFT:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("%s = %g %s\n")
+		    call pargstr (KY_MAXSHIFT)
+		    call pargr (apstatr (ap, MAXSHIFT))
+		    call pargstr (UN_CSCALEUNIT)
+	    } else {
+		call apsetr (ap, MAXSHIFT, rval)
+		if (stid > 1)
+		    call ap_rparam (out, KY_MAXSHIFT, rval, UN_CSCALEUNIT,
+			"maximum shift")
+		newfit = YES
+	    }
+
+	case CCMD_MINSNRATIO:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("%s = %g\n")
+		    call pargstr (KY_MINSNRATIO)
+		    call pargr (apstatr (ap, MINSNRATIO))
+	    } else {
+		call apsetr (ap, MINSNRATIO, rval)
+		if (stid > 1)
+		    call ap_rparam (out, KY_MINSNRATIO, rval,
+		        UN_CNUMBER, "minimum signal-to-noise ratio")
+		newfit = YES
+	    }
+
+	case CCMD_CMAXITER:
+	    call gargi (ival)
+	    if (nscan () == 1) {
+		call printf ("%s = %d\n")
+		    call pargstr (KY_CMAXITER)
+		    call pargi (apstati (ap, CMAXITER))
+	    } else {
+		call apseti (ap, CMAXITER, ival)
+		if (stid > 1)
+		    call ap_iparam (out, KY_CMAXITER, ival, UN_CNUMBER,
+			"maximum number of iterations")
+		newfit = YES
+	    }
+
+	case CCMD_CLEAN:
+	    call gargb (bval)
+	    if (nscan () == 1) {
+		call printf ("%s = %b\n")
+		    call pargstr (KY_CLEAN)
+		    call pargb (itob (apstati (ap, CLEAN)))
+	    } else {
+		call apseti (ap, CLEAN, btoi (bval))
+		if (stid > 1)
+		    call ap_bparam (out, KY_CLEAN, bval, UN_CSCALEUNIT,
+			"apply clean algorithm before centering")
+		newbuf = YES
+		newfit = YES
+	    }
+
+	case  CCMD_RCLEAN:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("%s = %g %s\n")
+		    call pargstr (KY_RCLEAN)
+		    call pargr (apstatr (ap, RCLEAN))
+		    call pargstr (UN_CSCALEUNIT)
+	    } else {
+		call apsetr (ap, RCLEAN, rval)
+		if (stid > 1)
+		    call ap_rparam (out, KY_RCLEAN, rval, UN_CSCALEUNIT,
+			"cleaning radius")
+		newbuf = YES
+		newfit = YES
+	    }
+
+	case CCMD_RCLIP:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("%s = %g %s\n")
+		    call pargstr (KY_RCLIP)
+		    call pargr (apstatr (ap, RCLIP))
+		    call pargstr (UN_CSCALEUNIT)
+	    } else {
+		call apsetr (ap, RCLIP, rval)
+		if (stid > 1)
+		    call ap_rparam (out, KY_RCLIP, rval, UN_CSCALEUNIT,
+			"clipping radius")
+		newbuf = YES
+		newfit = YES
+	    }
+
+	case CCMD_KCLEAN:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("%s = %g %s\n")
+		    call pargstr (KY_SIGMACLEAN)
+		    call pargr (apstatr (ap, SIGMACLEAN))
+		    call pargstr (UN_CSIGMA)
+	    } else {
+		call apsetr (ap, SIGMACLEAN, rval)
+		if (stid > 1)
+		    call ap_rparam (out, KY_SIGMACLEAN, rval,
+		        UN_CSIGMA, "k-sigma clean rejection criterion")
+		newbuf = YES
+		newfit = YES
+	    }
+
+	case CCMD_MKCENTER:
+	    call gargb (bval)
+	    if (nscan() == 1) {
+		call printf ("%s = %b\n")
+		    call pargstr (KY_MKCENTER)
+		    call pargb (itob (apstati (ap, MKCENTER)))
+	    } else
+		call apseti (ap, MKCENTER, btoi (bval))
+
+	default:
+	    call printf ("Unknown or ambiguous colon command\7\n")
+	}
+
+	call sfree (sp)
+end
+
+
+# AP_CIMCOLON --  Process colon commands for the center task that do
+# not affect the centering parameters
+
+procedure ap_cimcolon (ap, cmdstr)
+
+pointer	ap			# pointer to the apphot structure
+char	cmdstr[ARB]		# command string
+
+bool	bval
+int	ncmd
+pointer	sp, cmd
+bool	itob()
+int	strdic(), nscan(), apstati(), btoi()
+
+begin
+	call smark (sp)
+	call salloc (cmd, SZ_LINE, TY_CHAR)
+
+	# Get the command.
+	call sscan (cmdstr)
+	call gargwrd (Memc[cmd], SZ_LINE)
+	if (Memc[cmd] == EOS) {
+	    call sfree (sp)
+	    return
+	}
+
+	# Process the command.
+	ncmd = strdic (Memc[cmd], Memc[cmd], SZ_LINE, MISC)
+	switch (ncmd) {
+	case ACMD_SHOW:
+	    call gargwrd (Memc[cmd], SZ_LINE)
+	    ncmd = strdic (Memc[cmd], Memc[cmd], SZ_LINE, CSHOWARGS)
+	    switch (ncmd) {
+	    case CCMD_DATA:
+		call printf ("\n")
+		call ap_nshow (ap)
+		call printf ("\n")
+	    case CCMD_CENTER:
+		call printf ("\n")
+	        call ap_cpshow (ap)
+		call printf ("\n")
+	    default:
+		call printf ("\n")
+	        call ap_cshow (ap)
+		call printf ("\n")
+	    }
+	case ACMD_RADPLOTS:
+	    call gargb (bval)
+	    if (nscan () == 1) {
+		call printf ("%s = %b\n")
+		    call pargstr (KY_RADPLOTS)
+		    call pargb (itob (apstati (ap, RADPLOTS)))
+	    } else {
+		call apseti (ap, RADPLOTS, btoi (bval))
+	    }
+	default:
+	    call printf ("Unknown or ambiguous colon command\7\n")
+	}
+
+	call sfree (sp)
+end
diff --git a/noao/digiphot/apphot/center/apcenter.x b/noao/digiphot/apphot/center/apcenter.x
new file mode 100644
index 00000000..8e1b4c68
--- /dev/null
+++ b/noao/digiphot/apphot/center/apcenter.x
@@ -0,0 +1,357 @@
+include <ctype.h>
+include <gset.h>
+include <imhdr.h>
+include "../lib/apphot.h"
+include "../lib/display.h"
+include "../lib/center.h"
+
+define	HELPFILE	"apphot$center/center.key"
+
+# APCENTER -- Procedure to determine accurate centers for a list of objects
+# interactively.
+
+int procedure apcenter (ap, im, cl, gd, mgd, id, out, stid, interactive, cache)
+
+pointer ap			# pointer to the apphot structure
+pointer	im			# pointer to the IRAF image
+int	cl			# coordinate list file descriptor
+pointer	gd			# pointer to the graphics stream
+pointer	mgd			# pointer to the plot metacode stream
+pointer	id			# pointer to the image display stream
+int	out			# output file descriptor
+int	stid			# output file sequence number
+int	interactive		# interactive mode
+int	cache			# cache the input pixels
+
+real	wx, wy, xlist, ylist
+pointer	sp, cmd
+int	wcs, key, colonkey, newimage, newobject, newcenter, newlist, ier, oid
+int	ip, prev_num, req_num, ltid, buf_size, old_size, req_size, memstat
+
+real	apstatr()
+int	ctoi(), clgcur(), apgscur(), apfitcenter(), aprefitcenter(), apstati()
+int	apgqverify(), apgtverify(), apnew(), ap_memstat(), sizeof()
+
+define	endswitch_ 99
+
+begin
+	call smark (sp)
+	call salloc (cmd, SZ_LINE, TY_CHAR)
+
+	# Initialize cursor command.
+	key = ' '
+	Memc[cmd] = EOS
+
+	# Initialize fitting parameters.
+	newimage = NO
+	newobject = YES
+	newcenter = YES
+	ier = AP_OK
+
+	# Initialize sequencing.
+	newlist = NO
+	ltid = 0
+
+	# Loop over the cursor commands.
+	while (clgcur ("icommands", wx, wy, wcs, key, Memc[cmd], SZ_LINE) !=
+	    EOF) {
+
+	    # Store the current cursor coordinates.
+	    call ap_vtol (im, wx, wy, wx, wy, 1)
+	    call apsetr (ap, CWX, wx)
+	    call apsetr (ap, CWY, wy)
+
+	    # Test to see if the cursor has moved.
+	    if (apnew (ap, wx, wy, xlist, ylist, newlist) == YES) {
+		newobject = YES
+		newcenter = YES
+	    }
+
+	    # Process the colon commands.
+	    switch (key) {
+
+	    # Quit.
+	    case 'q':
+		if (interactive == YES) {
+		    if (apgqverify ("center", ap, key) == YES) {
+			call sfree (sp)
+			return (apgtverify (key))
+		    }
+		} else {
+		    call sfree (sp)
+		    return (NO)
+		}
+
+	    # Print error messages.
+	    case 'e':
+		if (interactive == YES)
+		    call ap_cerrors (ap, ier)
+
+	    # Get information on keystroke commands.
+	    case '?':
+		if ((id != NULL) && (id == gd))
+		    call gpagefile (id, HELPFILE, "")
+		else if (interactive == YES)
+		    call pagefile (HELPFILE, "[space=cmhelp,q=quit,?=help]")
+
+	    # Rewind the coord list.
+	    case 'r':
+		if (cl != NULL) {
+		    call seek (cl, BOF)
+		    ltid = 0
+		} else if (interactive == YES)
+		    call printf ("No coordinate list\n")
+
+	    # Draw a centered stellar radial profile plot.
+	    case 'd':
+		if (interactive == YES) {
+		    call ap_qrad (ap, im, wx, wy, gd)
+		    newobject = YES
+		    newcenter = YES
+		}
+
+	    # Interactively set the centering parameters.
+	    case 'i':
+		if (interactive == YES) {
+		    call ap_cradsetup (ap, im, wx, wy, gd, out, stid)
+		    newobject = YES
+		    newcenter = YES
+		}
+
+	    # Verify critical parameters.
+	    case 'v':
+		call ap_cconfirm (ap, out, stid)
+		newobject = YES
+		newcenter = YES
+
+	    # Measure the next object.
+	    case 'm', 'n':
+
+		# No coordinate file.
+		if (cl == NULL) {
+		    if (interactive == YES)
+		        call printf ("No coordinate list\n")
+		    goto endswitch_
+		}
+
+		# Need to rewind coordinate file.
+		prev_num = ltid
+		req_num = ltid + 1
+		if (apgscur (cl, id, xlist, ylist, prev_num, req_num, ltid) ==
+		    EOF) {
+		    if (interactive == YES)
+			call printf (
+			    "End of coordinate list, use r key to rewind\n")
+		    goto endswitch_
+		}
+
+		# Convert coordinates if necessary.
+		switch (apstati (ap, WCSIN)) {
+		case WCS_PHYSICAL, WCS_WORLD:
+		    call ap_itol (ap,  xlist, ylist, xlist, ylist, 1)
+		case WCS_TV:
+		    call ap_vtol (im, xlist, ylist, xlist, ylist, 1)
+		default:
+		    ;
+		}
+
+		# Move to next object.
+		newlist = YES
+		if (key == 'm') {
+		    newobject = YES
+		    newcenter = YES
+		    goto endswitch_
+		} 
+
+		# Measure next object.
+		ier = apfitcenter (ap, im, xlist, ylist)
+		if (id != NULL) {
+		    call apmark (ap, id, apstati (ap, MKCENTER), NO, NO)
+		    if (id == gd)
+		        call gflush (id)
+		    else
+		        call gframe (id)
+		}
+		call ap_cplot (ap, stid, gd, apstati (ap, RADPLOTS))
+		if (interactive == YES)
+		    call ap_qcenter (ap, ier)
+		if (stid == 1)
+		     call ap_param (ap, out, "center")
+		call ap_cplot (ap, stid, mgd, YES)
+		call ap_pcenter (ap, out, stid, ltid, ier)
+		stid = stid + 1
+		newobject = NO
+		newcenter = NO
+
+	    # Save centering parameters in the pset files.
+	    case 'w':
+		call ap_pcpars (ap)
+
+	    # Process apphot : commands.
+	    case ':':
+		for (ip = 1; IS_WHITE(Memc[cmd+ip-1]); ip = ip + 1)
+		    ;
+		colonkey = Memc[cmd+ip-1]
+
+		# Process the apphot :n command.
+		switch (colonkey) {
+		case 'm', 'n':
+
+		    # Show/set apphot commands
+		    if (Memc[cmd+ip] != EOS && Memc[cmd+ip] != ' ') {
+		        call apcentercolon (ap, im, cl, out, stid, ltid,
+			    Memc[cmd], newimage, newobject, newcenter)
+			goto endswitch_
+		    }
+
+		    # No coordinate list.
+		    if (cl == NULL) {
+		        if (interactive == YES)
+		            call printf ("No coordinate list\n")
+			goto endswitch_
+		    }
+
+		    # Get next object from list.
+		    ip = ip + 1
+		    prev_num = ltid
+		    if (ctoi (Memc[cmd], ip, req_num) <= 0)
+		        req_num = ltid + 1
+
+		    # Fetch the next object from the list.
+		    if (apgscur (cl, id, xlist, ylist, prev_num, req_num,
+			    ltid) == EOF) {
+		        if (interactive == YES)
+			    call printf (
+			    "End of coordinate list, use r key to rewind\n")
+			goto endswitch_
+		    }
+
+		    # Convert the coordinates.
+		    switch (apstati (ap, WCSIN)) {
+		    case WCS_PHYSICAL, WCS_WORLD:
+		        call ap_itol (ap,  xlist, ylist, xlist, ylist, 1)
+		    case WCS_TV:
+		        call ap_vtol (im, xlist, ylist, xlist, ylist, 1)
+		    default:
+		        ;
+		    }
+
+		    # Move to the next object.
+		    newlist = YES
+		    if (colonkey == 'm') {
+		    	newobject = YES
+			newcenter = YES
+			goto endswitch_
+		    }
+
+		    # Measure next object.
+		    ier = apfitcenter (ap, im, xlist, ylist)
+		    if (id != NULL) {
+		        call apmark (ap, id, apstati (ap, MKCENTER), NO, NO)
+		    	if (id == gd)
+			   call gflush (id)
+		    	else
+		           call gframe (id)
+		    }
+		    call ap_cplot (ap, stid, gd, apstati (ap, RADPLOTS))
+		    if (interactive == YES)
+		        call ap_qcenter (ap, ier)
+		    if (stid == 1)
+		        call ap_param (ap, out, "center")
+		    call ap_cplot (ap, stid, mgd, YES)
+		    call ap_pcenter (ap, out, stid, ltid, ier)
+		    stid = stid + 1
+		    newobject = NO
+		    newcenter = NO
+
+		# Show/set apphot parameters.
+		default:
+		    call apcentercolon (ap, im, cl, out, stid, ltid, Memc[cmd],
+		        newimage, newobject, newcenter)
+		}
+
+		if (newimage == YES) {
+		    if ((id != NULL) && (id != gd))
+		        call ap_gswv (id, Memc[cmd], im, 4)
+		    req_size = MEMFUDGE * IM_LEN(im,1) * IM_LEN(im,2) *
+			sizeof (IM_PIXTYPE(im))
+		    memstat = ap_memstat (cache, req_size, old_size)
+		    if (memstat == YES)
+		        call ap_pcache (im, INDEFI, buf_size)
+		}
+
+		newimage = NO
+
+	    # Fit the center, and save the results.
+	    case 'f',  ' ':
+
+		# Compute the center.
+		if (newobject == YES)
+		    ier = apfitcenter (ap, im, wx, wy)
+		else if (newcenter == YES)
+		    ier = aprefitcenter (ap, im, ier)
+		if (id != NULL) {
+		    call apmark (ap, id, apstati (ap, MKCENTER), NO, NO)
+		    if (id == gd)
+			call gflush (id)
+		    else
+		        call gframe (id)
+		}
+		call ap_cplot (ap, stid, gd, apstati (ap, RADPLOTS))
+		if (interactive == YES)
+		    call ap_qcenter (ap, ier)
+		newobject = NO; newcenter = NO
+
+		# Write the results to the output file(s).
+		if (key == ' ') {
+		    if (stid == 1)
+		        call ap_param (ap, out, "center")
+		    if (newlist == YES)
+		        call ap_pcenter (ap, out, stid, ltid, ier)
+		    else
+		        call ap_pcenter (ap, out, stid, 0, ier)
+		    call ap_cplot (ap, stid, mgd, YES)
+		    stid = stid + 1
+		}
+
+	    # Fit centers for the rest of the list.
+	    case 'l':
+
+		# No coordinate list.
+		if (cl == NULL) {
+		    if (interactive == YES)
+		        call printf ("No coordinate list\n")
+		    goto endswitch_
+		}
+
+		# Center rest of coordinate list.
+		ltid = ltid + 1
+		oid = stid
+		call apbcenter (ap, im, cl, out, stid, ltid, mgd, id, YES)
+		ltid = ltid + stid - oid + 1
+		if (id != NULL) {
+		    if (id == gd)
+		        call gflush (id)
+		    else
+		        call gframe (id)
+		}
+
+	    # Unknown command.
+	    default:
+		call printf ("Unknown or ambiguous keystroke command\n")
+	    }
+
+endswitch_
+	    # Setup for the next object by setting the default keystroke
+	    # command and storing the old cursor coordinates in the 
+	    # centering structure.
+
+	    key = ' '
+	    Memc[cmd] = EOS
+	    call apsetr (ap, WX, apstatr (ap, CWX))
+	    call apsetr (ap, WY, apstatr (ap, CWY))
+
+	}
+
+	call sfree (sp)
+end
diff --git a/noao/digiphot/apphot/center/apcerrors.x b/noao/digiphot/apphot/center/apcerrors.x
new file mode 100644
index 00000000..4f8c32bb
--- /dev/null
+++ b/noao/digiphot/apphot/center/apcerrors.x
@@ -0,0 +1,32 @@
+include "../lib/center.h"
+
+# AP_CERRORS -- Procedure to print the centering task error messages.
+
+procedure ap_cerrrors (ap, ier)
+
+pointer	ap		# apphot pointer (unused)
+int	ier		# centering error code
+
+begin
+	switch (ier) {
+	case AP_CTR_NOAREA:
+	    call printf ("There are no pixels in the centering box.\n")
+        case AP_CTR_OUTOFBOUNDS:
+	    call printf (
+	        "The centering region is partially outside the image.\n")
+	case AP_CTR_NTOO_SMALL:
+	    call printf ("The centering box has too few points.\n")
+	case AP_CTR_SINGULAR:
+	    call printf ("The centering solution is singular.\n")
+	case AP_CTR_NOCONVERGE:
+	    call printf ("The centering algorithm does not converge.\n")
+	case AP_CTR_BADSHIFT:
+	    call printf ("The center shift is large.\n")
+	case AP_CTR_LOWSNRATIO:
+	    call printf ("The signal to noise ratio is low.\n")
+	case AP_CTR_BADDATA:
+	    call printf ("Bad data in the centering subraster.\n")
+	default:
+	    call printf ("")
+	}
+end
diff --git a/noao/digiphot/apphot/center/apcfree.x b/noao/digiphot/apphot/center/apcfree.x
new file mode 100644
index 00000000..7fb219c7
--- /dev/null
+++ b/noao/digiphot/apphot/center/apcfree.x
@@ -0,0 +1,46 @@
+include "../lib/apphotdef.h"
+include "../lib/centerdef.h"
+
+# APCFREE -- Procedure to free the centering data structure.
+
+procedure apcfree (ap)
+
+pointer	ap		# pointer to the apphot structure
+
+begin
+	if (ap == NULL)
+	    return
+	if (AP_NOISE(ap) != NULL)
+	    call ap_noisecls (ap)
+	if (AP_PCENTER(ap) != NULL)
+	    call ap_ctrcls (ap)
+	if (AP_PDISPLAY(ap) != NULL)
+	    call ap_dispcls (ap)
+	if (AP_IMBUF(ap) != NULL)
+	    call mfree (AP_IMBUF(ap), TY_REAL)
+	if (AP_MW(ap) != NULL)
+	    call mw_close (AP_MW(ap))
+	call mfree (ap, TY_STRUCT)
+end
+
+
+# AP_CTRCLS -- Procedure to close up the centering structure arrays.
+
+procedure ap_ctrcls (ap)
+
+pointer	ap			# pointer to apphot structure
+
+pointer	ctr
+
+begin
+	ctr = AP_PCENTER(ap)
+	if (ctr == NULL)
+	    return
+	if (AP_CTRPIX(ctr) != NULL)
+	    call mfree (AP_CTRPIX(ctr), TY_REAL)
+	if (AP_XCTRPIX(ctr) != NULL)
+	    call mfree (AP_XCTRPIX(ctr), TY_REAL)
+	if (AP_YCTRPIX(ctr) != NULL)
+	    call mfree (AP_YCTRPIX(ctr), TY_REAL)
+	call mfree (AP_PCENTER(ap), TY_STRUCT)
+end
diff --git a/noao/digiphot/apphot/center/apcinit.x b/noao/digiphot/apphot/center/apcinit.x
new file mode 100644
index 00000000..a9f38fd7
--- /dev/null
+++ b/noao/digiphot/apphot/center/apcinit.x
@@ -0,0 +1,94 @@
+include "../lib/apphotdef.h"
+include "../lib/centerdef.h"
+include "../lib/center.h"
+
+# APCINIT - Procedure to initialize the centering structure.
+
+procedure apcinit (ap, function, cbox, fwhmpsf, noise)
+
+pointer	ap		# pointer to the apphot structure
+int	function	# centering algorithm
+real	cbox		# half width of centering box
+real	fwhmpsf		# FWHM of the PSF
+int	noise		# noise function
+
+begin
+	# Allocate space.
+	call malloc (ap, LEN_APSTRUCT, TY_STRUCT)
+
+	# Set up the global apphot package parameters.
+	call ap_defsetup (ap, fwhmpsf)
+
+	# Setup the noise structure.
+	call ap_noisesetup (ap, noise)
+	    
+	# Initialize the centering parameters.
+	call ap_ctrsetup (ap, function, cbox)
+
+	# Display the options.
+	call ap_dispsetup (ap)
+
+	# Unused structures are set to null.
+	AP_PSKY(ap) = NULL
+	AP_PPSF(ap) = NULL
+	AP_PPHOT(ap) = NULL
+	AP_POLY(ap) = NULL
+	AP_RPROF(ap) = NULL
+end
+
+
+# AP_CTRSETUP -- Procedure to setup the centering array parameters.
+
+procedure ap_ctrsetup (ap, function, cbox)
+
+pointer	ap		# pointer to apphot structure
+int	function	# centering function
+real	cbox		# radius of centering aperture
+
+pointer	ctr
+
+begin
+	# Allocate space for the centering structure.
+	call malloc (AP_PCENTER(ap), LEN_CENSTRUCT, TY_STRUCT)
+	ctr = AP_PCENTER(ap)
+	AP_CXCUR(ctr) = INDEFR
+	AP_CYCUR(ctr) = INDEFR
+
+	AP_CENTERFUNCTION(ctr) = function
+	switch (function) {
+	case AP_CENTROID1D:
+	    call strcpy ("centroid", AP_CSTRING(ctr), SZ_FNAME)
+	case AP_GAUSS1D:
+	    call strcpy ("gauss", AP_CSTRING(ctr), SZ_FNAME)
+	case AP_NONE:
+	    call strcpy ("none", AP_CSTRING(ctr), SZ_FNAME)
+	case AP_OFILT1D:
+	    call strcpy ("ofilter", AP_CSTRING(ctr), SZ_FNAME)
+	default:
+	    AP_CENTERFUNCTION(ctr) = DEF_CENTERFUNCTION
+	    call strcpy ("centroid", AP_CSTRING(ctr), SZ_FNAME)
+	}
+	AP_CAPERT(ctr) = cbox
+	AP_CTHRESHOLD(ctr) = DEF_CTHRESHOLD
+	AP_MINSNRATIO(ctr) = DEF_MINSNRATIO
+	AP_CMAXITER(ctr) = DEF_CMAXITER
+	AP_MAXSHIFT(ctr) = DEF_MAXSHIFT
+	AP_CLEAN(ctr) = DEF_CLEAN
+	AP_RCLEAN(ctr) = DEF_RCLEAN
+	AP_RCLIP(ctr) = DEF_RCLIP
+	AP_SIGMACLEAN(ctr) = DEF_CLEANSIGMA
+
+	AP_NCTRPIX(ctr) = 0
+	AP_CTRPIX(ctr) = NULL
+	AP_XCTRPIX(ctr) = NULL
+	AP_YCTRPIX(ctr) = NULL
+
+	AP_OXINIT(ctr) = INDEFR
+	AP_OYINIT(ctr) = INDEFR
+	AP_XCENTER(ctr) = INDEFR
+	AP_YCENTER(ctr) = INDEFR
+	AP_OXCENTER(ctr) = INDEFR
+	AP_OYCENTER(ctr) = INDEFR
+	AP_XERR(ctr) = INDEFR
+	AP_YERR(ctr) = INDEFR
+end
diff --git a/noao/digiphot/apphot/center/apclean.x b/noao/digiphot/apphot/center/apclean.x
new file mode 100644
index 00000000..2e3a827a
--- /dev/null
+++ b/noao/digiphot/apphot/center/apclean.x
@@ -0,0 +1,152 @@
+include "../lib/apphot.h"
+include "../lib/noise.h"
+include "../lib/center.h"
+
+# AP_CLEAN -- Procedure to clean a data array.
+
+procedure ap_clean (ap, pix, nx, ny, xc, yc)
+
+pointer	ap		# pointer to the apphot structure
+real	pix[nx,ny]	# data array
+int	nx, ny		# size of subarray
+real	xc, yc		# center of subarray
+
+int	apstati()
+real	apstatr()
+
+begin
+	switch (apstati (ap, NOISEFUNCTION)) {
+	case AP_NCONSTANT:
+	    call ap_cclean (pix, nx, ny, xc, yc, apstatr (ap, SCALE) *
+	        apstatr (ap, RCLIP), apstatr (ap, SIGMACLEAN),
+		apstatr (ap, SKYSIGMA), apstatr (ap, MAXSHIFT))
+	case AP_NPOISSON:
+	    call ap_pclean (pix, nx, ny, xc, yc, apstatr (ap, SCALE) *
+		apstatr (ap, RCLEAN), apstatr (ap, SCALE) * apstatr (ap,
+		RCLIP), apstatr (ap, SIGMACLEAN), apstatr (ap, SKYSIGMA),
+		apstatr (ap, EPADU), apstatr (ap, MAXSHIFT))
+	default:
+	    return
+	}
+end
+
+
+# AP_CCLEAN -- Procedure to clean the subraster assuming the noise is
+# due to a constant sky sigma.
+
+procedure ap_cclean (pix, nx, ny, cxc, cyc, rclip, kclean, skysigma,
+    maxshift)
+
+real	pix[nx, ny]		# array of pixels
+int	nx, ny			# dimensions of the subarray
+real	cxc, cyc		# initial center
+real	rclip			# cleaning and clipping radius
+real	kclean			# k-sigma clipping factor
+real	skysigma		# maxshift
+real	maxshift		# maximum shift
+
+int	i, ii, ixc, j, jj, jyc, ijunk, jjunk
+real	mindat, maxdat, rclip2, r2, ksigma 
+
+begin
+	# Return if indef valued sigma or sigma <= 0.
+	if (IS_INDEFR(skysigma) || (skysigma <= 0.0))
+	    return
+
+	# Find the maximum pixel in the subarray and treat this point as
+	# the center of symmetry if it is less than maxshift from the
+	# initial center.
+
+	call ap_2dalimr (pix, nx, ny, mindat, maxdat, ijunk, jjunk, ixc, jyc)
+	if (abs (cxc - ixc) > maxshift || abs (cyc - jyc) > maxshift) {
+	    ixc = int (cxc)
+	    jyc = int (cyc)
+	}
+
+	# Clip.
+	rclip2 = rclip ** 2
+	ksigma = kclean * skysigma
+	do j = 1, ny {
+	    jj = 2 * jyc - j
+	    if (jj < 1 || jj > ny)
+		next
+	    do i = 1, ixc {
+		ii = 2 * ixc - i
+		if (ii < 1 || ii > nx)
+		    next
+		r2 = (i - ixc) ** 2 + (j - jyc) ** 2
+		if (r2 > rclip2) {
+		    if (pix[ii,jj] > pix[i,j] + ksigma)
+			pix[ii,jj] = pix[i,j]
+		    else if (pix[i,j] > pix[ii,jj] + ksigma)
+			pix[i,j] = pix[ii,jj]
+		}
+	    }
+	}
+end
+
+
+# AP_PCLEAN -- Procedure to clean the subraster assuming that the noise is
+# due to a constant sky value plus poisson noise.
+
+procedure ap_pclean (pix, nx, ny, cxc, cyc, rclean, rclip, kclean, skysigma,
+    padu, maxshift)
+
+real	pix[nx, ny]		# array of pixels
+int	nx, ny			# dimensions of the subarray
+real	cxc, cyc		# initial center
+real	rclean, rclip		# cleaning and clipping radius
+real	kclean			# k-sigma clipping factor
+real	skysigma		# maxshift
+real	padu			# photons per adu
+real	maxshift		# maximum shift
+
+int	i, ii, ixc, j, jj, jyc, ijunk, jjunk
+real	mindat, maxdat, rclean2, rclip2, r2, ksigma, ksigma2 
+
+begin
+	# Return if indef-valued sigma.
+	if (IS_INDEFR(skysigma))
+	    return
+
+	# Find the maximum pixel in the subarray and treat this point as
+	# the center of symmetry if it is less than maxshift from the
+	# initial center.
+
+	call ap_2dalimr (pix, nx, ny, mindat, maxdat, ijunk, jjunk, ixc, jyc)
+	if (abs (cxc - ixc) > maxshift || abs (cyc - jyc) > maxshift) {
+	    ixc = int (cxc)
+	    jyc = int (cyc)
+	}
+
+	# Clip.
+	rclean2 = rclean ** 2
+	rclip2 = rclip ** 2
+	ksigma = kclean * skysigma
+	ksigma2 = ksigma ** 2
+	do j = 1, ny {
+	    jj = 2 * jyc - j
+	    if (jj < 1 || jj > ny)
+		next
+	    do i = 1, ixc {
+		ii = 2 * ixc - i
+		if (ii < 1 || ii > nx)
+		    next
+		r2 = (i - ixc) ** 2 + (j - jyc) ** 2
+		if (r2 > rclean2 && r2 <= rclip2) {
+		    if (pix[ii,jj] > (pix[i,j] + sqrt (pix[i,j] / padu +
+		        ksigma2)))
+			pix[ii,jj] = pix[i,j]
+		    else if (pix[i,j] > (pix[ii,jj] + sqrt (pix[ii,jj] / padu +
+		        ksigma2)))
+			pix[i,j] = pix[ii,jj]
+		}
+		if (r2 > rclip2) {
+		    if (pix[ii,jj] > pix[i,j] + ksigma)
+			pix[ii,jj] = pix[i,j]
+		    else if (pix[i,j] > pix[ii,jj] + ksigma)
+			pix[i,j] = pix[ii,jj]
+		}
+	    }
+	}
+end
diff --git a/noao/digiphot/apphot/center/apcplot.x b/noao/digiphot/apphot/center/apcplot.x
new file mode 100644
index 00000000..6e2845c3
--- /dev/null
+++ b/noao/digiphot/apphot/center/apcplot.x
@@ -0,0 +1,304 @@
+include <mach.h>
+include <gset.h>
+include <pkg/gtools.h>
+include "../lib/apphotdef.h"
+include "../lib/apphot.h"
+include "../lib/noise.h"
+include "../lib/centerdef.h"
+include "../lib/center.h"
+
+# AP_CPLOT -- Procedure to compute radial profile plots for the centering
+# routine.
+
+procedure ap_cplot (ap, sid, gd, makeplot)
+
+pointer	ap		# the pointer to the apphot structure
+int	sid		# the output sequence number of the star
+pointer	gd		# the graphics stream
+int	makeplot	# the make a plot ?
+
+int	nx, ny
+pointer	ctr, sp, str, r, gt
+real	xcenter, ycenter, xc, yc, rmin, rmax, imin, imax
+real	u1, u2, v1, v2, x1, x2, y1, y2
+int	apstati()
+pointer	ap_gtinit()
+real	apstatr()
+
+begin
+	# Check for enabled graphics stream.
+	if (gd == NULL || makeplot == NO)
+	    return
+
+	# Check for defined center.
+	xcenter = apstatr (ap, XCENTER)
+	ycenter = apstatr (ap, YCENTER)
+	if (IS_INDEFR(xcenter) || IS_INDEFR(ycenter))
+	    return
+
+	# Check for a centering algorithm.
+	if (apstati (ap, CENTERFUNCTION) == AP_NONE)
+	    return
+
+	# Get the pixel buffer parameters.
+	ctr = AP_PCENTER(ap)
+	nx = AP_CNX(ctr)
+	ny = AP_CNY(ctr)
+	if (AP_CTRPIX(ctr) == NULL || nx <= 0 || ny <= 0)
+	    return
+	xc = AP_CXC(ctr) + (xcenter - apstatr (ap, CXCUR))
+	yc = AP_CYC(ctr) + (ycenter - apstatr (ap, CYCUR))
+
+	# Allocate working space.
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+	call salloc (r, nx * ny, TY_REAL)
+
+	# Compute the radii and the plot limits.
+	call ap_ijtor2 (Memr[r], nx, ny, xc, yc)
+	call alimr (Memr[r], nx * ny, rmin, rmax)
+	call alimr (Memr[AP_CTRPIX(ctr)], nx * ny, imin, imax)
+
+	# Reactivate the work station.
+	call greactivate (gd, 0)
+
+	# Save old viewport and window coordinates.
+	call ggview (gd, u1, u2, v1, v2)
+	call ggwind (gd, x1, x2, y1, y2)
+
+	# Initialize the plot.
+	call apstats (ap, IMROOT, Memc[str], SZ_LINE)
+	call sprintf (Memc[str], SZ_LINE, "%s  Star %d")
+	    call pargstr (Memc[str])
+	    call pargi (sid)
+	gt = ap_gtinit (Memc[str], apstatr (ap,OXINIT), apstatr (ap, OYINIT))
+
+	# Make the plot.
+	call gclear (gd)
+	call ap_cpset (gd, gt, ap, rmin, rmax, imin, imax)
+	if (apstati (ap, POSITIVE) == YES)
+	    call ap_plotpts (gd, gt, Memr[r], Memr[AP_CTRPIX(ctr)], nx * ny,
+	        rmin, rmax, imin + EPSILONR, imax, "plus")
+	else
+	    call ap_plotpts (gd, gt, Memr[r], Memr[AP_CTRPIX(ctr)], nx * ny,
+	        rmin, rmax, imin, imax - EPSILONR, "plus")
+	call ap_cpreset (gd, gt, ap, rmin, rmax, imin, imax)
+	call ap_cpannotate (gd, ap)
+
+	# Restore the viewport and window coordinates.
+	call gsview (gd, u1, u2, v1, v2)
+	call gswind (gd, x1, x2, y1, y2)
+
+	# Free the space.
+	call ap_gtfree (gt)
+	call gdeactivate (gd, 0)
+	call sfree (sp)
+end
+
+
+# AP_CPSET -- Procedure to set up the parameters for the center radial profile
+# plot.
+
+procedure ap_cpset (gd, gt, ap, xmin, xmax, ymin, ymax)
+
+pointer	gd		# the graphics stream
+pointer	gt		# the gtools pointer
+pointer	ap		# the apphot pointer
+real	xmin, xmax	# the minimum and maximum radial distance
+real	ymin, ymax	# the minimum and maximum of y axis (ymin not used)
+
+int	fd
+pointer	sp, str, title
+real	scale, aspect, datalimit, skysigma, threshold, vx1, vx2, vy1, vy2
+real	cthreshold
+int	stropen(), apstati()
+real	apstatr(), gstatr()
+
+begin
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+	call salloc (title, SZ_LINE, TY_CHAR)
+
+	# Encode the parameter string.
+	fd = stropen (Memc[str], SZ_LINE, WRITE_ONLY)
+	call sysid (Memc[title], SZ_LINE)
+	call fprintf (fd, "%s\n")
+	    call pargstr (Memc[title])
+	call fprintf (fd,
+	    "Center: xc=%0.2f yc=%0.2f xerr=%0.2f yerr=%0.2f\n")
+	    call pargr (apstatr (ap, OXCENTER))
+	    call pargr (apstatr (ap, OYCENTER))
+	    call pargr (apstatr (ap, XERR))
+	    call pargr (apstatr (ap, YERR))
+	call gt_gets (gt, GTTITLE, Memc[title], SZ_LINE)
+	call fprintf (fd, "%s\n")
+	    call pargstr (Memc[title])
+	call strclose (fd)
+
+	# Store some default plotting parameters.
+	scale = apstatr (ap, SCALE)
+	aspect = gstatr (gd, G_ASPECT)
+	call gsetr (gd, G_ASPECT, 0.75)
+
+	# Set the labels and window.
+	datalimit = apstatr (ap, CDATALIMIT)
+	skysigma = apstatr (ap, SKYSIGMA)
+	cthreshold = apstatr (ap, CTHRESHOLD)
+	if (IS_INDEFR(skysigma) || IS_INDEFR(cthreshold))
+	    threshold = 0.0
+	else
+	    threshold = cthreshold * skysigma
+	if (apstati (ap, POSITIVE) == YES) {
+	    call gseti (gd, G_XDRAWAXES, 2)
+	    call gswind (gd, xmin / scale, xmax / scale, datalimit + threshold,
+		datalimit + threshold + ymax)
+	    call glabax (gd, Memc[str], "", "Intensity")
+	    call gseti (gd, G_YDRAWAXES, 0)
+	    call gseti (gd, G_XDRAWAXES, 1)
+	    call ggview (gd, vx1, vx2, vy1, vy2)
+	    call gswind (gd, vx1, vx2, vy1, vy2)
+	    call gswind (gd, xmin, xmax, datalimit + threshold,
+		datalimit + threshold + ymax)
+	    call glabax (gd, "",
+	        "Radial Distance (lower-pixels, upper-scale units)", "")
+	} else {
+	    call gseti (gd, G_XDRAWAXES, 2)
+	    call gswind (gd, xmin / scale, xmax / scale, datalimit -
+	        threshold - ymax, datalimit - threshold)
+	    call glabax (gd, Memc[str], "", "Intensity")
+	    call gseti (gd, G_YDRAWAXES, 0)
+	    call gseti (gd, G_XDRAWAXES, 1)
+	    call ggview (gd, vx1, vx2, vy1, vy2)
+	    call gswind (gd, vx1, vx2, vy1, vy2)
+	    call gswind (gd,  xmin,  xmax, datalimit - threshold - ymax,
+	        datalimit - threshold)
+	    call glabax (gd, "",
+	        "Radial Distance (lower-pixels, upper-scale units)", "")
+	}
+
+	# Set the window up for plotting the data.
+	call gseti (gd, G_YDRAWAXES, 3)
+	call gseti (gd, G_XDRAWAXES, 3)
+	call gsetr (gd, G_ASPECT, aspect)
+	call gt_sets (gt, GTTYPE, "mark")
+	call gt_setr (gt, GTXMIN, xmin)
+	call gt_setr (gt, GTXMAX, xmax)
+	if (apstati (ap, POSITIVE) == YES) {
+	    call gt_setr (gt, GTYMIN, ymin)
+	    call gt_setr (gt, GTYMAX, ymax)
+	} else {
+	    call gt_setr (gt, GTYMIN, ymin)
+	    call gt_setr (gt, GTYMAX, ymax)
+	}
+	call gt_swind (gd, gt)
+
+	call sfree (sp)
+end
+
+
+# AP_CPANNOTATE -- Procedure to annotate the radial plot in center.
+
+procedure ap_cpannotate (gd, ap)
+
+pointer	gd		# the graphics stream
+pointer	ap		# the apphot structure
+
+pointer	sp, str
+real	fwhmpsf, capert, datalimit, threshold, skysigma, cthreshold
+real	xmin, xmax, ymin, ymax
+int	apstati()
+real	apstatr()
+
+begin
+	# Allocate working space.
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+	call gseti (gd, G_PLTYPE, GL_DASHED)
+	call ggwind (gd, xmin, xmax, ymin, ymax)
+
+	fwhmpsf = 0.5 * apstatr (ap, FWHMPSF) * apstatr (ap, SCALE)
+	capert = 2.0 * fwhmpsf * apstatr (ap, CAPERT)
+	datalimit = apstatr (ap, CDATALIMIT)
+	skysigma = apstatr (ap, SKYSIGMA)
+	cthreshold = apstatr (ap, CTHRESHOLD)
+	if (IS_INDEFR(skysigma) || IS_INDEFR(cthreshold))
+	    threshold = 0.0
+	else
+	    threshold = cthreshold * skysigma
+	if (apstati (ap, POSITIVE) == YES)
+	    threshold = datalimit + threshold
+	else
+	    threshold = datalimit - threshold
+
+	# Plot the full width half maximum of the radial profile.
+	if (fwhmpsf >= xmin && fwhmpsf <= xmax) {
+	    call gamove (gd, fwhmpsf, ymin)
+	    call gadraw (gd, fwhmpsf, ymax)
+	    call sprintf (Memc[str], SZ_LINE, "hwhm = %0.2f")
+	        call pargr (fwhmpsf)
+	    call gtext (gd, fwhmpsf, ymax, Memc[str], "q=h;u=180;v=t;p=r")
+	}
+
+	# Mark the centering aperture.
+	if (capert >= xmin && capert <= xmax) {
+	    call gamove (gd, capert, ymin)
+	    call gadraw (gd, capert, ymax)
+	    call sprintf (Memc[str], SZ_LINE, "cbox half-width = %0.2f")
+	        call pargr (capert)
+	    call gtext (gd, capert, ymax, Memc[str], "q=h;u=180;v=t;p=r")
+	}
+
+	# Mark the threshold level for centering.
+        call sprintf (Memc[str], SZ_LINE, "threshold = %g")
+	    call pargr (threshold)
+	call gtext (gd, xmin, ymin, Memc[str], "q=h")
+
+	# Mark the sky sigma if defined.
+	if (! IS_INDEFR(skysigma) && skysigma >= ymin && skysigma <= ymax) {
+	    call gmark (gd, (xmin + xmax) / 2.0, (ymin + ymax) / 2.0, 
+		GM_VEBAR, -0.25, -skysigma)
+	    call sprintf (Memc[str], SZ_LINE, "sigma = %g")
+		call pargr (skysigma)
+	    call gtext (gd, (xmin + xmax) / 2.0, (ymin + ymax + skysigma) / 2.0,
+		Memc[str], "q=h;h=c")
+	}
+
+	call sfree (sp)
+end
+
+
+# AP_CPRSET -- Procedure to reset the plot window after the data points
+# have been plotted.
+
+procedure ap_cpreset (gd, gt, ap, xmin, xmax, ymin, ymax)
+
+pointer	gd		# the graphics stream
+pointer	gt		# the gtools pointer
+pointer	ap		# the apphot pointer
+real	xmin, xmax	# the minimum and maximum radial distance
+real	ymin, ymax	# the minimum and maximum of the y axis (ymin not used)
+
+real	datalimit, skysigma, threshold, cthreshold
+int	apstati()
+real	apstatr()
+
+begin
+	# Set the data window.
+	datalimit = apstatr (ap, CDATALIMIT)
+	skysigma = apstatr (ap, SKYSIGMA)
+	cthreshold = apstatr (ap, CTHRESHOLD)
+	if (IS_INDEFR(skysigma) || IS_INDEFR(cthreshold))
+	    threshold = 0.0
+	else
+	    threshold = cthreshold * skysigma
+	call gt_setr (gt, GTXMIN, xmin)
+	call gt_setr (gt, GTXMAX, xmax)
+	if (apstati (ap, POSITIVE) == YES) {
+	    call gt_setr (gt, GTYMIN, datalimit + threshold)
+	    call gt_setr (gt, GTYMAX, ymax + datalimit + threshold)
+	} else {
+	    call gt_setr (gt, GTYMIN, datalimit - ymax - threshold)
+	    call gt_setr (gt, GTYMAX, datalimit - threshold)
+	}
+	call gt_swind (gd, gt)
+end
diff --git a/noao/digiphot/apphot/center/apcpshow.x b/noao/digiphot/apphot/center/apcpshow.x
new file mode 100644
index 00000000..79fdf43e
--- /dev/null
+++ b/noao/digiphot/apphot/center/apcpshow.x
@@ -0,0 +1,66 @@
+include "../lib/center.h"
+include "../lib/display.h"
+
+# AP_CPSHOW -- Procedure to display the current centering algorithm
+# parameters.
+
+procedure ap_cpshow (ap)
+
+pointer	ap	# pointer to the apphot strucuture
+
+pointer	sp, str
+bool	itob()
+int	apstati()
+real	apstatr()
+
+begin
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	# Print the centering algorithm.
+	call printf ("Centering Parameters\n")
+	call apstats (ap, CSTRING, Memc[str], SZ_FNAME)
+	call printf ("    %s = %s %s\n")
+	    call pargstr (KY_CSTRING)
+	    call pargstr (Memc[str])
+	    call pargstr (UN_CALGORITHM)
+
+	# Print the rest of the centering parameters.
+	call printf ("    %s = %g %s    %s = %g %s    %s = %g\n")
+	    call pargstr (KY_CAPERT)
+	    call pargr (2.0 * apstatr (ap, CAPERT))
+	    call pargstr (UN_CSCALEUNIT)
+	    call pargstr (KY_CTHRESHOLD)
+	    call pargr (apstatr (ap, CTHRESHOLD))
+	    call pargstr (UN_CSIGMA)
+	    call pargstr (KY_MINSNRATIO)
+	    call pargr (apstatr (ap, MINSNRATIO))
+
+	call printf ("    %s = %d    %s = %g %s\n")
+	    call pargstr (KY_CMAXITER)
+	    call pargi (apstati (ap, CMAXITER))
+	    call pargstr (KY_MAXSHIFT)
+	    call pargr (apstatr (ap, MAXSHIFT))
+	    call pargstr (UN_CNUMBER)
+
+	call printf ("    %s = %b    %s = %g %s\n")
+	    call pargstr (KY_CLEAN)
+	    call pargb (itob (apstati (ap, CLEAN)))
+	    call pargstr (KY_SIGMACLEAN)
+	    call pargr (apstatr (ap, SIGMACLEAN))
+	    call pargstr (UN_CSIGMA)
+
+	call printf ("    %s = %g %s    %s = %g %s\n")
+	    call pargstr (KY_RCLEAN)
+	    call pargr (apstatr (ap, RCLEAN))
+	    call pargstr (UN_CSCALEUNIT)
+	    call pargstr (KY_RCLIP)
+	    call pargr (apstatr (ap, RCLIP))
+	    call pargstr (UN_CSCALEUNIT)
+
+	call printf ("    %s = %b\n")
+	    call pargstr (KY_MKCENTER)
+	    call pargb (itob (apstati (ap, MKCENTER)))
+
+	call sfree (sp)
+end
diff --git a/noao/digiphot/apphot/center/apcradsetup.x b/noao/digiphot/apphot/center/apcradsetup.x
new file mode 100644
index 00000000..8902c5ef
--- /dev/null
+++ b/noao/digiphot/apphot/center/apcradsetup.x
@@ -0,0 +1,99 @@
+include "../lib/display.h"
+
+define	HELPFILE	"apphot$center/icenter.key"
+
+# AP_CRADSETUP -- Procedure to set up apphot interactively
+
+procedure ap_cradsetup (ap, im, wx, wy, gd, out, stid)
+
+pointer	ap			# pointer to apphot structure
+pointer	im			# pointero to the IRAF image
+real	wx, wy			# cursor coordinates
+pointer	gd			# pointer to graphics stream
+int	out			# output file descriptor
+int	stid			# output file sequence number
+
+int	cier, wcs, key
+pointer	sp, cmd
+real	xcenter, ycenter, xc, yc, rmin, rmax, imin, imax, rval
+real	u1, u2, v1, v2, x1, x2, y1, y2
+
+int	ap_showplot(), apfitcenter(), apstati(), clgcur()
+real	ap_cfwhmpsf(), ap_ccapert(), ap_csigma()
+real	ap_cdatamin(), ap_cdatamax(), ap_crclean(), ap_crclip()
+
+begin
+	# Check for open graphics stream.
+	if (gd == NULL)
+	    return
+	call greactivate (gd, 0)
+
+	# Store old viewport and window coordinates.
+	call ggview (gd, u1, u2, v1, v2)
+	call ggwind (gd, x1, x2, y1, y2)
+
+	# Make the plot.
+	if (ap_showplot (ap, im, wx, wy, gd, xcenter, ycenter, rmin, rmax,
+	    imin, imax) == ERR) {
+	    call gdeactivate (gd, 0)
+	    return
+	}
+
+	# Allocate temporary space.
+	call smark (sp)
+	call salloc (cmd, SZ_LINE, TY_CHAR)
+
+	call printf (
+	"Waiting for setup menu command (?=help, v=default setup, q=quit):\n")
+
+	while (clgcur ("gcommands", xc, yc, wcs, key, Memc[cmd],
+	    SZ_LINE) != EOF) {
+
+	switch (key) {
+	    case 'q':
+	        break
+	    case '?':
+		call gpagefile (gd, HELPFILE, "")
+	    case 'f':
+	        rval = ap_cfwhmpsf (ap, gd, out, stid, rmin, rmax, imin, imax)
+	    case 's':
+	        rval = ap_csigma (ap, gd, out, stid, rmin, rmax, imin, imax)
+	    case 'l':
+	        rval = ap_cdatamin (ap, gd, out, stid, rmin, rmax, imin, imax)
+	    case 'u':
+	        rval = ap_cdatamax (ap, gd, out, stid, rmin, rmax, imin, imax)
+	    case 'c':
+		rval = ap_ccapert (ap, gd, out, stid, rmin, rmax, imin, imax)
+	    case 'n':
+		rval = ap_crclean (ap, gd, out, stid, rmin, rmax, imin, imax)
+	    case 'p':
+		rval = ap_crclip (ap, gd, out, stid, rmin, rmax, imin, imax)
+	    case 'v':
+	        rval = ap_cfwhmpsf (ap, gd, out, stid, rmin, rmax, imin, imax)
+		rval = ap_ccapert (ap, gd, out, stid, rmin, rmax, imin, imax)
+	        rval = ap_csigma (ap, gd, out, stid, rmin, rmax, imin, imax)
+	        #rval = ap_ccthresh (ap, gd, out, stid, rmin, rmax, imin, imax)
+	    default:
+		call printf ("Unknown or ambiguous keystroke command\007\n")
+	    }
+	    call printf (
+	 "Waiting for setup menu command (?=help, v=default setup, q=quit):\n")
+	}
+
+	# Set up interactively.
+	call printf (
+	    "Interactive setup is complete [Type w to save parameters].\n")
+
+	# Restore old window and viewport coordinates.
+	call gsview (gd, u1, u2, v1, v2)
+	call gswind (gd, x1, x2, y1, y2)
+
+	# Free plots and working space.
+	call sfree (sp)
+	call gdeactivate (gd, 0)
+
+	# Fit the center and plot the results.
+	cier = apfitcenter (ap, im, xcenter, ycenter)
+	call ap_cplot (ap, 0, gd, apstati (ap, RADPLOTS))
+	call ap_qcenter (ap, cier)
+end
diff --git a/noao/digiphot/apphot/center/apcshow.x b/noao/digiphot/apphot/center/apcshow.x
new file mode 100644
index 00000000..8cedbfc2
--- /dev/null
+++ b/noao/digiphot/apphot/center/apcshow.x
@@ -0,0 +1,19 @@
+include "../lib/display.h"
+
+# AP_CSHOW -- Procedure to display the current centering parameters.
+
+procedure ap_cshow (ap)
+
+pointer	ap	# pointer to the apphot strucuture
+
+bool	itob()
+int	apstati()
+
+begin
+	call ap_nshow (ap)
+	call printf ("\n")
+	call ap_cpshow (ap)
+	call printf ("    %s = %b\n")
+	    call pargstr (KY_RADPLOTS)
+	    call pargb (itob (apstati (ap, RADPLOTS)))
+end
diff --git a/noao/digiphot/apphot/center/apcsnratio.x b/noao/digiphot/apphot/center/apcsnratio.x
new file mode 100644
index 00000000..0fd0e20d
--- /dev/null
+++ b/noao/digiphot/apphot/center/apcsnratio.x
@@ -0,0 +1,92 @@
+include <mach.h>
+include "../lib/center.h"
+include "../lib/noise.h"
+
+# AP_CSNRATIO -- Procedure to estimate the signal to noise ratio for centering.
+
+real procedure ap_csnratio (ap, array, nx, ny, threshold)
+
+pointer	ap		# apphot structure
+real	array[nx,ny]	# data array
+int	nx, ny		# subraster dimensions
+real	threshold	# threshold value for snr computation
+
+int	apstati()
+real	apstatr(), ap_cratio(), ap_pratio()
+
+begin
+	switch (apstati (ap, NOISEFUNCTION)) {
+	case AP_NCONSTANT:
+	    return (ap_cratio (array, nx, ny, threshold, apstatr (ap,
+	        SKYSIGMA)))
+	case AP_NPOISSON:
+	    return (ap_pratio (array, nx, ny, threshold, apstatr (ap,
+	        SKYSIGMA), apstatr (ap, EPADU)))
+	default:
+	    return (MAX_REAL)
+	}
+end
+
+
+# AP_CRATIO -- Procedure to estimate the signal to noise ratio in the
+# centering aperture assuming that the noise is due to a constant sky
+# sigma. These computations are approximate only.
+
+real procedure ap_cratio (array, nx, ny, sky, noise)
+
+real	array[nx,ny]		# centering buffer
+int	nx, ny			# dimensions of the centering buffer
+real	sky			# mean sky value of the data in ADU
+real	noise			# estimate of sky noise in ADU
+
+int	npts
+real	signal, tnoise
+real	asumr()
+
+begin
+	npts = nx * ny
+	signal = asumr (array, npts) - npts * sky
+	if (IS_INDEFR(noise))
+	    tnoise = 0.0
+	else
+	    tnoise = sqrt (npts * noise ** 2)
+	if (signal <= 0.0)
+	    return (0.0)
+	else if (tnoise <= 0.0)
+	    return (MAX_REAL)
+	else
+	    return (signal / tnoise)
+end
+
+
+# AP_PRATIO -- Procedure to estimate the signal to noise ratio in the
+# centering aperture assuming the noise is due to a constant sky sigma
+# and poisson statistics in the image. These computations are
+# approximate only.
+
+real procedure ap_pratio (array, nx, ny, sky, noise, padu)
+
+real	array[nx,ny]		# centering buffer
+int	nx, ny			# dimensions of the centering buffer
+real	sky			# mean sky value of the data in ADU
+real	noise			# estimate of sky noise in ADU
+real	padu			# photons per adu
+
+int	npts
+real	signal, tnoise
+real	asumr()
+
+begin
+	npts = nx * ny
+	signal = asumr (array, npts) - npts * sky
+	if (IS_INDEFR(noise))
+	    tnoise = sqrt (abs (signal / padu))
+	else
+	    tnoise = sqrt (abs (signal / padu) + npts * noise ** 2)
+	if (signal <= 0.0)
+	    return (0.0)
+	else if (tnoise <= 0.0)
+	    return (MAX_REAL)
+	else
+	    return (signal / tnoise)
+end
diff --git a/noao/digiphot/apphot/center/apctr1d.x b/noao/digiphot/apphot/center/apctr1d.x
new file mode 100644
index 00000000..10ee109e
--- /dev/null
+++ b/noao/digiphot/apphot/center/apctr1d.x
@@ -0,0 +1,113 @@
+include <mach.h>
+include "../lib/center.h"
+
+# AP_CTR1D -- Procedure to compute the center from the 1D marginal
+# distributions.
+
+int procedure ap_ctr1d (ctrpix, nx, ny, norm, xc, yc, xerr, yerr)
+
+real	ctrpix[nx, ny]		# object to be centered
+int	nx, ny			# dimensions of subarray
+real	norm			# the normalization factor
+real	xc, yc			# computed centers
+real	xerr, yerr		# estimate of centering error
+
+pointer	sp, xm, ym
+
+begin
+	call smark (sp)
+	call salloc (xm, nx, TY_REAL)
+	call salloc (ym, ny, TY_REAL)
+
+	# Compute the marginal distributions.
+	call ap_mkmarg (ctrpix, Memr[xm], Memr[ym], nx, ny)
+	call adivkr (Memr[xm], real (ny), Memr[xm], nx)
+	call adivkr (Memr[ym], real (nx), Memr[ym], ny)
+
+	# Get the centers and errors.
+	call ap_cmarg (Memr[xm], nx, norm, xc, xerr)
+	call ap_cmarg (Memr[ym], ny, norm, yc, yerr)
+
+	call sfree (sp)
+	return (AP_OK)
+end
+
+
+# AP_CMARG -- Compute the center and estimate its error given the
+# marginal distribution and the number of points.
+
+procedure ap_cmarg (a, npts, norm, xc, err)
+
+real	a[npts]		# array
+int	npts		# number of points
+real	norm		# the normalization factor
+real	xc		# center value
+real	err		# error
+
+int	i, npos
+real	val, sumi, sumix, sumix2
+
+begin
+	# Initialize.
+	npos = 0
+	sumi = 0.0
+	sumix = 0.0
+	sumix2 = 0.0
+
+	# Accumulate the sums.
+	do i = 1, npts {
+	    val = a[i]
+	    if (val > 0.0)
+	        npos = npos + 1
+	    sumi = sumi + val
+	    sumix = sumix + val * i
+	    sumix2 = sumix2 + val * i ** 2
+	}
+
+	# Compute the position and the error.
+	if (npos <= 0) {
+	    xc =  (1.0 + npts) / 2.0
+	    err = INDEFR
+	} else {
+	    xc = sumix / sumi
+	    err = (sumix2 / sumi - xc * xc)
+	    if (err <= 0.0) {
+		err = 0.0
+	    } else {
+	        err = sqrt (err / (sumi * norm))
+		if (err > real (npts))
+		    err = INDEFR
+	    }
+	}
+end
+
+
+# AP_MKMARG -- Accumulate the marginal distributions.
+
+procedure ap_mkmarg (ctrpix, xm, ym, nx, ny)
+
+real	ctrpix[nx, ny]		# pixels
+real	xm[nx]			# x marginal distribution
+real	ym[ny]			# y marginal distribution
+int	nx, ny			# dimensions of array
+
+int	i, j
+real	sum
+
+begin
+	# Compute the x marginal.
+	do i = 1, nx {
+	    sum = 0.0
+	    do j = 1, ny
+		sum = sum + ctrpix[i,j]
+	    xm[i] = sum
+	}
+
+	# Compute the y marginal.
+	do j = 1, ny {
+	    sum = 0.0
+	    do i = 1, nx
+		sum = sum + ctrpix[i,j]
+	    ym[j] = sum
+	}
+end
diff --git a/noao/digiphot/apphot/center/apctrbuf.x b/noao/digiphot/apphot/center/apctrbuf.x
new file mode 100644
index 00000000..0a8053b8
--- /dev/null
+++ b/noao/digiphot/apphot/center/apctrbuf.x
@@ -0,0 +1,115 @@
+include <imhdr.h>
+include <math.h>
+include <mach.h>
+include "../lib/apphotdef.h"
+include "../lib/centerdef.h"
+include "../lib/center.h"
+
+# APCTRBUF -- Procedure to fetch the center pixels given the pointer to the
+# IRAF image, the coordinates of the initial center and the width of the APPHOT
+# centering box.
+
+int procedure apctrbuf (ap, im, wx, wy)
+
+pointer	ap		# pointer to apphot structure
+pointer	im		# pointer to the IRAF image
+real	wx, wy		# center coordinates
+
+int	icpix
+pointer	ctr
+real	cpix, gdatamin, gdatamax, datamin, datamax
+pointer	ap_ctrpix()
+
+begin
+	# Get pointer to centering structure.
+	ctr = AP_PCENTER(ap)
+
+	# Check for 0 sized aperture.
+	if (AP_CAPERT(ctr) <= 0.0)
+	    return (AP_CTR_NOAREA)
+
+	# Get the centering buffer of pixels.
+	cpix = max (1.0, AP_CAPERT(ctr) * AP_SCALE(ap))
+	icpix = 2 * int (cpix) + 1
+	if (AP_CTRPIX(ctr) != NULL)
+	    call mfree (AP_CTRPIX(ctr), TY_REAL)
+	AP_CTRPIX(ctr) = ap_ctrpix (im, wx, wy, icpix,
+	    AP_CXC(ctr), AP_CYC(ctr), AP_CNX(ctr), AP_CNY(ctr))
+	if (AP_CTRPIX(ctr) == NULL)
+	    return (AP_CTR_NOAREA)
+
+	# Compute the data limits.
+	if (IS_INDEFR(AP_DATAMIN(ap)))
+	    gdatamin = -MAX_REAL
+	else
+	    gdatamin = AP_DATAMIN(ap)
+	if (IS_INDEFR(AP_DATAMAX(ap)))
+	    gdatamax = MAX_REAL
+	else
+	    gdatamax = AP_DATAMAX(ap)
+	call alimr (Memr[AP_CTRPIX(ctr)], AP_CNX(ctr) * AP_CNY(ctr),
+	    datamin, datamax)
+
+	if (datamin < gdatamin || datamax > gdatamax)
+	    return (AP_CTR_BADDATA)
+	else if (AP_CNX(ctr) < icpix || AP_CNY(ctr) < icpix)
+	    return (AP_CTR_OUTOFBOUNDS)
+	else
+	    return (AP_OK)
+end
+
+
+# AP_CTRPIX -- Procedure to fetch the pixels to be used for centering.
+
+pointer procedure ap_ctrpix (im, wx, wy, capert, xc, yc, nx, ny)
+
+pointer	im		# pointer to IRAF image
+real	wx, wy		# center of subraster to be extracted
+int	capert		# width of subraster to be extracted
+real	xc, yc		# center of extracted subraster
+int	nx, ny		# dimensions of extracted subraster
+
+int	i, ncols, nlines, c1, c2, l1, l2, half_capert
+pointer	buf, lbuf
+real	xc1, xc2, xl1, xl2
+pointer	imgs2r()
+
+begin
+	# Check for nonsensical input.
+	half_capert = (capert - 1) / 2
+	if (half_capert <= 0)
+	    return (NULL)
+
+	# Test for out of bounds pixels
+	ncols = IM_LEN(im,1)
+	nlines = IM_LEN(im,2)
+	xc1 = wx - half_capert
+	xc2 = wx + half_capert
+	xl1 = wy - half_capert
+	xl2 = wy + half_capert
+	if (xc1 > real (ncols) || xc2 < 1.0 || xl1 > real (nlines) || xl2 < 1.0)
+	    return (NULL)
+
+	# Get column and line limits, dimensions and center of subraster.
+	c1 = max (1.0, min (real (ncols), xc1)) + 0.5
+	c2 = min (real (ncols), max (1.0, xc2)) + 0.5
+	l1 = max (1.0, min (real (nlines), xl1)) + 0.5
+	l2 = min (real (nlines), max (1.0, xl2)) + 0.5
+	nx = c2 - c1 + 1
+	ny = l2 - l1 + 1
+	xc = wx - c1 + 1
+	yc = wy - l1 + 1
+
+	# Get pixels and return.
+	if (nx < 1 && ny < 1)
+	    return (NULL)
+	else {
+	    call malloc (buf, nx * ny, TY_REAL)
+	    lbuf = buf
+	    do i = l1, l2 {
+	        call amovr (Memr[imgs2r (im, c1, c2, i, i)], Memr[lbuf], nx)
+	        lbuf = lbuf + nx
+	    }
+	    return (buf)
+	}
+end
diff --git a/noao/digiphot/apphot/center/apfitcen.x b/noao/digiphot/apphot/center/apfitcen.x
new file mode 100644
index 00000000..bd7224e8
--- /dev/null
+++ b/noao/digiphot/apphot/center/apfitcen.x
@@ -0,0 +1,291 @@
+include <mach.h>
+include "../lib/apphotdef.h"
+include "../lib/apphot.h"
+include "../lib/noisedef.h"
+include "../lib/centerdef.h"
+include "../lib/center.h"
+
+define	CONVERT		.424660900	# conversion from fwhmpsf to sigma
+
+# APFITCENTER -- Procedure to fit the centers using either 1) the intensity
+# weighted mean of the marginal distributions, 2) a Gaussian fit to the
+# marginals assuming a fixed sigma or 3) a simplified version of the optimal
+# filtering method using a Gaussian of fixed sigma to model the profile.
+
+int procedure apfitcenter (ap, im, wx, wy)
+
+pointer	ap		# pointer to the apphot structure
+pointer	im		# pointer to the IRAF image
+real	wx, wy		# object coordinates
+
+int	i, niter, ier, fier, lowsnr
+pointer	ctr, nse
+real	cthreshold, datamin, datamax, owx, owy, xshift, yshift, med
+int	apctrbuf(), ap_ctr1d(), ap_mctr1d(), ap_gctr1d(), ap_lgctr1d()
+real	ap_csnratio(), amedr()
+
+begin
+	ctr = AP_PCENTER(ap)
+	nse = AP_NOISE(ap)
+	ier = AP_OK
+
+	# Initialize.
+	AP_CXCUR(ctr) = wx
+	AP_CYCUR(ctr) = wy
+	AP_OXINIT(ctr) = INDEFR
+	AP_OYINIT(ctr) = INDEFR
+	AP_XCENTER(ctr) = INDEFR
+	AP_YCENTER(ctr) = INDEFR
+	AP_OXCENTER(ctr) = INDEFR
+	AP_OYCENTER(ctr) = INDEFR
+	AP_XSHIFT(ctr) = 0.0
+	AP_YSHIFT(ctr) = 0.0
+	AP_OXSHIFT(ctr) = 0.0
+	AP_OYSHIFT(ctr) = 0.0
+	AP_XERR(ctr) = INDEFR
+	AP_YERR(ctr) = INDEFR
+
+	# Return input coordinates if centering is disabled.
+	if (IS_INDEFR(wx) || IS_INDEFR(wy)) {
+	    return (AP_CTR_NOAREA)
+	} else if (AP_CENTERFUNCTION(ctr) == AP_NONE) {
+	    AP_XCENTER(ctr) = wx
+	    AP_YCENTER(ctr) = wy
+	    switch (AP_WCSOUT(ap)) {
+	    case WCS_WORLD, WCS_PHYSICAL:
+		call ap_ltoo (ap, wx, wy, AP_OXINIT(ctr), AP_OYINIT(ctr), 1)
+		call ap_ltoo (ap, wx, wy, AP_OXCENTER(ctr), AP_OYCENTER(ctr), 1)
+	    case WCS_TV:
+		call ap_ltov (im, wx, wy, AP_OXINIT(ctr), AP_OYINIT(ctr), 1)
+		call ap_ltov (im, wx, wy, AP_OXCENTER(ctr), AP_OYCENTER(ctr), 1)
+	    default:
+		AP_OXINIT(ctr) = wx
+		AP_OYINIT(ctr) = wy
+		AP_OXCENTER(ctr) = wx
+		AP_OYCENTER(ctr) = wy
+	    }
+	    return (AP_OK)
+	}
+
+	# Intialize.
+	owx = wx
+	owy = wy 
+	niter = 0
+	if (IS_INDEFR(AP_SKYSIGMA(nse)) || IS_INDEFR(AP_CTHRESHOLD(ctr)) ||
+	    AP_CTHRESHOLD(ctr) <= 0.0)
+	    cthreshold = 0.0
+	else
+	    cthreshold = AP_CTHRESHOLD(ctr) * AP_SKYSIGMA(nse) 
+
+	repeat {
+
+	    # Set initial cursor position.
+	    AP_CXCUR(ctr) = owx
+	    AP_CYCUR(ctr) = owy
+
+	    # Get centering pixels.
+	    ier = apctrbuf (ap, im, owx, owy)
+	    if (ier == AP_CTR_NOAREA) {
+	    	AP_XCENTER(ctr) = wx
+	    	AP_YCENTER(ctr) = wy
+	    	AP_XSHIFT(ctr) = 0.0
+	    	AP_YSHIFT(ctr) = 0.0
+	    	AP_XERR(ctr) = INDEFR
+	    	AP_YERR(ctr) = INDEFR
+	        switch (AP_WCSOUT(ap)) {
+	        case WCS_WORLD, WCS_PHYSICAL:
+		    call ap_ltoo (ap, wx, wy, AP_OXINIT(ctr), AP_OYINIT(ctr), 1)
+		    call ap_ltoo (ap, wx, wy, AP_OXCENTER(ctr),
+		        AP_OYCENTER(ctr), 1)
+	        case WCS_TV:
+		    call ap_ltov (im, wx, wy, AP_OXINIT(ctr), AP_OYINIT(ctr), 1)
+		    call ap_ltov (im, wx, wy, AP_OXCENTER(ctr),
+		        AP_OYCENTER(ctr), 1)
+	        default:
+		    AP_OXINIT(ctr) = wx
+		    AP_OYINIT(ctr) = wy
+		    AP_OXCENTER(ctr) = wx
+		    AP_OYCENTER(ctr) = wy
+	        }
+	    	AP_OXSHIFT(ctr) = 0.0
+	    	AP_OYSHIFT(ctr) = 0.0
+	    	return (ier)
+	    }
+
+	    # Clean the subraster.
+	    if (AP_CLEAN(ctr) == YES)
+	        call apclean (ap, Memr[AP_CTRPIX(ctr)], AP_CNX(ctr),
+		    AP_CNY(ctr), AP_CXC(ctr), AP_CYC(ctr))
+
+	    # Compute the datalimits.
+	    if (cthreshold <= 0.0)
+	        call alimr (Memr[AP_CTRPIX(ctr)], AP_CNX(ctr) * AP_CNY(ctr),
+		    datamin, datamax)
+	    else {
+		datamin = MAX_REAL
+		datamax = -MAX_REAL
+		do i = 1, AP_CNY(ctr) {
+		    med = amedr (Memr[AP_CTRPIX(ctr)+(i-1)*AP_CNX(ctr)],
+		        AP_CNX(ctr))
+		    if (med < datamin)
+			datamin = med
+		    if (med > datamax)
+			datamax = med
+		}
+	    }
+
+	    # Apply threshold and check for positive or negative features.
+	    if (AP_POSITIVE(ap) == YES) {
+		call apsetr (ap, CDATALIMIT, datamin)
+		call asubkr (Memr[AP_CTRPIX(ctr)], datamin +
+		    cthreshold, Memr[AP_CTRPIX(ctr)], AP_CNX(ctr) *
+		    AP_CNY(ctr))
+		call amaxkr (Memr[AP_CTRPIX(ctr)], 0.0, Memr[AP_CTRPIX(ctr)],
+		    AP_CNX(ctr) * AP_CNY(ctr))
+	    } else {
+		call apsetr (ap, CDATALIMIT, datamax)
+		call anegr (Memr[AP_CTRPIX(ctr)], Memr[AP_CTRPIX(ctr)],
+		    AP_CNX(ctr) * AP_CNY(ctr))
+		call aaddkr (Memr[AP_CTRPIX(ctr)], datamax -
+		    cthreshold, Memr[AP_CTRPIX(ctr)], AP_CNX(ctr) *
+		    AP_CNY(ctr))
+		call amaxkr (Memr[AP_CTRPIX(ctr)], 0.0,
+		    Memr[AP_CTRPIX(ctr)], AP_CNX(ctr) * AP_CNY(ctr))
+	    }
+
+	    # Test signal to noise ratio.
+	    if (ap_csnratio (ap, Memr[AP_CTRPIX(ctr)], AP_CNX(ctr),
+	        AP_CNY(ctr), 0.0) < AP_MINSNRATIO(ctr))
+	    	lowsnr = YES
+	    else
+		lowsnr = NO
+
+	    # Compute the x and y centers.
+	    switch (AP_CENTERFUNCTION(ctr)) {
+
+	    case AP_CENTROID1D:
+		if (AP_CTHRESHOLD(ctr) <= 0.0) {
+	            fier = ap_mctr1d (Memr[AP_CTRPIX(ctr)], AP_CNX(ctr),
+		        AP_CNY(ctr), AP_EPADU(nse), AP_XCENTER(ctr),
+			AP_YCENTER(ctr), AP_XERR(ctr), AP_YERR(ctr))
+		    if (IS_INDEFR (AP_XERR(ctr)))
+			AP_XCENTER(ctr) = AP_CXC(ctr)
+		    if (IS_INDEFR (AP_YERR(ctr)))
+			AP_YCENTER(ctr) = AP_CYC(ctr)
+		} else {
+	            fier = ap_ctr1d (Memr[AP_CTRPIX(ctr)], AP_CNX(ctr),
+		        AP_CNY(ctr), AP_EPADU(nse), AP_XCENTER(ctr),
+			AP_YCENTER(ctr), AP_XERR(ctr), AP_YERR(ctr))
+		    if (IS_INDEFR (AP_XERR(ctr)))
+			AP_XCENTER(ctr) = AP_CXC(ctr)
+		    if (IS_INDEFR (AP_YERR(ctr)))
+			AP_YCENTER(ctr) = AP_CYC(ctr)
+		}
+
+	    case AP_GAUSS1D:
+
+	        fier = ap_gctr1d (Memr[AP_CTRPIX(ctr)], AP_CNX(ctr),
+		    AP_CNY(ctr), CONVERT * AP_FWHMPSF(ap) * AP_SCALE(ap),
+		    AP_CMAXITER(ctr), AP_XCENTER(ctr), AP_YCENTER(ctr),
+		    AP_XERR(ctr), AP_YERR(ctr))
+
+	    case AP_OFILT1D:
+
+	        fier = ap_lgctr1d (Memr[AP_CTRPIX(ctr)], AP_CNX(ctr),
+		    AP_CNY(ctr), AP_CXC(ctr), AP_CYC(ctr), CONVERT *
+		    AP_FWHMPSF(ap) * AP_SCALE(ap), AP_CMAXITER(ctr),
+		    AP_EPADU(nse), AP_SKYSIGMA(nse), AP_XCENTER(ctr),
+		    AP_YCENTER(ctr), AP_XERR(ctr), AP_YERR(ctr))
+
+	    default:
+		# do nothing gracefully
+	    }
+
+	    # Confine the next x center to the data box.
+	    AP_XCENTER(ctr) = max (0.5, min (AP_CNX(ctr) + 0.5,
+	        AP_XCENTER(ctr)))
+	    xshift = AP_XCENTER(ctr) - AP_CXC(ctr)
+	    AP_XCENTER(ctr) = xshift + owx
+	    AP_XSHIFT(ctr) = AP_XCENTER(ctr) - wx
+
+	    # Confine the next y center to the data box.
+	    AP_YCENTER(ctr) = max (0.5, min (AP_CNY(ctr) + 0.5,
+	        AP_YCENTER(ctr)))
+	    yshift = AP_YCENTER(ctr) - AP_CYC(ctr)
+	    AP_YCENTER(ctr) = yshift + owy
+	    AP_YSHIFT(ctr) = AP_YCENTER(ctr) - wy
+
+	    switch (AP_WCSOUT(ap)) {
+	    case WCS_PHYSICAL, WCS_WORLD:
+		call ap_ltoo (ap, AP_XCENTER(ctr), AP_YCENTER(ctr),
+		    AP_OXCENTER(ctr), AP_OYCENTER(ctr), 1)
+		call ap_ltoo (ap, AP_XCENTER(ctr) - AP_XSHIFT(ctr),
+		    AP_YCENTER(ctr) - AP_YSHIFT(ctr), AP_OXINIT(ctr),
+		    AP_OYINIT(ctr), 1)
+		AP_OXSHIFT(ctr) = AP_OXCENTER(ctr) - AP_OXINIT(ctr)
+		AP_OYSHIFT(ctr) = AP_OYCENTER(ctr) - AP_OYINIT(ctr)
+	    case WCS_TV:
+		call ap_ltov (im, AP_XCENTER(ctr), AP_YCENTER(ctr),
+		    AP_OXCENTER(ctr), AP_OYCENTER(ctr), 1)
+		call ap_ltov (im, AP_XCENTER(ctr) - AP_XSHIFT(ctr),
+		    AP_YCENTER(ctr) - AP_YSHIFT(ctr), AP_OXINIT(ctr),
+		    AP_OYINIT(ctr), 1)
+		AP_OXSHIFT(ctr) = AP_OXCENTER(ctr) - AP_OXINIT(ctr)
+		AP_OYSHIFT(ctr) = AP_OYCENTER(ctr) - AP_OYINIT(ctr)
+	    default:
+		AP_OXINIT(ctr) = AP_XCENTER(ctr) - AP_XSHIFT(ctr)
+		AP_OYINIT(ctr) = AP_YCENTER(ctr) - AP_YSHIFT(ctr)
+		AP_OXCENTER(ctr) = AP_XCENTER(ctr)
+		AP_OYCENTER(ctr) = AP_YCENTER(ctr)
+		AP_OXSHIFT(ctr) = AP_XSHIFT(ctr)
+		AP_OYSHIFT(ctr) = AP_YSHIFT(ctr)
+		
+	    }
+
+	    # Setup for next iteration.
+	    niter = niter + 1
+	    owx = AP_XCENTER(ctr)
+	    owy = AP_YCENTER(ctr)
+
+	} until ((fier != AP_OK && fier != AP_CTR_NOCONVERGE) ||
+	    (niter >= AP_CMAXITER(ctr)) || (abs (xshift) < 1.0 &&
+	    abs (yshift) < 1.0))
+
+	# Return appropriate error code.
+	if (fier != AP_OK) {
+	    AP_XCENTER(ctr) = wx
+	    AP_YCENTER(ctr) = wy
+	    AP_XSHIFT(ctr) = 0.0
+	    AP_YSHIFT(ctr) = 0.0
+	    AP_XERR(ctr) = INDEFR
+	    AP_YERR(ctr) = INDEFR
+	    switch (AP_WCSOUT(ap)) {
+	    case WCS_WORLD, WCS_PHYSICAL:
+		call ap_ltoo (ap, wx, wy, AP_OXINIT(ctr), AP_OYINIT(ctr), 1)
+		call ap_ltoo (ap, wx, wy, AP_OXCENTER(ctr), AP_OYCENTER(ctr), 1)
+	    case WCS_TV:
+		call ap_ltov (im, wx, wy, AP_OXINIT(ctr), AP_OYINIT(ctr), 1)
+		call ap_ltov (im, wx, wy, AP_OXCENTER(ctr), AP_OYCENTER(ctr), 1)
+	    default:
+		AP_OXINIT(ctr) = wx
+		AP_OYINIT(ctr) = wy
+		AP_OXCENTER(ctr) = wx
+		AP_OYCENTER(ctr) = wy
+	    }
+	    AP_OXSHIFT(ctr) = 0.0
+	    AP_OYSHIFT(ctr) = 0.0
+	    return (fier)
+	} else if (ier == AP_CTR_BADDATA) {
+	    return (AP_CTR_BADDATA)
+	} else if (lowsnr == YES) {
+	    return (AP_CTR_LOWSNRATIO)
+	} else if (abs (AP_XSHIFT(ctr)) > (AP_MAXSHIFT(ctr) * AP_SCALE(ap))) {
+	    return (AP_CTR_BADSHIFT)
+	} else if (abs (AP_YSHIFT(ctr)) > (AP_MAXSHIFT(ctr) * AP_SCALE(ap))) {
+	    return (AP_CTR_BADSHIFT)
+	} else if (ier == AP_CTR_OUTOFBOUNDS) {
+	    return (AP_CTR_OUTOFBOUNDS)
+	} else {
+	    return (AP_OK)
+	}
+end
diff --git a/noao/digiphot/apphot/center/apgcpars.x b/noao/digiphot/apphot/center/apgcpars.x
new file mode 100644
index 00000000..e40db12e
--- /dev/null
+++ b/noao/digiphot/apphot/center/apgcpars.x
@@ -0,0 +1,27 @@
+include "../lib/display.h"
+include "../lib/noise.h"
+include "../lib/center.h"
+
+# AP_GCPARS -- Procedure to read in the centering parameters from the
+# appropriate parameters files.
+
+procedure ap_gcpars (ap)
+
+pointer	ap			# pointer to the apphot structure
+
+bool	clgetb()
+int	btoi()
+
+begin
+	# Initialize the structure.
+	call apcinit (ap, AP_CENTROID1D, 2.5, 2.0, AP_NPOISSON)
+
+	# Get the data dependent parameters.
+	call ap_gdapars (ap)
+
+	# Get the centering algorithm parameters.
+	call ap_gcepars (ap)
+
+	# Make radial plots on stdgraph.
+	call apseti (ap, RADPLOTS, btoi (clgetb ("radplots")))
+end
diff --git a/noao/digiphot/apphot/center/apgctr1d.x b/noao/digiphot/apphot/center/apgctr1d.x
new file mode 100644
index 00000000..e0ac8077
--- /dev/null
+++ b/noao/digiphot/apphot/center/apgctr1d.x
@@ -0,0 +1,105 @@
+include <math/nlfit.h>
+include "../lib/center.h"
+
+define	NPARS	4		# the total number of parameters
+define	NAPARS	3		# the total number of active parameters
+define	TOL	0.001		# the tolerance for convergence
+
+# AP_GCTR1D -- Procedure to compute the x and y centers from the 1D marginal
+# distributions using 1D Gaussian fits. Three parameters are fit for each
+# marginal, the amplitude, the center of the Gaussian function itself
+# and a constant background value. The sigma is set by the user and is
+# assumed to be fixed.
+
+int procedure ap_gctr1d (ctrpix, nx, ny, sigma, maxiter, xc, yc, xerr, yerr)
+
+real	ctrpix[nx, ny]		# data subarray
+int	nx, ny			# dimensions of data subarray
+real	sigma			# sigma of PSF
+int	maxiter			# maximum number of iterations
+real	xc, yc			# computed centers
+real	xerr, yerr		# estimate of centering error
+
+extern	cgauss1d, cdgauss1d
+int	i, minel, maxel, xier, yier, npar, npts
+pointer	sp, x, xm, ym, w, fit, list, nl
+real	chisqr, variance, p[NPARS], dp[NPARS]
+int	locpr()
+
+begin
+	# Check the number of points.
+	if (nx < NAPARS || ny < NAPARS)
+	    return (AP_CTR_NTOO_SMALL)
+	npts = max (nx, ny)
+
+	call smark (sp)
+	call salloc (list, NAPARS, TY_INT)
+	call salloc (xm, nx, TY_REAL)
+	call salloc (ym, ny, TY_REAL)
+	call salloc (x, npts, TY_REAL)
+	call salloc (w, npts, TY_REAL)
+	call salloc (fit, npts, TY_REAL)
+	do i = 1, npts
+	    Memr[x+i-1] = i
+
+	# Compute the marginal distributions.
+	call ap_mkmarg (ctrpix, Memr[xm], Memr[ym], nx, ny)
+	call adivkr (Memr[xm], real (nx), Memr[xm], nx)
+	call adivkr (Memr[ym], real (ny), Memr[ym], ny)
+
+	# Specify which parameters are to be fit.
+	Memi[list] = 1
+	Memi[list+1] = 2
+	Memi[list+2] = 4
+
+	# Initialize the x fit parameters.
+	call ap_alimr (Memr[xm], nx, p[4], p[1], minel, maxel)
+	p[1] = p[1] - p[4]
+	p[2] = maxel
+	p[3] = sigma ** 2
+
+	# Compute the x center and error.
+	call nlinitr (nl, locpr (cgauss1d), locpr (cdgauss1d), p, dp, NPARS,
+	    Memi[list], NAPARS, TOL, maxiter)
+	call nlfitr (nl, Memr[x], Memr[xm], Memr[w], nx, 1, WTS_UNIFORM, xier)
+	call nlvectorr (nl, Memr[x], Memr[fit], nx, 1)
+	call nlpgetr (nl, p, npar)
+	call nlerrorsr (nl, Memr[xm], Memr[fit], Memr[w], nx, variance,
+	    chisqr, dp)
+	call nlfreer (nl)
+	xc = p[2]
+	xerr = dp[2] / sqrt (real (nx))
+	if (xerr > real (nx))
+	    xerr = INDEFR
+
+	# Initialize the y fit parameters.
+	call ap_alimr (Memr[ym], ny, p[4], p[1], minel, maxel)
+	p[1] = p[1] - p[4]
+	p[2] = maxel
+	p[3] = sigma ** 2
+
+	# Fit the y marginal.
+	call nlinitr (nl, locpr (cgauss1d), locpr (cdgauss1d), p, dp, NPARS,
+	    Memi[list], NAPARS, TOL, maxiter)
+	call nlfitr (nl, Memr[x], Memr[ym], Memr[w], ny, 1, WTS_UNIFORM, yier)
+	call nlvectorr (nl, Memr[x], Memr[fit], ny, 1)
+	call nlpgetr (nl, p, npar)
+	call nlerrorsr (nl, Memr[ym], Memr[fit], Memr[w], ny, variance,
+	    chisqr, dp)
+	call nlfreer (nl)
+	yc = p[2]
+	yerr = dp[2] / sqrt (real (ny))
+	if (yerr > real (ny))
+	    yerr = INDEFR
+
+	# Return the appropriate error code.
+	call sfree (sp)
+	if (xier == NO_DEG_FREEDOM || yier == NO_DEG_FREEDOM)
+	    return (AP_CTR_NTOO_SMALL)
+	else if (xier == SINGULAR || yier == SINGULAR)
+	    return (AP_CTR_SINGULAR)
+	else if (xier == NOT_DONE || yier == NOT_DONE)
+	    return (AP_CTR_NOCONVERGE)
+	else
+	    return (AP_OK)
+end
diff --git a/noao/digiphot/apphot/center/apictr.x b/noao/digiphot/apphot/center/apictr.x
new file mode 100644
index 00000000..3a3bd31e
--- /dev/null
+++ b/noao/digiphot/apphot/center/apictr.x
@@ -0,0 +1,48 @@
+# AP_ICTR -- Given a subraster of pixels and an initial center compute
+# the centroid of the pixels.
+
+procedure ap_ictr (im, wx, wy, radius, emission, xcenter, ycenter)
+
+pointer	im		# pointer to the iraf image
+real	wx		# initial x coordinate
+real	wy		# initial y coordinate
+int	radius		# half width of centering box
+int	emission	# emission feature
+real	xcenter		# fitted x coordinate
+real	ycenter		# fitted y coordinate
+
+int	nx, ny, ier
+pointer	cbuf
+real	xc, yc, datamin, datamax, junk
+int	ap_mctr1d()
+pointer	ap_ctrpix()
+
+begin
+	# Get the pixels.
+	cbuf = ap_ctrpix (im, wx, wy, radius, xc, yc, nx, ny)
+	if (cbuf == NULL) {
+	    xcenter = wx
+	    ycenter = wy
+	    return
+	}
+
+	# Fit the center of the subraster.
+	call alimr (Memr[cbuf], nx * ny, datamin, datamax)
+	if (emission == YES)
+	    call asubkr (Memr[cbuf], datamin, Memr[cbuf], nx * ny)
+	else {
+	    call anegr (Memr[cbuf], Memr[cbuf], nx * ny)
+	    call aaddkr (Memr[cbuf], datamax, Memr[cbuf], nx * ny)
+	}
+	ier = ap_mctr1d (Memr[cbuf], nx, ny, 1.0, xcenter, ycenter, junk, junk)
+
+	# Compute the center in image coordinates.
+	if (IS_INDEFR(xcenter))
+	    xcenter = wx
+	else
+	    xcenter = xcenter + wx - xc
+	if (IS_INDEFR(ycenter))
+	    ycenter = wy
+	else
+	    ycenter = ycenter + wy - yc
+end
diff --git a/noao/digiphot/apphot/center/aplgctr1d.x b/noao/digiphot/apphot/center/aplgctr1d.x
new file mode 100644
index 00000000..ce6d7c53
--- /dev/null
+++ b/noao/digiphot/apphot/center/aplgctr1d.x
@@ -0,0 +1,87 @@
+include <math.h>
+include "../lib/center.h"
+
+define	TOL	0.001		# tolerance for fitting algorithm
+
+# AP_LGCTR1D -- Procedure to compute the center from the 1D marginal
+# distributions using a simplified version of the optimal filtering
+# technique and addopting a Gaussian model for the fit. The method
+# is streamlined by replacing the Gaussian with a simple triangle
+# following L.Goad.
+
+int procedure ap_lgctr1d (ctrpix, nx, ny, cx, cy, sigma, maxiter, norm,
+	skysigma, xc, yc, xerr, yerr)
+
+real	ctrpix[nx, ny]		# object to be centered
+int	nx, ny			# dimensions of subarray
+real	cx, cy			# center in subraster coordinates
+real	sigma			# sigma of PSF
+int	maxiter			# maximum number of iterations
+real	norm			# the normalization factor
+real	skysigma		# standard deviation of the pixels
+real	xc, yc			# computed centers
+real	xerr, yerr		# first guess at errors
+
+int	nxiter, nyiter
+pointer	sp, xm, ym
+real	ratio, constant
+
+int	aptopt()
+real	asumr()
+
+begin
+	# Allocate working space.
+	call smark (sp)
+	call salloc (xm, nx, TY_REAL)
+	call salloc (ym, ny, TY_REAL)
+
+	# Compute the marginal distributions.
+	call ap_mkmarg (ctrpix, Memr[xm], Memr[ym], nx, ny)
+	xerr = asumr (Memr[xm], nx)
+	yerr = asumr (Memr[ym], ny)
+	call adivkr (Memr[xm], real (nx), Memr[xm], nx)
+	call adivkr (Memr[ym], real (ny), Memr[ym], ny)
+
+	# Compute the x center and error.
+	xc = cx
+	nxiter = aptopt (Memr[xm], nx, xc, sigma, TOL, maxiter, YES)
+	if (xerr <= 0.0)
+	    xerr = INDEFR
+	else {
+	    if (IS_INDEFR(skysigma))
+	        constant = 0.0
+	    else
+	        constant = 4.0 * SQRTOFPI * sigma * skysigma ** 2 
+	    ratio = constant / xerr
+	    xerr = sigma ** 2 / (xerr * norm)
+	    xerr = sqrt (max (xerr, ratio * xerr))
+	    if (xerr > real (nx))
+		xerr = INDEFR
+	}
+
+	# Compute the y center and error.
+	yc = cy
+	nyiter = aptopt (Memr[ym], ny, yc, sigma, TOL, maxiter, YES)
+	if (yerr <= 0.0)
+	    yerr = INDEFR
+	else {
+	    if (IS_INDEFR(skysigma))
+	        constant = 0.0
+	    else
+	        constant = 4.0 * SQRTOFPI * sigma * skysigma ** 2 
+	    ratio = constant / yerr
+	    yerr = sigma ** 2 / (yerr * norm)
+	    yerr = sqrt (max (yerr, ratio * yerr))
+	    if (yerr > real (ny))
+		yerr = INDEFR
+	}
+
+	# Return appropriate error code.
+	call sfree (sp)
+	if (nxiter < 0 || nyiter < 0)
+	    return (AP_CTR_SINGULAR)
+	else if (nxiter > maxiter || nyiter > maxiter)
+	    return (AP_CTR_NOCONVERGE)
+	else
+	    return (AP_OK)
+end
diff --git a/noao/digiphot/apphot/center/apmctr1d.x b/noao/digiphot/apphot/center/apmctr1d.x
new file mode 100644
index 00000000..5dc47f60
--- /dev/null
+++ b/noao/digiphot/apphot/center/apmctr1d.x
@@ -0,0 +1,88 @@
+include <mach.h>
+include "../lib/center.h"
+
+# AP_MCTR1D -- Procedure to compute the center from the 1D marginal
+# distributions.
+
+int procedure ap_mctr1d (ctrpix, nx, ny, norm, xc, yc, xerr, yerr)
+
+real	ctrpix[nx, ny]		# object to be centered
+int	nx, ny			# dimensions of subarray
+real	norm			# the normalization factor
+real	xc, yc			# computed centers
+real	xerr, yerr		# estimate of centering error
+
+pointer	sp, xm, ym
+
+begin
+	call smark (sp)
+	call salloc (xm, nx, TY_REAL)
+	call salloc (ym, ny, TY_REAL)
+
+	# Compute the marginal distributions.
+	call ap_mkmarg (ctrpix, Memr[xm], Memr[ym], nx, ny)
+	call adivkr (Memr[xm], real (ny), Memr[xm], nx)
+	call adivkr (Memr[ym], real (nx), Memr[ym], ny)
+
+	# Get the centers and errors.
+	call ap_cmmarg (Memr[xm], nx, norm, xc, xerr)
+	call ap_cmmarg (Memr[ym], ny, norm, yc, yerr)
+
+	call sfree (sp)
+	return (AP_OK)
+end
+
+
+# AP_CMMARG -- Compute the center and estimate its error given the
+# marginal distribution and the number of points. The center is estimated
+# using a centroid which is an optimal solution for a gaussian star with
+# no background noise. The centering error is computing using a variance
+# estimate based on a gaussian star with no background noise.
+
+procedure ap_cmmarg (a, npts, norm, xc, err)
+
+real	a[npts]		# array
+int	npts		# number of points
+real	norm		# the normalization factor
+real	xc		# center value
+real	err		# error
+
+int	i, npos
+real	sumi, sumix, sumix2, mean, val
+real	asumr()
+
+begin
+	# Initialize.
+	mean = asumr (a, npts) / npts
+	npos = 0
+	sumi = 0.0
+	sumix = 0.0
+	sumix2 = 0.0
+
+	# Accumulate the sums.
+	do i = 1, npts {
+	    val = (a[i] - mean)
+	    if (val > 0.0) {
+		npos = npos + 1
+	        sumi = sumi + val
+	        sumix = sumix + val * i
+	        sumix2 = sumix2 + val * i ** 2
+	    }
+	}
+
+	# Compute the position and the error.
+	if (npos <= 0) {
+	    xc =  (1.0 + npts) / 2.0
+	    err = INDEFR
+	} else {
+	    xc = sumix / sumi
+	    err = (sumix2 / sumi - xc * xc)
+	    if (err <= 0.0) {
+		err = 0.0
+	    } else {
+	        err = sqrt (err / ((sumi + mean * npts) * norm))
+		if (err > real (npts))
+		    err = INDEFR
+	    }
+	}
+end
diff --git a/noao/digiphot/apphot/center/appcenter.x b/noao/digiphot/apphot/center/appcenter.x
new file mode 100644
index 00000000..c2590fa3
--- /dev/null
+++ b/noao/digiphot/apphot/center/appcenter.x
@@ -0,0 +1,84 @@
+include "../lib/apphot.h"
+include "../lib/center.h"
+
+# AP_PCENTER -- Procedure to write the output of centering task to a file.
+
+procedure ap_pcenter (ap, fd, id, lid, ier)
+
+pointer	ap	# pointer to apphot structure
+int	fd	# output file descriptor
+int	id	# id of the star
+int	lid	# list number
+int	ier	# comment string
+
+real	xpos, ypos
+real	apstatr()
+
+begin
+	if (fd == NULL)
+	    return
+
+	# Print description of object.
+	xpos = apstatr (ap, OXINIT)
+	ypos = apstatr (ap, OYINIT)
+	call ap_wid (ap, fd, xpos, ypos, id, lid, '\\')
+
+	# Print the computed centers.
+	call ap_wcres (ap, fd, ier, ' ')
+end
+
+
+# AP_QCENTER -- Procedure to print a quick summary of the center task output
+# on the standard output.
+
+procedure ap_qcenter (ap, ier)
+
+pointer	ap	# pointer to apphot structure
+int	ier	# comment string
+
+real	owx, owy, wx, wy
+pointer	sp, imname
+real	apstatr()
+
+begin
+	call smark (sp)
+	call salloc (imname, SZ_FNAME, TY_CHAR)
+
+	owx = apstatr (ap, OXINIT)
+	owy = apstatr (ap, OYINIT)
+	wx = apstatr (ap, OXCENTER)
+	wy = apstatr (ap, OYCENTER)
+
+	call apstats (ap, IMROOT, Memc[imname], SZ_FNAME)
+	call printf ( "%s  %7.2f %7.2f  %7.2f %7.2f  ")
+	    call pargstr (Memc[imname])
+	    call pargr (owx)
+	    call pargr (owy)
+	    call pargr (wx)
+	    call pargr (wy)
+	call printf ("%6.2f %6.2f  %s\n")
+	    call pargr (apstatr (ap, XERR))
+	    call pargr (apstatr (ap, YERR))
+	    if (ier != AP_OK)
+		call pargstr ("err")
+	    else
+		call pargstr ("ok")
+
+	call sfree (sp)
+end
+
+
+# AP_CPHDR -- Procedure to write the banner for the center task to the
+# output file.
+
+procedure ap_cphdr (ap, fd)
+
+pointer	ap		# apphot descriptor
+int	fd		# output file descriptor
+
+begin
+	if (fd == NULL)
+	    return
+	call ap_idhdr (ap, fd)
+	call ap_chdr (ap, fd)
+end
diff --git a/noao/digiphot/apphot/center/appcpars.x b/noao/digiphot/apphot/center/appcpars.x
new file mode 100644
index 00000000..272fcc07
--- /dev/null
+++ b/noao/digiphot/apphot/center/appcpars.x
@@ -0,0 +1,21 @@
+include "../lib/display.h"
+
+# AP_PCPARS -- Procedure to write out the current centering parameters
+# to the current parameter files.
+
+procedure ap_pcpars (ap)
+
+pointer	ap		# pointer to apphot structure
+bool	itob()
+int	apstati()
+
+begin
+	# Write the data dependent parameters.
+	call ap_dapars (ap)
+
+	# Write the centering parameters.
+	call ap_cepars (ap)
+
+	# Set the plotting command.
+	call clputb ("radplots", itob (apstati (ap, RADPLOTS)))
+end
diff --git a/noao/digiphot/apphot/center/aprefitcen.x b/noao/digiphot/apphot/center/aprefitcen.x
new file mode 100644
index 00000000..c50e2efb
--- /dev/null
+++ b/noao/digiphot/apphot/center/aprefitcen.x
@@ -0,0 +1,193 @@
+include <mach.h>
+include "../lib/apphotdef.h"
+include "../lib/apphot.h"
+include "../lib/noisedef.h"
+include "../lib/centerdef.h"
+include "../lib/center.h"
+
+define	CONVERT	.424660900	# conversion factor from fwhmpsf to sigma
+
+# APREFITCENTER -- Procedure to refit the centers assuming that the appropriate
+# pixel buffer is in memory. See apfitcenter for further information.
+
+int procedure aprefitcenter (ap, im, ier)
+
+pointer	ap		# pointer to the apphot structure
+pointer	im		# the input image descriptor
+int	ier		# previous error code
+
+int	fier
+pointer	cen, nse
+int	ap_ctr1d(), ap_mctr1d(), ap_gctr1d(), ap_lgctr1d()
+
+begin
+	cen = AP_PCENTER(ap)
+	nse = AP_NOISE(ap)
+
+	# Initialize
+        AP_XCENTER(cen) = AP_CXCUR(cen)
+	AP_YCENTER(cen) = AP_CYCUR(cen)
+	AP_XSHIFT(cen) = 0.0
+	AP_YSHIFT(cen) = 0.0
+	AP_OXSHIFT(cen) = 0.0
+	AP_OYSHIFT(cen) = 0.0
+	AP_XERR(cen) = INDEFR
+	AP_YERR(cen) = INDEFR
+
+	# Return if the center is undefined.
+	if (IS_INDEFR(AP_CXCUR(cen)) || IS_INDEFR(AP_CYCUR(cen))) {
+	    AP_OXINIT(cen) = INDEFR
+	    AP_OYINIT(cen) = INDEFR
+            AP_OXCENTER(cen) = INDEFR
+	    AP_OYCENTER(cen) = INDEFR
+	    return (AP_CTR_NOAREA)
+	}
+
+	# Convert the coordinates.
+	switch (AP_WCSOUT(ap)) {
+        case WCS_WORLD, WCS_PHYSICAL:
+            call ap_ltoo (ap, AP_CXCUR(cen), AP_CYCUR(cen), AP_OXINIT(cen),
+		AP_OYINIT(cen), 1)
+            call ap_ltoo (ap, AP_CXCUR(cen), AP_CYCUR(cen), AP_OXCENTER(cen),
+		AP_OYCENTER(cen), 1)
+        case WCS_TV:
+            call ap_ltov (im, AP_CXCUR(cen), AP_CYCUR(cen), AP_OXINIT(cen),
+		AP_OYINIT(cen), 1)
+            call ap_ltov (im, AP_CXCUR(cen), AP_CYCUR(cen), AP_OXCENTER(cen),
+		AP_OYCENTER(cen), 1)
+        default:
+	    AP_OXINIT(cen) = AP_CXCUR(cen)
+	    AP_OYINIT(cen) = AP_CYCUR(cen)
+            AP_OXCENTER(cen) = AP_CXCUR(cen)
+	    AP_OYCENTER(cen) = AP_CYCUR(cen)
+        }
+
+	# Return input coordinates if no center fitting.
+	if (AP_CENTERFUNCTION(cen) == AP_NONE)
+	    return (AP_OK)
+
+	# Choose the centering algorithm.
+	switch (AP_CENTERFUNCTION(cen)) {
+
+	case AP_CENTROID1D:
+
+	    # Compute the x and y centroids.
+	    if (AP_CTHRESHOLD(cen) <= 0.0) {
+	        fier = ap_mctr1d (Memr[AP_CTRPIX(cen)], AP_CNX(cen),
+		    AP_CNY(cen), AP_EPADU(nse), AP_XCENTER(cen),
+		    AP_YCENTER(cen), AP_XERR(cen), AP_YERR(cen))
+		if (IS_INDEFR(AP_XERR(cen)))
+		    AP_XCENTER(cen) = AP_CXC(cen)
+		if (IS_INDEFR(AP_YERR(cen)))
+		    AP_YCENTER(cen) = AP_CYC(cen)
+	    } else {
+	        fier = ap_ctr1d (Memr[AP_CTRPIX(cen)], AP_CNX(cen),
+		    AP_CNY(cen), AP_EPADU(nse), AP_XCENTER(cen),
+		    AP_YCENTER(cen), AP_XERR(cen), AP_YERR(cen))
+		if (IS_INDEFR(AP_XERR(cen)))
+		    AP_XCENTER(cen) = AP_CXC(cen)
+		if (IS_INDEFR(AP_YERR(cen)))
+		    AP_YCENTER(cen) = AP_CYC(cen)
+	    }
+	    AP_XCENTER(cen) = AP_XCENTER(cen) + AP_CXCUR(cen) - AP_CXC(cen)
+	    AP_YCENTER(cen) = AP_YCENTER(cen) + AP_CYCUR(cen) - AP_CYC(cen)
+	    AP_XSHIFT(cen) = AP_XCENTER(cen) - AP_CXCUR(cen)
+	    AP_YSHIFT(cen) = AP_YCENTER(cen) - AP_CYCUR(cen)
+
+	case AP_GAUSS1D:
+
+	    # Compute the 1D Gaussian x and y centers.
+	    fier = ap_gctr1d (Memr[AP_CTRPIX(cen)], AP_CNX(cen), AP_CNY(cen),
+		CONVERT * AP_FWHMPSF(ap) * AP_SCALE(ap), AP_CMAXITER(cen),
+		AP_XCENTER(cen), AP_YCENTER(cen), AP_XERR(cen), AP_YERR(cen))
+	    AP_XCENTER(cen) = AP_XCENTER(cen) + AP_CXCUR(cen) - AP_CXC(cen)
+	    AP_YCENTER(cen) = AP_YCENTER(cen) + AP_CYCUR(cen) - AP_CYC(cen)
+	    AP_XSHIFT(cen) = AP_XCENTER(cen) - AP_CXCUR(cen)
+	    AP_YSHIFT(cen) = AP_YCENTER(cen) - AP_CYCUR(cen)
+
+	case AP_OFILT1D:
+
+	    # Compute the Goad 1D x and y centers.
+	    fier = ap_lgctr1d (Memr[AP_CTRPIX(cen)], AP_CNX(cen), AP_CNY(cen),
+		 AP_CXC(cen), AP_CYC(cen), CONVERT * AP_FWHMPSF(ap) *
+		 AP_SCALE(ap), AP_CMAXITER(cen), AP_EPADU(nse),
+		 AP_SKYSIGMA(nse), AP_XCENTER(cen), AP_YCENTER(cen),
+		 AP_XERR(cen), AP_YERR(cen))
+	    AP_XCENTER(cen) = AP_XCENTER(cen) + AP_CXCUR(cen) - AP_CXC(cen)
+	    AP_YCENTER(cen) = AP_YCENTER(cen) + AP_CYCUR(cen) - AP_CYC(cen)
+	    AP_XSHIFT(cen) = AP_XCENTER(cen) - AP_CXCUR(cen)
+	    AP_YSHIFT(cen) = AP_YCENTER(cen) - AP_CYCUR(cen)
+
+	default:
+
+	    # do nothing gracefully
+        }
+
+	# Convert the coordinates.
+        switch (AP_WCSOUT(ap)) {
+        case WCS_WORLD, WCS_PHYSICAL:
+            call ap_ltoo (ap, AP_XCENTER(cen), AP_YCENTER(cen),
+                AP_OXCENTER(cen), AP_OYCENTER(cen), 1)
+            call ap_ltoo (ap, AP_XCENTER(cen) - AP_XSHIFT(cen),
+                AP_YCENTER(cen) - AP_YSHIFT(cen), AP_OXINIT(cen),
+                AP_OYINIT(cen), 1)
+            AP_OXSHIFT(cen) = AP_OXCENTER(cen) - AP_OXINIT(cen)
+            AP_OYSHIFT(cen) = AP_OYCENTER(cen) - AP_OYINIT(cen)
+        case WCS_TV:
+            call ap_ltov (im, AP_XCENTER(cen), AP_YCENTER(cen),
+                AP_OXCENTER(cen), AP_OYCENTER(cen), 1)
+            call ap_ltov (im, AP_XCENTER(cen) - AP_XSHIFT(cen),
+                AP_YCENTER(cen) - AP_YSHIFT(cen), AP_OXINIT(cen),
+                AP_OYINIT(cen), 1)
+            AP_OXSHIFT(cen) = AP_OXCENTER(cen) - AP_OXINIT(cen)
+            AP_OYSHIFT(cen) = AP_OYCENTER(cen) - AP_OYINIT(cen)
+        default:
+	    AP_OXCENTER(cen) = AP_XCENTER(cen)
+	    AP_OYCENTER(cen) = AP_YCENTER(cen)
+	    AP_OXINIT(cen) = AP_XCENTER(cen) - AP_XSHIFT(cen)
+	    AP_OYINIT(cen) = AP_YCENTER(cen) - AP_YSHIFT(cen)
+	    AP_OXSHIFT(cen) = AP_XSHIFT(cen)
+	    AP_OYSHIFT(cen) = AP_YSHIFT(cen)
+        }
+
+	# Return appropriate error code.
+	if (fier != AP_OK) {
+	    AP_XCENTER(cen) = AP_CXCUR(cen)
+	    AP_YCENTER(cen) = AP_CYCUR(cen)
+	    AP_XSHIFT(cen) = 0.0
+	    AP_YSHIFT(cen) = 0.0
+	    AP_XERR(cen) = INDEFR
+	    AP_YERR(cen) = INDEFR
+	    switch (AP_WCSOUT(ap)) {
+            case WCS_WORLD, WCS_PHYSICAL:
+                call ap_ltoo (ap, AP_CXCUR(cen), AP_CYCUR(cen), AP_OXINIT(cen),
+		    AP_OYINIT(cen), 1)
+                call ap_ltoo (ap, AP_CXCUR(cen), AP_CYCUR(cen),
+		    AP_OXCENTER(cen), AP_OYCENTER(cen), 1)
+            case WCS_TV:
+                call ap_ltov (im, AP_CXCUR(cen), AP_CYCUR(cen), AP_OXINIT(cen),
+		    AP_OYINIT(cen), 1)
+                call ap_ltov (im, AP_CXCUR(cen), AP_CYCUR(cen),
+		    AP_OXCENTER(cen), AP_OYCENTER(cen), 1)
+            default:
+	        AP_OXCENTER(cen) = AP_CXCUR(cen)
+	        AP_OYCENTER(cen) = AP_CYCUR(cen)
+	        AP_OXINIT(cen) = AP_CXCUR(cen)
+	        AP_OYINIT(cen) = AP_CYCUR(cen)
+            }
+            AP_OXSHIFT(cen) = 0.0
+            AP_OYSHIFT(cen) = 0.0
+	    return (fier)
+	} else if (ier == AP_CTR_BADDATA) {
+	    return (AP_CTR_BADDATA)
+	} else if (ier == AP_CTR_LOWSNRATIO) {
+	    return (AP_CTR_LOWSNRATIO)
+	} else if (abs (AP_XSHIFT(cen)) > (AP_MAXSHIFT(cen) * AP_SCALE(ap))) {
+	    return (AP_CTR_BADSHIFT)
+	} else if (abs (AP_YSHIFT(cen)) > (AP_MAXSHIFT(cen) * AP_SCALE(ap))) {
+	    return (AP_CTR_BADSHIFT)
+	} else if (ier == AP_CTR_OUTOFBOUNDS) {
+	    return (AP_CTR_OUTOFBOUNDS)
+	} else
+	    return (AP_OK)
+end
diff --git a/noao/digiphot/apphot/center/center.key b/noao/digiphot/apphot/center/center.key
new file mode 100644
index 00000000..f71cd8d3
--- /dev/null
+++ b/noao/digiphot/apphot/center/center.key
@@ -0,0 +1,78 @@
+	Interactive Keystroke Commands
+
+?	Print help
+:	Colon commands
+v	Verify the critical parameters
+w	Save the current parameters
+d	Plot radial profile of current star
+i	Interactively set parameters using current star
+f	Fit center of current star
+spbar	Fit center of current star, output results
+m	Move to next star in coordinate list
+n	Center next star in coordinate list, output results
+l	Center remaining stars in coordinate list, output results
+e	Print error messages
+r	Rewind the coordinate list
+q	Exit task
+
+
+	Colon Commands
+
+:show	[data/center]	List the parameters
+:m [n]	Move to next [nth] star in coordinate list
+:n [n]	Center next [nth] star in coordinate list, output results
+
+
+	Colon Parameter Editing Commands
+
+# Image and file name parameters
+
+:image		[string]	Image name
+:coords		[string]	Coordinate file name
+:output 	[string]	Output file name
+
+# Data dependent parameters
+
+:scale		[value]		Image scale (units per pixel)
+:fwhmpsf	[value]		Full-width half-maximum of PSF (scale units)
+:emission	[y/n]		Emission feature (y), absorption (n)
+:sigma		[value]		Standard deviation of sky (counts)
+:datamin	[value]		Minimum good data value (counts)
+:datamax	[value]		Maximum good data value (counts)
+
+# Noise parameters
+
+:noise 		[string]	Noise model (constant|poisson)
+:gain		[string]	Gain image header keyword
+:ccdread	[string]	Readout noise image header keyword
+:epadu		[value]		Gain (electrons per adu)
+:readnoise	[value]		Readout noise (electrons)
+
+# Observations parameters
+
+:exposure	[string]	Exposure time image header keyword
+:airmass	[string]	Airmass image header keyword
+:filter		[string]	Filter image header keyword
+:obstime	[string]	Time of observation image header keyword
+:itime		[value]		Exposure time (time units)
+:xairmass	[value]		Airmass value (number)
+:ifilter	[string]	Filter id string
+:otime		[string]	Time of observation (time units)
+
+# Centering parameters 
+
+:calgorithm	[string]	Centering algorithm
+:cbox		[value]		Width of centering box (scale units)
+:cthreshold	[value]		Centering intensity threshold (sigma)
+:cmaxiter	[value]		Maximum number of iterations
+:maxshift	[value]		Maximum center shift (scale units)
+:minsnratio	[value]		Minimum signal to noise for centering
+:clean		[y/n]		Clean subraster before centering
+:rclean		[value]		Cleaning radius (scale units)
+:rclip		[value]		Clipping radius (scale units)
+:kclean		[value]		Clean K-sigma rejection limit (sigma)
+
+# Plotting and marking parameters
+
+:mkcenter	[y/n]		Mark computed centers on the display
+:radplot	[y/n]		Plot radial profile of object
diff --git a/noao/digiphot/apphot/center/icenter.key b/noao/digiphot/apphot/center/icenter.key
new file mode 100644
index 00000000..0e5379f0
--- /dev/null
+++ b/noao/digiphot/apphot/center/icenter.key
@@ -0,0 +1,12 @@
+                    Interactive Center Setup Menu
+
+	v	Mark and verify the critical center parameters (f,s,c)
+
+	f	Mark and verify the full-width half-maximum of the psf
+	s	Mark and verify the standard deviation of the background
+	l	Mark and verify the minimum good data value
+	u	Mark and verify the maximum good data value
+
+	c	Mark and verify the centering box half-width
+	n	Mark and verify the cleaning radius
+	p	Mark and verify the clipping radius
diff --git a/noao/digiphot/apphot/center/mkpkg b/noao/digiphot/apphot/center/mkpkg
new file mode 100644
index 00000000..dde81929
--- /dev/null
+++ b/noao/digiphot/apphot/center/mkpkg
@@ -0,0 +1,54 @@
+# CENTER task
+
+$checkout libpkg.a ".."
+$update libpkg.a
+$checkin libpkg.a ".."
+$exit
+
+libpkg.a:
+	apbcenter.x	../lib/apphot.h		../lib/display.h	\
+			<fset.h>
+	apcconfirm.x	../lib/apphot.h		../lib/noise.h          \
+			../lib/center.h
+	apcencolon.x	../lib/apphot.h		../lib/center.h         \
+			../lib/display.h	../lib/noise.h
+	apcenter.x	../lib/apphot.h		../lib/center.h		\
+			../lib/display.h	<ctype.h>		\
+			<gset.h>                <imhdr.h>
+	apcerrors.x	../lib/center.h
+	apcfree.x	../lib/apphotdef.h	../lib/centerdef.h
+	apictr.x
+	apcinit.x	../lib/apphotdef.h	../lib/center.h		\
+			../lib/centerdef.h
+	apclean.x	../lib/apphot.h		../lib/center.h		\
+			../lib/noise.h
+	appcpars.x	../lib/display.h
+	apcplot.x	../lib/apphot.h		../lib/apphotdef.h	\
+			../lib/center.h		../lib/centerdef.h	\
+			../lib/noise.h		<gset.h>		\
+			<pkg/gtools.h>		<mach.h>
+	apcpshow.x	../lib/center.h		../lib/display.h
+	apcshow.x	../lib/display.h
+	apcsnratio.x	../lib/center.h		../lib/noise.h          \
+			<mach.h>
+	apctr1d.x	../lib/center.h		<mach.h>
+	apctrbuf.x	../lib/apphotdef.h	../lib/center.h         \
+			../lib/centerdef.h	<imhdr.h> <math.h>      \
+			<mach.h>
+	apmctr1d.x	../lib/center.h		<mach.h>
+	apfitcen.x	../lib/apphotdef.h	../lib/center.h         \
+			../lib/centerdef.h	../lib/noisedef.h	\
+			../lib/apphot.h         <mach.h>
+	apgcpars.x	../lib/center.h		../lib/display.h        \
+			../lib/noise.h
+	apgctr1d.x	../lib/center.h		<math/nlfit.h>
+	aplgctr1d.x	../lib/center.h		<math.h>
+	appcenter.x	../lib/apphot.h		../lib/center.h
+	apcradsetup.x	../lib/display.h
+	aprefitcen.x	../lib/apphotdef.h	../lib/center.h         \
+			../lib/centerdef.h	../lib/noisedef.h	\
+			../lib/apphot.h         <mach.h>
+	t_center.x	../lib/apphot.h		../lib/noise.h          \
+			<fset.h>		<gset.h>                \
+			<imhdr.h>
+	;
diff --git a/noao/digiphot/apphot/center/t_center.x b/noao/digiphot/apphot/center/t_center.x
new file mode 100644
index 00000000..456b93d2
--- /dev/null
+++ b/noao/digiphot/apphot/center/t_center.x
@@ -0,0 +1,306 @@
+include <fset.h>
+include <gset.h>
+include <imhdr.h>
+include "../lib/apphot.h"
+include "../lib/noise.h"
+
+# T_CENTER -- Procedure to compute accurate centers for a list of objects
+# given a list of starting centers, a set of starting parameters and a list
+# of images.
+
+procedure t_center ()
+
+pointer	image	               	# pointer to name of the image
+pointer	output	 		# pointer to output file name
+pointer	coords			# pointer to coordinate file
+pointer plotfile		# name of plot metacode file
+pointer	graphics		# pointer to graphics device
+pointer display			# pointer to display device
+int	interactive		# interactive mode
+int	cache			# cache image buffer in memory
+int	verify			# verify parameters
+int	update			# update parameters
+int	verbose			# verbose mode
+
+pointer	sp, str, cname, outfname, ap, im, id, gd, mgd
+int	limlist, lclist, lolist, sid, lid, cl, pfd, out, root, stat, memstat
+int	imlist, clist, olist, wcs, req_size, old_size, buf_size
+
+pointer	immap(), gopen()
+int	imtlen(), imtgetim(), clplen(), clgfil(), btoi(), fnldir(), strncmp()
+int	open(), strlen(), apcenter(), imtopenp(), clpopnu(), clgwrd()
+int	ap_memstat(), sizeof()
+bool	clgetb(), streq()
+errchk	gopen
+
+begin
+	# Allocate working space.
+	call smark (sp)
+	call salloc (image, SZ_FNAME, TY_CHAR)
+	call salloc (output, SZ_FNAME, TY_CHAR)
+	call salloc (coords, SZ_FNAME, TY_CHAR)
+	call salloc (plotfile, SZ_FNAME, TY_CHAR)
+	call salloc (graphics, SZ_FNAME, TY_CHAR)
+	call salloc (display, SZ_FNAME, TY_CHAR)
+	call salloc (outfname, SZ_FNAME, TY_CHAR)
+	call salloc (cname, SZ_FNAME, TY_CHAR)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	# Set standard output to flush on newline.
+	call fseti (STDOUT, F_FLUSHNL, YES)
+
+	# Get input and output file names.
+	imlist = imtopenp ("image")
+	limlist = imtlen (imlist)
+	clist = clpopnu ("coords") 
+	lclist = clplen (clist)
+	olist = clpopnu ("output")
+	lolist = clplen (olist)
+
+	# Check that image and coordinate list lengths match.
+	if (limlist < 1 || (lclist > 1 && lclist != limlist)) {
+	    call imtclose (imlist)
+	    call clpcls (clist)
+	    call clpcls (olist)
+	    call error (0, "Imcompatible image and coordinate list lengths")
+	}
+
+	# Check that image and output list lengths match.
+	if (lolist > 1 && lolist != limlist) {
+	    call imtclose (imlist)
+	    call clpcls (clist)
+	    call clpcls (olist)
+	    call error (0, "Imcompatible image and output list lengths")
+	}
+
+	call clgstr ("icommands.p_filename", Memc[cname], SZ_FNAME)
+	if (Memc[cname] != EOS)
+	    interactive = NO
+	#else if (lclist == 0)
+	    #interactive = YES
+	else
+	    interactive = btoi (clgetb ("interactive"))
+	cache = btoi (clgetb ("cache"))
+	verbose = btoi (clgetb ("verbose"))
+	verify = btoi (clgetb ("verify"))
+	update = btoi (clgetb ("update"))
+
+	# Get the centering parameters.
+	call ap_gcpars (ap)
+
+	# confirm the centering algorithm parameters.
+	if (verify == YES && interactive == NO) {
+	    call ap_cconfirm (ap, NULL, 1)
+	    if (update == YES)
+		call ap_pcpars (ap)
+	}
+
+	# Get the wcs information.
+        wcs = clgwrd ("wcsin", Memc[str], SZ_LINE, WCSINSTR)
+        if (wcs <= 0) {
+	    call eprintf (
+	        "Warning: Setting the input coordinate system to logical\n")
+            wcs = WCS_LOGICAL
+	}
+        call apseti (ap, WCSIN, wcs)
+        wcs = clgwrd ("wcsout", Memc[str], SZ_LINE, WCSOUTSTR)
+        if (wcs <= 0) {
+	    call eprintf (
+	        "Warning: Setting the output coordinate system to logical\n")
+            wcs = WCS_LOGICAL
+	}
+        call apseti (ap, WCSOUT, wcs)
+	
+	# Get the graphics and display devices.
+	call clgstr ("graphics", Memc[graphics], SZ_FNAME)
+	call clgstr ("display", Memc[display], SZ_FNAME)
+
+	# Open the graphics and image deisplay devices.
+	if (interactive == YES) {
+	    if (Memc[graphics] == EOS)
+		gd = NULL
+	    else {
+		iferr {
+		    gd =  gopen (Memc[graphics], APPEND+AW_DEFER, STDGRAPH)
+		} then {
+		    call eprintf (
+			"Warning: Error opening the graphics device.\n")
+		    gd = NULL
+		}
+	    }
+	    if (Memc[display] == EOS)
+		id = NULL
+	    else if (streq (Memc[graphics], Memc[display]))
+		id = gd
+	    else {
+		iferr {
+		    id = gopen (Memc[display], APPEND, STDIMAGE)
+		} then {
+		    call eprintf (
+		"Warning: Graphics overlay not available for display device.\n")
+		    id = NULL
+		}
+	    }
+	} else {
+	    gd = NULL
+	    id = NULL
+	}
+
+	# Open the plot metacode file.
+	call clgstr ("plotfile", Memc[plotfile], SZ_FNAME)
+	if (Memc[plotfile] == EOS)
+	    pfd = NULL
+	else
+	    pfd = open (Memc[plotfile], APPEND, BINARY_FILE)
+	if (pfd != NULL)
+	    mgd = gopen (Memc[graphics], NEW_FILE, pfd)
+	else
+	    mgd = NULL
+
+	# Begin looping over the image list.
+	sid = 1
+	while (imtgetim (imlist, Memc[image], SZ_FNAME) != EOF) {
+
+	    # Open image.
+	    im = immap (Memc[image], READ_ONLY, 0)
+	    call apimkeys (ap, im, Memc[image])
+
+	    # Set the image display viewport. 
+	    if ((id != NULL) && (id != gd))
+		call ap_gswv (id, Memc[image], im, 4)
+
+	    # Cache the input image pixels.
+	    req_size = MEMFUDGE * IM_LEN(im,1) * IM_LEN(im,2) *
+	        sizeof (IM_PIXTYPE(im))
+	    memstat = ap_memstat (cache, req_size, old_size)
+	    if (memstat == YES)
+		call ap_pcache (im, INDEFI, buf_size)
+
+	    # Open the coordinate file; where coords is assumed to be a simple
+	    # text file in which the x and y positions are in columns 1 and 2
+	    # respectively and all remaining fields are ignored.
+
+	    if (lclist <= 0) {
+		cl = NULL
+		call strcpy ("", Memc[outfname], SZ_FNAME)
+	    } else {
+		stat = clgfil (clist, Memc[coords], SZ_FNAME)
+		root = fnldir (Memc[coords], Memc[outfname], SZ_FNAME)
+		if (strncmp ("default", Memc[coords+root], 7) == 0 || root ==
+		    strlen (Memc[coords])) {
+		    call ap_inname (Memc[image], Memc[outfname], "coo",
+		        Memc[outfname], SZ_FNAME)
+		    lclist = limlist
+	            cl = open (Memc[outfname], READ_ONLY, TEXT_FILE)
+		} else if (stat != EOF) {
+		    call strcpy (Memc[coords], Memc[outfname], SZ_FNAME)
+	            cl = open (Memc[outfname], READ_ONLY, TEXT_FILE)
+		} else { 
+		    call apstats (ap, CLNAME, Memc[outfname], SZ_FNAME)
+		    call seek (cl, BOF)
+		}
+	    }
+	    call apsets (ap, CLNAME, Memc[outfname])
+	    call apfroot (Memc[outfname], Memc[str], SZ_LINE)
+	    call apsets (ap, CLROOT, Memc[str])
+
+	    # Open output text file; if output is "default", dir$default  or
+	    # a directory specification then the extension "ctr" is added to
+	    # the root image name and a suitable version number is appended to
+	    # the output name. If the output string is null then no output
+	    # file is created.
+
+	    if (lolist == 0) {
+		out = NULL
+		call strcpy ("", Memc[outfname], SZ_FNAME)
+	    } else {
+	        stat = clgfil (olist, Memc[output], SZ_FNAME)
+		root = fnldir (Memc[output], Memc[outfname], SZ_FNAME)
+		if (strncmp ("default", Memc[output+root], 7) == 0 || root == 
+		    strlen (Memc[output])) {
+	            call apoutname (Memc[image], Memc[outfname], "ctr",
+		        Memc[outfname], SZ_FNAME)
+	            out = open (Memc[outfname], NEW_FILE, TEXT_FILE)
+		    lolist = limlist
+		} else if (stat != EOF) {
+		    call strcpy (Memc[output], Memc[outfname], SZ_FNAME)
+	            out = open (Memc[outfname], NEW_FILE, TEXT_FILE)
+		} else
+		    call apstats (ap, OUTNAME, Memc[outfname], SZ_FNAME)
+	    }
+	    call apsets (ap, OUTNAME, Memc[outfname])
+
+	    # Fit the centers.
+	    if (interactive == NO) {
+	        if (Memc[cname] != EOS)
+		    stat = apcenter (ap, im, cl, NULL, mgd, NULL, out, sid,
+		        NO, cache)
+	        else if (cl != NULL) {
+		    lid = 1
+	            call apbcenter (ap, im, cl, out, sid, lid, mgd, id,
+		        verbose)
+		    stat = NO
+		} else
+		    stat = NO
+	    } else
+	        stat = apcenter (ap, im, cl, gd, mgd, id, out, sid, YES, cache)
+
+	    # Close up image, coordinates, and output file.
+	    call imunmap (im)
+	    if (cl != NULL) {
+		if (lclist > 1)
+	            call close (cl)
+	    }
+	    if (out != NULL && lolist != 1) {
+		call close (out)
+		if (sid <= 1) {
+		    call apstats (ap, OUTNAME, Memc[outfname], SZ_FNAME)
+		    call delete (Memc[outfname])
+		}
+		sid = 1
+	    }
+
+	    # Uncache memory.
+	    call fixmem (old_size)
+
+	    if (stat == YES)
+		break
+
+	}
+
+	# If only one coordinate file for a list of images close it.
+	if (cl != NULL && lclist == 1)
+	    call close (cl)
+
+	# If only one output file defined for a list of images close it.
+	if (out != NULL && lolist == 1) {
+	    call close (out)
+	    if (sid <= 1) {
+		call apstats (ap, OUTNAME, Memc[outfname], SZ_FNAME)
+		call delete (Memc[outfname])
+	    }
+	}
+
+	# Close up graphics and image display streams and the plot files.
+	if (id == gd && id != NULL)
+	    call gclose (gd)
+	else {
+	    if (gd != NULL)
+	        call gclose (gd)
+	    if (id != NULL)
+		call gclose (id)
+	}
+	if (mgd != NULL)
+	    call gclose (mgd)
+	if (pfd != NULL)
+	    call close (pfd)
+
+	# Free centering structure.
+	call apcfree (ap)
+
+	# Close up the file lists.
+	call imtclose (imlist)
+	call clpcls (clist)
+	call clpcls (olist)
+	call sfree (sp)
+end