Initial commit

17 files changed, 2213 insertions, 0 deletions

diff --git a/noao/digiphot/apphot/aputil/apbsmooth.x b/noao/digiphot/apphot/aputil/apbsmooth.x
new file mode 100644
index 00000000..28c53c00
--- /dev/null
+++ b/noao/digiphot/apphot/aputil/apbsmooth.x
@@ -0,0 +1,33 @@
+# AP_BSMOOTH - Box car smooth a histogram 1, 2 or 3 times.
+
+procedure ap_bsmooth (hgm, shgm, nbins, nker, iter)
+
+real	hgm[ARB]		# the original histogram
+real	shgm[ARB]		# the smoothed histogram
+int	nbins			# length of the histogram
+int	nker			# half width of box kernel
+int	iter			# number of iterations
+
+pointer	sp, work1, work2
+
+begin
+	# Smooth the histogram
+	switch (iter) {
+	case 1:
+	    call ap_sboxr (hgm, shgm, nbins, nker)
+	case 2:
+	    call smark (sp)
+	    call salloc (work1, nbins, TY_REAL)
+	    call ap_sboxr (hgm, Memr[work1], nbins, nker)
+	    call ap_sboxr (Memr[work1], shgm, nbins, nker)
+	    call sfree (sp)
+	default:
+	    call smark (sp)
+	    call salloc (work1, nbins, TY_REAL)
+	    call salloc (work2, nbins, TY_REAL)
+	    call ap_sboxr (hgm, Memr[work1], nbins, nker)
+	    call ap_sboxr (Memr[work1], Memr[work2], nbins, nker)
+	    call ap_sboxr (Memr[work2], shgm, nbins, nker)
+	    call sfree (sp)
+	}
+end
diff --git a/noao/digiphot/apphot/aputil/apclip.x b/noao/digiphot/apphot/aputil/apclip.x
new file mode 100644
index 00000000..62a52c30
--- /dev/null
+++ b/noao/digiphot/apphot/aputil/apclip.x
@@ -0,0 +1,23 @@
+# AP_CLIP -- Clip the ends of a sorted pixel distribution by a certain
+# percent.
+
+int procedure ap_clip (skypix, index, npix, loclip, hiclip, loindex, hiindex)
+
+real	skypix[ARB]		# the unsorted array of sky pixels
+int	index[ARB]		# the array of sorted indices
+int	npix			# the number of sky pixels
+real	loclip, hiclip		# the clipping factors in percent
+int	loindex, hiindex	# the clipping indices
+
+begin
+	# Sort the pixels.
+	call apqsort (skypix, index, index, npix)
+
+	# Determine the clipping factors.
+	loindex = nint (0.01 * loclip * npix) + 1
+	hiindex = npix - nint (0.01 * hiclip * npix)
+	if ((hiindex - loindex + 1) <= 0)
+	    return (npix)
+	else
+	    return (loindex - 1 + npix - hiindex)
+end
diff --git a/noao/digiphot/apphot/aputil/apdate.x b/noao/digiphot/apphot/aputil/apdate.x
new file mode 100644
index 00000000..21faea88
--- /dev/null
+++ b/noao/digiphot/apphot/aputil/apdate.x
@@ -0,0 +1,29 @@
+include <time.h>
+
+# APDATE -- Fetch the date and time strings required for the apphot output
+# files.
+
+procedure apdate (date, time, maxch)
+
+char	date[ARB]	# the date string
+char	time[ARB]	# the time string
+int	maxch		# the maximum number of character in the string
+
+int	tm[LEN_TMSTRUCT]
+long	ctime
+long	clktime()
+#long	lsttogmt()
+
+begin
+	ctime = clktime (long(0))
+	#ctime = lsttogmt (ctime)
+	call brktime (ctime, tm)
+	call sprintf (date, maxch, "%04d-%02d-%02d")
+	    call pargi (TM_YEAR(tm))
+	    call pargi (TM_MONTH(tm))
+	    call pargi (TM_MDAY(tm))
+	call sprintf (time, maxch, "%02d:%02d:%02d")
+	    call pargi (TM_HOUR(tm))
+	    call pargi (TM_MIN(tm))
+	    call pargi (TM_SEC(tm))
+end
diff --git a/noao/digiphot/apphot/aputil/apdiverr.x b/noao/digiphot/apphot/aputil/apdiverr.x
new file mode 100644
index 00000000..d3c36444
--- /dev/null
+++ b/noao/digiphot/apphot/aputil/apdiverr.x
@@ -0,0 +1,7 @@
+# AP_DIVERR -- Error condition procedure which returns a value of unity.
+
+real procedure ap_diverr ()
+
+begin
+	return (1.0)
+end
diff --git a/noao/digiphot/apphot/aputil/apfnames.x b/noao/digiphot/apphot/aputil/apfnames.x
new file mode 100644
index 00000000..95e56167
--- /dev/null
+++ b/noao/digiphot/apphot/aputil/apfnames.x
@@ -0,0 +1,280 @@
+
+# APINNAME -- Construct an apphot input file name.
+# If input is null or a directory, a name is constructed from the root
+# of the image name and the extension. The disk is searched to avoid
+# name collisions.
+
+procedure apinname (image, input, ext, name, maxch)
+
+char	image[ARB]		# image name
+char	input[ARB]		# input directory or name
+char	ext[ARB]		# extension
+char	name[ARB]		# input name
+int	maxch			# maximum size of name
+
+int	ndir, nimdir, clindex, clsize
+pointer	sp, root, str
+int	fnldir(), strlen()
+
+begin
+	call smark (sp)
+	call salloc (root, SZ_FNAME, TY_CHAR)
+	call salloc (str, SZ_FNAME, TY_CHAR)
+
+	ndir = fnldir (input, name, maxch)
+	if (strlen (input) == ndir) {
+            call imparse (image, Memc[root], SZ_FNAME, Memc[str], SZ_FNAME,
+                Memc[str], SZ_FNAME, clindex, clsize)
+            nimdir = fnldir (Memc[root], Memc[str], SZ_FNAME)
+            if (clindex >= 0) {
+                call sprintf (name[ndir+1], maxch, "%s%d.%s.*")
+                    call pargstr (Memc[root+nimdir])
+                    call pargi (clindex)
+                    call pargstr (ext)
+            } else {
+                call sprintf (name[ndir+1], maxch, "%s.%s.*")
+                    call pargstr (Memc[root+nimdir])
+                    call pargstr (ext)
+            }
+	    call apiversion (name, name, maxch)
+	} else
+	    call strcpy (input, name, maxch)
+
+	call sfree (sp)
+end
+
+
+# APOUTNAME -- Construct an apphot output file name.
+# If output is null or a directory, a name is constructed from the root
+# of the image name and the extension. The disk is searched to avoid
+# name collisions.
+
+procedure apoutname (image, output, ext, name, maxch)
+
+char	image[ARB]		# image name
+char	output[ARB]		# output directory or name
+char	ext[ARB]		# extension
+char	name[ARB]		# output name
+int	maxch			# maximum size of name
+
+int	ndir, nimdir, clindex, clsize
+pointer	sp, root, str
+int	fnldir(), strlen()
+
+begin
+	call smark (sp)
+	call salloc (root, SZ_FNAME, TY_CHAR)
+	call salloc (str, SZ_FNAME, TY_CHAR)
+
+	ndir = fnldir (output, name, maxch)
+	if (strlen (output) == ndir) {
+            call imparse (image, Memc[root], SZ_FNAME, Memc[str], SZ_FNAME,
+                Memc[str], SZ_FNAME, clindex, clsize)
+            nimdir = fnldir (Memc[root], Memc[str], SZ_FNAME)
+            if (clindex >= 0) {
+                call sprintf (name[ndir+1], maxch, "%s%d.%s.*")
+                    call pargstr (Memc[root+nimdir])
+                    call pargi (clindex)
+                    call pargstr (ext)
+            } else {
+                call sprintf (name[ndir+1], maxch, "%s.%s.*")
+                    call pargstr (Memc[root+nimdir])
+                    call pargstr (ext)
+            }
+	    call apoversion (name, name, maxch)
+	} else
+	    call strcpy (output, name, maxch)
+
+	call sfree (sp)
+end
+
+
+# APTMPIMAGE -- Generate a temporary image name either by calling a system
+# routine or by appending the image name to a user specified prefix.
+
+int procedure aptmpimage (image, prefix, tmp, name, maxch)
+
+char	image[ARB]		# image name
+char	prefix[ARB]		# user supplied prefix
+char	tmp[ARB]		# user supplied temporary root
+char	name[ARB]		# output name
+int	maxch			# max number of chars
+
+int	npref, ndir, nimdir, clindex, clsize
+pointer sp, root, str
+int	fnldir(), strlen()
+
+begin
+	npref = strlen (prefix)
+	ndir = fnldir (prefix, name, maxch)
+	if (npref == ndir) {
+	    call mktemp (tmp, name[ndir+1], maxch)
+	    return (NO)
+	} else {
+	    call smark (sp)
+	    call salloc (root, SZ_FNAME, TY_CHAR)
+	    call salloc (str, SZ_FNAME, TY_CHAR)
+	    call strcpy (prefix, name, npref)
+	    call imparse (image, Memc[root], SZ_FNAME, Memc[str], SZ_FNAME,
+	        Memc[str], SZ_FNAME, clindex, clsize)
+		nimdir = fnldir (Memc[root], Memc[str], SZ_FNAME)
+	    if (clindex >= 0) {
+	        call sprintf (name[npref+1], maxch, "%s%d")
+		    call pargstr (Memc[root+nimdir])
+		    call pargi (clindex)
+	    } else {
+	        call sprintf (name[npref+1], maxch, "%s")
+		    call pargstr (Memc[root+nimdir])
+	    }
+	    call sfree (sp)
+	    return (YES)
+	}
+end
+
+
+# APIMROOT -- Fetch the root image name minus the directory specification
+# and the section notation.
+
+procedure apimroot (image, root, maxch)
+
+char	image[ARB]		# image specification
+char	root[ARB]		# output root name
+int	maxch			# maximum number of characters
+
+pointer	sp, imroot, kernel, section, str
+int	clindex, clsize, nchars
+int	fnldir()
+
+begin
+	call smark (sp)
+	call salloc (imroot, SZ_PATHNAME, TY_CHAR)
+	call salloc (kernel, SZ_FNAME, TY_CHAR)
+	call salloc (section, SZ_FNAME, TY_CHAR)
+	call salloc (str, SZ_PATHNAME, TY_CHAR)
+
+        call imparse (image, Memc[imroot], SZ_PATHNAME, Memc[kernel], SZ_FNAME,
+            Memc[section], SZ_FNAME, clindex, clsize)
+	nchars = fnldir (Memc[imroot], Memc[str], SZ_PATHNAME)
+	if (clindex >= 0) {
+	    call sprintf (root, maxch, "%s[%d]%s%s")
+	        call pargstr (Memc[imroot+nchars])
+		call pargi (clindex)
+	        call pargstr (Memc[kernel])
+	        call pargstr (Memc[section])
+	} else {
+	    call sprintf (root, maxch, "%s%s%s")
+	        call pargstr (Memc[imroot+nchars])
+	        call pargstr (Memc[kernel])
+	        call pargstr (Memc[section])
+	}
+
+	call sfree (sp)
+end
+
+
+
+# APFROOT -- Fetch the file name minus the directory specification,
+
+procedure apfroot (filename, root, maxch)
+
+char	filename[ARB]		# input file name 
+char	root[ARB]		# output root file name
+int	maxch			# maximum number of characters
+
+pointer	sp, str
+int	nchars
+int	fnldir()
+
+begin
+	call smark (sp)
+	call salloc (str, SZ_PATHNAME, TY_CHAR)
+
+	nchars = fnldir (filename, Memc[str], SZ_PATHNAME)
+	call strcpy (filename[nchars+1], root, maxch)
+
+	call sfree (sp)
+end
+
+
+# APOVERSION -- Compute the next available version number of a given file
+# name template and output the new file name.
+
+procedure apoversion (template, filename, maxch)
+
+char	template[ARB]			# name template
+char	filename[ARB]			# output name
+int	maxch				# maximum number of characters
+
+char	period
+int	newversion, version, len
+pointer	sp, list, name
+int	fntgfnb() strldx(), ctoi(), fntopnb()
+
+begin
+	# Allocate temporary space
+	call smark (sp)
+	call salloc (name, maxch, TY_CHAR)
+	period = '.'
+	list = fntopnb (template, NO)
+
+	# Loop over the names in the list searchng for the highest version.
+	newversion = 0
+	while (fntgfnb (list, Memc[name], maxch) != EOF) {
+	    len = strldx (period, Memc[name])
+	    len = len + 1
+	    if (ctoi (Memc[name], len, version) <= 0)
+		next
+	    newversion = max (newversion, version)
+	}
+
+	# Make new output file name.
+	len = strldx (period, template)
+	call strcpy (template, filename, len)
+	call sprintf (filename[len+1], maxch, "%d")
+	    call pargi (newversion + 1)
+
+	call fntclsb (list)
+	call sfree (sp)
+end
+
+
+# APIVERSION -- Compute the highest available version number of a given file
+# name template and output the file name.
+
+procedure apiversion (template, filename, maxch)
+
+char	template[ARB]			# name template
+char	filename[ARB]			# output name
+int	maxch				# maximum number of characters
+
+char	period
+int	newversion, version, len
+pointer	sp, list, name
+int	fntgfnb() strldx(), ctoi(), fntopnb()
+
+begin
+	# Allocate temporary space
+	call smark (sp)
+	call salloc (name, maxch, TY_CHAR)
+	period = '.'
+	list = fntopnb (template, NO)
+
+	# Loop over the names in the list searchng for the highest version.
+	newversion = 0
+	while (fntgfnb (list, Memc[name], maxch) != EOF) {
+	    len = strldx (period, Memc[name])
+	    len = len + 1
+	    if (ctoi (Memc[name], len, version) <= 0)
+		next
+	    newversion = max (newversion, version)
+	}
+
+	# Make new output file name.
+	len = strldx (period, template)
+	call strcpy (template, filename, len)
+	call sprintf (filename[len+1], maxch, "%d")
+	    call pargi (newversion)
+
+	call fntclsb (list)
+	call sfree (sp)
+end
diff --git a/noao/digiphot/apphot/aputil/apgscur.x b/noao/digiphot/apphot/aputil/apgscur.x
new file mode 100644
index 00000000..5cc8a6df
--- /dev/null
+++ b/noao/digiphot/apphot/aputil/apgscur.x
@@ -0,0 +1,88 @@
+include <gset.h>
+include <fset.h>
+
+# APGSCUR -- Read the x and y coordinates of an object from a file and move
+# the cursor to those coordinates.
+
+int procedure apgscur (sl, gd, xcur, ycur, prev_num, req_num, num)
+
+int	sl		# coordinate file descriptor
+pointer	gd		# pointer to graphics stream
+real	xcur, ycur	# x cur and y cur
+int	prev_num	# previous number
+int	req_num		# requested number
+int	num		# list number
+
+int	stdin, nskip, ncount
+pointer	sp, fname
+int	fscan(), nscan(), strncmp()
+errchk	greactivate, gdeactivate, gscur
+
+begin
+	if (sl == NULL)
+	    return (EOF)
+
+	call smark (sp)
+	call salloc (fname, SZ_FNAME, TY_CHAR)
+
+	# Find the number of objects to be skipped.
+	call fstats (sl, F_FILENAME, Memc[fname], SZ_FNAME)
+	if (strncmp ("STDIN", Memc[fname], 5) == 0) {
+	    stdin = YES
+	    nskip = 1
+	} else {
+	    stdin = NO
+	    if (req_num <= prev_num) {
+		call seek (sl, BOF)
+		nskip = req_num
+	    } else
+		nskip = req_num - prev_num
+	}
+
+	# Initialize.
+	ncount = 0
+	num = prev_num
+
+	# Read the coordinates and write the cursor.
+	repeat {
+
+	    # Print the prompt if file is STDIN.
+	    if (stdin == YES) {
+		call printf ("Type object x and y coordinates: ")
+		call flush (STDOUT)
+	    }
+
+	    # Fetch the coordinates.
+	    if (fscan (sl) != EOF) {
+		call gargr (xcur)
+		call gargr (ycur)
+		if (nscan () == 2) {
+		    ncount = ncount + 1
+		    num = num + 1
+		}
+	    } else 
+		ncount = EOF
+
+	    # Move the cursor.
+	    if (gd != NULL && (ncount == nskip || ncount == EOF)) {
+		iferr {
+		    call greactivate (gd, 0)
+	            call gscur (gd, xcur, ycur)
+		    call gdeactivate (gd, 0)
+		} then
+		    ;
+	    }
+
+	} until (ncount == EOF || ncount == nskip)
+
+	call sfree (sp)
+
+	if (ncount == EOF) {
+	    return (EOF)
+	} else if (nskip == req_num) {
+	    num = ncount
+	    return (ncount)
+	} else {
+	    return (num)
+	}
+end
diff --git a/noao/digiphot/apphot/aputil/apgtools.x b/noao/digiphot/apphot/aputil/apgtools.x
new file mode 100644
index 00000000..e35cb543
--- /dev/null
+++ b/noao/digiphot/apphot/aputil/apgtools.x
@@ -0,0 +1,141 @@
+include <gset.h>
+include <pkg/gtools.h>
+
+# AP_GTINIT -- Initialize the gtools package for the apphot routines.
+
+pointer procedure ap_gtinit (image, wx, wy)
+
+char	image[ARB]	# the image name
+real	wx, wy		# center of sky subraster
+
+pointer	sp, gt, str
+pointer	gt_init()
+
+begin
+	# Allocate working space.
+    	gt = gt_init ()
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	# Set the plot title.
+	call sprintf (Memc[str], SZ_LINE, "Image: %s:  %.2f %.2f\n")
+	    call pargstr (image)
+	    call pargr (wx)
+	    call pargr (wy)
+	call gt_sets (gt, GTTITLE, Memc[str])
+
+	call sfree (sp)
+	return (gt)
+end
+
+
+# AP_GTFREE -- Free the gtools package.
+
+procedure ap_gtfree (gt)
+
+pointer	gt		# pointer to gtools structure
+
+begin
+	call gt_free (gt)
+end
+
+
+# AP_PLOTRAD -- Plot the radial profile of a list of pixels.
+
+procedure ap_plotrad (gd, gt, r, i, npts, polymark) 
+
+pointer	gd		# pointer to graphics stream
+pointer	gt		# the GTOOLS pointer
+real	r[ARB]		# the radii array
+real	i[ARB]		# the intensity array
+int	npts		# number of points
+char	polymark[ARB]	# polyline type
+
+begin
+	call gt_sets (gt, GTTYPE, "mark")
+	call gt_sets (gt, GTMARK, polymark)
+	call gt_plot (gd, gt, r, i, npts)
+end
+
+
+# AP_PLOTPTS -- Plot the radial profile of a list of pixels excluding points
+# that are outside the plotting window altogether.
+
+procedure ap_plotpts (gd, gt, r, i, npts, xmin, xmax, ymin, ymax, polymark) 
+
+pointer	gd		# pointer to graphics stream
+pointer	gt		# the GTOOLS pointer
+real	r[ARB]		# the radii array
+real	i[ARB]		# the intensity array
+int	npts		# number of points
+real	xmin, xmax	# the x plot limits
+real	ymin, ymax	# the y plot limits
+char	polymark[ARB]	# polyline type
+
+int	j
+
+begin
+	call gt_sets (gt, GTTYPE, "mark")
+	call gt_sets (gt, GTMARK, polymark)
+	do j = 1, npts {
+	    if (r[j] < xmin || r[j] > xmax)
+		next
+	    if (i[j] < ymin || i[j] > ymax)
+		next
+	    call gt_plot (gd, gt, r[j], i[j], 1)
+	    
+	}
+end
+
+
+# AP_RSET -- Set up the parameters for the radial profile plot.
+
+procedure ap_rset (gd, gt, xmin, xmax, ymin, ymax, xscale)
+
+pointer	gd		# pointer to GRAPHICS stream
+pointer	gt		# pointer to GTOOLS structure
+real	xmin, xmax	# min and max of x vector
+real	ymin, ymax	# min and max of y vector
+real	xscale		# image scale
+
+pointer	sp, str, title
+real	vx1, vx2, vy1, vy2, aspect
+real	gstatr()
+
+begin
+	call smark (sp)
+	call salloc (title, 3 * SZ_LINE, TY_CHAR)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	# Reset the aspect ratio.
+	aspect = gstatr (gd, G_ASPECT)
+	call gsetr (gd, G_ASPECT, 0.75)
+
+	# Construct the title.
+	call sysid (Memc[title], 3 * SZ_LINE)
+	call strcat ("\n", Memc[title], 3 * SZ_LINE)
+	call gt_gets (gt, GTTITLE, Memc[str], SZ_LINE) 
+	call strcat (Memc[str], Memc[title], 3 * SZ_LINE)
+	call strcat ("\n\n", Memc[title], 3 * SZ_LINE)
+
+	# Draw three axes.
+	call gseti (gd, G_XDRAWAXES, 2)
+	call gswind (gd, xmin / xscale, xmax / xscale, ymin, ymax)
+	call glabax (gd, Memc[title], "", "Intensity") 
+
+	# Draw the bottom x axis 
+	call gseti (gd, G_YDRAWAXES, 0)
+	call gseti (gd, G_XDRAWAXES, 1)
+	call ggview (gd, vx1, vx2, vy1, vy2)
+	call gsview (gd, vx1, vx2, vy1, vy2)
+	call gswind (gd, xmin, xmax, ymin, ymax)
+	call glabax (gd, "",
+	    "Radial Distance (lower-pixels, upper-scale units)", "")
+
+	# Reset to standard gio parameters.
+	call gseti (gd, G_YDRAWAXES, 3)
+	call gseti (gd, G_XDRAWAXES, 3)
+	call gsetr (gd, G_ASPECT, aspect)
+
+	call sfree (sp)
+end
diff --git a/noao/digiphot/apphot/aputil/aplucy.x b/noao/digiphot/apphot/aputil/aplucy.x
new file mode 100644
index 00000000..f70a39ed
--- /dev/null
+++ b/noao/digiphot/apphot/aputil/aplucy.x
@@ -0,0 +1,47 @@
+# AP_LUCY_SMOOTH - Lucy smooth a histogram with a box shaped function.
+
+procedure ap_lucy_smooth (hgm, shgm, nbins, nker, iter)
+
+real	hgm[ARB]		# the original histogram
+real	shgm[ARB]		# smoothed histogram
+int	nbins			# length of the histogram
+int	nker			# length of box kernel
+int	iter			# number of iterations
+
+int	i, j, nsmooth
+pointer	sp, work1, work2
+real	ap_diverr()
+
+begin
+	# Allocate space and clear the work arrays.
+	call smark (sp)
+	call salloc (work1, nbins, TY_REAL)
+	call salloc (work2, nbins, TY_REAL)
+	call aclrr (Memr[work1], nbins)
+	call aclrr (Memr[work2], nbins)
+	call aclrr (shgm, nbins)
+
+	# Compute the smoothing length and initialize output array.
+	nsmooth = nbins - nker + 1
+	#call ap_aboxr (hgm, shgm[1+nker/2], nsmooth, nker)
+	call ap_sboxr (hgm, shgm, nbins, nker)
+
+	# Iterate on the solution.
+	do i = 1, iter {
+	    #call ap_aboxr (shgm, Memr[work1+nker/2], nsmooth, nker)
+	    call ap_sboxr (shgm, Memr[work1], nbins, nker)
+	    #do j = 1 + nker / 2, nsmooth + nker / 2 {
+	    do j = 1, nbins {
+		if (Memr[work1+j-1] == 0.0)
+		    Memr[work1+j-1] = ap_diverr ()
+		else
+		    #Memr[work1+j-1] = hgm[j] / Memr[work1+j-1]
+		    Memr[work1+j-1] = shgm[j] / Memr[work1+j-1]
+	    }
+	    #call ap_aboxr (Memr[work1], Memr[work2+nker/2], nsmooth, nker)
+	    call ap_sboxr (Memr[work1], Memr[work2], nbins, nker)
+	    call amulr (shgm, Memr[work2], shgm, nbins)
+	}
+
+	call sfree (sp)
+end
diff --git a/noao/digiphot/apphot/aputil/apmapr.x b/noao/digiphot/apphot/aputil/apmapr.x
new file mode 100644
index 00000000..0d4159ab
--- /dev/null
+++ b/noao/digiphot/apphot/aputil/apmapr.x
@@ -0,0 +1,21 @@
+# APMAPR -- Vector linear transformation.  Map the range of pixel values
+# a1, a2 from a into the range b1, b2 into b.  It is assumed that a1 < a2
+# and b1 < b2.
+
+real procedure apmapr (a, a1, a2, b1, b2)
+
+real	a		# the value to be mapped
+real	a1, a2		# the numbers specifying the input data range
+real	b1, b2		# the numbers specifying the output data range
+
+real	minout, maxout, aoff, boff
+real	scalar
+
+begin
+	scalar = (real (b2) - real (b1)) / (real (a2) - real (a1))
+	minout = min (b1, b2)
+	maxout = max (b1, b2)
+	aoff = a1
+	boff = b1
+	return (max(minout, min(maxout, real((a - aoff) * scalar) + boff)))
+end
diff --git a/noao/digiphot/apphot/aputil/apmeds.x b/noao/digiphot/apphot/aputil/apmeds.x
new file mode 100644
index 00000000..0a206c8a
--- /dev/null
+++ b/noao/digiphot/apphot/aputil/apmeds.x
@@ -0,0 +1,154 @@
+# APSMED -- Compute the averaged median given a sorted array and an averaging
+# half width.
+
+real procedure apsmed (pix, index, npts, medcut)
+
+real	pix[ARB]	# array of sky pixels
+int	index[ARB]	# sorted index array
+int	npts		# number of pixels
+int	medcut		# averaging half width
+
+int	med, j, nmed, medlo, medhi
+real	sumed
+
+begin
+	med = (npts + 1) / 2
+	if (mod (med, 2) == 1) {
+	    medlo = max (1, med - medcut)
+	    medhi = min (npts, med + medcut)
+	} else {
+	    medlo = max (1, med - medcut)
+	    medhi = min (npts, med + medcut + 1)
+	}
+
+	sumed = 0.0
+	nmed = 0
+	do j = medlo, medhi {
+	    sumed = sumed + pix[index[j]]
+	    nmed = nmed + 1
+	}
+
+	return (sumed / nmed)
+end
+
+
+# APIMED -- Compute the index of new median value. Weight is an arbitrary
+# weight array which is assumed to be zero if the pixels has been rejected
+# and is positive otherwise.
+
+int procedure apimed (weight, index, lo, hi, nmed)
+
+real	weight[ARB]		# array of weights
+int	index[ARB]		# array of sorted indices
+int	lo, hi			# ranges of weights
+int	nmed			# number of good sky pixels
+
+int	npts, med
+
+begin
+	npts = 0
+	for (med = lo; med <= hi && npts < nmed; med = med + 1) {
+	    if (weight[index[med]] > 0.0)
+		npts = npts + 1
+	}
+
+	if (npts == 0)
+	    return (0)
+	else
+	    return (med)
+end
+
+
+# APWSMED -- Compute the new averaged median given a sorted input array,
+# an averaging half-width, and assuming that there has been pixel rejection.
+
+real procedure apwsmed (pix, index, weight, npix, med, medcut)
+
+real	pix[ARB]		# pixel values array
+int	index[ARB]		# sorted indices array
+real	weight[ARB]		# the weights array
+int	npix			# number of pixels
+int	med			# index of median value
+int	medcut			# of median cut
+
+int	j, nmed, maxmed
+real	sumed
+
+begin
+	sumed = pix[index[med]]
+	if (mod (med, 2) == 1)
+	    maxmed = 2 * medcut + 1 
+	else
+	    maxmed = 2 * medcut + 2 
+
+	nmed = 1
+	for (j = med - 1; j >= 1; j = j - 1) {
+	    if (nmed >=  medcut + 1)
+		break
+	    if (weight[index[j]] > 0.0) {
+		sumed = sumed + pix[index[j]]
+		nmed = nmed + 1
+	    }
+	}
+	for (j = med + 1; j <= npix; j = j + 1) {
+	    if (nmed >= maxmed)
+		break
+	    if (weight[index[j]] > 0.0) {
+		sumed = sumed + pix[index[j]]
+		nmed = nmed + 1
+	    }
+	}
+
+	return (sumed / nmed)
+end
+
+
+# APMEDR -- Vector median selection. The selection is carried out in a temporary
+# array, leaving the input vector unmodified.  Especially demanding applications
+# may wish to call the asok routine directory to avoid the call to the memory
+# allocator.
+
+real procedure apmedr (a, index, npix)
+
+real	a[ARB]				# input array of values
+int	index[ARB]			# sorted index array
+int	npix				# number of pixels
+
+int	i
+pointer	sp, aa
+real	median
+real	asokr()		# select the Kth smallest element from A
+
+begin
+	switch (npix) {
+	case 1, 2:
+	    return (a[1])
+
+	case 3:
+	    if (a[1] < a[2]) {
+		if (a[2] < a[3])
+		    return (a[2])
+		else if (a[1] < a[3])
+		    return (a[3])
+		else
+		    return (a[1])
+	    } else {
+		if (a[2] > a[3])
+		    return (a[2])
+		else if (a[1] < a[3])
+		    return (a[1])
+		else
+		    return (a[3])
+	    }
+
+	default:
+	    call smark (sp)
+	    call salloc (aa, npix, TY_REAL)
+	    do i = 1, npix
+		Memr[aa+i-1] = a[index[i]]
+	    median = asokr (Memr[aa], npix, (npix + 1) / 2)
+	    call sfree (sp)
+
+	    return (median)
+	}
+end
diff --git a/noao/digiphot/apphot/aputil/apmoments.x b/noao/digiphot/apphot/aputil/apmoments.x
new file mode 100644
index 00000000..5a1bfbd4
--- /dev/null
+++ b/noao/digiphot/apphot/aputil/apmoments.x
@@ -0,0 +1,138 @@
+# APMOMENTS -- Procedure to compute the first three moments of a
+# distribution given the appropriate sums.
+
+procedure apmoments (sumpx, sumsqpx, sumcbpx, npix, sky_zero, mean, sigma, skew)
+
+double	sumpx		# sum of the pixels
+double	sumsqpx		# sum of pixels squared
+double	sumcbpx		# sum of cubes of pixels
+int	npix		# number of pixels
+real	sky_zero	# sky zero point for moment analysis
+real	mean		# mean of pixels
+real	sigma		# sigma of pixels
+real	skew		# skew of pixels
+
+double	dmean, dsigma, dskew
+
+begin
+	# Recompute the moments.
+	dmean = sumpx / npix
+	dsigma = sumsqpx / npix - dmean ** 2
+	if (dsigma <= 0.0d0) {
+	    sigma = 0.0
+	    skew = 0.0
+	} else {
+	    dskew = sumcbpx / npix - 3.0d0 * dmean * dsigma - dmean ** 3
+	    sigma = sqrt (dsigma)
+	    if (dskew < 0.0d0)
+	        skew = - (abs (dskew) ** (1.0d0 / 3.0d0))
+	    else if (dskew > 0.0d0)
+	        skew = dskew ** (1.0d0 / 3.0d0)
+	    else
+		skew = 0.0
+	}
+	mean = sky_zero + dmean
+end
+
+
+# APFMOMENTS -- Procedure to compute the first three moments of a distribution
+# given the data. The sums are returned as well.
+
+procedure apfmoments (pix, npix, sky_zero, sumpx, sumsqpx, sumcbpx, mean,
+	sigma, skew)
+
+real	pix[npix]	# array of pixels
+int	npix		# number of pixels
+real	sky_zero	# sky zero point for moment analysis
+double	sumpx		# sum of the pixels
+double	sumsqpx		# sum of pixels squared
+double	sumcbpx		# sum of cubes of pixels
+real	mean		# mean of pixels
+real	sigma		# sigma of pixels
+real	skew		# skew of pixels
+
+double	dpix, dmean, dsigma, dskew
+int	i
+
+begin
+	# Zero and accumulate the sums.
+	sumpx = 0.0d0
+	sumsqpx = 0.0d0
+	sumcbpx = 0.0d0
+	do i = 1, npix {
+	    dpix = pix[i] - sky_zero
+	    sumpx = sumpx + dpix
+	    sumsqpx = sumsqpx + dpix * dpix
+	    sumcbpx = sumcbpx + dpix * dpix * dpix
+	}
+
+	# Compute the moments.
+	dmean = sumpx / npix
+	dsigma = sumsqpx / npix - dmean ** 2
+	if (dsigma <= 0.0) {
+	    sigma = 0.0
+	    skew = 0.0
+	} else {
+	    dskew = sumcbpx / npix - 3.0d0 * dmean * dsigma - dmean ** 3
+	    sigma = sqrt (dsigma)
+	    if (dskew < 0.0d0)
+	        skew = - (abs (dskew) ** (1.0d0 / 3.0d0))
+	    else if (dskew > 0.0d0)
+	        skew = dskew ** (1.0d0 / 3.0d0)
+	    else
+		skew = 0.0
+	}
+	mean = sky_zero + dmean
+end
+
+
+# APFIMOMENTS -- Procedure to compute the first three moments of a distribution
+# given the data. The sums are returned as well.
+
+procedure apfimoments (pix, index, npix, sky_zero, sumpx, sumsqpx, sumcbpx,
+	mean, sigma, skew)
+
+real	pix[ARB]	# array of pixels
+int	index[ARB]	# the index array
+int	npix		# number of pixels
+real	sky_zero	# sky zero point for moment analysis
+double	sumpx		# sum of the pixels
+double	sumsqpx		# sum of pixels squared
+double	sumcbpx		# sum of cubes of pixels
+real	mean		# mean of pixels
+real	sigma		# sigma of pixels
+real	skew		# skew of pixels
+
+double	dpix, dmean, dsigma, dskew
+int	i
+
+begin
+	# Zero and accumulate the sums.
+	sumpx = 0.0d0
+	sumsqpx = 0.0d0
+	sumcbpx = 0.0d0
+	do i = 1, npix {
+	    dpix = pix[index[i]] - sky_zero
+	    sumpx = sumpx + dpix
+	    sumsqpx = sumsqpx + dpix * dpix
+	    sumcbpx = sumcbpx + dpix * dpix * dpix
+	}
+
+	# Compute the moments.
+	dmean = sumpx / npix
+	dsigma = sumsqpx / npix - dmean ** 2
+	if (dsigma <= 0.0d0) {
+	    sigma = 0.0
+	    skew = 0.0
+	} else {
+	    dskew = sumcbpx / npix - 3.0d0 * dmean * dsigma - dmean ** 3
+	    sigma = sqrt (dsigma)
+	    if (dskew < 0.0d0)
+	        skew = - (abs (dskew) ** (1.0d0 / 3.0d0))
+	    else if (dskew > 0.0d0)
+	        skew = dskew ** (1.0d0 / 3.0d0)
+	    else
+		skew = 0.0
+	}
+	mean = dmean + sky_zero
+end
diff --git a/noao/digiphot/apphot/aputil/apnlfuncs.x b/noao/digiphot/apphot/aputil/apnlfuncs.x
new file mode 100644
index 00000000..290b9e90
--- /dev/null
+++ b/noao/digiphot/apphot/aputil/apnlfuncs.x
@@ -0,0 +1,375 @@
+
+include	<math.h>
+include	<mach.h>
+
+# CGAUSS1D - Compute the value of a 1-D Gaussian function  on a constant
+# background.
+
+procedure cgauss1d (x, nvars, p, np, z)
+
+real	x[ARB]		# variables, x[1] = position coordinate
+int	nvars		# the number of variables, not used
+real	p[ARB]		# p[1]=amplitude p[2]=center p[3]=variance p[4]=sky
+int	np		# number of parameters np = 4
+real	z		# function return
+
+real	r2
+
+begin
+	if (p[3] == 0.)
+	    r2 = 36.0
+	else
+	    r2 = (x[1] - p[2]) ** 2 / (2. * p[3])
+	if (abs (r2) > 25.0)
+	    z = p[4]
+	else
+	    z = p[1] * exp (-r2) + p[4]
+end
+
+
+# CDGAUSS1D -- Compute the value a 1-D Gaussian profile on a constant
+# background and its derivatives.
+
+procedure cdgauss1d (x, nvars p, dp, np, z, der)
+
+real	x[ARB]		# variables, x[1] = position coordinate
+int	nvars		# the number of variables, not used
+real	p[ARB]		# p[1]=amplitude, p[2]=center, p[3]=sky p[4]=variance
+real	dp[ARB]		# parameter derivatives
+int	np		# number of parameters np=4
+real	z		# function value
+real	der[ARB]	# derivatives
+
+real	dx, r2
+
+begin
+	dx = x[1] - p[2]
+	if (p[3] == 0.)
+	    r2 = 36.0
+	else
+	    r2 = dx * dx / (2.0 * p[3])
+	if (abs (r2) > 25.0) {
+	    z = p[4]
+	    der[1] = 0.0
+	    der[2] = 0.0
+	    der[3] = 0.0
+	    der[4] = 1.0
+	} else {
+	    der[1] = exp (-r2)
+	    z = p[1] * der[1]
+	    der[2] = z * dx / p[3]
+	    der[3] = z * r2 / p[3]
+	    der[4] = 1.0
+	    z = z + p[4]
+	}
+end
+
+
+# GAUSSR -- Compute the value of a 2-D radially symmetric Gaussian profile
+# which is assumed to be sitting on a constant background.
+
+# Parameter Allocation:
+# 1	Amplitude
+# 2	X-center
+# 3	Y-center
+# 4	Variance
+# 5	Sky
+
+procedure gaussr (x, nvars, p, np, z)
+
+real	x[ARB]		# the input variables
+int	nvars		# the number of variables
+real	p[np]		# parameter vector
+int	np		# number of parameters
+real	z		# function return
+
+real	dx, dy, r2
+
+begin
+	dx = x[1] - p[2]
+	dy = x[2] - p[3]
+	if (p[4] == 0.)
+	    r2 = 36.0
+	else
+	    r2 = (dx * dx + dy * dy) / (2.0 * p[4])
+	if (abs (r2) > 25.0)
+	    z = p[5]
+	else
+	    z = p[1] * exp (- r2) + p[5]
+end
+
+
+# DGAUSSR -- Compute the value of a 2-D Gaussian profile and its derivatives
+# which assumed to be sitting on top of a constant background.
+
+procedure dgaussr (x, nvars, p, dp, np, z, der)
+
+real	x[ARB]		# the input variables
+int	nvars		# the number of variables
+real	p[np]		# parameter vector
+real	dp[np]		# dummy array of parameter increments
+int	np		# number of parameters
+real	z		# function return
+real	der[np]		# derivatives
+
+real	dx, dy, r2
+
+begin
+	dx = x[1] - p[2]
+	dy = x[2] - p[3]
+	if (p[4] == 0.)
+	    r2 = 36.0
+	else
+	    r2 = (dx * dx + dy * dy) / (2.0 * p[4])
+
+	if (abs (r2) > 25.0) {
+	    z = p[5]
+	    der[1] = 0.0
+	    der[2] = 0.0
+	    der[3] = 0.0
+	    der[4] = 0.0
+	    der[5] = 1.0
+	} else {
+	    der[1] = exp (-r2)
+	    z = p[1] * der[1]
+	    der[2] = z * dx / p[4] 
+	    der[3] = z * dy / p[4]
+	    der[4] = z * r2 / p[4]
+	    z = z + p[5]
+	    der[5] = 1.0
+	}
+end
+
+
+# ELGAUSS -- Compute the value of a 2-D elliptical Gaussian function which
+# is assumed to be sitting on top of a constant background.
+
+# Parameter Allocation:
+# 1	Amplitude
+# 2	X-center
+# 3	Y-center
+# 4	Variance-x
+# 5	Variance-y
+# 6	Theta-rotation
+# 7	Sky
+
+procedure elgauss (x, nvars, p, np, z)
+
+real	x[ARB]		# input variables, x[1] = x, x[2] = y
+int	nvars		# number of variables, not used
+real	p[np]		# parameter vector
+int	np		# number of parameters
+real	z		# function return
+
+real	dx, dy, crot, srot, xt, yt, r2
+
+begin
+	dx = x[1] - p[2]
+	dy = x[2] - p[3]
+	crot = cos (p[6])
+	srot = sin (p[6])
+	xt = (dx * crot + dy * srot)
+	yt = (-dx * srot + dy * crot)
+	if (p[4] == 0. || p[5] == 0.)
+	    r2 = 36.0
+	else
+	    r2 = (xt ** 2 / p[4] + yt ** 2 / p[5]) / 2.0
+	if (abs (r2) > 25.0)
+	    z = p[7]
+	else
+	    z = p[1] * exp (-r2) + p[7]
+end
+
+
+# DELGAUSS -- Compute the value of a 2-D elliptical Gaussian assumed to
+# sitting on top of a constant background and its derivatives.
+
+procedure delgauss (x, nvars, p, dp, np, z, der)
+
+real	x[ARB]		# input variables, x[1] = x, x[2] = y
+int	nvars		# number of variables, not used
+real	p[np]		# parameter vector
+real	dp[np]		# delta of parameters
+int	np		# number of parameters
+real	z		# function value
+real	der[np]		# function return
+
+real	crot, srot, crot2, srot2, sigx2, sigy2, a, b, c
+real	dx, dy, dx2, dy2, r2
+
+begin
+	crot = cos (p[6])
+	srot = sin (p[6])
+	crot2 = crot ** 2
+	srot2 = srot ** 2
+	sigx2 = p[4]
+	sigy2 = p[5]
+	if (sigx2 == 0. || sigy2 == 0.)
+	    r2 = 36.0
+	else {
+	    a = (crot2 / sigx2 + srot2 / sigy2)
+	    b = 2.0 * crot * srot * (1.0 / sigx2 - 1.0 /sigy2)
+	    c = (srot2 / sigx2 + crot2 / sigy2)
+
+	    dx = x[1] - p[2]
+	    dy = x[2] - p[3]
+	    dx2 = dx ** 2
+	    dy2 = dy ** 2
+	    r2 = 0.5 * (a * dx2 + b * dx * dy + c * dy2)
+	}
+
+	if (abs (r2) > 25.0) {
+	    z = p[7]
+	    der[1] = 0.0
+	    der[2] = 0.0
+	    der[3] = 0.0
+	    der[4] = 0.0
+	    der[5] = 0.0
+	    der[6] = 0.0
+	    der[7] = 1.0
+	} else {
+	    der[1] = exp (-r2)
+	    z = p[1] * der[1]
+	    der[2] = z * (2.0 * a * dx + b * dy)
+	    der[3] = z * (b * dx + 2.0 * c * dy)
+	    der[4] = z * (crot2 * dx2 + 2.0 * crot * srot * dx * dy +
+	        srot2 * dy2) / (2.0 * sigx2 * sigx2)
+	    der[5] = z * (srot2 * dx2 - 2.0 * crot * srot * dx * dy +
+	        crot2 * dy2) / (2.0 * sigy2 * sigy2)
+	    der[6] = z * (b * dx2 + 2.0 * (c - a) * dx * dy - b * dy2)  
+	    z = z + p[7]
+	    der[7] = 1.0
+	}
+end
+
+
+
+# GAUSS1D -- Compute the profile of a 1d Gaussian with a background value
+# of zero.
+
+procedure gauss1d (x, nvars, p, np, z)
+
+real	x[ARB]		# list of variables, x[1] = position coordinate
+int	nvars		# number of variables
+real	p[ARB]		# p[1]=amplitude p[2]=center p[3]=sigma
+int	np		# number of parameters == 3
+real	z		# function return
+
+real	r
+
+begin
+	if (p[3] == 0.)
+	    r = 6.0
+	else
+	    r = (x[1] - p[2]) / (p[3] * SQRTOF2)
+	if (abs (r) > 5.0)
+	    z = 0.0
+	else
+	    z = p[1] * exp (- r ** 2)
+end
+
+
+# DGAUSS1D -- Compute the function value and derivatives of a 1-D Gaussian
+# function with a background value of zero.
+
+procedure dgauss1d (x, nvars, p, dp, np, z, der)
+
+real	x[ARB]		# list of variables, x[1] = position coordinate
+int	nvars		# number of variables
+real	p[ARB]		# p[1]=amplitude, p[2]=center, p[3]=sigma
+real	dp[ARB]		# parameter derivatives
+int	np		# number of parameters
+real	z		# function value
+real	der[ARB]	# derivatives
+
+real	r
+
+begin
+	if (p[3] == 0.)
+	    r = 6.0
+	else
+	    r = (x[1] - p[2]) / (SQRTOF2 * p[3])
+	if (abs (r) > 5.0) {
+	    z = 0.0
+	    der[1] = 0.0
+	    der[2] = 0.0
+	    der[3] = 0.0
+	} else {
+	    der[1] = exp (- r ** 2)
+	    z = der[1] * p[1]
+	    der[2] = z * r  * SQRTOF2 / p[3]
+	    der[3] = der[2] * SQRTOF2 * r
+	}
+end
+
+
+# GAUSSKEW - Compute the value of a 1-D skewed Gaussian profile.
+# The background value is assumed to be zero.
+
+procedure gausskew (x, nvars, p, np, z)
+
+real	x[ARB]		# list of variables, x[1] = position coordinate
+int	nvars		# number of variables, not used
+real	p[ARB]		# p[1]=amplitude p[2]=center p[3]=variance p[4]=skew
+int	np		# number of parameters == 3
+real	z		# function return
+
+real	dx, r2, r3
+
+begin
+	dx = (x[1] - p[2])
+	if (p[3] == 0.)
+	    r2 = 36.0
+	else {
+	    r2 = dx ** 2 / (2.0 * p[3])
+	    r3 = r2 * dx / sqrt (2.0 * abs (p[3])) 
+	}
+	if (abs (r2) > 25.0)
+	    z = 0.0
+	else
+	    z = (1.0 + p[4] * r3) * p[1] * exp (-r2)
+end
+
+
+# DGAUSSKEW -- Compute the value of a 1-D skewed Gaussian and its derivatives. 
+# The background value is assumed to be zero.
+
+procedure dgausskew (x, nvars, p, dp, np, z, der)
+
+real	x[ARB]		# list of variables, x[1] = position coordinate
+int	nvars		# number of variables, not used
+real	p[ARB]		# p[1]=amplitude, p[2]=center, p[3]=variance, p[4]=skew
+real	dp[ARB]		# parameter derivatives
+int	np		# number of parameters
+real	z		# function value
+real	der[ARB]	# derivatives
+
+real	dx, d1, d2, d3, r, r2, r3, rint
+
+begin
+	dx = x[1] - p[2]
+	if (p[3] == 0.)
+	    r2 = 36.0
+	else
+	    r2 = dx ** 2 / (2.0 * p[3])
+	if (abs (r2) > 25.0) {
+	    z = 0.0
+	    der[1] = 0.0
+	    der[2] = 0.0
+	    der[3] = 0.0
+	    der[4] = 0.0
+	} else {
+	    r = dx / sqrt (2.0 * abs (p[3]))
+	    r3 = r2 * r
+	    d1 = exp (-r2)
+	    z = d1 * p[1]
+	    d2 = z * dx / p[3]
+	    d3 = z * r2 / p[3]
+	    rint = 1.0 + p[4] * r3
+	    der[1] = d1 * rint
+	    der[2] = d2 * (rint - 1.5 * p[4] * r)
+	    der[3] = d3 * (rint - 1.5 * p[4] * r)
+	    der[4] = z * r3
+	    z = z * rint
+	}
+end
diff --git a/noao/digiphot/apphot/aputil/appcache.x b/noao/digiphot/apphot/aputil/appcache.x
new file mode 100644
index 00000000..d1059583
--- /dev/null
+++ b/noao/digiphot/apphot/aputil/appcache.x
@@ -0,0 +1,83 @@
+include <imhdr.h>
+include <imset.h>
+
+# AP_MEMSTAT -- Figure out if there is  enough memory to cache the image
+# pixels. If it is necessary to request more memory and the memory is
+# avalilable return YES otherwise return NO.
+
+int procedure ap_memstat (cache, req_size, old_size)
+
+int	cache			#I cache memory ?
+int	req_size		#I the requested working set size in chars 
+int	old_size		#O the original working set size in chars 
+
+int	cur_size, max_size
+int	begmem()
+
+begin
+        # Find the default working set size.
+        cur_size = begmem (0, old_size, max_size)
+
+	# If cacheing is disabled return NO regardless of the working set size.
+	if (cache == NO)
+	    return (NO)
+
+	# If the requested working set size is less than the current working
+	# set size return YES.
+	if (req_size <= cur_size)
+	    return (YES)
+
+	# Reset the current working set size.
+	cur_size = begmem (req_size, old_size, max_size)
+	if (req_size <= cur_size) {
+	    return (YES)
+	} else {
+	    return (NO)
+	}
+end
+
+
+# AP_PCACHE -- Cache the image pixels im memory by resetting the  default image
+# buffer size. If req_size is INDEF the size of the image is used to determine
+# the size of the image i/o buffers.
+
+procedure ap_pcache (im, req_size, buf_size)
+
+pointer im                      #I the input image point
+int     req_size                #I the requested working set size in chars
+int	buf_size		#O the new image buffer size
+
+int     def_size, new_imbufsize
+int     sizeof(), imstati()
+
+begin
+	# Find the default buffer size.
+	def_size = imstati (im, IM_BUFSIZE)
+
+        # Return if the image is not 2-dimensional.
+        if (IM_NDIM(im) != 2) {
+	    buf_size = def_size
+            return
+	}
+
+        # Compute the new required image i/o buffer size in chars.
+        if (IS_INDEFI(req_size)) {
+            new_imbufsize = IM_LEN(im,1) * IM_LEN(im,2) *
+                sizeof (IM_PIXTYPE(im))
+        } else {
+            new_imbufsize = req_size
+        }
+
+	# If the default image i/o buffer size is already bigger than
+	# the requested size do nothing.
+	if (def_size >= new_imbufsize) {
+	    buf_size = def_size
+	    return
+	}
+
+        # Reset the image i/o buffer.
+        call imseti (im, IM_BUFSIZE, new_imbufsize)
+        call imseti (im, IM_BUFFRAC, 0)
+	buf_size = new_imbufsize
+        return 
+end
diff --git a/noao/digiphot/apphot/aputil/aprmwhite.x b/noao/digiphot/apphot/aputil/aprmwhite.x
new file mode 100644
index 00000000..071ea7ed
--- /dev/null
+++ b/noao/digiphot/apphot/aputil/aprmwhite.x
@@ -0,0 +1,22 @@
+include <ctype.h>
+
+# AP_RMWHITE -- Remove whitespace from a string.
+
+procedure ap_rmwhite (instr, outstr, maxch)
+
+char	instr[ARB]		# the input string
+char	outstr[ARB]		# the output string, may be the same as instr
+int	maxch			# maximum number of characters in outstr
+
+int	ip, op
+
+begin
+	op = 1
+	for (ip = 1; (instr[ip] != EOS) && (op <= maxch); ip = ip + 1) {
+	    if (IS_WHITE(instr[ip]))
+		next
+	    outstr[op] = instr[ip]
+	    op = op + 1
+	}
+	outstr[op] = EOS
+end
diff --git a/noao/digiphot/apphot/aputil/aptopt.x b/noao/digiphot/apphot/aputil/aptopt.x
new file mode 100644
index 00000000..0e5bc212
--- /dev/null
+++ b/noao/digiphot/apphot/aputil/aptopt.x
@@ -0,0 +1,232 @@
+include	<mach.h>
+include	<math.h>
+
+# APTOPT - One-dimensional centering routine using repeated convolutions to
+# locate image center.
+
+define	MAX_SEARCH	3	# Max initial search steps
+
+int procedure aptopt (data, npix, center, sigma, tol, maxiter, ortho)
+
+real	data[ARB]		# initial data
+int	npix			# number of pixels
+real	center			# initial guess at center
+real	sigma			# sigma of Gaussian
+real	tol			# gap tolerance for sigma
+int	maxiter			# maximum number of iterations
+int	ortho			# orthogonalize weighting vector
+
+int	i, iter
+pointer	sp, wgt
+real	newx, x[3], news, s[3], delx
+real	adotr(), apqzero()
+
+begin
+	if (sigma <= 0.0)
+	    return (-1)
+
+	# Allocate working space.
+	call smark (sp)
+	call salloc (wgt, npix, TY_REAL)
+
+	# Initialize.
+	x[1] = center
+	call mkt_prof_derv (Memr[wgt], npix, x[1], sigma, ortho)
+	s[1] = adotr (Memr[wgt], data, npix)
+	#if (abs (s[1]) <= EPSILONR) {
+	if (s[1] == 0.0) {
+	    center = x[1]
+	    call sfree (sp)
+	    return (0)
+	} else
+	    s[3] = s[1]
+
+	# Search for the correct interval.
+	for (i = 1; (s[3] * s[1] >= 0.0) && (i <= MAX_SEARCH); i = i + 1) {
+	    s[3] = s[1]
+	    x[3] = x[1]
+	    x[1] = x[3] + sign (sigma, s[3])
+	    call mkt_prof_derv (Memr[wgt], npix, x[1], sigma, ortho)
+	    s[1] = adotr (Memr[wgt], data, npix)
+	    #if (abs (s[1]) <= EPSILONR) {
+	    if (s[1] == 0.0) {
+		center = x[1]
+		call sfree (sp)
+		return (0)
+	    }
+	}
+
+	# Location not bracketed.
+	if (s[3] * s[1] > 0.0) {
+	    call sfree (sp)
+	    return (-1)
+	}
+
+	# Intialize the quadratic search.
+	delx = x[1] - x[3]
+	x[2] = x[3] - s[3] * delx / (s[1] - s[3])
+	call mkt_prof_derv (Memr[wgt], npix, x[2], sigma, ortho)
+	s[2] = adotr (Memr[wgt], data, npix)
+	#if (abs (s[2]) <= EPSILONR) {
+	if (s[2] == 0.0) {
+	    center = x[2]
+	    call sfree (sp)
+	    return (1)
+	}
+
+	# Search quadratically.
+	for (iter = 2; iter <= maxiter; iter = iter + 1)  {
+
+	    # Check for completion.
+	    #if (abs (s[2]) <= EPSILONR)
+	    if (s[2] == 0.0)
+		break
+	    if (abs (x[2] - x[1]) <= tol)
+		break
+	    if (abs (x[3] - x[2]) <= tol)
+		break
+
+	    # Compute new intermediate value.
+	    newx = x[1] + apqzero (x, s)
+	    call mkt_prof_derv (Memr[wgt], npix, newx, sigma, ortho)
+	    news = adotr (Memr[wgt], data, npix)
+
+	    if (s[1] * s[2] > 0.0) {
+		s[1] = s[2]
+		x[1] = x[2]
+		s[2] = news
+		x[2] = newx
+	    } else {
+		s[3] = s[2]
+		x[3] = x[2]
+		s[2] = news
+		x[2] = newx
+	    }
+	}
+
+	# Evaluate the center.
+	center = x[2]
+	call sfree (sp)
+	return (iter)
+end
+
+
+# AP_TPROFDER -- Procedure to estimate the approximating triangle function
+# and its derivatives.
+
+procedure ap_tprofder (data, der, npix, center, sigma, ampl)
+
+real	data[ARB]		# input data
+real	der[ARB]		# derivatives
+int	npix			# number of pixels
+real	center			# center of input Gaussian function
+real	sigma			# sigma of input Gaussian function
+real	ampl			# amplitude
+
+int	i
+real	x, xabs, width
+
+begin
+	width = sigma * 2.35
+	do i = 1, npix {
+	    x = (i - center) / width
+	    xabs = abs (x)
+	    if (xabs <= 1.0) {
+		data[i] = ampl * (1.0 - xabs)
+		der[i] = x * data[i]
+	    } else {
+		data[i] = 0.0
+		der[i] = 0.0
+	    }
+	}
+end
+
+
+# MKT_PROF_DERV - Make orthogonal profile derivative vector.
+
+procedure mkt_prof_derv (weight, npix, center, sigma, norm)
+
+real	weight[ARB]	# input weight
+int	npix		# number of pixels
+real	center		# center
+real	sigma		# center
+int	norm		# orthogonalise weight
+
+pointer	sp, der
+real	coef
+real	asumr(),  adotr()
+
+begin
+	call smark (sp)
+	call salloc (der, npix, TY_REAL)
+
+	# Fetch the weighting function and derivatives.
+	call ap_tprofder (Memr[der], weight, npix, center, sigma, 1.0)
+
+	if (norm == YES) {
+
+	    # Make orthogonal to level background.
+	    coef = -asumr (weight, npix) / npix
+	    call aaddkr (weight, coef, weight, npix)
+	    coef = -asumr (Memr[der], npix) / npix
+	    call aaddkr (Memr[der], coef, Memr[der], npix)
+
+	    # Make orthogonal to profile vector.
+	    coef = adotr (Memr[der], Memr[der], npix)
+	    if (coef <= 0.0)
+		coef = 1.0
+	    else
+	        coef = adotr (weight, Memr[der], npix) / coef
+	    call amulkr (Memr[der], coef, Memr[der], npix)
+	    call asubr (weight, Memr[der], weight, npix)
+
+	    # Normalize the final vector.
+	    coef = adotr (weight, weight, npix)
+	    if (coef <= 0.0)
+		coef = 1.0
+	    else
+	        coef = sqrt (1.0 / coef)
+	    call amulkr (weight, coef, weight, npix)
+	}
+
+	call sfree (sp)
+end
+
+define	QTOL	.125
+
+# APQZERO - Solve for the root of a quadratic function defined by three
+# points.
+
+real procedure apqzero (x, y)
+
+real	x[3]
+real	y[3]
+
+real	a, b, c, det, dx
+real	x2, x3, y2, y3
+
+begin
+	# Compute the determinant.
+	x2 = x[2] - x[1]
+	x3 = x[3] - x[1]
+	y2 = y[2] - y[1]
+	y3 = y[3] - y[1]
+	det = x2 * x3 * (x2 - x3)
+
+	# Compute the shift in x.
+	#if (abs (det) > 100.0 * EPSILONR) {
+	if (abs (det) > 0.0) {
+	    a = (x3 * y2 - x2 * y3) / det
+	    b = - (x3 * x3 * y2 - x2 * x2 * y3) / det
+	    c =  a * y[1] / (b * b)
+	    if (abs (c) > QTOL)
+		dx = (-b / (2.0 * a)) * (1.0 - sqrt (1.0 - 4.0 * c))
+	    else
+		dx = - (y[1] / b) * (1.0 + c)
+	    return (dx)
+	#} else if (abs (y3) > EPSILONR)
+	} else if (abs (y3) > 0.0)
+	    return (-y[1] * x3 / y3)
+	else
+	    return (0.0)
+end
diff --git a/noao/digiphot/apphot/aputil/apvectors.x b/noao/digiphot/apphot/aputil/apvectors.x
new file mode 100644
index 00000000..7c9dbcfd
--- /dev/null
+++ b/noao/digiphot/apphot/aputil/apvectors.x
@@ -0,0 +1,515 @@
+
+# AP_ABOXR -- Vector boxcar smooth. The input vector is convolved with a
+# uniform kernel of length knpix. Boundary extension is assumed to have
+# been performed on the input array before entering the routine.
+
+procedure ap_aboxr (in, out, npix, knpix)
+
+real	in[npix+knpix-1]	# input array with boundary exntension
+real	out[npix]		# output array
+int	npix			# number of pixels
+int	knpix			# length of the kernel
+
+int	i
+real	sum
+
+begin
+	sum = 0.0
+	do i = 1, knpix - 1
+	    sum = sum + in[i]
+
+	do i = 1, npix {
+	    sum = sum + in[i+knpix-1] 
+	    out[i] = sum / knpix
+	    sum = sum - in[i]
+	}
+
+end
+
+
+# AP_SBOXR -- Boxcar smooth for a vector that has not been boundary
+# extended.
+
+procedure ap_sboxr (in, out, npix, nsmooth)
+
+real	in[npix]		# input array 
+real	out[npix]		# output array
+int	npix			# number of pixels
+int	nsmooth			# half width of smoothing box
+
+int	i, j, ib, ie, ns
+real	sum
+
+begin
+	ns = 2 * nsmooth + 1
+	do i = 1, npix {
+	    ib = max (i - nsmooth, 1)
+	    ie = min (i + nsmooth, npix)
+	    sum = 0.0
+	    do j = ib, ie
+	        sum = sum + in[j]
+	    out[i] = sum / ns
+	}
+end
+
+
+# AP_ALIMR -- Compute the maximum and minimum data values and indices of a
+# 1D array.
+
+procedure ap_alimr (data, npts, mindat, maxdat, imin, imax)
+
+real	data[npts]	# data array
+int	npts		# number of points
+real	mindat, maxdat	# min and max data value
+int	imin, imax	# indices of min and max data values
+
+int	i
+
+begin
+	imin = 1
+	imax = 1
+	mindat = data[1]
+	maxdat = data[1]
+
+	do i = 2, npts {
+	    if (data[i] > maxdat) {
+		imax = i
+		maxdat = data[i]
+	    }
+	    if (data[i] < mindat) {
+		imin = i
+		mindat = data[i]
+	    }
+	}
+end
+
+
+# AP_IALIMR -- Compute the maximum and minimum data values of a 1D indexed
+# array
+
+procedure ap_ialimr (data, index, npts, mindat, maxdat)
+
+real	data[npts]	# data array
+int	index[npts]	# index array
+int	npts		# number of points
+real	mindat, maxdat	# min and max data value
+
+int	i
+
+begin
+	mindat = data[index[1]]
+	maxdat = data[index[1]]
+	do i = 2, npts {
+	    if (data[index[i]] > maxdat)
+		maxdat = data[index[i]]
+	    if (data[index[i]] < mindat)
+		mindat = data[index[i]]
+	}
+end
+
+
+# AP_ASUMR - Compute the sum of an index sorted array
+
+real procedure ap_asumr (data, index, npts)
+
+real	data[npts]	# data array
+int	index[npts]	# index array
+int	npts		# number of points
+
+double	sum
+int	i
+
+begin
+	sum = 0.0d0
+	do i = 1, npts
+	    sum = sum + data[index[i]]
+	return (real (sum))
+end
+
+
+# AP_AMAXEL -- Find the maximum value and its index of a 1D array.
+
+procedure ap_amaxel (a, npts, maxdat, imax)
+
+real	a[ARB]		# the data array
+int	npts		# number of points
+real	maxdat		# maximum value
+int	imax		# imdex of max value
+
+int	i
+
+begin
+	maxdat = a[1]
+	imax = 1
+	do i = 2, npts {
+	    if (a[i] > maxdat) {
+		maxdat = a[i]
+		imax = i
+	    }
+	}
+end
+
+
+# AHGMR -- Accumulate the histogram of the input vector.  The output vector
+# hgm (the histogram) should be cleared prior to the first call. The procedure
+# returns the number of data values it could not include in the histogram.
+
+int procedure aphgmr (data, wgt, npix, hgm, nbins, z1, z2)
+
+real 	data[ARB]		# data vector
+real	wgt[ARB]		# weights vector
+int	npix			# number of pixels
+real	hgm[ARB]		# output histogram
+int	nbins			# number of bins in histogram
+real	z1, z2			# greyscale values of first and last bins
+
+real	dz
+int	bin, i, nreject
+
+begin
+	if (nbins < 2)
+	    return (0)
+
+	nreject = 0
+	dz = real (nbins - 1) / real (z2 - z1)
+	do i = 1, npix {
+	    if (data[i] < z1 || data[i] > z2) {
+		nreject = nreject + 1
+		wgt[i] = 0.0
+		next
+	    }
+	    bin = int ((data[i] - z1) * dz) + 1
+	    hgm[bin] = hgm[bin] + 1.0
+	}
+
+	return (nreject)
+end
+
+
+# APHIGMR -- Accumulate the histogram of the input vector.  The output vector
+# hgm (the histogram) should be cleared prior to the first call. The procedure
+# returns the number of data values it could not include in the histogram.
+
+int procedure aphigmr (data, wgt, index, npix, hgm, nbins, z1, z2)
+
+real 	data[ARB]		# data vector
+real	wgt[ARB]		# weights vector
+int	index[ARB]		# index vector
+int	npix			# number of pixels
+real	hgm[ARB]		# output histogram
+int	nbins			# number of bins in histogram
+real	z1, z2			# greyscale values of first and last bins
+
+real	dz
+int	bin, i, nreject
+
+begin
+	if (nbins < 2)
+	    return (0)
+
+	nreject = 0
+	dz = real (nbins - 1) / real (z2 - z1)
+	do i = 1, npix {
+	    if (data[index[i]] < z1 || data[index[i]] > z2) {
+		nreject = nreject + 1
+		wgt[index[i]] = 0.0
+		next
+	    }
+	    bin = int ((data[index[i]] - z1) * dz) + 1
+	    hgm[bin] = hgm[bin] + 1.0
+	}
+
+	return (nreject)
+end
+
+
+# AP_IJTOR -- Compute radius values given the center coordinates, the line
+# number in a the subraster and the length of the subraster x axis.
+
+procedure ap_ijtor (r, nr, line, xc, yc)
+
+real	r[nr]		# array of output r values
+int	nr		# length of r array
+int	line		# line number
+real	xc, yc		# subraster center
+
+int	i
+real	temp
+
+begin
+	temp = (line - yc ) ** 2
+	do i = 1, nr
+	    r[i] = sqrt ((i - xc) ** 2 + temp)
+end
+
+
+# AP_IJTOR2 -- Compute radius values given the center coordinates and the size
+# of the subraster.
+
+procedure ap_ijtor2 (r, nx, ny, xc, yc)
+
+real	r[nx,ARB]	# array of output r values
+int	nx		# x dimension of output array
+int	ny		# Y dimension of output array
+real	xc, yc		# subraster center
+
+int	i, j
+real	temp
+
+begin
+	do j = 1, ny {
+	    temp = (j - yc) ** 2
+	    do i = 1, nx
+	        r[i,j] = sqrt ((i - xc) ** 2 + temp)
+	}
+end
+
+
+# AP_2DALIMR -- Compute min and max values of a 2D array along with
+# the indices of the minimum and maximum pixel.
+
+procedure ap_2dalimr (data, nx, ny, mindat, maxdat, imin, jmin, imax, jmax)
+
+real	data[nx, ny]	# data
+int	nx, ny		# array dimensions
+real	mindat, maxdat	# min, max data values
+int	imin, jmin	# indices of min element
+int	imax, jmax	# indices of max element
+
+int	i, j
+
+begin
+	imin = 1
+	jmin = 1
+	imax = 1
+	jmax = 1
+	mindat = data[1,1]
+	maxdat = data[1,1]
+
+	do j = 2, ny {
+	    do i = 1, nx {
+		if (data[i,j] < mindat) {
+		    imin = i
+		    jmin = j
+		    mindat = data[i,j]
+		} else if (data[i,j] > maxdat) {
+		    imax = i
+		    jmax = j
+		    maxdat = data[i,j]
+		}
+	    }
+	}
+end
+
+
+define	LOGPTR		20			# log2(maxpts) (1e6)
+
+# APQSORT -- Vector Quicksort. In this version the index array is
+# sorted.
+
+procedure apqsort (data, a, b, npix)
+
+real	data[ARB]		# data array
+int	a[ARB], b[ARB]		# index array
+int	npix			# number of pixels
+
+int	i, j, lv[LOGPTR], p, uv[LOGPTR], temp
+real	pivot
+
+begin
+	# Initialize the indices for an inplace sort.
+	do i = 1, npix
+	    a[i] = i
+	call amovi (a, b, npix)
+
+	p = 1
+	lv[1] = 1
+	uv[1] = npix
+	while (p > 0) {
+
+	    # If only one elem in subset pop stack otherwise pivot line.
+	    if (lv[p] >= uv[p])
+		p = p - 1
+	    else {
+		i = lv[p] - 1
+		j = uv[p]
+		pivot = data[b[j]]
+
+		while (i < j) {
+		    for (i=i+1;  data[b[i]] < pivot;  i=i+1)
+			;
+		    for (j=j-1;  j > i;  j=j-1)
+			if (data[b[j]] <= pivot)
+			    break
+		    if (i < j) {		# out of order pair
+			temp = b[j]		# interchange elements
+			b[j] = b[i]
+			b[i] = temp
+		    }
+		}
+
+		j = uv[p]			# move pivot to position i
+		temp = b[j]			# interchange elements
+		b[j] = b[i]
+		b[i] = temp
+
+		if (i-lv[p] < uv[p] - i) {	# stack so shorter done first
+		    lv[p+1] = lv[p]
+		    uv[p+1] = i - 1
+		    lv[p] = i + 1
+		} else {
+		    lv[p+1] = i + 1
+		    uv[p+1] = uv[p]
+		    uv[p] = i - 1
+		}
+
+		p = p + 1			# push onto stack
+	    }
+	}
+end
+
+
+
+# APRECIPROCAL -- Compute the reciprocal of the absolute value of a vector.
+
+procedure apreciprocal (a, b, npts, value)
+
+real	a[ARB]		# the input vector
+real	b[ARB]		# the output vector
+int	npts		# the number of data points
+real	value		# the value to be assigned to b[i] if a[i] = 0
+
+int	i
+
+begin
+	do i = 1, npts {
+	    if (a[i] > 0.0)
+		b[i] = 1.0 / a[i]
+	    else if (a[i] < 0.0)
+		b[i] = - 1.0 / a[i]
+	    else
+		b[i] = value
+	}
+end
+
+
+# AP_WLIMR -- Set the weights of all data points outside a given minimum and
+# maximum values to zero. Compute the minimum and maximum values and their
+# indices of the remaining data.
+
+procedure ap_wlimr (pix, w, npts, datamin, datamax, dmin, dmax, imin, imax)
+
+real	pix[ARB]		# input pixel array
+real	w[ARB]			# weight array
+int	npts			# number of points
+real	datamin			# minimum good data point
+real	datamax			# maximum good data point
+real	dmin			# output data minimum
+real	dmax			# output data maximum
+int	imin			# index of the data minimum
+int	imax			# index of the data maximum
+
+int	i
+real	value
+
+begin
+	dmin = datamax
+	dmax = datamin
+	imin = 1
+	imax = 1
+
+	do i = 1, npts {
+	    value = pix[i]
+	    if ((value < datamin) || (value > datamax)) {
+		w[i] = 0.0
+		next
+	    }
+	    if (value  < dmin) {
+		dmin = value
+		imin = i
+	    } else if (value  > dmax) {
+		dmax = value
+		imax = i
+	    }
+	}
+end
+
+
+# APWSSQR -- Compute the weighted sum of the squares of a vector.
+
+real procedure apwssqr (a, wgt, npts)
+
+real	a[ARB]		# data array
+real	wgt[ARB]	# array of weights
+int	npts		# number of points
+
+int	i
+real	sum
+
+begin
+	sum = 0.0
+	do i = 1, npts
+	    sum = sum + wgt[i] * a[i] * a[i]
+	return (sum)
+end
+
+
+# AP_XYTOR -- Change the single integer coord of the sky pixels to a radial
+# distance value. The integer coordinate is equal to coord = (i - xc + 1) +
+# blklen * (j - yc).
+
+procedure ap_xytor (coords, index, r, nskypix, xc, yc, blklen)
+
+int	coords[ARB]		# coordinate array
+int	index[ARB]		# the index array
+real	r[ARB]			# radial coordinates
+int	nskypix			# number of sky pixels
+real	xc, yc			# center of sky subraster
+int	blklen			# x dimension of sky subraster
+
+int	i
+real	x, y
+
+begin
+	do i = 1, nskypix {
+	    x = real (mod (coords[index[i]], blklen))
+	    if (x == 0)
+		x = real (blklen)
+	    y = (coords[index[i]] - x) / blklen + 1 
+	    r[i] = sqrt ((x - xc) ** 2 + (y - yc) ** 2)
+	}
+end
+
+
+# AP_W1SUR -- Linearly combine two vectors where the weight for the first
+# vectors is 1.0 and is k2 for the second vector
+
+procedure ap_w1sur (a, b, c, npts, k2)
+
+real	a[ARB]		# the first input vector
+real	b[ARB]		# the second input vector
+real	c[ARB]		# the output vector
+int	npts		# number of points
+real	k2		# the weigting factor for the second vector
+
+int	i
+
+begin
+	do i = 1, npts
+	    c[i] = a[i] + k2 * b[i]
+end
+
+
+# AP_INDEX -- Define an index array.
+
+procedure ap_index (index, npix)
+
+int	index[ARB]		# the index array
+int	npix			# the number of pixels
+
+int	i
+
+begin
+	do i = 1, npix
+	    index[i] = i
+end
diff --git a/noao/digiphot/apphot/aputil/mkpkg b/noao/digiphot/apphot/aputil/mkpkg
new file mode 100644
index 00000000..eefeed9d
--- /dev/null
+++ b/noao/digiphot/apphot/aputil/mkpkg
@@ -0,0 +1,25 @@
+# Miscellaneous Routines used by the APPHOT Package
+
+$checkout libpkg.a ".."
+$update libpkg.a
+$checkin libpkg.a ".."
+$exit
+
+libpkg.a:
+	apbsmooth.x
+	apclip.x
+	apdate.x		<time.h>
+	apdiverr.x
+	apfnames.x
+	apgscur.x		<gset.h>		<fset.h>
+	apgtools.x		<pkg/gtools.h>		<gset.h>
+	aplucy.x
+	apnlfuncs.x		<mach.h>		<math.h>
+	apmapr.x
+	apmeds.x
+	apmoments.x
+	appcache.x		<imhdr.h>		<imset.h>
+	aprmwhite.x		<ctype.h>
+	aptopt.x		<math.h> <mach.h>
+	apvectors.x
+	;