Initial commit

34 files changed, 5468 insertions, 0 deletions

diff --git a/noao/digiphot/apphot/fitsky/apavsky.x b/noao/digiphot/apphot/fitsky/apavsky.x
new file mode 100644
index 00000000..c7e38c0e
--- /dev/null
+++ b/noao/digiphot/apphot/fitsky/apavsky.x
@@ -0,0 +1,107 @@
+include "../lib/display.h"
+include "../lib/fitsky.h"
+
+# AP_AVSKY -- Compute an estimate of the sky value by averaging several
+# individual sky values measured in different parts of the frame.
+
+int procedure ap_avsky (ap, im, stid, sd, id, gd, interactive)
+
+pointer	ap			# the pointer to the main apphot data structure
+pointer	im			# the pointer to the input image
+int	stid			# the current sequence number
+int	sd			# the sky file descriptor
+pointer	id			# the display stream descriptor
+pointer	gd			# the graphics stream descriptor
+int	interactive		# interactive mode
+
+int	wcs, key, nmsky, nmsigma, nmskew, nsky, nrej, sier, ier
+pointer	sp, cmd
+real	wx, wy, msky, msigma, mskew
+int	clgcur(), apfitsky(), apstati()
+real	apstatr()
+
+begin
+	call smark (sp)
+	call salloc (cmd, SZ_LINE, TY_CHAR)
+
+	# Initialize the skyvalues.
+	msky = 0.0
+	nmsky = 0
+	msigma = 0.0
+	nmsigma = 0
+	mskew = 0.0
+	nmskew = 0
+	nsky = 0
+	nrej = 0
+	ier = AP_OK
+
+	call printf (
+	"\nMeasure sky around several cursor positions (t=measure, q=quit)\n")
+	while (clgcur ("icommands", wx, wy, wcs, key, Memc[cmd], SZ_LINE) !=
+	    EOF) {
+
+	    switch (key) {
+	    case 'q':
+		break
+
+	    case 't':
+		sier = apfitsky (ap, im, wx, wy, sd, gd)
+		if (sier != AP_OK)
+		    ier = sier
+
+		if (id != NULL) {
+		    call apmark (ap, id, NO, apstati (ap, MKSKY), NO)
+		    if (id == gd)
+			call gflush (id)
+		    else
+		        call gframe (id)
+		}
+		call ap_splot (ap, stid, gd, apstati (ap, RADPLOTS))
+		if (interactive == YES)
+		    call ap_qspsky (ap, sier)
+
+		if (! IS_INDEFR(apstatr (ap, SKY_MODE))) {
+		    msky = msky + apstatr (ap, SKY_MODE)
+		    nmsky = nmsky + 1
+		}
+		if (! IS_INDEFR(apstatr (ap, SKY_SIGMA))) {
+		    msigma = msigma + apstatr (ap, SKY_SIGMA)
+		    nmsigma = nmsigma + 1
+		}
+		if (! IS_INDEFR(apstatr (ap, SKY_SKEW))) {
+		    mskew = mskew + apstatr (ap, SKY_SKEW)
+		    nmskew = nmskew + 1
+		}
+		nsky = nsky + apstati (ap, NSKY)
+		nrej = nrej + apstati (ap, NSKY_REJECT)
+
+	    default:
+		;
+	    }
+	}
+
+	# Compute the average values.
+	if (nmsky > 0)
+	    msky = msky / nmsky
+	else
+	    msky = INDEFR
+	if (nmsigma > 0)
+	    msigma = msigma / nmsigma
+	else
+	    msigma = INDEFR
+	if (nmskew > 0)
+	    mskew = mskew / nmskew
+	else
+	    mskew = INDEFR
+
+	# Store the average values.
+	call apsetr (ap, SKY_MODE, msky)
+	call apsetr (ap, SKY_SIGMA, msigma)
+	call apsetr (ap, SKY_SKEW, mskew)
+	call apseti (ap, NSKY, nsky)
+	call apseti (ap, NSKY_REJECT, nrej)
+
+	call sfree (sp)
+
+	return (ier)
+end
diff --git a/noao/digiphot/apphot/fitsky/apbsky.x b/noao/digiphot/apphot/fitsky/apbsky.x
new file mode 100644
index 00000000..28e96288
--- /dev/null
+++ b/noao/digiphot/apphot/fitsky/apbsky.x
@@ -0,0 +1,98 @@
+include <fset.h>
+include "../lib/apphot.h"
+include "../lib/display.h"
+
+# APBSKY -- Procedure to determine the sky statistics for a list of objects
+# in batch mode using a simple coordinate list.
+
+procedure apbsky (ap, im, cl, sd, out, id, ld, gd, mgd, gid, interactive)
+
+pointer ap			# pointer to apphot structure
+pointer	im			# pointer to IRAF image
+int	cl			# starlist file descriptor
+int	sd			# sky file descriptor
+int	out			# output file descriptor
+int	id, ld			# sequence and list numbers
+int	gd			# pointer to stdgraph stream
+pointer	mgd			# pointer to graphics metacode file
+pointer	gid			# pointer to image display stream
+int	interactive		# interactive mode
+
+int	stdin, ier, ild
+pointer	sp, str
+real	wx, wy
+int	apfitsky(), 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 the 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) {
+
+	    # Fetch 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 sky value.
+	    ier = apfitsky (ap, im, wx, wy, sd, gd)
+
+	    # Print the sky values.
+	    if (interactive == YES) {
+		call ap_qspsky (ap, ier)
+		if (gid != NULL)
+		    call apmark (ap, gid, NO, apstati (ap, MKSKY), NO)
+	    }
+	    call ap_splot(ap, id, mgd, YES)
+	    if (id == 1)
+	        call ap_param (ap, out, "fitsky")
+	    call ap_pssky (ap, out, id, ild, ier)
+
+	    # Set up for the next object.
+	    id = id + 1
+	    ild = ild + 1
+	    call apsetr (ap, WX, wx)
+	    call apsetr (ap, WY, wy)
+
+	    # print query
+	    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/fitsky/apcentroid.x b/noao/digiphot/apphot/fitsky/apcentroid.x
new file mode 100644
index 00000000..310db85d
--- /dev/null
+++ b/noao/digiphot/apphot/fitsky/apcentroid.x
@@ -0,0 +1,244 @@
+include <mach.h>
+include "../lib/fitsky.h"
+
+# AP_CENTROID -- Procedure to find the mode, width and skew of the sky
+# distribution by computing the moments of the histogram. 
+
+int procedure ap_centroid (skypix, coords, wgt, index, nskypix, snx, sny, k1,
+	hwidth, binsize, smooth, losigma, hisigma, rgrow, maxiter, sky_mode,
+	sky_sigma, sky_skew, nsky, nsky_reject)
+
+real	skypix[ARB]		# array of sky pixels
+int	coords[ARB]		# array of coordinates for regions growing
+real	wgt[ARB]		# array of weights for rejection
+int	index[ARB]		# array of sorted indices
+int	nskypix			# the number of sky pixels
+int	snx, sny		# the maximum dimensions of sky raster
+real	k1			# extent of the histogram in skysigma
+real	hwidth			# width of the histogram
+real	binsize			# the size of the histogram in sky sigma
+int	smooth			# smooth the histogram before fitting
+real	losigma, hisigma	# upper and lower k-sigma rejection limits
+real	rgrow			# region growing radius in pixels
+int	maxiter			# maximum number of rejection cycles
+real	sky_mode		# computed sky value
+real	sky_sigma		# computed sigma of the sky pixels
+real	sky_skew		# skew of sky pixels
+int	nsky			# number of sky pixels used in fit
+int	nsky_reject		# number of sky pixels rejected
+
+double	dsky, sumpx, sumsqpx, sumcbpx
+int	nreject, nbins, nker, ier, i, j
+pointer	sp, hgm, shgm
+real	dmin, dmax, sky_mean, hmin, hmax, dh, locut, hicut, cut
+real	sky_zero
+int	ap_grow_hist2(), aphigmr()
+real	ap_asumr(), apmedr()
+
+begin
+	# Intialize.
+	nsky = nskypix
+	nsky_reject = 0
+	sky_mode = INDEFR
+	sky_sigma = INDEFR
+	sky_skew = INDEFR
+	if (nskypix <= 0)
+	    return (AP_NOSKYAREA)
+
+	# Compute the histogram width and binsize.
+	sky_zero = ap_asumr (skypix, index, nskypix) / nskypix
+	call ap_ialimr (skypix, index, nskypix, dmin, dmax)
+	call apfimoments (skypix, index, nskypix, sky_zero, sumpx, sumsqpx,
+	    sumcbpx, sky_mean, sky_sigma, sky_skew)
+	sky_mean = apmedr (skypix, index, nskypix)
+	sky_mean = max (dmin, min (sky_mean, dmax))
+	if (! IS_INDEFR(hwidth) && hwidth > 0.0) {
+	    hmin = sky_mean - k1 * hwidth
+	    hmax = sky_mean + k1 * hwidth
+	    dh = binsize * hwidth
+	} else {
+	    cut = min (sky_mean - dmin, dmax - sky_mean, k1 * sky_sigma)
+	    hmin = sky_mean - cut
+	    hmax = sky_mean + cut
+	    dh = binsize * cut / k1
+	}
+
+	# Compute the number of histogram bins and the histogram resolution.
+	if (dh <= 0.0) {
+	    nbins = 1
+	    dh = 0.0
+	} else {
+	    nbins = 2 * nint ((hmax - sky_mean) / dh) + 1
+	    dh = (hmax - hmin) / (nbins - 1)
+	}
+
+	# Check for a valid histogram.
+	if (nbins < 2 || k1 <= 0.0 || sky_sigma < dh || dh <= 0.0 ||
+	    sky_sigma <= 0.0) {
+	    sky_mode = sky_mean
+	    sky_sigma = 0.0
+	    sky_skew = 0.0
+	    return (AP_NOHISTOGRAM)
+	}
+
+	# Allocate temporary space.
+	call smark (sp)
+	call salloc (hgm, nbins, TY_REAL)
+	call salloc (shgm, nbins, TY_REAL)
+
+	# Accumulate the histogram.
+	call aclrr (Memr[hgm], nbins)
+	nsky_reject = nsky_reject + aphigmr (skypix, wgt, index, nskypix,
+	    Memr[hgm], nbins, hmin, hmax)
+	nsky = nskypix - nsky_reject
+
+	# Do the initial rejection.
+	if (nsky_reject > 0) {
+	    do i = 1, nskypix {
+		if (wgt[index[i]] <= 0.0) {
+		    dsky = skypix[index[i]] - sky_zero
+		    sumpx = sumpx - dsky
+		    sumsqpx = sumsqpx - dsky ** 2
+		    sumcbpx = sumcbpx - dsky ** 3
+		}
+	    }
+	    call apmoments (sumpx, sumsqpx, sumcbpx, nsky, sky_zero, sky_mean,
+	        sky_sigma, sky_skew)
+	}
+
+	# Fit the mode, sigma an skew of the histogram.
+	if (smooth == YES) {
+	    nker = max (1, nint (sky_sigma / dh))
+	    #call ap_lucy_smooth (Memr[hgm], Memr[shgm], nbins, nker, 2)
+	    call ap_bsmooth (Memr[hgm], Memr[shgm], nbins, nker, 2)
+	    call ap_imode (Memr[shgm], nbins, hmin, hmax, YES, sky_mode, ier)
+	} else
+	    call ap_imode (Memr[hgm], nbins, hmin, hmax, NO, sky_mode, ier)
+	sky_mode = max (dmin, min (sky_mode, dmax))
+	if ((IS_INDEFR(losigma) && IS_INDEFR(hisigma)) || (sky_sigma < dh) ||
+	    maxiter < 1) {
+	    call sfree (sp)
+	    return (AP_OK)
+	}
+
+	# Fit histogram with pixel rejection and optional region growing.
+	do i = 1, maxiter {
+
+	    # Compute new histogram limits.
+	    if (IS_INDEFR(losigma))
+		locut = -MAX_REAL
+	    else
+	        locut = sky_mode - losigma * sky_sigma
+	    if (IS_INDEFR(hisigma))
+		hicut = MAX_REAL
+	    else
+	        hicut = sky_mode + hisigma * sky_sigma
+
+	    # Reject pixels.
+	    nreject = 0
+	    do j = 1, nskypix {
+		if (skypix[index[j]] >= locut && skypix[index[j]] <= hicut)
+		    next
+		if (rgrow > 0.0)
+		    nreject = nreject + ap_grow_hist2 (skypix, coords,
+			wgt, nskypix, sky_zero, index[j], snx, sny, Memr[hgm],
+			nbins, hmin, hmax, rgrow, sumpx, sumsqpx, sumcbpx)
+		else if (wgt[index[j]] > 0.0) {
+		    call ap_hgmsub2 (Memr[hgm], nbins, hmin, hmax,
+		        skypix[index[j]], sky_zero, sumpx, sumsqpx, sumcbpx)
+		    wgt[index[j]] = 0.0
+		    nreject = nreject + 1
+		}
+	    }
+	    if (nreject == 0)
+		break
+	    nsky_reject = nsky_reject + nreject
+	    nsky = nskypix - nsky_reject 
+	    if (nsky <= 0)
+		break
+
+	    # Recompute the moments.
+	    call apmoments (sumpx, sumsqpx, sumcbpx, nsky, sky_zero, sky_mean,
+	        sky_sigma, sky_skew)
+
+	    # Recompute the histogram peak.
+	    if (smooth == YES) {
+		nker = max (1, nint (sky_sigma / dh))
+	        #call ap_lucy_smooth (Memr[hgm], Memr[shgm], nbins, nker, 2)
+	        call ap_bsmooth (Memr[hgm], Memr[shgm], nbins, nker, 2)
+	        call ap_imode (Memr[shgm], nbins, hmin, hmax, YES, sky_mode,
+		    ier)
+	    } else
+	        call ap_imode (Memr[hgm], nbins, hmin, hmax, NO, sky_mode, ier)
+	    sky_mode = max (dmin, min (sky_mode, dmax))
+	    if (sky_sigma <= dh || ier != AP_OK)
+		break
+	}
+
+	# Return the error codes.
+	call sfree (sp)
+	if (nsky == 0 || nsky_reject == nskypix || ier == AP_NOHISTOGRAM) {
+	    nsky = 0
+	    nsky_reject = nskypix
+	    sky_mode = INDEFR
+	    sky_sigma = INDEFR
+	    sky_skew = INDEFR
+	    return (AP_NSKY_TOO_SMALL)
+	} else
+	    return (AP_OK)
+end
+
+
+# AP_IMODE -- Procedure to compute the 1st, 2nd and third moments of the
+# histogram of sky pixels.
+
+procedure ap_imode (hgm, nbins, z1, z2, smooth, sky_mode, ier)
+
+real	hgm[ARB]		# histogram
+int	nbins			# number of bins
+real	z1, z2			# min and max of the histogram
+int	smooth			# is the histogram smoothed
+real	sky_mode		# mode of histogram
+int	ier			# error code
+
+int	i, noccup
+double	sumi, sumix, x, dh
+real	hmean, dz
+real	asumr()
+
+begin
+	# Initialize the sums.
+	sumi = 0.0
+	sumix = 0.0
+
+	# Compute a continuum level.
+	if (smooth == NO) 
+	    hmean = asumr (hgm, nbins) / nbins
+	else {
+	    call alimr (hgm, nbins, dz, hmean)
+	    hmean = 2.0 * hmean / 3.0
+	}
+
+	# Accumulate the sums.
+	noccup = 1
+	dz = (z2 - z1) / (nbins - 1)
+	x = z1 + 0.5 * dz
+	do i = 1, nbins {
+	    dh = hgm[i] - hmean
+	    if (dh > 0.0d0) {
+	        sumi = sumi + dh
+	        sumix = sumix + dh * x
+	        noccup = noccup + 1
+	    }
+	    x = x + dz
+	}
+
+	# Compute the sky mode, sigma and skew.
+	if (sumi > 0.0) {
+	    sky_mode = sumix / sumi
+	    ier = AP_OK
+	} else {
+	    sky_mode = INDEFR
+	    ier = AP_NOHISTOGRAM
+	}
+end
diff --git a/noao/digiphot/apphot/fitsky/apcrosscor.x b/noao/digiphot/apphot/fitsky/apcrosscor.x
new file mode 100644
index 00000000..757addcc
--- /dev/null
+++ b/noao/digiphot/apphot/fitsky/apcrosscor.x
@@ -0,0 +1,288 @@
+include <mach.h>
+include "../lib/fitsky.h"
+
+# AP_CROSSCOR -- Procedure to compute the sky value by calculating the
+# cross-correlation function of the histogram of the sky pixels and
+# a Gaussian function with the same sigma as the sky distribution.
+# The peak of the cross-correlation function is found by parabolic
+# interpolation.
+
+int procedure ap_crosscor (skypix, coords, wgt, index, nskypix, snx, sny, k1,
+	hwidth, binsize, smooth, losigma, hisigma, rgrow, maxiter, sky_mode,
+	sky_sigma, sky_skew, nsky, nsky_reject)
+
+real	skypix[ARB]		# array of sky pixels
+int	coords[ARB]		# array of sky coordinates for region growing
+real	wgt[ARB]		# array of weights for rejection
+int	index[ARB]		# array of sort indices
+int	nskypix			# the number of sky pixels
+int	snx, sny		# the maximum dimensions of sky raster
+real	k1			# half-width of the histogram in sky sigma
+real	hwidth			# the input sky sigma
+real	binsize			# the size of the histogram in sky sigma
+int	smooth			# smooth the histogram before fitting (not used)
+real	losigma, hisigma	# upper and lower rejection limits
+real	rgrow			# region growing radius in pixels
+int	maxiter			# maximum number of rejection cycles
+real	sky_mode		# computed sky value
+real	sky_sigma		# computed standard deviation of the sky pixels
+real	sky_skew		# computed skew of sky pixels
+int	nsky			# number of sky pixels used in fit
+int	nsky_reject		# number of sky pixels rejected
+
+double 	dsky, sumpx, sumsqpx, sumcbpx
+int	nreject, nbins, nker, nsmooth, ier, i, j
+pointer	sp, x, hgm, shgm, kernel
+real	dmin, dmax, hmin, hmax, dh, kmin, kmax, locut, hicut, sky_mean, cut
+real	sky_zero
+int	ap_grow_hist2(), aphigmr()
+real	ap_asumr(), apmedr()
+
+begin
+	# Initialize.
+	nsky = nskypix
+	nsky_reject = 0
+	sky_mode = INDEFR
+	sky_sigma = INDEFR
+	sky_skew = INDEFR
+	if (nskypix <= 0)
+	    return (AP_NOSKYAREA)
+
+	# Set up initial guess at sky mean, sigma and skew.
+	sky_zero = ap_asumr (skypix, index, nskypix) / nskypix
+	call ap_ialimr (skypix, index, nskypix, dmin, dmax)
+	call apfimoments (skypix, index, nskypix, sky_zero,
+	    sumpx, sumsqpx, sumcbpx, sky_mean, sky_sigma, sky_skew)
+	sky_mean = apmedr (skypix, index, nskypix)
+	sky_mean = max (dmin, min (sky_mean, dmax))
+
+	# Compute the width and bin size of the histogram.
+	if (! IS_INDEFR(hwidth) && hwidth > 0.0) {
+	    hmin = sky_mean - k1 * hwidth
+	    hmax = sky_mean + k1 * hwidth
+	    dh = binsize * hwidth
+	} else {
+	    cut = min (sky_mean - dmin, dmax - sky_mean, k1 * sky_sigma)
+	    hmin = sky_mean - cut
+	    hmax = sky_mean + cut
+	    dh = binsize * cut / k1
+	} 
+
+	# Compute the number of bins in and the width of the kernel.
+	if (dh <= 0.0) {
+	    nbins = 1
+	    nker = 1
+	    nsmooth = 1
+	    dh = 0.0
+	} else {
+    	    nbins = 2 * nint ((hmax - sky_mean) / dh) + 1
+	    nker = 2 * nint (2.0 * (hmax - sky_mean) / (dh * 3.0)) + 1
+	    nsmooth = nbins - nker + 1
+	    dh = (hmax - hmin) / (nbins - 1)
+	}
+	kmin = - dh * (nker / 2 + 0.5)
+	kmax = dh * (nker / 2 + 0.5)
+
+	# Test for a valid histogram.
+	if (nbins < 2 || k1 <= 0.0 || sky_sigma <= 0.0 || dh <= 0.0 ||
+	    sky_sigma <= dh) {
+	    sky_mode = sky_mean
+	    sky_sigma = 0.0
+	    sky_skew = 0.0
+	    return (AP_NOHISTOGRAM)
+	}
+
+	# Allocate space for the histogram and kernel.
+	call smark (sp)
+	call salloc (x, nbins, TY_REAL)
+	call salloc (hgm, nbins, TY_REAL)
+	call salloc (shgm, nbins, TY_REAL)
+	call salloc (kernel, nker, TY_REAL)
+
+	# Set up x array.
+	do i = 1, nbins
+	    Memr[x+i-1] = i
+	call amapr (Memr[x], Memr[x], nbins, 1.0, real (nbins),
+	    hmin + 0.5 * dh, hmax + 0.5 * dh)
+
+	# Accumulate the histogram.
+	call aclrr (Memr[hgm], nbins)
+	call aclrr (Memr[shgm], nbins)
+	nsky_reject = nsky_reject + aphigmr (skypix, wgt, index, nskypix,
+	    Memr[hgm], nbins, hmin, hmax)
+	nsky = nskypix - nsky_reject
+
+	# Perform the initial rejection cycle.
+	if (nsky_reject > 0.0) {
+	    do i = 1, nskypix {
+	        if (wgt[index[i]] <= 0.0) {
+		    dsky = skypix[index[i]] - sky_zero
+		    sumpx = sumpx - dsky
+		    sumsqpx = sumsqpx - dsky ** 2
+		    sumcbpx = sumcbpx - dsky ** 3
+		}
+	    }
+	    call apmoments (sumpx, sumsqpx, sumcbpx, nsky, sky_zero, sky_mean,
+	        sky_sigma, sky_skew)
+	}
+
+	# Construct kernel and convolve with histogram.
+	if (sky_sigma > 0.0) {
+	    call ap_gauss_kernel (Memr[kernel], nker, kmin, kmax, sky_sigma)
+	    call acnvr (Memr[hgm], Memr[shgm+nker/2], nsmooth, Memr[kernel],
+	        nker)
+	} else
+	    call amovr (Memr[hgm], Memr[shgm], nbins)
+	call ap_corfit (Memr[x], Memr[shgm], nbins, sky_mode, ier)
+	sky_mode = max (dmin, min (sky_mode, dmax))
+	if (ier != OK) {
+	    call sfree (sp)
+	    return (ier)
+	}
+	if ((IS_INDEFR(losigma) && IS_INDEFR(hisigma)) || (sky_sigma <= dh) ||
+	    maxiter < 1) {
+	    call sfree (sp)
+	    return (AP_OK)
+	}
+
+	# Fit histogram with pixel rejection and optional region growing.
+	do i = 1, maxiter {
+
+	    # Compute new histogram limits.
+	    if (IS_INDEFR(losigma))
+		locut = -MAX_REAL
+	    else
+	        locut = sky_mode - losigma * sky_sigma
+	    if (IS_INDEFR(hisigma))
+		hicut = MAX_REAL
+	    else
+	        hicut = sky_mode + hisigma * sky_sigma
+
+	    # Detect rejected pixels.
+	    nreject = 0
+	    do j = 1, nskypix {
+		if (skypix[index[j]] >= locut && skypix[index[j]] <= hicut)
+		    next
+		if (rgrow > 0.0)
+		    nreject = nreject + ap_grow_hist2 (skypix, coords,
+			wgt, nskypix, sky_zero, index[j], snx, sny,
+			Memr[hgm], nbins, hmin, hmax, rgrow, sumpx, sumsqpx,
+			sumcbpx)
+		else if (wgt[index[j]] > 0.0) {
+		    call ap_hgmsub2 (Memr[hgm], nbins, hmin, hmax,
+			skypix[index[j]], sky_zero, sumpx, sumsqpx, sumcbpx)
+		    wgt[index[j]] = 0.0
+		    nreject = nreject + 1
+		}
+	    }
+	    if (nreject == 0)
+		break
+
+	    # Update the sky parameters.
+	    nsky_reject = nsky_reject + nreject
+	    nsky = nskypix - nsky_reject 
+	    if (nsky <= 0)
+		break
+	    call apmoments (sumpx, sumsqpx, sumcbpx, nsky, sky_zero, sky_mean,
+	        sky_sigma, sky_skew)
+	    if (sky_sigma <= dh)
+		break
+
+	    # Recompute the peak of the histogram.
+	    call ap_gauss_kernel (Memr[kernel], nker, kmin, kmax, sky_sigma)
+	    call aclrr (Memr[shgm], nbins)
+	    call acnvr (Memr[hgm], Memr[shgm+nker/2], nsmooth, Memr[kernel],
+	        nker)
+	    call ap_corfit (Memr[x], Memr[shgm], nbins, sky_mode, ier)
+	    sky_mode = max (dmin, min (sky_mode, dmax))
+	    if (ier != AP_OK)
+		 break
+	}
+
+	# Return the appropriate error code.
+	call sfree (sp)
+	if (nsky == 0 || nsky_reject == nskypix) {
+	    nsky = 0
+	    nsky_reject = nskypix
+	    sky_mode = INDEFR
+	    sky_sigma = INDEFR
+	    sky_skew = INDEFR
+	    return (AP_NSKY_TOO_SMALL)
+	} else if (ier != AP_OK) {
+	    sky_mode = sky_mean
+	    sky_sigma = 0.0
+	    sky_skew = 0.0 
+	    return (ier)
+	} else
+	    return (AP_OK)
+end
+
+
+# AP_GAUSS_KERNEL -- Procedure to compute a Gaussian kernel of given length
+# and sigma.
+
+procedure ap_gauss_kernel (kernel, nker, kmin, kmax, sky_sigma)
+
+real	kernel[ARB]		# kernel 
+int	nker			# length of kernel
+real	kmin, kmax		# limits of the kernel
+real	sky_sigma		# sigma of the sky
+
+int	i
+real	dk, x, sumx
+
+begin
+	# Return 1 if unit sized kernel.
+	if (nker == 1) {
+	    kernel[1] = 1.0
+	    return
+	}
+
+	# Intialize.
+	sumx = 0.0
+	x = kmin
+	dk = (kmax - kmin ) / (nker - 1)
+
+ 	# Compute and normalize the kernel.
+	do i = 1, nker {
+	    kernel[i] = exp (- (x ** 2) / (2. * sky_sigma ** 2))
+	    sumx = sumx + kernel[i]
+	    x = x + dk
+	}
+        do i = 1, nker
+            kernel[i] = kernel[i] / sumx
+end
+
+
+# AP_CORFIT -- Procedure to compute the peak of the cross-correlation
+# function using parabolic interpolation.
+
+procedure ap_corfit (x, shgm, nbins, sky_mode, ier)
+
+real	x[ARB]				# x coordinates of histogram
+real	shgm[ARB]			# convolved histogram
+int	nbins				# number of histogram bins
+real	sky_mode			# computed sky_mode
+int	ier				# error code
+
+int	bin
+real	max, xo, dh1, dh2
+
+begin
+	call ap_amaxel (shgm, nbins, max, bin)
+	if (max <= 0) {
+	    ier = AP_FLAT_HIST
+	} else if (bin == 1) {
+	    sky_mode = x[1]
+	    ier = AP_OK
+	} else if (bin == nbins) {
+	    sky_mode = x[nbins]
+	    ier = AP_OK
+	} else {
+	    xo = 0.5 * (x[bin] + x[bin-1])
+	    dh1 = shgm[bin] - shgm[bin-1] 
+	    dh2 = shgm[bin] - shgm[bin+1]
+	    sky_mode = xo + (x[bin] - x[bin-1]) * (dh1 / (dh1 + dh2))
+	    ier = AP_OK
+	}
+end
diff --git a/noao/digiphot/apphot/fitsky/apfitsky.x b/noao/digiphot/apphot/fitsky/apfitsky.x
new file mode 100644
index 00000000..1c105b79
--- /dev/null
+++ b/noao/digiphot/apphot/fitsky/apfitsky.x
@@ -0,0 +1,452 @@
+include "../lib/apphotdef.h"
+include "../lib/apphot.h"
+include "../lib/noisedef.h"
+include "../lib/fitskydef.h"
+include "../lib/fitsky.h"
+
+# APFITSKY -- Procedure to the compute the sky value in an annular region
+# around a given position in the IRAF image.
+
+int procedure apfitsky (ap, im, wx, wy, sd, gd)
+
+pointer	ap		# pointer to the apphot structure
+pointer	im		# pointer to the IRAF image
+real	wx		# object x coordinate
+real	wy		# object y coordinate
+int	sd		# pointer to input text file containing sky values
+pointer	gd		# pointer to graphics stream
+
+int	ier, nclip, nsky, ilo, ihi
+pointer	sky, nse, gt
+real	x, y
+int	apskybuf(), ap_mode(), ap_centroid(), ap_histplot(), ap_readsky()
+int	ap_median(), ap_radplot(), ap_gauss(), ap_lgsky(), ap_crosscor()
+int	ap_mean(), ap_clip()
+pointer	ap_gtinit()
+
+begin
+	# Initialize.
+	sky = AP_PSKY(ap)
+	nse = AP_NOISE(ap)
+	AP_SXCUR(sky) = wx
+	AP_SYCUR(sky) = wy
+        if (IS_INDEFR(wx) || IS_INDEFR(wy)) {
+            AP_OSXCUR(sky) = INDEFR
+            AP_OSYCUR(sky) = INDEFR
+        } else {
+            switch (AP_WCSOUT(ap)) {
+            case WCS_WORLD, WCS_PHYSICAL:
+                call ap_ltoo (ap, wx, wy, AP_OSXCUR(sky), AP_OSYCUR(sky), 1)
+            case WCS_TV:
+                call ap_ltov (im, wx, wy, AP_OSXCUR(sky), AP_OSYCUR(sky), 1)
+            default:
+                AP_OSXCUR(sky) = wx
+                AP_OSYCUR(sky) = wy
+            }
+        }
+
+	AP_SKY_MODE(sky) = INDEFR
+	AP_SKY_SIG(sky) = INDEFR
+	AP_SKY_SKEW(sky) = INDEFR
+	AP_NSKY(sky) = 0
+	AP_NSKY_REJECT(sky) = 0
+	if (IS_INDEFR(wx) || IS_INDEFR(wy))
+	    return (AP_NOSKYAREA)
+
+	switch (AP_SKYFUNCTION(sky)) {
+
+	case AP_MEAN:
+
+	    # Fetch  the sky pixels.
+	    ier = apskybuf (ap, im, wx, wy)
+	    if (ier != AP_OK)
+		return (ier)
+
+	    # Initialze the weights.
+	    call amovkr (1.0, Memr[AP_SWGT(sky)], AP_NSKYPIX(sky))
+
+	    # Clip the data.
+	    if (AP_SLOCLIP(sky) > 0.0 || AP_SHICLIP(sky) > 0.0) {
+		nclip = ap_clip (Memr[AP_SKYPIX(sky)], Memi[AP_INDEX(sky)],
+		    AP_NSKYPIX(sky), AP_SLOCLIP(sky), AP_SHICLIP(sky), ilo,
+		    ihi)
+		if (nclip >= AP_NSKYPIX(sky))
+		    return (AP_NSKY_TOO_SMALL)
+		nsky = AP_NSKYPIX(sky) - nclip
+	    } else {
+		nclip = 0
+		call ap_index (Memi[AP_INDEX(sky)], AP_NSKYPIX(sky))
+		ilo = 1
+		nsky = AP_NSKYPIX(sky) 
+	    }
+
+	    # Compute the mean of the sky pixel distribution with pixel
+	    # rejection and region growing.
+	    ier = ap_mean (Memr[AP_SKYPIX(sky)], Memi[AP_COORDS(sky)],
+	        Memr[AP_SWGT(sky)], Memi[AP_INDEX(sky)+ilo-1],
+		nsky, AP_SNX(sky), AP_SNY(sky), AP_SLOREJECT(sky),
+		AP_SHIREJECT(sky), AP_RGROW(sky) * AP_SCALE(ap),
+		AP_SNREJECT(sky), AP_SKY_MODE(sky), AP_SKY_SIG(sky),
+		AP_SKY_SKEW(sky), AP_NSKY(sky), AP_NSKY_REJECT(sky))
+		AP_NSKY_REJECT(sky) = AP_NBADSKYPIX(sky) + nclip +
+		    AP_NSKY_REJECT(sky)
+
+	    return (ier)
+
+	case AP_MEDIAN:
+
+	    # Fetch  the sky pixels.
+	    ier = apskybuf (ap, im, wx, wy)
+	    if (ier != AP_OK)
+		return (ier)
+
+	    # Initialze the weights.
+	    call amovkr (1.0, Memr[AP_SWGT(sky)], AP_NSKYPIX(sky))
+
+	    # Clip the data.
+	    if (AP_SLOCLIP(sky) > 0.0 || AP_SHICLIP(sky) > 0.0) {
+		nclip = ap_clip (Memr[AP_SKYPIX(sky)], Memi[AP_INDEX(sky)],
+		    AP_NSKYPIX(sky), AP_SLOCLIP(sky), AP_SHICLIP(sky), ilo,
+		    ihi)
+		if (nclip >= AP_NSKYPIX(sky))
+		    return (AP_NSKY_TOO_SMALL)
+		nsky = AP_NSKYPIX(sky) - nclip
+	    } else {
+		nclip = 0
+		call apqsort (Memr[AP_SKYPIX(sky)], Memi[AP_INDEX(sky)],
+		    Memi[AP_INDEX(sky)], AP_NSKYPIX(sky))
+		ilo = 1
+		nsky = AP_NSKYPIX(sky) 
+	    }
+
+	    # Compute the median of the sky pixel distribution with pixel
+	    # rejection and region growing.
+	    ier = ap_median (Memr[AP_SKYPIX(sky)], Memi[AP_COORDS(sky)],
+	        Memr[AP_SWGT(sky)], Memi[AP_INDEX(sky)+ilo-1], nsky,
+		AP_SNX(sky), AP_SNY(sky), AP_SLOREJECT(sky),
+		AP_SHIREJECT(sky), AP_RGROW(sky) * AP_SCALE(ap),
+		AP_SNREJECT(sky), AP_SKY_MODE(sky), AP_SKY_SIG(sky),
+		AP_SKY_SKEW(sky), AP_NSKY(sky), AP_NSKY_REJECT(sky))
+	    AP_NSKY_REJECT(sky) = AP_NBADSKYPIX(sky) + nclip +
+	        AP_NSKY_REJECT(sky)
+
+	    return (ier)
+
+	case AP_MODE:
+
+	    # Fetch the sky pixels.
+	    ier = apskybuf (ap, im, wx, wy)
+	    if (ier != AP_OK)
+		return (ier)
+
+	    # Initialze the weights.
+	    call amovkr (1.0, Memr[AP_SWGT(sky)], AP_NSKYPIX(sky))
+
+	    # Clip the data.
+	    if (AP_SLOCLIP(sky) > 0.0 || AP_SHICLIP(sky) > 0.0) {
+		nclip = ap_clip (Memr[AP_SKYPIX(sky)], Memi[AP_INDEX(sky)],
+		    AP_NSKYPIX(sky), AP_SLOCLIP(sky), AP_SHICLIP(sky), ilo,
+		    ihi)
+		if (nclip >= AP_NSKYPIX(sky))
+		    return (AP_NSKY_TOO_SMALL)
+		nsky = AP_NSKYPIX(sky) - nclip
+	    } else {
+		nclip = 0
+		call apqsort (Memr[AP_SKYPIX(sky)], Memi[AP_INDEX(sky)],
+		    Memi[AP_INDEX(sky)], AP_NSKYPIX(sky))
+		ilo = 1
+		nsky = AP_NSKYPIX(sky) 
+	    }
+
+	    # Compute the median of the sky pixel distribution with pixel
+	    # rejection and region growing.
+	    ier = ap_mode (Memr[AP_SKYPIX(sky)], Memi[AP_COORDS(sky)],
+	        Memr[AP_SWGT(sky)], Memi[AP_INDEX(sky)+ilo-1], nsky,
+		AP_SNX(sky), AP_SNY(sky), AP_SLOREJECT(sky),
+		AP_SHIREJECT(sky), AP_RGROW(sky) * AP_SCALE(ap),
+		AP_SNREJECT(sky), AP_SKY_MODE(sky), AP_SKY_SIG(sky),
+		AP_SKY_SKEW(sky), AP_NSKY(sky), AP_NSKY_REJECT(sky))
+	    AP_NSKY_REJECT(sky) = AP_NBADSKYPIX(sky) + nclip +
+	        AP_NSKY_REJECT(sky)
+
+	    return (ier)
+
+	case AP_CENTROID:
+
+	    # Fetch the sky pixels.
+	    ier = apskybuf (ap, im, wx, wy)
+	    if (ier != AP_OK)
+		return (ier)
+
+	    # Initialze the weights.
+	    call amovkr (1.0, Memr[AP_SWGT(sky)], AP_NSKYPIX(sky))
+
+	    # Clip the data.
+	    if (AP_SLOCLIP(sky) > 0.0 || AP_SHICLIP(sky) > 0.0) {
+		nclip = ap_clip (Memr[AP_SKYPIX(sky)], Memi[AP_INDEX(sky)],
+		    AP_NSKYPIX(sky), AP_SLOCLIP(sky), AP_SHICLIP(sky), ilo,
+		    ihi)
+		if (nclip >= AP_NSKYPIX(sky))
+		    return (AP_NSKY_TOO_SMALL)
+		nsky = AP_NSKYPIX(sky) - nclip
+	    } else {
+		nclip = 0
+		call ap_index (Memi[AP_INDEX(sky)], AP_NSKYPIX(sky))
+		ilo = 1
+		nsky = AP_NSKYPIX(sky) 
+	    }
+
+	    # Compute the sky value by performing a moment analysis of the
+	    # sky pixel histogram.
+	    ier = ap_centroid (Memr[AP_SKYPIX(sky)], Memi[AP_COORDS(sky)],
+	        Memr[AP_SWGT(sky)], Memi[AP_INDEX(sky)+ilo-1], nsky,
+		AP_SNX(sky), AP_SNY(sky), AP_K1(sky), INDEFR,
+		AP_BINSIZE(sky), AP_SMOOTH(sky), AP_SLOREJECT(sky),
+		AP_SHIREJECT(sky), AP_RGROW(sky) * AP_SCALE(ap),
+		AP_SMAXITER(sky), AP_SKY_MODE(sky), AP_SKY_SIG(sky),
+		AP_SKY_SKEW(sky), AP_NSKY(sky), AP_NSKY_REJECT(sky))
+	    AP_NSKY_REJECT(sky) = AP_NBADSKYPIX(sky) + nclip +
+	        AP_NSKY_REJECT(sky)
+
+	    return (ier)
+
+	case AP_CONSTANT:
+
+	    # Set the sky value to a constant.
+	    AP_SKY_MODE(sky) = AP_SKYBACKGROUND(sky)
+	    AP_SKY_SIG(sky) = AP_SKYSIGMA(nse)
+	    AP_SKY_SKEW(sky) = INDEFR
+	    AP_NSKY(sky) = 0
+	    AP_NSKY_REJECT(sky) = 0
+	    return (AP_OK)
+
+	case AP_SKYFILE:
+
+	    # Read the sky values from a file.
+	    if (sd == NULL)
+		return (AP_NOSKYFILE)
+	    ier = ap_readsky (sd, x, y, AP_SKY_MODE(sky), AP_SKY_SIG(sky),
+		AP_SKY_SKEW(sky), AP_NSKY(sky), AP_NSKY_REJECT(sky))
+	    if (ier == EOF)
+		return (AP_EOFSKYFILE)
+	    else if (ier != 7)
+		return (AP_BADSKYSCAN)
+	    else if (AP_NSKY(sky) <= 0)
+		return (AP_NOSKYAREA)
+	    else
+		return (AP_OK)
+
+	case AP_RADPLOT:
+
+	    # Check the status of the graphics stream.
+	    if (gd == NULL)
+		return (AP_NOGRAPHICS)
+
+	    # Fetch the sky pixels.
+	    ier = apskybuf (ap, im, wx, wy)
+	    if (ier != AP_OK)
+		return (ier)
+
+	    # Clip the data.
+	    if (AP_SLOCLIP(sky) > 0.0 || AP_SHICLIP(sky) > 0.0) {
+		nclip = ap_clip (Memr[AP_SKYPIX(sky)], Memi[AP_INDEX(sky)],
+		    AP_NSKYPIX(sky), AP_SLOCLIP(sky), AP_SHICLIP(sky), ilo,
+		    ihi)
+		if (nclip >= AP_NSKYPIX(sky))
+		    return (AP_NSKY_TOO_SMALL)
+		nsky = AP_NSKYPIX(sky) - nclip
+	    } else {
+		nclip = 0
+		call ap_index (Memi[AP_INDEX(sky)], AP_NSKYPIX(sky))
+		ilo = 1
+		nsky = AP_NSKYPIX(sky) 
+	    }
+
+	    # Mark the sky level on the radial profile plot.
+	    call gactivate (gd, 0)
+	    gt = ap_gtinit (AP_IMROOT(ap), wx, wy)
+	    ier = ap_radplot (gd, gt, Memr[AP_SKYPIX(sky)],
+	        Memi[AP_COORDS(sky)], Memi[AP_INDEX(sky)+ilo-1], nsky,
+		AP_SXC(sky), AP_SYC(sky), AP_SNX(sky), AP_SNY(sky),
+		AP_SCALE(ap), AP_SKY_MODE(sky), AP_SKY_SKEW(sky),
+		AP_SKY_SIG(sky), AP_NSKY(sky), AP_NSKY_REJECT(sky))
+	    AP_NSKY_REJECT(sky) = AP_NBADSKYPIX(sky) + nclip +
+	        AP_NSKY_REJECT(sky)
+	    call ap_gtfree (gt)
+	    call gdeactivate (gd, 0)
+
+	    return (ier)
+
+	case AP_HISTPLOT:
+
+	    # Check the status of the graphics stream.
+	    if (gd == NULL)
+		return (AP_NOGRAPHICS)
+
+	    # Fetch the sky pixels.
+	    ier = apskybuf (ap, im, wx, wy)
+	    if (ier != AP_OK)
+		return (ier)
+
+	    # Initialze the weights.
+	    call amovkr (1.0, Memr[AP_SWGT(sky)], AP_NSKYPIX(sky))
+
+	    # Clip the data.
+	    if (AP_SLOCLIP(sky) > 0.0 || AP_SHICLIP(sky) > 0.0) {
+		nclip = ap_clip (Memr[AP_SKYPIX(sky)], Memi[AP_INDEX(sky)],
+		    AP_NSKYPIX(sky), AP_SLOCLIP(sky), AP_SHICLIP(sky),
+		    ilo, ihi)
+		if (nclip >= AP_NSKYPIX(sky))
+		    return (AP_NSKY_TOO_SMALL)
+		nsky = AP_NSKYPIX(sky) - nclip
+	    } else {
+		nclip = 0
+		call ap_index (Memi[AP_INDEX(sky)], AP_NSKYPIX(sky))
+		ilo = 1
+		nsky = AP_NSKYPIX(sky) 
+	    }
+
+	    # Mark the peak of the histogram on the histogram plot.
+	    #call gactivate (gd, 0)
+	    gt = ap_gtinit (AP_IMROOT(ap), wx, wy)
+	    ier = ap_histplot (gd, gt, Memr[AP_SKYPIX(sky)],
+	        Memr[AP_SWGT(sky)], Memi[AP_INDEX(sky)+ilo-1], nsky,
+		AP_K1(sky), INDEFR, AP_BINSIZE(sky), AP_SMOOTH(sky),
+		AP_SKY_MODE(sky), AP_SKY_SIG(sky), AP_SKY_SKEW(sky),
+		AP_NSKY(sky), AP_NSKY_REJECT(sky))	
+	    AP_NSKY_REJECT(sky) = AP_NBADSKYPIX(sky) + nclip +
+	        AP_NSKY_REJECT(sky)
+	    call ap_gtfree (gt)
+	    #call gdeactivate (gd, 0)
+
+	    return (ier)
+
+	case AP_OFILT:
+
+	    # Fetch the sky pixels.
+	    ier = apskybuf (ap, im, wx, wy)
+	    if (ier != AP_OK)
+		return (ier)
+
+	    # Initialze the weights.
+	    call amovkr (1.0, Memr[AP_SWGT(sky)], AP_NSKYPIX(sky))
+
+	    # Clip the data.
+	    if (AP_SLOCLIP(sky) > 0.0 || AP_SHICLIP(sky) > 0.0) {
+		nclip = ap_clip (Memr[AP_SKYPIX(sky)], Memi[AP_INDEX(sky)],
+		    AP_NSKYPIX(sky), AP_SLOCLIP(sky), AP_SHICLIP(sky),
+		    ilo, ihi)
+		if (nclip >= AP_NSKYPIX(sky))
+		    return (AP_NSKY_TOO_SMALL)
+		nsky = AP_NSKYPIX(sky) - nclip
+	    } else {
+		nclip = 0
+		call ap_index (Memi[AP_INDEX(sky)], AP_NSKYPIX(sky))
+		ilo = 1
+		nsky = AP_NSKYPIX(sky) 
+	    }
+
+	    # Compute the sky value using the histogram of the sky pixels
+	    # and a variation of the optimal filtering technique.
+	    ier = ap_lgsky (Memr[AP_SKYPIX(sky)], Memi[AP_COORDS(sky)],
+	        Memr[AP_SWGT(sky)], Memi[AP_INDEX(sky)+ilo-1], nsky,
+		AP_SNX(sky), AP_SNY(sky), AP_K1(sky), INDEFR,
+		AP_BINSIZE(sky), AP_SMOOTH(sky), AP_SLOREJECT(sky),
+		AP_SHIREJECT(sky), AP_RGROW(sky) * AP_SCALE(ap),
+		AP_SMAXITER(sky), AP_SKY_MODE(sky), AP_SKY_SIG(sky),
+		AP_SKY_SKEW(sky), AP_NSKY(sky), AP_NSKY_REJECT(sky))
+	    AP_NSKY_REJECT(sky) = AP_NBADSKYPIX(sky) + nclip +
+	        AP_NSKY_REJECT(sky)
+
+	    return (ier)
+
+	case AP_GAUSS:
+
+	    # Fetch the sky pixels.
+	    ier = apskybuf (ap, im, wx, wy)
+	    if (ier != AP_OK)
+		return (ier)
+
+	    # Initialze the weights.
+	    call amovkr (1.0, Memr[AP_SWGT(sky)], AP_NSKYPIX(sky))
+
+	    # Clip the data.
+	    if (AP_SLOCLIP(sky) > 0.0 || AP_SHICLIP(sky) > 0.0) {
+		nclip = ap_clip (Memr[AP_SKYPIX(sky)], Memi[AP_INDEX(sky)],
+		    AP_NSKYPIX(sky), AP_SLOCLIP(sky), AP_SHICLIP(sky),
+		    ilo, ihi)
+		if (nclip >= AP_NSKYPIX(sky))
+		    return (AP_NSKY_TOO_SMALL)
+		nsky = AP_NSKYPIX(sky) - nclip
+	    } else {
+		nclip = 0
+		call ap_index (Memi[AP_INDEX(sky)], AP_NSKYPIX(sky))
+		ilo = 1
+		nsky = AP_NSKYPIX(sky) 
+	    }
+
+	    # Compute the sky value by a fitting a skewed Gaussian function
+	    # to the sky pixel histogram.
+	    ier = ap_gauss (Memr[AP_SKYPIX(sky)], Memi[AP_COORDS(sky)],
+	        Memr[AP_SWGT(sky)], Memi[AP_INDEX(sky)+ilo-1], nsky,
+		AP_SNX(sky), AP_SNY(sky), AP_SMAXITER(sky), AP_K1(sky),
+		INDEFR, AP_BINSIZE(sky), AP_SMOOTH(sky),
+		AP_SLOREJECT(sky), AP_SHIREJECT(sky), AP_RGROW(sky) *
+		AP_SCALE(ap), AP_SNREJECT(sky), AP_SKY_MODE(sky),
+		AP_SKY_SIG(sky), AP_SKY_SKEW(sky), AP_NSKY(sky),
+		AP_NSKY_REJECT(sky))
+	    AP_NSKY_REJECT(sky) = AP_NBADSKYPIX(sky) + nclip +
+	        AP_NSKY_REJECT(sky)
+
+	    return (ier)
+
+	case AP_CROSSCOR:
+
+	    # Fetch the sky pixels.
+	    ier = apskybuf (ap, im, wx, wy)
+	    if (ier != AP_OK)
+		return (ier)
+
+	    # Initialze the weights.
+	    call amovkr (1.0, Memr[AP_SWGT(sky)], AP_NSKYPIX(sky))
+
+	    # Clip the data.
+	    if (AP_SLOCLIP(sky) > 0.0 || AP_SHICLIP(sky) > 0.0) {
+		nclip = ap_clip (Memr[AP_SKYPIX(sky)], Memi[AP_INDEX(sky)],
+		    AP_NSKYPIX(sky), AP_SLOCLIP(sky), AP_SHICLIP(sky),
+		    ilo, ihi)
+		if (nclip >= AP_NSKYPIX(sky))
+		    return (AP_NSKY_TOO_SMALL)
+		nsky = AP_NSKYPIX(sky) - nclip
+	    } else {
+		nclip = 0
+		call ap_index (Memi[AP_INDEX(sky)], AP_NSKYPIX(sky))
+		ilo = 1
+		nsky = AP_NSKYPIX(sky) 
+	    }
+
+	    # Fit the sky value  by computing the cross-correlation
+	    # function of the histogram and an estimate of the noise
+	    # distribution.
+	    ier = ap_crosscor (Memr[AP_SKYPIX(sky)], Memi[AP_COORDS(sky)],
+	        Memr[AP_SWGT(sky)], Memi[AP_INDEX(sky)+ilo-1], nsky,
+		AP_SNX(sky), AP_SNY(sky), AP_K1(sky), INDEFR,
+		AP_BINSIZE(sky), AP_SMOOTH(sky), AP_SLOREJECT(sky),
+		AP_SHIREJECT(sky), AP_RGROW(sky) * AP_SCALE(ap),
+		AP_SMAXITER(sky), AP_SKY_MODE(sky), AP_SKY_SIG(sky),
+		AP_SKY_SKEW(sky), AP_NSKY(sky), AP_NSKY_REJECT(sky))
+	    AP_NSKY_REJECT(sky) = AP_NBADSKYPIX(sky) + nclip +
+	        AP_NSKY_REJECT(sky)
+
+	    return (ier)
+
+	default:
+
+	    AP_SKY_MODE(sky) = INDEFR
+	    AP_SKY_SIG(sky) = INDEFR
+	    AP_SKY_SKEW(sky) = INDEFR
+	    AP_NSKY(sky) = AP_NSKYPIX(sky)
+	    AP_NSKY_REJECT(sky) = 0
+	    return (AP_OK)
+	}
+end
diff --git a/noao/digiphot/apphot/fitsky/apgauss.x b/noao/digiphot/apphot/fitsky/apgauss.x
new file mode 100644
index 00000000..860301dc
--- /dev/null
+++ b/noao/digiphot/apphot/fitsky/apgauss.x
@@ -0,0 +1,296 @@
+include <mach.h>
+include <math/nlfit.h>
+include "../lib/fitsky.h"
+
+define	TOL		.001		# Fitting tolerance
+
+# AP_GAUSS -- Procedure to compute the peak, width and skew of the histogram
+# by fitting a skewed Gaussian function to the histogram.
+
+int procedure ap_gauss (skypix, coords, wgt, index, nskypix, snx, sny, maxfit,
+	k1, hwidth, binsize, smooth, losigma, hisigma, rgrow, maxiter,
+	sky_mode, sky_sigma, sky_skew, nsky, nsky_reject)
+
+real	skypix[ARB]		# array of unsorted sky pixels
+int	coords[ARB]		# array of coordinates for region growing
+real	wgt[ARB]		# weights for pixel rejection
+int	index[ARB]		# array of sort indices
+int	nskypix			# the number of sky pixels
+int	snx, sny		# the maximum dimensions of sky raster
+int	maxfit			# maximum number of iterations per fit
+real	k1			# extent of the histogram in skysigma
+real	hwidth			# width of the histogram
+real	binsize			# the size of the histogram in sky sigma
+int	smooth			# smooth the histogram before fitting
+real	losigma, hisigma	# upper and lower sigma rejection criterion
+real	rgrow			# region growing radius in pixels
+int	maxiter			# maximum number of rejection cycles
+real	sky_mode		# computed sky value
+real	sky_sigma		# computed sigma of the sky pixels
+real	sky_skew		# skew of sky pixels
+int	nsky			# number of sky pixels used in fit
+int	nsky_reject		# number of sky pixels rejected
+
+double	dsky, sumpx, sumsqpx, sumcbpx
+int	i, j, nreject, nbins, nker, ier
+pointer	sp, x, hgm, shgm, w
+real	sky_mean, sky_msigma, dmin, dmax, hmin, hmax, dh, locut, hicut, cut
+real	sky_zero
+int	ap_grow_hist2(), aphigmr()
+real	ap_asumr(), apmedr()
+
+begin
+	# Initialize.
+	nsky = nskypix
+	nsky_reject = 0
+	sky_mode = INDEFR
+	sky_sigma = INDEFR
+	sky_skew = INDEFR
+	if (nskypix <= 0)
+	    return (AP_NOSKYAREA)
+
+	# Compute initial guess for sky statistics.
+	sky_zero = ap_asumr (skypix, index, nskypix) / nskypix
+	call ap_ialimr (skypix, index, nskypix, dmin, dmax)
+	call apfimoments (skypix, index, nskypix, sky_zero, sumpx, sumsqpx,
+	    sumcbpx, sky_mean, sky_msigma, sky_skew)
+	sky_mean = apmedr (skypix, index,  nskypix)
+	sky_mean = max (dmin, min (sky_mean, dmax))
+
+	# Compute the width and bin size of the sky histogram.
+	if (! IS_INDEFR(hwidth) && hwidth > 0.0) {
+	    hmin = sky_mean - k1 * hwidth
+	    hmax = sky_mean + k1 * hwidth
+	    dh = binsize * hwidth
+	} else {
+	    cut = min (sky_mean - dmin, dmax - sky_mean, k1 * sky_msigma)
+	    hmin = sky_mean - cut
+	    hmax = sky_mean + cut
+	    dh = binsize * cut / k1
+	}
+
+	# Compute the number of histogram bins and width of smoothing kernel.
+	if (dh <= 0.0) {
+	    nbins = 1
+	    dh = 0.0
+	} else {
+	    nbins = 2 * nint ((hmax - sky_mean) / dh) + 1
+	    dh = (hmax - hmin) / (nbins - 1)
+	}
+
+	# Test for a valid histogram.
+	if (sky_msigma <= dh || dh <= 0.0 ||  k1 <= 0.0 || sky_msigma <= 0.0 ||
+	    nbins < 4) {
+	    sky_mode = sky_mean
+	    sky_sigma = 0.0
+	    sky_skew = 0.0
+	    return (AP_NOHISTOGRAM)
+	}
+
+	# Allocate temporary working space.
+	call smark (sp)
+	call salloc (x, nbins, TY_REAL)
+	call salloc (hgm, nbins, TY_REAL)
+	call salloc (shgm, nbins, TY_REAL)
+	call salloc (w, nbins, TY_REAL)
+
+	# Compute the x array.
+	do i = 1, nbins
+	    Memr[x+i-1] = i
+	call amapr (Memr[x], Memr[x], nbins, 1.0, real (nbins),
+	    hmin + 0.5 * dh, hmax + 0.5 * dh)
+
+	# Accumulate the histogram.
+	call aclrr (Memr[hgm], nbins)
+	nsky_reject = nsky_reject + aphigmr (skypix, wgt, index, nskypix,
+	    Memr[hgm], nbins, hmin, hmax)
+	nsky = nskypix - nsky_reject
+
+	# Do the initial rejection.
+	if (nsky_reject > 0) {
+	    do i = 1, nskypix {
+		if (wgt[index[i]] <= 0.0) {
+		    dsky = skypix[index[i]] - sky_zero
+		    sumpx = sumpx - dsky
+		    sumsqpx = sumsqpx - dsky ** 2
+		    sumcbpx = sumcbpx - dsky ** 3
+		}
+	    }
+	    call apmoments (sumpx, sumsqpx, sumcbpx, nsky, sky_zero,
+	        sky_mean, sky_msigma, sky_skew)
+	}
+
+	# Find the mode, sigma and skew of the histogram.
+	if (smooth == YES) {
+	    nker = max (1, nint (sky_msigma / dh))
+	    #call ap_lucy_smooth (Memr[hgm], Memr[shgm], nbins, nker, 2)
+	    call ap_bsmooth (Memr[hgm], Memr[shgm], nbins, nker, 2)
+	    call ap_hist_mode (Memr[x], Memr[shgm], Memr[w], nbins,
+		sky_mode, sky_sigma, sky_skew, maxfit, TOL, ier)
+	} else
+	    call ap_hist_mode (Memr[x], Memr[hgm], Memr[w], nbins,
+		sky_mode, sky_sigma, sky_skew, maxfit, TOL, ier)
+	sky_mode = max (dmin, min (sky_mode, dmax))
+	if (ier != AP_OK) {
+	    call sfree (sp)
+	    return (ier)
+	}
+	if ((IS_INDEFR(losigma) && IS_INDEFR(hisigma)) || sky_sigma <= dh ||
+	    maxiter < 1) {
+	    call sfree (sp)
+	    return (AP_OK)
+	}
+
+	# Fit histogram with pixel rejection and optional region growing.
+	do i = 1, maxiter {
+
+	    # Compute the new rejection limits.
+	    if (IS_INDEFR(losigma))
+		locut = -MAX_REAL
+	    else
+	        locut = sky_mode - losigma * sky_msigma
+	    if (IS_INDEFR(hisigma))
+		hicut = MAX_REAL
+	    else
+	        hicut = sky_mode + hisigma * sky_msigma
+
+	    # Detect and reject pixels.
+	    nreject = 0
+	    do j = 1, nskypix {
+		if (skypix[index[j]] >= locut && skypix[index[j]] <= hicut)
+		    next
+		if (rgrow > 0.0)
+		    nreject = nreject + ap_grow_hist2 (skypix, coords,
+			wgt, nskypix, sky_zero, index[j], snx, sny, Memr[hgm],
+			nbins, hmin, hmax, rgrow, sumpx, sumsqpx, sumcbpx)
+		else if (wgt[index[j]] > 0.0) {
+		    call ap_hgmsub2 (Memr[hgm], nbins, hmin, hmax,
+		        skypix[index[j]], sky_zero, sumpx, sumsqpx, sumcbpx)
+		    wgt[index[j]] = 0.0
+		    nreject = nreject + 1
+		}
+	    }
+	    if (nreject == 0)
+		break
+	    nsky_reject = nsky_reject + nreject
+	    nsky = nskypix - nsky_reject 
+	    if (nsky <= 0)
+		break
+	    call apmoments (sumpx, sumsqpx, sumcbpx, nsky, sky_zero,
+	        sky_mean, sky_msigma, sky_skew)
+
+	    # Recompute mean, mode, sigma and skew.
+	    if (smooth == YES) {
+		nker = max (1, nint (sky_msigma / dh))
+	        #call ap_lucy_smooth (Memr[hgm], Memr[shgm], nbins, nker, 2)
+	        call ap_bsmooth (Memr[hgm], Memr[shgm], nbins, nker, 2)
+	    	call ap_hist_mode (Memr[x], Memr[shgm], Memr[w], nbins,
+		    sky_mode, sky_sigma, sky_skew, maxfit, TOL, ier)
+	    } else
+	    	call ap_hist_mode (Memr[x], Memr[hgm], Memr[w], nbins,
+		    sky_mode, sky_sigma, sky_skew, maxfit, TOL, ier)
+	    sky_mode = max (dmin, min (sky_mode, dmax))
+	    if (ier != AP_OK)
+		break
+	    if (sky_sigma <= dh)
+		break
+	}
+
+	# Return appropriate error code.
+	call sfree (sp)
+	if (nsky == 0 || nsky_reject == nskypix) {
+	    nsky = 0
+	    nsky_reject = nskypix
+	    sky_mode = INDEFR
+	    sky_sigma = INDEFR
+	    sky_skew = INDEFR
+	    return (AP_NSKY_TOO_SMALL)
+	} else if (ier == AP_NSKY_TOO_SMALL) {
+	    return (AP_NSKY_TOO_SMALL)
+	} else if (ier != AP_OK) {
+	    return (ier)
+	} else
+	    return (AP_OK)
+end
+
+
+# AP_HIST_MODE -- Procedure to fit the skewed Gaussian function to the histogram
+# of the sky pixels.
+
+procedure ap_hist_mode (x, hgm, w, nbins, sky_mode, sky_sigma, sky_skew,
+    maxiter, tol, ier)
+
+real	x[ARB]		# x array
+real	hgm[ARB]	# histogram
+real	w[ARB]		# weights
+int	nbins		# number of bins
+real	sky_mode	# sky value
+real	sky_sigma	# sky sigma
+real	sky_skew	# sky skew
+int	maxiter		# max number of iterations
+real	tol		# tolerances
+int	ier		# error code
+
+extern  gausskew, dgausskew
+int	i, imin, imax, np, fier
+pointer	sp, list, fit, nl
+real	p[4], dp[4], dummy1
+int	locpr()
+
+begin
+	# Allocate working memory.
+	call smark (sp)
+	call salloc (list, 4, TY_INT)
+	call salloc (fit, nbins, TY_REAL)
+
+	# Initialize.
+	do i = 1, 4
+	    Memi[list+i-1] = i
+
+	# Compute initial guesses for the parameters.
+	call ap_alimr (hgm, nbins, dummy1, p[1], imin, imax)
+	p[2] = x[imax]
+	#p[3] = max (sky_sigma ** 2, abs (x[2] - x[1]) ** 2)
+	p[3] = abs ((x[nbins] - x[1]) / 6.0) ** 2
+	p[4] = 0.0
+	np = 4
+
+	# Fit the histogram.
+	call nlinitr (nl, locpr (gausskew), locpr (dgausskew), p, dp, np,
+	    Memi[list], 4, tol, maxiter)
+	call nlfitr (nl, x, hgm, w, nbins, 1, WTS_UNIFORM, fier)
+	call nlvectorr (nl, x, Memr[fit], nbins, 1)
+	call nlpgetr (nl, p, np)
+	call nlfreer (nl)
+
+	call sfree (sp)
+
+	# Return the appropriate error code.
+	ier = AP_OK
+	if (fier == NO_DEG_FREEDOM) {
+	    sky_mode = INDEFR
+	    sky_sigma = INDEFR
+	    sky_skew = INDEFR
+	    ier = AP_NSKY_TOO_SMALL
+	} else {
+	    if (fier == SINGULAR)
+		ier = AP_SKY_SINGULAR
+	    else if (fier == NOT_DONE)
+		ier = AP_SKY_NOCONVERGE
+	    if (p[2] < x[1] || p[2] > x[nbins]) {
+		sky_mode = x[imax]
+		ier = AP_BADPARAMS
+	    } else
+	        sky_mode = p[2]
+	    if (p[3] <= 0.0) {
+		sky_sigma = 0.0
+		sky_skew = 0.0
+		ier = AP_BADPARAMS
+	    } else {
+	        sky_sigma = sqrt (p[3])
+	        sky_skew = 1.743875281 * abs (p[4]) ** (1.0 / 3.0) * sky_sigma
+	        if (p[4] < 0.0)
+		    sky_skew = - sky_skew
+	    }
+	}
+end
diff --git a/noao/digiphot/apphot/fitsky/apgrowhist.x b/noao/digiphot/apphot/fitsky/apgrowhist.x
new file mode 100644
index 00000000..6f01a789
--- /dev/null
+++ b/noao/digiphot/apphot/fitsky/apgrowhist.x
@@ -0,0 +1,136 @@
+# AP_GROW_HIST -- Procedure to reject pixels with region growing.
+
+int procedure ap_grow_hist (skypix, coords, wgt, nskypix, index, snx, sny, hgm,
+    nbins, z1, z2, rgrow)
+
+real	skypix[ARB]		# sky pixels
+int	coords[ARB]		# sky pixel coordinates
+real	wgt[ARB]		# array of weights
+int	nskypix			# number of sky pixels
+int	index			# index of pixel to be rejected
+int	snx, sny		# size of sky subraster
+real	hgm[ARB]		# histogram
+int	nbins			# number of bins
+real	z1, z2			# value of first and last histogram bin
+real	rgrow			# region growing radius
+
+int	j, k, ixc, iyc, ymin, ymax, xmin, xmax, nreject, cstart, c, bin
+real	dh, r2, rgrow2, d
+
+begin
+	# Find the x and y coordinates of the pixel to be rejected.
+	ixc = mod (coords[index], snx)
+	if (ixc == 0)
+	    ixc = snx
+	iyc = (coords[index] - ixc) / snx + 1
+
+	# Find the coordinate space to be used for regions growing.
+	ymin = max (1, int (iyc - rgrow))
+	ymax = min (sny, int (iyc + rgrow))
+	xmin = max (1, int (ixc - rgrow))
+	xmax = min (snx, int (ixc + rgrow))
+	if (ymin <= iyc)
+	    cstart = min (nskypix, max (1, index - int (rgrow) + snx *
+	        (ymin - iyc)))
+	else
+	    cstart = index
+
+	# Perform the region growing.
+	dh = real (nbins - 1) / (z2 - z1)
+	rgrow2 = rgrow ** 2
+	nreject = 0
+	do j = ymin, ymax {
+	    d = rgrow2 - (j - iyc) ** 2
+	    if (d <= 0.0)
+		d = 0.0
+	    else
+		d = sqrt (d)
+	    do k = max (xmin, int (ixc - d)), min (xmax, int (ixc + d)) {
+		c = k + (j - 1) * snx
+		while (coords[cstart] < c && cstart < nskypix)
+		    cstart = cstart + 1
+		r2 = (k - ixc) ** 2 + (j - iyc) ** 2
+		if (r2 <= rgrow2 && c == coords[cstart] && wgt[cstart] > 0.0) {
+		    nreject = nreject + 1
+		    wgt[cstart] = 0.0
+		    if (skypix[cstart] >= z1 && skypix[cstart] <= z2) {
+		        bin = int ((skypix[cstart] - z1) * dh) + 1
+		        hgm[bin] = hgm[bin] - 1.0
+		    }
+		}
+	    }
+	}
+
+	return (nreject)
+end
+
+
+# AP_GROW_HIST2 -- Procedure to reject pixels with region growing.
+
+int procedure ap_grow_hist2 (skypix, coords, wgt, nskypix, sky_zero, index,
+    snx, sny, hgm, nbins, z1, z2, rgrow, sumpx, sumsqpx, sumcbpx)
+
+real	skypix[ARB]		# sky pixels
+int	coords[ARB]		# coordinates
+real	wgt[ARB]		# array of weights
+int	nskypix			# number of sky pixels
+real	sky_zero		# the sky zero point for moment analysis
+int	index			# index of pixel to be rejected
+int	snx, sny		# size of sky subraster
+real	hgm[ARB]		# histogram
+int	nbins			# number of bins
+real	z1, z2			# value of first and last histogram bin
+real	rgrow			# region growing radius
+double	sumpx			# sum of sky values
+double	sumsqpx			# sum of sky values squared
+double	sumcbpx			# sum of sky values cubed
+
+double	dsky
+int	j, k, ixc, iyc, ymin, ymax, xmin, xmax, nreject, cstart, c, bin
+real	dh, r2, rgrow2, d
+
+begin
+	# Find the coordinates of the region growing center.
+	ixc = mod (coords[index], snx)
+	if (ixc == 0)
+	    ixc = snx
+	iyc = (coords[index] - ixc) / snx + 1
+
+	ymin = max (1, int (iyc - rgrow))
+	ymax = min (sny, int (iyc + rgrow))
+	xmin = max (1, int (ixc - rgrow))
+	xmax = min (snx, int (ixc + rgrow))
+	dh = real (nbins - 1) / (z2 - z1)
+	if (ymin <= iyc)
+	    cstart = min (nskypix, max (1, index - int (rgrow) + snx *
+	        (ymin - iyc)))
+	else
+	    cstart = index
+
+	# Perform the region growing.
+	nreject = 0
+	rgrow2 = rgrow ** 2
+	do j = ymin, ymax {
+	    d = sqrt (rgrow2 - (j - iyc) ** 2)
+	    do k = max (xmin, int (ixc - d)), min (xmax, int (ixc + d)) {
+		c = k + (j - 1) * snx
+		while (coords[cstart] < c && cstart < nskypix)
+		    cstart = cstart + 1
+		r2 = (k - ixc) ** 2 + (j - iyc) ** 2
+		if (r2 <= rgrow2 && c == coords[cstart] && wgt[cstart] > 0.0) {
+		    nreject = nreject + 1
+		    wgt[cstart] = 0.0
+		    dsky = skypix[cstart] - sky_zero
+		    sumpx = sumpx - dsky
+		    sumsqpx = sumsqpx - dsky ** 2
+		    sumcbpx = sumcbpx - dsky ** 3
+		    if (skypix[cstart] >= z1 && skypix[cstart] <= z2) {
+		        bin = int ((skypix[cstart] - z1) * dh) + 1
+		        hgm[bin] = hgm[bin] - 1.0
+		    }
+		}
+	    }
+	}
+
+	return (nreject)
+end
diff --git a/noao/digiphot/apphot/fitsky/apgspars.x b/noao/digiphot/apphot/fitsky/apgspars.x
new file mode 100644
index 00000000..8ca7f52c
--- /dev/null
+++ b/noao/digiphot/apphot/fitsky/apgspars.x
@@ -0,0 +1,26 @@
+include "../lib/display.h"
+include "../lib/noise.h"
+include "../lib/fitsky.h"
+
+# AP_SGPARS -- Procedure to fetch the parameters for the fitsky task.
+
+procedure ap_sgpars (ap)
+
+pointer	ap		# pointer to apphot structure
+
+bool	clgetb()
+int	btoi()
+
+begin
+	# Open the apphot structure.
+	call apsinit (ap, AP_MODE, 10.0, 10.0, 2.0, AP_NPOISSON)
+
+	# Get the data dependent parameters.
+	call ap_gdapars (ap)
+
+	# Get the sky fitting parameters.
+	call ap_gsapars (ap)
+
+	# Get radial plots.
+	call apseti (ap, RADPLOTS, btoi (clgetb ("radplots")))
+end
diff --git a/noao/digiphot/apphot/fitsky/aphgmsub.x b/noao/digiphot/apphot/fitsky/aphgmsub.x
new file mode 100644
index 00000000..ae969bba
--- /dev/null
+++ b/noao/digiphot/apphot/fitsky/aphgmsub.x
@@ -0,0 +1,51 @@
+# AP_HGMSUB -- Procedure to subtract a point from an existing  histogram.
+
+procedure ap_hgmsub (hgm, nbins, z1, z2, skypix)
+
+real	hgm[ARB]		# histogram
+int	nbins			# number of bins
+real	z1, z2			# range of histogram
+real	skypix	               	# sky value
+
+int	bin
+real	dh
+
+begin
+	if (skypix < z1 || skypix > z2)
+	    return
+	dh = real (nbins - 1) / (z2 - z1)
+	bin = int ((skypix - z1) * dh) + 1
+	hgm[bin] = hgm[bin] - 1.0
+end
+
+
+# AP_HGMSUB2 -- Procedure to subract points from the accumulated sums
+# and the existing histogram.
+
+procedure ap_hgmsub2 (hgm, nbins, z1, z2, skypix, sky_zero, sumpx, sumsqpx,
+	sumcbpx)
+
+real	hgm[ARB]		# histogram
+int	nbins			# number of bins
+real	z1, z2			# range of histogram
+real	skypix	               	# sky value
+real	sky_zero		# sky zero point for moment analysis
+double	sumpx			# sum of the sky pixel values
+double	sumsqpx			# sum of the squares of the sky pixel values
+double	sumcbpx			# sum of the cubes of the sky pixel values
+
+double	dsky
+int	bin
+real	dh
+
+begin
+	if (skypix < z1 || skypix > z2)
+	    return
+	dsky = skypix - sky_zero
+	sumpx = sumpx - dsky
+	sumsqpx = sumsqpx - dsky ** 2
+	sumcbpx = sumcbpx - dsky ** 3
+	dh = real (nbins - 1) / (z2 - z1)
+	bin = int ((skypix - z1) * dh) + 1
+	hgm[bin] = hgm[bin] - 1.0
+end
diff --git a/noao/digiphot/apphot/fitsky/aphistplot.x b/noao/digiphot/apphot/fitsky/aphistplot.x
new file mode 100644
index 00000000..49c00071
--- /dev/null
+++ b/noao/digiphot/apphot/fitsky/aphistplot.x
@@ -0,0 +1,233 @@
+include <gset.h>
+include <pkg/gtools.h>
+include "../lib/fitsky.h"
+
+# AP_HISTPLOT -- Procedure to the sky value using a plot of the sky histogram
+# and cursor readback.
+
+int procedure ap_histplot (gd, gt, skypix, wgt, index, nskypix, k1, hwidth,
+	binsize, smooth, sky_mode, sky_sigma, sky_skew, nsky, nsky_reject)
+
+pointer	gd			# pointer to graphics stream
+pointer	gt			# pointer to GTOOLS structure
+real	skypix[ARB]		# array of unsorted sky pixels
+real	wgt[ARB]		# array of weights for rejection
+int	index[ARB]		# array of sort indices
+int	nskypix			# number of sky pixels
+real	k1			# rejection criterion
+real	hwidth			# half width of histogram in k1 units
+real	binsize			# histogram binsize in units of sigma
+int	smooth			# smooth the histogram
+real	sky_mode		# sky value
+real	sky_sigma		# sigma of sky pixels
+real	sky_skew		# skew of sky pixels
+int	nsky			# number of sky pixels
+int	nsky_reject		# number of rejected sky pixels
+
+double	dsky, sumpx, sumsqpx, sumcbpx
+int	i, nbins, nker, wcs, key
+pointer	sp, x, hgm, shgm, cmd
+real	dmin, dmax, hmin, hmax, dh, ymin, ymax, symin, symax, wx, wy
+real	u1, u2, v1, v2, x1, x2, y1, y2, cut, sky_mean, sky_zero
+int	clgcur(), aphigmr()
+real	apmedr(), ap_asumr()
+
+begin
+	# Check for valid graphics stream.
+	if (gd == NULL)
+	    return (AP_NOGRAPHICS)
+
+	# Initialize.
+	nsky = nskypix
+	nsky_reject = 0
+	sky_mode = INDEFR
+	sky_sigma = INDEFR
+	sky_skew = INDEFR
+	if (nskypix <= 0)
+	    return (AP_SKY_OUTOFBOUNDS)
+
+	# Compute an initial guess at the sky distribution.
+	sky_zero = ap_asumr (skypix, index, nskypix) / nskypix
+	call ap_ialimr (skypix, index, nskypix, dmin, dmax)
+	call apfimoments (skypix, index, nskypix, sky_zero, sumpx, sumsqpx,
+	    sumcbpx, sky_mean, sky_sigma, sky_skew)
+	sky_mode = apmedr (skypix, index, nskypix)
+	sky_mode = max (dmin, min (sky_mode, dmax))
+
+	# Compute histogram width and binsize.
+	if (! IS_INDEFR(hwidth) && hwidth > 0.0) {
+	    hmin = sky_mode - k1 * hwidth
+	    hmax = sky_mode + k1 * hwidth
+	    dh = binsize * hwidth
+	} else {
+	    cut = min (sky_mode - dmin, dmax - sky_mode, k1 * sky_sigma)
+	    hmin = sky_mode - cut
+	    hmax = sky_mode + cut
+	    dh = binsize * cut / k1
+	}
+
+	# Compute the number of histgram bins and the histogram resolution.
+	if (dh <= 0.0) {
+	    nbins = 1
+	    dh = 0.0
+	} else {
+	    nbins = 2 * nint ((hmax - sky_mode) / dh) + 1
+	    dh = (hmax - hmin) / (nbins - 1)
+	}
+
+	# Test for a valid histogram.
+	if (dh <= 0.0 || k1 <= 0.0 || sky_sigma <= 0.0 || sky_sigma <= dh ||
+	    nbins < 2) {
+	    sky_mode = sky_mode
+	    sky_sigma = 0.0
+	    sky_skew = 0.0
+	    return (AP_NOHISTOGRAM)
+	}
+
+	# Allocate temporary space.
+	call smark (sp)
+	call salloc (x, nbins, TY_REAL)
+	call salloc (hgm, nbins, TY_REAL)
+	call salloc (shgm, nbins, TY_REAL)
+	call salloc (cmd, SZ_LINE, TY_CHAR)
+
+	# Compute the x array and accumulate the histogram.
+	do i = 1, nbins
+	    Memr[x+i-1] = i
+	call amapr (Memr[x], Memr[x], nbins, 1.0, real (nbins),
+	    hmin + 0.5 * dh, hmax + 0.5 * dh) 
+	call aclrr (Memr[hgm], nbins)
+	nsky_reject = nsky_reject + aphigmr (skypix, wgt, index, nskypix,
+	    Memr[hgm], nbins, hmin, hmax)
+	nsky = nskypix - nsky_reject
+
+	# Subtract rejected pixels and recompute the moments.
+	if (nsky_reject > 0) {
+	    do i = 1, nskypix {
+	        if (wgt[index[i]] <= 0.0) {
+		    dsky = skypix[index[i]] - sky_zero
+		    sumpx = sumpx - dsky
+		    sumsqpx = sumsqpx - dsky ** 2
+		    sumcbpx = sumcbpx - dsky ** 3
+	        }
+	    }
+	}
+	call apmoments (sumpx, sumsqpx, sumcbpx, nsky, sky_zero, sky_mean,
+	    sky_sigma, sky_skew)
+
+	# Smooth the histogram and compute the histogram plot limits.
+	if (smooth == YES) {
+	    nker = max (1, nint (sky_sigma / dh))
+	    #call ap_lucy_smooth (Memr[hgm], Memr[shgm], nbins, nker, 2)
+	    call ap_bsmooth (Memr[hgm], Memr[shgm], nbins, nker, 2)
+	    call alimr (Memr[hgm], nbins, ymin, ymax)
+	    call alimr (Memr[shgm], nbins, symin, symax)
+	    ymin = min (ymin, symin)
+	    ymax = max (ymax, symax)
+	} else
+	    call alimr (Memr[hgm], nbins, ymin, ymax)
+
+	# Store the old viewport and window coordinates.
+	call greactivate (gd, 0)
+	call ggview (gd, u1, u2, v1, v2)
+	call ggwind (gd, x1, x2, y1, y2)
+
+	# Plot the raw and smoothed histograms.
+	call ap_set_hplot (gd, gt, hmin, hmax, ymin, ymax, nbins, smooth)
+	if (smooth == YES) {
+	    call ap_plothist (gd, gt, Memr[x], Memr[hgm], nbins, "histogram",
+	        GL_SOLID)
+	    call ap_plothist (gd, gt, Memr[x], Memr[shgm], nbins, "line",
+	        GL_SOLID)
+	} else
+	    call ap_plothist (gd, gt, Memr[x], Memr[hgm], nbins, "histogram",
+	        GL_SOLID)
+
+	# Mark the peak of the histogram with the cursor.
+	call printf (
+	    "Mark histogram peak (%g) [space=mark,q=quit,:.help=help:]")
+	    call pargr (sky_mode)
+	call gscur (gd, sky_mode, (ymin + ymax) / 2.0)
+	while (clgcur ("gcommands", wx, wy, wcs, key, Memc[cmd], SZ_LINE) !=
+	    EOF) {
+		if (key == 'q')
+		    break
+		else
+		    sky_mode = max (hmin, min (wx, hmax))
+	    call printf (
+	        "Mark histogram peak (%g) [space=mark,q=quit,:.help=help:]")
+	        call pargr (sky_mode)
+	}
+
+	# Store the old viewport and window coordinates.
+	call gsview (gd, u1, u2, v1, v2)
+	call gswind (gd, x1, x2, y1, y2)
+
+	# Close up.
+	call gflush (gd)
+	call gdeactivate (gd, 0)
+	call sfree (sp)
+	return (AP_OK)
+end
+
+
+# AP_PLOTHIST -- Procedure plot the histogram.
+
+procedure ap_plothist (gd, gt, x, hgm, nbins, plottype, polytype) 
+
+pointer	gd		# pointer to graphics stream
+pointer	gt		# GTOOLS pointer
+real	x[ARB]		# the histogram bin values
+real	hgm[ARB]	# histogram
+int	nbins		# number of bins
+char	plottype[ARB]	# the plot type "histogram" or "line"
+int	polytype	# polyline type
+
+begin
+	call gt_sets (gt, GTTYPE, plottype)
+	call gt_seti (gt, GTLINE, polytype)
+	call gt_plot (gd, gt, x, hgm, nbins)
+	call gflush (gd)
+end
+
+
+# AP_SET_HPLOT -- Procedure to set up the histogram plot.
+
+procedure ap_set_hplot (gd, gt, xmin, xmax, ymin, ymax, nbins, smooth)
+
+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
+int	nbins		# number of bins
+int	smooth		# smooth histogram
+
+pointer	sp, str
+
+begin
+	# Initialize
+	call gclear (gd)
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	# Set up the gtools parameter string.
+	call sprintf (Memc[str], SZ_LINE,
+	    "Sky Histogram: nbins = %d hmin = %g hmax = %g smooth=%b")
+	    call pargi (nbins)
+	    call pargr (xmin)
+	    call pargr (xmax)
+	    call pargi (smooth)
+	call gt_sets (gt, GTPARAMS, Memc[str])
+
+	# Set up the plot axes.
+	call gt_sets (gt, GTXLABEL, "Sky Values")
+	call gt_sets (gt, GTYLABEL, "Number of Pixels")
+	call gt_setr (gt, GTXMIN, xmin)
+	call gt_setr (gt, GTXMAX, xmax)
+	call gt_setr (gt, GTYMIN, ymin)
+	call gt_setr (gt, GTYMAX, ymax)
+	call gt_swind (gd, gt)
+	call gt_labax (gd, gt)
+
+	call sfree (sp)
+end
diff --git a/noao/digiphot/apphot/fitsky/aplgsky.x b/noao/digiphot/apphot/fitsky/aplgsky.x
new file mode 100644
index 00000000..db5b4a19
--- /dev/null
+++ b/noao/digiphot/apphot/fitsky/aplgsky.x
@@ -0,0 +1,216 @@
+include <mach.h>
+include "../lib/fitsky.h"
+
+define	TOL		0.001		# Fitting tolerance
+
+# AP_LGSKY -- Procedure to fit the peak and width of the histogram using
+# repeated convolutions and a triangle function.
+
+int procedure ap_lgsky (skypix, coords, wgt, index, nskypix, snx, sny, k1,
+	hwidth, binsize, smooth, losigma, hisigma, rgrow, maxiter, sky_mode,
+	sky_sigma, sky_skew, nsky, nsky_reject)
+
+real	skypix[ARB]		# array of sky pixels
+int	coords[ARB]		# array of coordinates of region growing
+real	wgt[ARB]		# array of weights for rejection
+int	index[ARB]		# array of sort indices
+int	nskypix			# the number of sky pixels
+int	snx, sny		# the maximum dimensions of sky raster
+real	k1			# extent of the histogram in skysigma
+real	hwidth			# width of histogram
+real	binsize			# the size of the histogram in sky sigma
+int	smooth			# smooth the histogram before fitting
+real	losigma, hisigma	# upper and lower sigma rejection limits
+real	rgrow			# region growing radius in pixels
+int	maxiter			# maximum number of rejection cycles
+real	sky_mode		# computed sky value
+real	sky_sigma		# computed sigma of the sky pixels
+real	sky_skew		# skew of sky pixels
+int	nsky			# number of sky pixels used in fit
+int	nsky_reject		# number of sky pixels rejected
+
+double 	dsky, sumpx, sumsqpx, sumcbpx
+int	nreject, nbins, nker, i, j, iter
+pointer	sp, hgm, shgm
+real	dmin, dmax, hmin, hmax, dh, locut, hicut, sky_mean, center, cut
+real	sky_zero
+int	ap_grow_hist2(), aphigmr(), aptopt()
+real	ap_asumr(), apmedr(), apmapr()
+
+begin
+	# Initialize.
+	nsky = nskypix
+	nsky_reject = 0
+	sky_mode = INDEFR
+	sky_sigma = INDEFR
+	sky_skew = INDEFR
+	if (nskypix <= 0)
+	    return (AP_NOSKYAREA)
+
+	# Compute a first guess for the parameters.
+	sky_zero = ap_asumr (skypix, index, nskypix) / nskypix
+	call ap_ialimr (skypix, index, nskypix, dmin, dmax)
+	call apfimoments (skypix, index, nskypix, sky_zero, sumpx, sumsqpx,
+	    sumcbpx, sky_mean, sky_sigma, sky_skew)
+	sky_mean = apmedr (skypix, index, nskypix)
+	sky_mean = max (dmin, min (sky_mean, dmax))
+
+	# Compute the width and bin size of histogram.
+	if (! IS_INDEFR(hwidth) && hwidth > 0.0) {
+	    hmin = sky_mean - k1 * hwidth
+	    hmax = sky_mean + k1 * hwidth
+	    dh = binsize * hwidth
+	} else {
+	    cut = min (sky_mean - dmin, dmax - sky_mean, k1 * sky_sigma)
+	    hmin = sky_mean - cut
+	    hmax = sky_mean + cut
+	    dh = binsize * cut / k1
+	}
+
+	# Compute the number of histogram bins and the resolution.
+	# filter.
+	if (dh <= 0.0) {
+	    nbins = 1
+	    dh = 0.0
+	} else {
+	    nbins = 2 * nint ((hmax - sky_mean) / dh) + 1
+	    dh  = (hmax - hmin) / (nbins - 1)
+	}
+
+	# Test for a valid histogram.
+	if (nbins < 2 || k1 <= 0.0 || sky_sigma <= 0.0 || dh <= 0.0 ||
+	    sky_sigma <= dh) {
+	    sky_mode = sky_mean
+	    sky_sigma = 0.0
+	    sky_skew = 0.0
+	    return (AP_NOHISTOGRAM)
+	}
+
+	# Allocate temporary space.
+	call smark (sp)
+	call salloc (hgm, nbins, TY_REAL)
+	call salloc (shgm, nbins, TY_REAL)
+
+	# Accumulate the histogram.
+	call aclrr (Memr[hgm], nbins)
+	nsky_reject = nsky_reject + aphigmr (skypix, wgt, index, nskypix,
+	    Memr[hgm], nbins, hmin, hmax)
+	nsky = nskypix - nsky_reject
+
+	# Perform the initial rejection.
+	if (nsky_reject > 0) {
+	    do i = 1, nskypix {
+	        if (wgt[index[i]] <= 0.0) {
+		    dsky = skypix[index[i]] - sky_zero
+		    sumpx = sumpx - dsky
+		    sumsqpx = sumsqpx - dsky ** 2
+		    sumcbpx = sumcbpx - dsky ** 3
+		}
+	    }
+	    call apmoments (sumpx, sumsqpx, sumcbpx, nsky, sky_zero,
+	        sky_mean, sky_sigma, sky_skew)
+	}
+
+	# Fit the peak of the histogram.
+	center = apmapr ((hmin + hmax) / 2.0, hmin + 0.5 * dh,
+	    hmax + 0.5 * dh, 1.0, real (nbins))
+	if (smooth == YES) {
+	    nker = max (1, nint (sky_sigma / dh))
+	    #call ap_lucy_smooth (Memr[hgm], Memr[shgm], nbins, nker, 2)
+	    call ap_bsmooth (Memr[hgm], Memr[shgm], nbins, nker, 2)
+	    iter =  aptopt (Memr[shgm], nbins, center, sky_sigma / dh,
+	        TOL, maxiter, NO)
+	} else
+	    iter = aptopt (Memr[hgm], nbins, center, sky_sigma / dh, TOL,
+	        maxiter, NO)
+	sky_mode = apmapr (center, 1.0, real (nbins), hmin + 0.5 * dh,
+	    hmax + 0.5 * dh)
+	sky_mode = max (dmin, min (sky_mode, dmax))
+	if (iter < 0) {
+	    call sfree (sp)
+	    return (AP_SKY_NOCONVERGE)
+	}
+	if ((IS_INDEFR(losigma) && IS_INDEFR(hisigma)) || (sky_sigma <= dh) ||
+	    (maxiter < 1)) {
+	    call sfree (sp)
+	    return (AP_OK)
+	}
+
+	# Fit the histogram with pixel rejection and optional region growing.
+	do i = 1, maxiter {
+
+	    # Compute new histogram limits.
+	    if (IS_INDEFR(losigma))
+		locut = -MAX_REAL
+	    else
+	        locut = sky_mode - losigma * sky_sigma
+	    if (IS_INDEFR(hisigma))
+		hicut = MAX_REAL
+	    else
+	        hicut = sky_mode + hisigma * sky_sigma
+
+	    # Detect and reject the pixels.
+	    nreject = 0
+	    do j = 1, nskypix {
+		if (skypix[index[j]] >= locut && skypix[index[j]] <= hicut)
+		    next
+		if (rgrow > 0.0)
+		    nreject = nreject + ap_grow_hist2 (skypix, coords,
+			wgt, nskypix, sky_zero, index[j], snx, sny, Memr[hgm],
+			nbins, hmin, hmax, rgrow, sumpx, sumsqpx, sumcbpx)
+		else if (wgt[index[j]] > 0.0) {
+		    call ap_hgmsub2 (Memr[hgm], nbins, hmin, hmax,
+		        skypix[index[j]], sky_zero, sumpx, sumsqpx, sumcbpx)
+		    wgt[index[j]] = 0.0
+		    nreject = nreject + 1
+		}
+	    }
+	    if (nreject == 0)
+		break
+
+	    # Recompute the data limits.
+	    nsky_reject = nsky_reject + nreject
+	    nsky = nskypix - nsky_reject 
+	    if (nsky <= 0)
+		break
+	    call apmoments (sumpx, sumsqpx, sumcbpx, nsky, sky_zero,
+	        sky_mean, sky_sigma, sky_skew)
+	    if (sky_sigma <= dh)
+		break
+
+	    # Refit the sky.
+	    if (smooth == YES) {
+		nker = max (1, nint (sky_sigma / dh))
+	        #call ap_lucy_smooth (Memr[hgm], Memr[shgm], nbins, nker, 2)
+	        call ap_bsmooth (Memr[hgm], Memr[shgm], nbins, nker, 2)
+	    	iter = aptopt (Memr[shgm], nbins, center, sky_sigma / dh,
+		    TOL, maxiter, NO)
+	    } else
+	    	iter = aptopt (Memr[hgm], nbins, center, sky_sigma / dh,
+		    TOL, maxiter, NO)
+	    sky_mode = apmapr (center, 1.0, real (nbins), hmin + 0.5 * dh,
+	        hmax + 0.5 * dh)
+	    sky_mode = max (dmin, min (sky_mode, dmax))
+	    if (iter < 0) 
+		break
+	}
+
+	# Return an appropriate error code.
+	call sfree (sp)
+	if (nsky == 0 || nsky_reject == nskypix) {
+	    nsky = 0
+	    nsky_reject = nskypix
+	    sky_mode = INDEFR
+	    sky_sigma = INDEFR
+	    sky_skew = INDEFR
+	    return (AP_NSKY_TOO_SMALL)
+	} else if (sky_sigma <= 0.0) {
+	    sky_sigma = 0.0
+	    sky_skew = 0.0
+	    return (AP_OK)
+	} else if (iter < 0) {
+	    return (AP_SKY_NOCONVERGE)
+	} else {
+	    return (AP_OK)
+	}
+end
diff --git a/noao/digiphot/apphot/fitsky/apmean.x b/noao/digiphot/apphot/fitsky/apmean.x
new file mode 100644
index 00000000..f134f02c
--- /dev/null
+++ b/noao/digiphot/apphot/fitsky/apmean.x
@@ -0,0 +1,120 @@
+include <mach.h>
+include "../lib/fitsky.h"
+
+# AP_MEAN -- Procedure to calculate the mean of the sky pixel array.
+
+int procedure ap_mean (skypix, coords, wgt, index, nskypix, snx, sny, losigma,
+	hisigma, rgrow, maxiter, sky_mean, sky_sigma, sky_skew, nsky,
+	nsky_reject)
+
+real	skypix[ARB]		# unsorted array of skypixels
+int	coords[ARB]		# coordinate array for region growing
+real	wgt[ARB]		# the weight array for rejection
+int	index[ARB]		# sorted array of indices
+int	nskypix			# total number of sky pixels
+int	snx, sny		# dimensions of the sky subraster
+real	losigma, hisigma	# number of sky_sigma for rejection
+real	rgrow			# radius of region growing
+int	maxiter			# maximum number of cycles of rejection
+real	sky_mean		# computed sky value
+real	sky_sigma		# the computed sky sigma
+real	sky_skew		# skewness of sky distribution
+int	nsky			# the number of sky pixels used
+int	nsky_reject		# the number of sky pixels rejected
+
+double  dsky, sumpx, sumsqpx, sumcbpx
+int	i, j, nreject
+real	sky_zero, dmin, dmax, locut, hicut
+int	ap_grow_regions()
+real	ap_asumr()
+
+begin
+	# Intialize.
+	nsky = nskypix
+	nsky_reject = 0
+	sky_mean = INDEFR
+	sky_sigma = INDEFR
+	sky_skew = INDEFR
+	if (nskypix <= 0)
+	    return (AP_NOSKYAREA)
+
+	# Compute the mean, sigma and skew.
+	sky_zero = ap_asumr (skypix, index, nskypix) / nskypix
+	call apfimoments (skypix, index, nskypix, sky_zero, sumpx, sumsqpx,
+	    sumcbpx, sky_mean, sky_sigma, sky_skew)
+	call ap_ialimr (skypix, index, nskypix, dmin, dmax)
+	sky_mean = max (dmin, min (sky_mean, dmax))
+
+	# Decide whether to do the rejection cycle.
+	if (maxiter < 1 || (IS_INDEFR(losigma) && IS_INDEFR(hisigma)) ||
+	    sky_sigma <= 0.0)
+	    return (AP_OK)
+
+	# Reject points within k1 * sky_sigma of the median.
+	do i = 1, maxiter {
+
+	    # Compute the new rejection limits.
+	    if (IS_INDEFR(losigma))
+		locut = max (-MAX_REAL, dmin)
+	    else if (i == 1)
+		locut = sky_mean -  min (sky_mean - dmin, dmax - sky_mean,
+		    losigma * sky_sigma)
+	    else
+	        locut = sky_mean - losigma * sky_sigma
+	    if (IS_INDEFR(hisigma))
+		hicut = min (MAX_REAL, dmax)
+	    else if (i == 1)
+		hicut = sky_mean + min (dmax - sky_mean, sky_mean - dmin,
+		    hisigma * sky_sigma)
+	    else
+	        hicut = sky_mean + hisigma * sky_sigma
+
+	    nreject = 0
+	    do j = 1, nskypix {
+	        if (wgt[index[j]] <= 0.0)
+		    next
+	        if (skypix[index[j]] < locut || skypix[index[j]] > hicut) {
+		    if (rgrow > 0.0) {
+		        nreject = ap_grow_regions (skypix, coords, wgt,
+			    nskypix, sky_zero, index[j], snx, sny, rgrow,
+			    sumpx, sumsqpx, sumcbpx)
+		    } else {
+			dsky = (skypix[index[j]] - sky_zero)
+		        sumpx = sumpx - dsky
+		        sumsqpx = sumsqpx - dsky ** 2
+		        sumcbpx = sumcbpx - dsky ** 3
+		        wgt[index[j]] = 0.0 
+		        nreject = nreject + 1
+		    }
+	        }
+	    }
+
+	    # Test the number of rejected pixels.
+	    if (nreject <= 0)
+		break
+	    nsky_reject = nsky_reject + nreject
+
+	    # Test that some pixels actually remain.
+	    nsky = nskypix - nsky_reject
+	    if (nsky <= 0)
+		break
+
+	    # Recompute the mean, sigma and skew.
+	    call apmoments (sumpx, sumsqpx, sumcbpx, nsky, sky_zero,
+	        sky_mean, sky_sigma, sky_skew)
+	    if (sky_sigma <= 0.0)
+		break
+	    sky_mean = max (dmin, min (sky_mean, dmax))
+	}
+
+	# Return an appropriate error code.
+	if (nsky == 0 || nsky_reject == nskypix) {
+	    nsky = 0
+	    nsky_reject = nskypix
+	    sky_mean = INDEFR
+	    sky_sigma = INDEFR
+	    sky_skew = INDEFR
+	    return (AP_NSKY_TOO_SMALL)
+	} else
+	    return (AP_OK)
+end
diff --git a/noao/digiphot/apphot/fitsky/apmedian.x b/noao/digiphot/apphot/fitsky/apmedian.x
new file mode 100644
index 00000000..5998492c
--- /dev/null
+++ b/noao/digiphot/apphot/fitsky/apmedian.x
@@ -0,0 +1,184 @@
+include <mach.h>
+include "../lib/fitsky.h"
+
+define	MEDCUT	0.025
+
+# AP_MEDIAN -- Procedure to calculate the median of the sky pixel array.
+
+int procedure ap_median (skypix, coords, wgt, index, nskypix, snx, sny, losigma,
+	hisigma, rgrow, maxiter, sky_med, sky_sigma, sky_skew, nsky,
+	nsky_reject)
+
+real	skypix[ARB]		# unsorted array of sky pixels
+int	coords[ARB]		# coordinate array for regions growing
+real	wgt[ARB]		# array of weights for rejections
+int	index[ARB]		# sky pixel indices in sort order
+int	nskypix			# total number of sky pixels
+int	snx, sny		# dimensions of the sky subraster
+real	losigma, hisigma	# upper and lower sigma for rejection
+real	rgrow			# radius of region growing
+int	maxiter			# maximum number of cycles of rejection
+real	sky_med			# computed sky value
+real	sky_sigma		# the computed sky sigma
+real	sky_skew		# skewness of sky distribution
+int	nsky			# the number of sky pixels used
+int	nsky_reject		# the number of sky pixels rejected
+
+double  dsky, sumpx, sumsqpx, sumcbpx
+int	i, j, ilo, ihi, il, ih, med, medcut
+real	sky_zero, sky_mean, locut, hicut, dmin, dmax
+int	ap_grow_regions(), apimed()
+real	apsmed(), apwsmed(), ap_asumr()
+
+begin
+	# Intialize.
+	nsky = nskypix
+	nsky_reject = 0
+	sky_med = INDEFR
+	sky_sigma = INDEFR
+	sky_skew = INDEFR
+	if (nskypix <= 0)
+	    return (AP_NOSKYAREA)
+
+	# Sort the sky pixels and compute the median, mean, sigma and skew.
+	# MEDCUT tries to correct for quantization effects
+	call ap_ialimr (skypix, index, nskypix, dmin, dmax)
+	sky_zero = ap_asumr (skypix, index, nskypix) / nskypix
+	medcut = nint (MEDCUT * real (nskypix))
+	sky_med = apsmed (skypix, index, nskypix, medcut)
+	sky_med = max (dmin, min (sky_med, dmax))
+	call apfimoments (skypix, index, nskypix, sky_zero, sumpx, sumsqpx,
+	    sumcbpx, sky_mean, sky_sigma, sky_skew)
+	sky_mean = max (dmin, min (sky_mean, dmax))
+	if (maxiter < 1 || (IS_INDEFR(losigma) && IS_INDEFR(hisigma)) ||
+	    sky_sigma <= 0.0)
+	    return (AP_OK)
+	#call printf (
+	    #"mean=%g med=%g sigma=%g nsky=%d nrej=%d dmin=%g dmax=%g\n")
+	    #call pargr (sky_mean)
+	    #call pargr (sky_med)
+	    #call pargr (sky_sigma)
+	    #call pargi (nsky)
+	    #call pargi (nsky_reject)
+	    #call pargr (dmin)
+	    #call pargr (dmax)
+
+	# Reject points outside losigma * sky_sigma and hisigma * sky_sigma
+	# of the median.
+	ilo = 1
+	ihi = nskypix
+	do i = 1, maxiter {
+
+	    # Compute the new rejection limits.
+	    if (IS_INDEFR(losigma))
+		locut = max (-MAX_REAL, dmin)
+	    else if (i == 1)
+		locut = sky_med - min (sky_med - dmin, dmax - sky_med,
+		    losigma * sky_sigma)
+	    else
+	        locut = sky_med - losigma * sky_sigma
+	    if (IS_INDEFR(hisigma))
+		hicut = min (MAX_REAL, dmax)
+	    else if (i == 1)
+		hicut = sky_med + min (dmax - sky_med, sky_med - dmin,
+		    hisigma * sky_sigma)
+	    else
+	        hicut = sky_med + hisigma * sky_sigma
+	    #call printf ("    locut=%g hicut=%g\n")
+		#call pargr (locut)
+		#call pargr (hicut)
+
+	    # Detect pixels to be rejected.
+	    for (il = ilo; il <= nskypix; il = il + 1) {
+		if (skypix[index[il]] >= locut)
+		    break
+	    }
+	    for (ih = ihi; ih >= 1; ih = ih - 1) {
+		if (skypix[index[ih]] <= hicut)
+		    break
+	    }
+	    if (il == ilo && ih == ihi)
+		break
+
+	    # Reject pixels with optional region growing.
+	    if (rgrow > 0.0) {
+
+		# Reject low side pixels with region growing.
+	        do j = ilo, il - 1 
+		    nsky_reject = nsky_reject + ap_grow_regions (skypix, coords,
+		        wgt, nskypix, sky_zero, index[j], snx, sny, rgrow,
+			sumpx, sumsqpx, sumcbpx)
+
+		# Reject high side pixels with region growing.
+		do j = ih + 1, ihi
+		    nsky_reject = nsky_reject + ap_grow_regions (skypix, coords,
+		        wgt, nskypix, sky_zero, index[j], snx, sny, rgrow,
+			sumpx, sumsqpx, sumcbpx)
+
+		# Recompute the median.
+		# and the median pixel
+		nsky = nskypix - nsky_reject
+		med = apimed (wgt, index, il, ih, (nsky + 1) / 2)
+
+
+	    } else {
+
+		# Recompute the number of sky pixels, the number of rejected
+		# pixels and median pixel
+		nsky_reject = nsky_reject + (il - ilo) + (ihi - ih)
+		nsky = nskypix - nsky_reject
+		med = (ih + il) / 2
+
+		# Reject pixels on the low side.
+		do j = ilo, il - 1 {
+		    dsky = skypix[index[j]] - sky_zero
+		    sumpx = sumpx - dsky
+		    sumsqpx = sumsqpx - dsky ** 2
+		    sumcbpx = sumcbpx - dsky ** 3
+		}
+
+		# Reject pixels on the high side.
+		do j = ih + 1, ihi {
+		    dsky = skypix[index[j]] - sky_zero
+		    sumpx = sumpx - dsky
+		    sumsqpx = sumsqpx - dsky ** 2
+		    sumcbpx = sumcbpx - dsky ** 3
+		}
+	    }
+	    if (nsky <= 0)
+		break
+
+	    # Recompute the median, sigma and skew.
+	    medcut = nint (MEDCUT * real (nsky))
+	    sky_med = apwsmed (skypix, index, wgt, nskypix, med, medcut)
+	    sky_med = max (dmin, min (sky_med, dmax))
+	    call apmoments (sumpx, sumsqpx, sumcbpx, nsky, sky_zero,
+	        sky_mean, sky_sigma, sky_skew)
+	    sky_mean = max (dmin, min (sky_mean, dmax))
+	    #call printf (
+	    #"    mean=%g med=%g sigma=%g nsky=%d nrej=%d dmin=%g dmax=%g\n")
+	        #call pargr (sky_mean)
+	        #call pargr (sky_med)
+	        #call pargr (sky_sigma)
+	        #call pargi (nsky)
+	        #call pargi (nsky_reject)
+	        #call pargr (dmin)
+	        #call pargr (dmax)
+
+	    if (sky_sigma <= 0.0)
+		break
+	    ilo = il
+	    ihi = ih
+	}
+
+	# Return an appropriate error code.
+	if (nsky == 0 || nsky_reject == nskypix) {
+	    nsky = 0
+	    nsky_reject = nskypix
+	    sky_med = INDEFR
+	    sky_sigma = INDEFR
+	    sky_skew = INDEFR
+	    return (AP_NSKY_TOO_SMALL)
+	} else
+	    return (AP_OK)
+end
diff --git a/noao/digiphot/apphot/fitsky/apmode.x b/noao/digiphot/apphot/fitsky/apmode.x
new file mode 100644
index 00000000..fc32ba43
--- /dev/null
+++ b/noao/digiphot/apphot/fitsky/apmode.x
@@ -0,0 +1,185 @@
+include <mach.h>
+include "../lib/fitsky.h"
+
+define	MEDCUT		0.025
+
+# AP_MODE -- Procedure to calculate the mode of the sky pixel array.
+
+int procedure ap_mode (skypix, coords, wgt, index, nskypix, snx, sny, losigma,
+	hisigma, rgrow, maxiter, sky_mode, sky_sigma, sky_skew, nsky,
+	nsky_reject)
+
+real	skypix[ARB]		# unsorted array of skypixels
+int	coords[ARB]		# coordinate array for regions growing
+real	wgt[ARB]		# array of weights for rejection
+int	index[ARB]		# indices in sort order
+int	nskypix			# total number of sky pixels
+int	snx, sny		# dimensions of the sky subraster
+real	losigma, hisigma	# number of sky_sigma for rejection
+real	rgrow			# radius of region growing
+int	maxiter			# maximum number of cycles of rejection
+real	sky_mode		# computed sky value
+real	sky_sigma		# the computed sky sigma
+real	sky_skew		# skew of sky pixels
+int	nsky			# the number of sky pixels used
+int	nsky_reject		# the number of sky pixels rejected
+
+double	dsky, sumpx, sumsqpx, sumcbpx
+int	i, j, ilo, ihi, il, ih, med, medcut
+real	dmin, dmax, locut, hicut, sky_zero, sky_mean, sky_med
+int	ap_grow_regions(), apimed()
+real	apsmed(), apwsmed(), ap_asumr()
+
+begin
+	# Initialize.
+	nsky = nskypix
+	nsky_reject = 0
+	sky_mode = INDEFR
+	sky_sigma = INDEFR
+	sky_skew = INDEFR
+	if (nskypix <= 0)
+	    return (AP_NOSKYAREA)
+
+	# Compute the median, sigma, skew and sky mode.
+	sky_zero = ap_asumr (skypix, index, nskypix) / nskypix
+	call ap_ialimr (skypix, index, nskypix, dmin, dmax)
+	medcut = nint (MEDCUT * real (nskypix))
+	sky_med = apsmed (skypix, index, nskypix, medcut)
+	sky_med = max (dmin, min (sky_med, dmax))
+	call apfimoments (skypix, index, nskypix, sky_zero, sumpx, sumsqpx,
+	    sumcbpx, sky_mean, sky_sigma, sky_skew)
+	sky_mean = max (dmin, min (sky_mean, dmax))
+	if (sky_mean < sky_med)
+	    sky_mode = sky_mean
+	else
+	    sky_mode = 3.0 * sky_med - 2.0 * sky_mean
+	sky_mode = max (dmin, min (sky_mode, dmax))
+	if (maxiter < 1 || (IS_INDEFR(losigma) && IS_INDEFR(hisigma)) ||
+	    sky_sigma <= 0.0)
+	    return (AP_OK)
+
+	# Reject points outside losigma * sky_sigma and hisigma * sky_sigma
+	# of the mode.
+	ilo = 1
+	ihi = nskypix
+	do i = 1, maxiter {
+
+	    # Compute the new rejection limits.
+	    if (i == 1) {
+		if (IS_INDEFR(losigma))
+		    locut = max (-MAX_REAL, dmin)
+		else
+	            locut = sky_med - min (sky_med - dmin, dmax - sky_med,
+		        losigma * sky_sigma)
+		if (IS_INDEFR(hisigma))
+		    hicut = min (MAX_REAL, dmax)
+		else
+	            hicut = sky_med + min (sky_med - dmin, dmax - sky_med,
+		        hisigma * sky_sigma)
+	    } else {
+		if (IS_INDEFR(losigma))
+		    locut = max (-MAX_REAL, dmin)
+		else
+	            locut = sky_mode - losigma * sky_sigma
+		if (IS_INDEFR(hisigma))
+		    hicut = min (MAX_REAL, dmax)
+		else
+	            hicut = sky_mode + hisigma * sky_sigma
+	    }
+	    #call eprintf ("i=%d mean=%g median=%g mode=%g locut=%g hicut=%g\n")
+		#call pargi (i)
+		#call pargr (sky_mean)
+		#call pargr (sky_med)
+		#call pargr (sky_mode)
+		#call pargr (locut)
+		#call pargr (hicut)
+
+	    # Perform lower bound pixel rejection.
+	    for (il = ilo; il <= nskypix; il = il + 1) {
+		if (skypix[index[il]] >= locut)
+		    break
+	    }
+
+	    # Perform upper bound pixel rejection.
+	    for (ih = ihi; ih >= 1; ih = ih - 1) {
+		if (skypix[index[ih]] <= hicut)
+		    break
+	    }
+	    if (il == ilo && ih == ihi)
+		break
+
+	    # Compute number of rejected pixels with optional region growing.
+	    if (rgrow > 0.0) {
+
+		# Reject lower bound pixels with region growing.
+	        do j = ilo, il - 1 
+		    nsky_reject = nsky_reject + ap_grow_regions (skypix, coords,
+		        wgt, nskypix, sky_zero, index[j], snx, sny, rgrow,
+			sumpx, sumsqpx, sumcbpx)
+
+		# Reject upper bound pixels with region growing.
+		do j = ih + 1, ihi
+		    nsky_reject = nsky_reject + ap_grow_regions (skypix, coords,
+		        wgt, nskypix, sky_zero, index[j], snx, sny, rgrow,
+			sumpx, sumsqpx, sumcbpx)
+
+		# Compute the new median.
+		nsky = nskypix - nsky_reject
+		med = apimed (wgt, index, il, ihi, (nsky + 1) / 2)
+
+	    } else {
+
+		# Recompute the number of sky pixels, number of rejected
+		# pixels and the median.
+		nsky_reject = nsky_reject + (il - ilo) + (ihi - ih)
+	        nsky = nskypix - nsky_reject
+		med = (ih + il) / 2
+
+		# Reject number of lower bound pixels.
+		do j = ilo, il - 1 {
+		    dsky = skypix[index[j]] - sky_zero
+		    sumpx = sumpx - dsky
+		    sumsqpx = sumsqpx - dsky ** 2
+		    sumcbpx = sumcbpx - dsky ** 3
+		}
+
+		# Reject number of upper bound pixels.
+		do j = ih + 1, ihi {
+		    dsky = skypix[index[j]] - sky_zero
+		    sumpx = sumpx - dsky
+		    sumsqpx = sumsqpx - dsky ** 2
+		    sumcbpx = sumcbpx - dsky ** 3
+		}
+	    }
+	    if (nsky <= 0)
+		break
+
+	    # Recompute mean, median, mode, sigma and skew.
+	    medcut = nint (MEDCUT * real (nsky))
+	    sky_med = apwsmed (skypix, index, wgt, nskypix, med, medcut)
+	    sky_med = max (dmin, min (sky_med, dmax))
+	    call apmoments (sumpx, sumsqpx, sumcbpx, nsky, sky_zero,
+	        sky_mean, sky_sigma, sky_skew)
+	    sky_mean = max (dmin, min (sky_mean, dmax))
+	    if (sky_mean < sky_med)
+		sky_mode = sky_mean
+	    else
+	        sky_mode = 3.0 * sky_med - 2.0 * sky_mean
+	    sky_mode = max (dmin, min (sky_mode, dmax))
+	    if (sky_sigma <= 0.0)
+		break
+	    ilo = il
+	    ihi = ih
+	}
+
+	# Return an appropriate error code.
+	if (nsky == 0 || nsky_reject == nskypix) {
+	    nsky = 0
+	    nsky_reject = nskypix
+	    sky_mode = INDEFR
+	    sky_sigma = INDEFR
+	    sky_skew = INDEFR
+	    return (AP_NSKY_TOO_SMALL)
+	} else
+	    return (AP_OK)
+end
diff --git a/noao/digiphot/apphot/fitsky/appsky.x b/noao/digiphot/apphot/fitsky/appsky.x
new file mode 100644
index 00000000..8da9773e
--- /dev/null
+++ b/noao/digiphot/apphot/fitsky/appsky.x
@@ -0,0 +1,103 @@
+include "../lib/apphot.h"
+include "../lib/fitsky.h"
+
+# AP_PSSKY -- Procedure to write the results of the fitsky task to
+# the output file.
+
+procedure ap_pssky (ap, fd, id, ld, ier)
+
+pointer	ap		# pointer to apphot structure
+int	fd		# output file descriptor
+int	id		# sequence number of star
+int	ld		# list number of star
+int	ier		# error code
+
+real	apstatr()
+
+begin
+	# Return if NULL file descriptor.
+	if (fd == NULL)
+	    return
+
+	# Print the object id and computed sky values.
+	call ap_wid (ap, fd, apstatr (ap, OSXCUR), apstatr (ap, OSYCUR), id,
+	    ld, '\\')
+	call ap_wsres (ap, fd, ier, ' ')
+end
+
+
+# AP_QSPSKY -- Procedure to print a quick summary of the fitsky task on the
+# standard output.
+
+procedure ap_qspsky (ap, ier)
+
+pointer	ap		# pointer to apphot structure
+int	ier		# error code
+
+pointer	sp, imname
+int	apstati()
+real	apstatr()
+
+begin
+	# Print out the results on the standard output.
+	call smark (sp)
+	call salloc (imname, SZ_FNAME, TY_CHAR)
+	call apstats (ap, IMROOT, Memc[imname], SZ_FNAME)
+	call printf ( "%s  %8.2f %8.2f  %8g %8g ") 
+	    call pargstr (Memc[imname])
+	    call pargr (apstatr (ap, OSXCUR))
+	    call pargr (apstatr (ap, OSYCUR))
+	    call pargr (apstatr (ap, SKY_MODE))
+	    call pargr (apstatr (ap, SKY_SIGMA))
+	call printf ("%8g  %5d %5d  %s\n")
+	    call pargr (apstatr (ap, SKY_SKEW))
+	    call pargi (apstati (ap, NSKY))
+	    call pargi (apstati (ap, NSKY_REJECT))
+	    if (ier != AP_OK)
+		call pargstr ("err")
+	    else
+		call pargstr ("ok")
+	call sfree (sp)
+end
+
+
+# AP_QASPSKY -- Procedure to print a quick summary of the fitsky task on the
+# standard output.
+
+procedure ap_qaspsky (ap, ier)
+
+pointer	ap		# pointer to apphot structure
+int	ier		# error code
+
+int	apstati()
+real	apstatr()
+
+begin
+	# Print out the results on the standard output.
+	call printf ( "    Averages  %8g %8g ") 
+	    call pargr (apstatr (ap, SKY_MODE))
+	    call pargr (apstatr (ap, SKY_SIGMA))
+	call printf ("%8g  %5d %5d  %s\n\n")
+	    call pargr (apstatr (ap, SKY_SKEW))
+	    call pargi (apstati (ap, NSKY))
+	    call pargi (apstati (ap, NSKY_REJECT))
+	    if (ier != AP_OK)
+		call pargstr ("err")
+	    else
+		call pargstr ("ok")
+end
+
+# AP_SPHDR -- Procedure to write the banner for the fitsky task to the
+# output file.
+
+procedure ap_sphdr (ap, fd)
+
+pointer	ap		# pointer to apphot structure
+int	fd		# output file descriptor
+
+begin
+	if (fd == NULL)
+	    return
+	call ap_idhdr (ap, fd)
+	call ap_shdr (ap, fd)
+end
diff --git a/noao/digiphot/apphot/fitsky/appspars.x b/noao/digiphot/apphot/fitsky/appspars.x
new file mode 100644
index 00000000..50da0f84
--- /dev/null
+++ b/noao/digiphot/apphot/fitsky/appspars.x
@@ -0,0 +1,22 @@
+include "../lib/display.h"
+
+# AP_PSPARS -- Procedure to write out the current sky fitting parameters
+# to the parameter files.
+
+procedure ap_pspars (ap)
+
+pointer	ap		# pointer to apphot structure
+
+bool	itob()
+int	apstati()
+
+begin
+	# Write the data dependent parameters.
+	call ap_dapars (ap)
+
+	# Write the sky fitting parameters.
+	call ap_sapars (ap)
+
+	# Radial profile plots
+	call clputb ("radplots", itob (apstati (ap, RADPLOTS)))
+end
diff --git a/noao/digiphot/apphot/fitsky/apradplot.x b/noao/digiphot/apphot/fitsky/apradplot.x
new file mode 100644
index 00000000..d9d19992
--- /dev/null
+++ b/noao/digiphot/apphot/fitsky/apradplot.x
@@ -0,0 +1,91 @@
+include "../lib/fitsky.h"
+
+# AP_RADPLOT -- Procedure to compute the mode, sigma, and skew of the sky by
+# eye using a radial profile plot of the sky pixels and cursor readback.
+
+int procedure ap_radplot (gd, gt, skypix, coords, index, nskypix, sxc, syc,
+	snx, sny, scale, sky_mode, sky_skew, sky_sigma, nsky,
+	nsky_reject)
+
+pointer	gd			# pointer to graphics stream
+pointer	gt			# pointer to gtools structure
+real	skypix[ARB]		# array of sky pixels
+int	coords[ARB]		# array of sky coordinates
+int	index[ARB]		# the index array
+int	nskypix			# number of sky pixels
+real	sxc, syc		# sky subraster center
+int	snx, sny		# sky subraster size
+real	scale			# the image scale
+real	sky_mode		# computed sky value
+real	sky_sigma		# computed sigma of sky pixels
+real	sky_skew		# computed skew of sky pixels
+int	nsky			# number of sky pixels used in fit
+int	nsky_reject		# number of rejected sky pixels
+
+double	sumpx, sumsqpx, sumcbpx
+int	wcs, key
+pointer	sp, r, cmd
+real	wx, wy, xmin, xmax, ymin, ymax, u1, u2, v1, v2, x1, x2, y1, y2
+real	sky_zero
+int	clgcur()
+real	ap_asumr()
+
+begin
+	if (gd == NULL)
+	    return (AP_NOGRAPHICS)
+
+	# Initialize.
+	nsky = nskypix
+	nsky_reject = 0
+	sky_mode = INDEFR
+	sky_sigma = INDEFR
+	sky_skew = INDEFR
+	if (nskypix <= 0)
+	    return (AP_SKY_OUTOFBOUNDS)
+
+	call smark (sp)
+	call salloc (r, nskypix, TY_REAL)
+	call salloc (cmd, SZ_LINE, TY_CHAR)
+
+	# Compute an initial guess at the data characteristics.
+	nsky = nskypix
+	sky_zero = ap_asumr (skypix, index, nskypix) / nskypix
+	call apfimoments (skypix, index, nskypix, sky_zero, sumpx, sumsqpx,
+	    sumcbpx, sky_mode, sky_sigma, sky_skew)
+
+	# Store the old window and viewport coordinates.
+	call ggview (gd, u1, u2, v1, v2)
+	call ggwind (gd, x1, x2, y1, y2)
+
+	# Compute the radial profile
+	call ap_xytor (coords, index, Memr[r], nskypix, sxc, syc, snx)
+	call alimr (Memr[r], nskypix, xmin, xmax)
+	call alimr (skypix, nskypix, ymin, ymax)
+
+	# Plot the radial profile.
+	call gclear (gd)
+	call ap_rset (gd, gt, xmin, xmax, ymin, ymax, scale)
+	call ap_plotrad (gd, gt, Memr[r], skypix, nskypix, "plus")
+
+	# Mark the sky level with the cursor.
+	call printf ("Mark sky level (%g) [space=mark,q=quit,:.help=help]:")
+	    call pargr (sky_mode)
+	call gscur (gd, (xmin + xmax) / 2.0, sky_mode)
+	while (clgcur ("gcommands", wx, wy, wcs, key, Memc[cmd], SZ_LINE) !=
+	    EOF) {
+	    if (key == 'q')
+		break
+	    else
+		sky_mode = wy
+	    call printf ("Mark sky level (%g) [space=mark,q=quit,:.help=help]:")
+	        call pargr (sky_mode)
+	    call gscur (gd, (xmin + xmax) / 2.0, sky_mode)
+	}
+
+	# Store the old window and viewport coordinates.
+	call gsview (gd, u1, u2, v1, v2)
+	call gswind (gd, x1, x2, y1, y2)
+
+	call sfree (sp)
+	return (AP_OK)
+end
diff --git a/noao/digiphot/apphot/fitsky/apreadsky.x b/noao/digiphot/apphot/fitsky/apreadsky.x
new file mode 100644
index 00000000..3db9a183
--- /dev/null
+++ b/noao/digiphot/apphot/fitsky/apreadsky.x
@@ -0,0 +1,40 @@
+# AP_READSKY -- Procedure to read sky values, sigma and skew from a file
+# The x and y positions, sky mode, sigma, and skew values, number of sky
+# pixels and number of rejected sky pixels are assumed to be columns
+# 1 to 7 respectively.
+
+int procedure ap_readsky (fd, x, y, sky_mode, sky_sigma, sky_skew, nsky,
+    nsky_reject)
+
+int	fd		# sky file descriptor
+real	x, y		# center of sky annulus
+real	sky_mode	# sky valye
+real	sky_sigma	# sky sigma
+real	sky_skew	# skew of sky pixels
+int	nsky		# number of sky pixels
+int	nsky_reject	# number of rejected pixesl
+
+int	stat
+int	fscan(), nscan()
+
+begin
+	# Initialize.
+	sky_mode = INDEFR
+	sky_sigma = INDEFR
+	sky_skew = INDEFR
+	nsky = 0
+	nsky_reject = 0
+
+	# Read in and decode a sky file text line.
+	stat = fscan (fd)
+	if (stat == EOF)
+	    return (EOF)
+	call gargr (x)
+	call gargr (y)
+	call gargr (sky_mode)
+	call gargr (sky_sigma)
+	call gargr (sky_skew)
+	call gargi (nsky)
+	call gargi (nsky_reject)
+	return (nscan ())
+end
diff --git a/noao/digiphot/apphot/fitsky/aprefitsky.x b/noao/digiphot/apphot/fitsky/aprefitsky.x
new file mode 100644
index 00000000..bc4ed6c2
--- /dev/null
+++ b/noao/digiphot/apphot/fitsky/aprefitsky.x
@@ -0,0 +1,371 @@
+include "../lib/apphotdef.h"
+include "../lib/apphot.h"
+include "../lib/noisedef.h"
+include "../lib/fitskydef.h"
+include "../lib/fitsky.h"
+
+# APREFITSKY -- Procedure to fit the sky using the pixels currently stored in
+# the sky fitting buffers.
+
+int procedure aprefitsky (ap, im, gd)
+
+pointer	ap		# pointer to the apphot structure
+pointer	im		# the input image descriptor
+pointer	gd		# pointer to graphics stream
+
+int	ier, nclip, nsky, ilo, ihi
+pointer	sky, nse, gt
+int	ap_mode(), ap_centroid(), ap_histplot(), ap_median()
+int	ap_radplot(), ap_gauss(), ap_lgsky(), ap_crosscor()
+int	ap_mean(), ap_clip()
+pointer	ap_gtinit()
+
+begin
+	# Initialize.
+	sky = AP_PSKY(ap)
+	nse = AP_NOISE(ap)
+	AP_SKY_MODE(sky) = INDEFR
+	AP_SKY_SIG(sky) = INDEFR
+	AP_SKY_SKEW(sky) = INDEFR
+	AP_NSKY(sky) = 0
+	AP_NSKY_REJECT(sky) = 0
+	if (IS_INDEFR(AP_SXCUR(sky)) || IS_INDEFR(AP_SYCUR(sky))) {
+            AP_OSXCUR(sky) = AP_SXCUR(sky)
+            AP_OSYCUR(sky) = AP_SYCUR(sky)
+	} else {
+            switch (AP_WCSOUT(ap)) {
+            case WCS_WORLD, WCS_PHYSICAL:
+                call ap_ltoo (ap, AP_SXCUR(sky), AP_SYCUR(sky), AP_OSXCUR(sky),
+		    AP_OSYCUR(sky), 1)
+            case WCS_TV:
+                call ap_ltov (im, AP_SXCUR(sky), AP_SYCUR(sky), AP_OSXCUR(sky),
+		    AP_OSYCUR(sky), 1)
+            default:
+                AP_OSXCUR(sky) = AP_SXCUR(sky)
+                AP_OSYCUR(sky) = AP_SYCUR(sky)
+            }
+	}
+
+	if (IS_INDEFR(AP_SXCUR(sky)) || IS_INDEFR(AP_SYCUR(sky)))
+	    return (AP_NOSKYAREA)
+
+	switch (AP_SKYFUNCTION(sky)) {
+
+	case AP_MEAN:
+
+	    # Initialize the weights.
+	    call amovkr (1.0, Memr[AP_SWGT(sky)], AP_NSKYPIX(sky))
+
+	    # Clip the data.
+	    if (AP_SLOCLIP(sky) > 0.0 || AP_SHICLIP(sky) > 0.0) {
+		nclip = ap_clip (Memr[AP_SKYPIX(sky)], Memi[AP_INDEX(sky)],
+		    AP_NSKYPIX(sky), AP_SLOCLIP(sky), AP_SHICLIP(sky), ilo,
+		    ihi)
+		if (nclip >= AP_NSKYPIX(sky))
+		    return (AP_NSKY_TOO_SMALL)
+		nsky = AP_NSKYPIX(sky) - nclip
+	    } else {
+		nclip = 0
+		call ap_index (Memi[AP_INDEX(sky)], AP_NSKYPIX(sky))
+		ilo = 1
+		nsky = AP_NSKYPIX(sky) 
+	    }
+
+	    ier = ap_mean (Memr[AP_SKYPIX(sky)], Memi[AP_COORDS(sky)],
+	        Memr[AP_SWGT(sky)], Memi[AP_INDEX(sky)+ilo-1], nsky,
+		AP_SNX(sky), AP_SNY(sky), AP_SLOREJECT(sky),
+		AP_SHIREJECT(sky), AP_RGROW(sky) * AP_SCALE(ap),
+		AP_SNREJECT(sky), AP_SKY_MODE(sky), AP_SKY_SIG(sky),
+		AP_SKY_SKEW(sky), AP_NSKY(sky), AP_NSKY_REJECT(sky))
+	    AP_NSKY_REJECT(sky) = AP_NBADSKYPIX(sky) + nclip +
+	        AP_NSKY_REJECT(sky) 
+
+	    return (ier)
+
+	case AP_MEDIAN:
+
+	    # Initialize the weights.
+	    call amovkr (1.0, Memr[AP_SWGT(sky)], AP_NSKYPIX(sky))
+
+	    # Clip the data.
+	    if (AP_SLOCLIP(sky) > 0.0 || AP_SHICLIP(sky) > 0.0) {
+		nclip = ap_clip (Memr[AP_SKYPIX(sky)], Memi[AP_INDEX(sky)],
+		    AP_NSKYPIX(sky), AP_SLOCLIP(sky), AP_SHICLIP(sky), ilo,
+		    ihi)
+		if (nclip >= AP_NSKYPIX(sky))
+		    return (AP_NSKY_TOO_SMALL)
+		nsky = AP_NSKYPIX(sky) - nclip
+	    } else {
+		nclip = 0
+		call apqsort (Memr[AP_SKYPIX(sky)], Memi[AP_INDEX(sky)],
+		    Memi[AP_INDEX(sky)], AP_NSKYPIX(sky))
+		ilo = 1
+		nsky = AP_NSKYPIX(sky) 
+	    }
+
+	    ier = ap_median (Memr[AP_SKYPIX(sky)], Memi[AP_COORDS(sky)],
+	        Memr[AP_SWGT(sky)], Memi[AP_INDEX(sky)+ilo-1], nsky,
+		AP_SNX(sky), AP_SNY(sky), AP_SLOREJECT(sky),
+		AP_SHIREJECT(sky), AP_RGROW(sky) * AP_SCALE(ap),
+		AP_SNREJECT(sky), AP_SKY_MODE(sky), AP_SKY_SIG(sky),
+		AP_SKY_SKEW(sky), AP_NSKY(sky), AP_NSKY_REJECT(sky))
+	    AP_NSKY_REJECT(sky) = AP_NBADSKYPIX(sky) + nclip +
+	        AP_NSKY_REJECT(sky) 
+	    return (ier)
+
+	case AP_MODE:
+
+	    # Initialize the weights.
+	    call amovkr (1.0, Memr[AP_SWGT(sky)], AP_NSKYPIX(sky))
+
+	    # Clip the data.
+	    if (AP_SLOCLIP(sky) > 0.0 || AP_SHICLIP(sky) > 0.0) {
+		nclip = ap_clip (Memr[AP_SKYPIX(sky)], Memi[AP_INDEX(sky)],
+		    AP_NSKYPIX(sky), AP_SLOCLIP(sky), AP_SHICLIP(sky), ilo,
+		    ihi)
+		if (nclip >= AP_NSKYPIX(sky))
+		    return (AP_NSKY_TOO_SMALL)
+		nsky = AP_NSKYPIX(sky) - nclip
+	    } else {
+		nclip = 0
+		call apqsort (Memr[AP_SKYPIX(sky)], Memi[AP_INDEX(sky)],
+		    Memi[AP_INDEX(sky)], AP_NSKYPIX(sky))
+		ilo = 1
+		nsky = AP_NSKYPIX(sky) 
+	    }
+
+	    ier = ap_mode (Memr[AP_SKYPIX(sky)], Memi[AP_COORDS(sky)],
+	        Memr[AP_SWGT(sky)], Memi[AP_INDEX(sky)+ilo-1], nsky,
+		AP_SNX(sky), AP_SNY(sky), AP_SLOREJECT(sky),
+		AP_SHIREJECT(sky), AP_RGROW(sky) * AP_SCALE(ap),
+		AP_SNREJECT(sky), AP_SKY_MODE(sky), AP_SKY_SIG(sky),
+		AP_SKY_SKEW(sky), AP_NSKY(sky), AP_NSKY_REJECT(sky))
+	    AP_NSKY_REJECT(sky) = AP_NBADSKYPIX(sky) + nclip +
+	        AP_NSKY_REJECT(sky) 
+
+	    return (ier)
+
+	case AP_CENTROID:
+
+	    # Initialize the weights.
+	    call amovkr (1.0, Memr[AP_SWGT(sky)], AP_NSKYPIX(sky))
+
+	    # Clip the data.
+	    if (AP_SLOCLIP(sky) > 0.0 || AP_SHICLIP(sky) > 0.0) {
+		nclip = ap_clip (Memr[AP_SKYPIX(sky)], Memi[AP_INDEX(sky)],
+		    AP_NSKYPIX(sky), AP_SLOCLIP(sky), AP_SHICLIP(sky), ilo,
+		    ihi)
+		if (nclip >= AP_NSKYPIX(sky))
+		    return (AP_NSKY_TOO_SMALL)
+		nsky = AP_NSKYPIX(sky) - nclip
+	    } else {
+		nclip = 0
+		call ap_index (Memi[AP_INDEX(sky)], AP_NSKYPIX(sky))
+		ilo = 1
+		nsky = AP_NSKYPIX(sky) 
+	    }
+
+	    ier = ap_centroid (Memr[AP_SKYPIX(sky)], Memi[AP_COORDS(sky)],
+	        Memr[AP_SWGT(sky)+ilo-1], Memi[AP_INDEX(sky)+ilo-1], nsky,
+		AP_SNX(sky), AP_SNY(sky), AP_K1(sky), INDEFR,
+		AP_BINSIZE(sky), AP_SMOOTH(sky), AP_SLOREJECT(sky),
+		AP_SHIREJECT(sky), AP_RGROW(sky) * AP_SCALE(ap),
+		AP_SMAXITER(sky), AP_SKY_MODE(sky), AP_SKY_SIG(sky),
+		AP_SKY_SKEW(sky), AP_NSKY(sky), AP_NSKY_REJECT(sky))
+	    AP_NSKY_REJECT(sky) = AP_NBADSKYPIX(sky) + nclip +
+	        AP_NSKY_REJECT(sky) 
+
+	    return (ier)
+
+	case AP_CONSTANT:
+
+	    AP_SKY_MODE(sky) = AP_SKYBACKGROUND(sky)
+	    AP_SKY_SIG(sky) = AP_SKYSIGMA(nse)
+	    AP_SKY_SKEW(sky) = INDEFR
+	    AP_NSKY(sky) = 0
+	    AP_NSKY_REJECT(sky) = 0
+	    return (AP_OK)
+
+	case AP_SKYFILE:
+
+	    return (AP_OK)
+
+	case AP_RADPLOT:
+
+	    # Check the status of the graphics stream.
+	    if (gd == NULL)
+		return (AP_NOGRAPHICS)
+
+	    # Clip the data.
+	    if (AP_SLOCLIP(sky) > 0.0 || AP_SHICLIP(sky) > 0.0) {
+		nclip = ap_clip (Memr[AP_SKYPIX(sky)], Memi[AP_INDEX(sky)],
+		    AP_NSKYPIX(sky), AP_SLOCLIP(sky), AP_SHICLIP(sky), ilo,
+		    ihi)
+		if (nclip >= AP_NSKYPIX(sky))
+		    return (AP_NSKY_TOO_SMALL)
+		nsky = AP_NSKYPIX(sky) - nclip
+	    } else {
+		nclip = 0
+		call ap_index (Memi[AP_INDEX(sky)], AP_NSKYPIX(sky))
+		ilo = 1
+		nsky = AP_NSKYPIX(sky) 
+	    }
+
+	    call gactivate (gd, 0)
+	    gt = ap_gtinit (AP_IMROOT(ap), AP_SXCUR(sky), AP_SYCUR(sky))
+	    ier = ap_radplot (gd, gt, Memr[AP_SKYPIX(sky)],
+	        Memi[AP_COORDS(sky)], Memi[AP_INDEX(sky)+ilo-1], nsky,
+		AP_SXC(sky), AP_SYC(sky), AP_SNX(sky), AP_SNY(sky),
+		AP_SCALE(ap), AP_SKY_MODE(sky), AP_SKY_SIG(sky),
+		AP_SKY_SKEW(sky), AP_NSKY(sky), AP_NSKY_REJECT(sky))
+	    AP_NSKY_REJECT(sky) = AP_NBADSKYPIX(sky) + nclip +
+	        AP_NSKY_REJECT(sky) 
+	    call ap_gtfree (gt)
+	    call gdeactivate (gd, 0)
+
+	    return (ier)
+
+	case AP_HISTPLOT:
+
+	    # Check the status of the graphics stream.
+	    if (gd == NULL)
+		return (AP_NOGRAPHICS)
+
+	    # Initialize the weights.
+	    call amovkr (1.0, Memr[AP_SWGT(sky)], AP_NSKYPIX(sky))
+
+	    # Clip the data.
+	    if (AP_SLOCLIP(sky) > 0.0 || AP_SHICLIP(sky) > 0.0) {
+		nclip = ap_clip (Memr[AP_SKYPIX(sky)], Memi[AP_INDEX(sky)],
+		    AP_NSKYPIX(sky), AP_SLOCLIP(sky), AP_SHICLIP(sky), ilo,
+		    ihi)
+		if (nclip >= AP_NSKYPIX(sky))
+		    return (AP_NSKY_TOO_SMALL)
+		nsky = AP_NSKYPIX(sky) - nclip
+	    } else {
+		nclip = 0
+		call ap_index (Memi[AP_INDEX(sky)], AP_NSKYPIX(sky))
+		ilo = 1
+		nsky = AP_NSKYPIX(sky) 
+	    }
+
+	    #call gactivate (gd, 0)
+	    gt = ap_gtinit (AP_IMROOT(ap), AP_SXCUR(sky), AP_SYCUR(sky))
+	    ier = ap_histplot (gd, gt, Memr[AP_SKYPIX(sky)],
+	        Memr[AP_SWGT(sky)], Memi[AP_INDEX(sky)+ilo-1], nsky,
+		AP_K1(sky), INDEFR, AP_BINSIZE(sky), AP_SMOOTH(sky),
+		AP_SKY_MODE(sky), AP_SKY_SIG(sky), AP_SKY_SKEW(sky),
+		AP_NSKY(sky), AP_NSKY_REJECT(sky))	
+	    AP_NSKY_REJECT(sky) = AP_NBADSKYPIX(sky) + nclip +
+	        AP_NSKY_REJECT(sky) 
+	    call ap_gtfree (gt)
+
+	    return (ier)
+
+	case AP_OFILT:
+
+	    # Initialize the weights.
+	    call amovkr (1.0, Memr[AP_SWGT(sky)], AP_NSKYPIX(sky))
+
+	    # Clip the data.
+	    if (AP_SLOCLIP(sky) > 0.0 || AP_SHICLIP(sky) > 0.0) {
+		nclip = ap_clip (Memr[AP_SKYPIX(sky)], Memi[AP_INDEX(sky)],
+		    AP_NSKYPIX(sky), AP_SLOCLIP(sky), AP_SHICLIP(sky), ilo,
+		    ihi)
+		if (nclip >= AP_NSKYPIX(sky))
+		    return (AP_NSKY_TOO_SMALL)
+		nsky = AP_NSKYPIX(sky) - nclip
+	    } else {
+		nclip = 0
+		call ap_index (Memi[AP_INDEX(sky)], AP_NSKYPIX(sky))
+		ilo = 1
+		nsky = AP_NSKYPIX(sky) 
+	    }
+
+	    ier = ap_lgsky (Memr[AP_SKYPIX(sky)], Memi[AP_COORDS(sky)],
+	        Memr[AP_SWGT(sky)], Memi[AP_INDEX(sky)+ilo-1], nsky,
+		AP_SNX(sky), AP_SNY(sky), AP_K1(sky), INDEFR,
+		AP_BINSIZE(sky), AP_SMOOTH(sky), AP_SLOREJECT(sky),
+		AP_SHIREJECT(sky), AP_RGROW(sky) * AP_SCALE(ap),
+		AP_SMAXITER(sky), AP_SKY_MODE(sky), AP_SKY_SIG(sky),
+		AP_SKY_SKEW(sky), AP_NSKY(sky), AP_NSKY_REJECT(sky))
+	    AP_NSKY_REJECT(sky) = AP_NBADSKYPIX(sky) + nclip +
+	        AP_NSKY_REJECT(sky) 
+
+	    return (ier)
+
+	case AP_GAUSS:
+
+	    # Initialize the weights.
+	    call amovkr (1.0, Memr[AP_SWGT(sky)], AP_NSKYPIX(sky))
+
+	    # Clip the data.
+	    if (AP_SLOCLIP(sky) > 0.0 || AP_SHICLIP(sky) > 0.0) {
+		nclip = ap_clip (Memr[AP_SKYPIX(sky)], Memi[AP_INDEX(sky)],
+		    AP_NSKYPIX(sky), AP_SLOCLIP(sky), AP_SHICLIP(sky), ilo,
+		    ihi)
+		if (nclip >= AP_NSKYPIX(sky))
+		    return (AP_NSKY_TOO_SMALL)
+		nsky = AP_NSKYPIX(sky) - nclip
+	    } else {
+		nclip = 0
+		call ap_index (Memi[AP_INDEX(sky)], AP_NSKYPIX(sky))
+		ilo = 1
+		nsky = AP_NSKYPIX(sky) 
+	    }
+
+	    ier = ap_gauss (Memr[AP_SKYPIX(sky)], Memi[AP_COORDS(sky)],
+	        Memr[AP_SWGT(sky)], Memi[AP_INDEX(sky)+ilo-1], nsky,
+		AP_SNX(sky), AP_SNY(sky), AP_SMAXITER(sky),
+		AP_K1(sky), INDEFR, AP_BINSIZE(sky), AP_SMOOTH(sky),
+		AP_SLOREJECT(sky), AP_SHIREJECT(sky), AP_RGROW(sky) *
+		AP_SCALE(ap), AP_SNREJECT(sky), AP_SKY_MODE(sky),
+		AP_SKY_SIG(sky), AP_SKY_SKEW(sky), AP_NSKY(sky),
+		AP_NSKY_REJECT(sky))
+	    AP_NSKY_REJECT(sky) = AP_NBADSKYPIX(sky) + nclip +
+	        AP_NSKY_REJECT(sky) 
+
+	    return (ier)
+
+	case AP_CROSSCOR:
+
+	    # Initialize the weights.
+	    call amovkr (1.0, Memr[AP_SWGT(sky)], AP_NSKYPIX(sky))
+
+	    # Clip the data.
+	    if (AP_SLOCLIP(sky) > 0.0 || AP_SHICLIP(sky) > 0.0) {
+		nclip = ap_clip (Memr[AP_SKYPIX(sky)], Memi[AP_INDEX(sky)],
+		    AP_NSKYPIX(sky), AP_SLOCLIP(sky), AP_SHICLIP(sky), ilo,
+		    ihi)
+		if (nclip >= AP_NSKYPIX(sky))
+		    return (AP_NSKY_TOO_SMALL)
+		nsky = AP_NSKYPIX(sky) - nclip
+	    } else {
+		nclip = 0
+		call ap_index (Memi[AP_INDEX(sky)], AP_NSKYPIX(sky))
+		ilo = 1
+		nsky = AP_NSKYPIX(sky) 
+	    }
+
+	    ier = ap_crosscor (Memr[AP_SKYPIX(sky)], Memi[AP_COORDS(sky)],
+	        Memr[AP_SWGT(sky)], Memi[AP_INDEX(sky)+ilo-1], nsky,
+		AP_SNX(sky), AP_SNY(sky), AP_K1(sky), INDEFR,
+		AP_BINSIZE(sky), AP_SMOOTH(sky), AP_SLOREJECT(sky),
+		AP_SHIREJECT(sky), AP_RGROW(sky) * AP_SCALE(ap),
+		AP_SMAXITER(sky), AP_SKY_MODE(sky), AP_SKY_SIG(sky),
+		AP_SKY_SKEW(sky), AP_NSKY(sky), AP_NSKY_REJECT(sky))
+	    AP_NSKY_REJECT(sky) = AP_NBADSKYPIX(sky) + nclip +
+	        AP_NSKY_REJECT(sky) 
+
+	    return (ier)
+
+	default:
+
+	    AP_SKY_MODE(sky) = INDEFR
+	    AP_SKY_SIG(sky) = INDEFR
+	    AP_SKY_SKEW(sky) = INDEFR
+	    AP_NSKY(sky) = AP_NSKYPIX(sky)
+	    AP_NSKY_REJECT(sky) = 0
+	    return (AP_OK)
+	}
+end
diff --git a/noao/digiphot/apphot/fitsky/aprgrow.x b/noao/digiphot/apphot/fitsky/aprgrow.x
new file mode 100644
index 00000000..902b91d6
--- /dev/null
+++ b/noao/digiphot/apphot/fitsky/aprgrow.x
@@ -0,0 +1,64 @@
+# AP_GROW_REGIONS -- Perform region growing around a rejected pixel.
+
+int procedure ap_grow_regions (skypix, coords, wgt, nskypix, sky_zero,
+	index, snx, sny, rgrow, sumpx, sumsqpx, sumcbpx)
+
+real	skypix[ARB]		# sky pixels
+int	coords[ARB]		# sky cordinates
+real	wgt[ARB]		# weights
+int	nskypix			# total number of sky pixels
+real	sky_zero		# sky zero point for moment analysis
+int	index			# index of pixel to be rejected
+int	snx, sny		# size of sky subraster
+real	rgrow			# region growing radius
+double	sumpx, sumsqpx, sumcbpx	# sum and sum of squares of sky pixels
+
+double	dsky
+int	j, k, ixc, iyc, xmin, xmax, ymin, ymax, cstart, c, nreject
+real	rgrow2, r2, d
+
+begin
+	# Find the center of the region to be rejected.
+	ixc = mod (coords[index], snx)
+	if (ixc == 0)
+	    ixc = snx
+	iyc = (coords[index] - ixc) / snx + 1
+
+	# Define the region to be searched.
+	rgrow2 = rgrow ** 2
+	ymin = max (1, int (iyc - rgrow)) 
+	ymax = min (sny, int (iyc + rgrow))
+	xmin = max (1, int (ixc - rgrow))
+	xmax = min (snx, int (ixc + rgrow)) 
+	if (ymin <= iyc)
+	    cstart = min (nskypix, max (1, index - int (rgrow) + snx *
+	    (ymin - iyc)))
+	else
+	    cstart = index
+
+	# Reject the pixels.
+	nreject = 0
+	do j = ymin, ymax {
+	    d = rgrow2 - (j - iyc) ** 2
+	    if (d <= 0.0)
+		d = 0.0
+	    else
+		d = sqrt (d)
+	    do k = max (xmin, int (ixc - d)), min (xmax, int (ixc + d)) {
+		c = k + (j - 1) * snx
+	        while (coords[cstart] < c && cstart < nskypix)
+	    	    cstart = cstart + 1
+		r2 = (k - ixc) ** 2 + (j - iyc) ** 2
+		if (r2 <= rgrow2 && c == coords[cstart] && wgt[cstart] > 0.0) {
+		    dsky = skypix[cstart] - sky_zero
+		    sumpx = sumpx - dsky
+		    sumsqpx = sumsqpx - dsky ** 2
+		    sumcbpx = sumcbpx - dsky ** 3
+		    nreject = nreject + 1
+		    wgt[cstart] = 0.0
+		}
+	    }
+	}
+
+	return (nreject)
+end
diff --git a/noao/digiphot/apphot/fitsky/apsconfirm.x b/noao/digiphot/apphot/fitsky/apsconfirm.x
new file mode 100644
index 00000000..43c0b168
--- /dev/null
+++ b/noao/digiphot/apphot/fitsky/apsconfirm.x
@@ -0,0 +1,74 @@
+include "../lib/apphot.h"
+include "../lib/noise.h"
+include "../lib/fitsky.h"
+
+# AP_SCONFIRM -- Procedure to confirm the critical fitsky parameters.
+
+procedure ap_sconfirm (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	annulus, dannulus, skysigma, datamin, datamax
+int	apstati()
+real	apstatr(), ap_vannulus(), ap_vdannulus(), ap_vsigma()
+real	ap_vdatamin(), ap_vdatamax()
+
+begin
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	call printf ("\n")
+
+	# Confirm the sky fitting algorithm.
+	call ap_vsstring (ap, Memc[str], SZ_FNAME)
+
+	# Confirm the remaining parameters.
+	if (apstati (ap, SKYFUNCTION) != AP_CONSTANT &&
+	    apstati (ap, SKYFUNCTION) != AP_SKYFILE) {
+
+	    # Confirm the sky annulus parameter.
+	    annulus = ap_vannulus (ap)
+
+	    # Confirm the width of the sky annulus.
+	    dannulus = ap_vdannulus (ap)
+
+	} else {
+
+	    annulus = apstatr (ap, ANNULUS)
+	    dannulus = apstatr (ap, DANNULUS)
+
+	}
+
+	# Confirm the sky sigma parameter.
+	if (apstati (ap, SKYFUNCTION) != AP_SKYFILE)
+	    skysigma = ap_vsigma (ap)
+	else
+	    skysigma = apstatr (ap, SKYSIGMA)
+
+	# 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_SSTRING, Memc[str], UN_SALGORITHM,
+		"sky fitting algorithm")
+	    call ap_rparam (out, KY_ANNULUS, annulus, UN_SSCALEUNIT,
+	        "inner radius of sky annulus")
+	    call ap_rparam (out, KY_DANNULUS, dannulus, UN_SSCALEUNIT,
+	        "width of the sky annulus")
+	    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/fitsky/apserrors.x b/noao/digiphot/apphot/fitsky/apserrors.x
new file mode 100644
index 00000000..80c8f0a7
--- /dev/null
+++ b/noao/digiphot/apphot/fitsky/apserrors.x
@@ -0,0 +1,42 @@
+include "../lib/fitsky.h"
+
+# AP_SERRORS -- Program to print out detailed fitsky error messages when the
+# program is run in interactive mode.
+
+procedure ap_serrors (ap, ier)
+
+pointer	ap		# pointer to apphot structure (not used)
+int	ier		# integer error code
+
+begin
+	switch (ier) {
+	case AP_NOSKYAREA:
+	    call printf ("The are no pixels in the sky annulus.\n")
+	case AP_SKY_OUTOFBOUNDS:
+	    call printf ("The sky annulus is outside of the image.\n")
+	case AP_NOHISTOGRAM:
+	    call printf ("The sky histogram has no width.\n")
+	case AP_FLAT_HIST:
+	    call printf ("The sky histogram is flat or concave.\n")
+	case AP_NSKY_TOO_SMALL:
+	    call printf ("The number of sky points is too small.\n")
+	case AP_SKY_SINGULAR:
+	    call printf ("The sky fit is singular.\n")
+	case AP_SKY_NOCONVERGE:
+	    call printf ("The sky fit did not converge.\n")
+	case AP_NOGRAPHICS:
+	    call printf ("Interactive graphics are not available.\n")
+	case AP_NOSKYFILE:
+	    call printf (
+	        "The text file containing sky values does not exist.\n")
+	case AP_EOFSKYFILE:
+	    call printf ("The sky file is at EOF.\n")
+	case AP_BADSKYSCAN:
+	    call printf (
+	    "An error occurred in decoding the current line in the sky file.\n")
+	case AP_BADPARAMS:
+	    call printf ("Out of range mode or -ve sigma in Gaussian fit.\n")
+	default:
+	    call printf ("")
+	}
+end
diff --git a/noao/digiphot/apphot/fitsky/apsfree.x b/noao/digiphot/apphot/fitsky/apsfree.x
new file mode 100644
index 00000000..63d41090
--- /dev/null
+++ b/noao/digiphot/apphot/fitsky/apsfree.x
@@ -0,0 +1,48 @@
+include "../lib/apphotdef.h"
+include "../lib/fitskydef.h"
+
+# APSFREE -- Procedure to free the sky fitting structure.
+
+procedure apsfree (ap)
+
+pointer	ap		# pointer to the apphot structure
+
+begin
+	if (ap == NULL)
+	    return
+	if (AP_NOISE(ap) != NULL)
+	    call ap_noisecls (ap)
+	if (AP_PDISPLAY(ap) != NULL)
+	    call ap_dispcls (ap)
+	if (AP_PSKY(ap) != NULL)
+	    call ap_skycls (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_SKYCLS -- Procedure to close up the sky fitting arrays.
+
+procedure ap_skycls (ap)
+
+pointer	ap		# pointer to the apphot structure
+
+pointer	sky
+
+begin
+	sky = AP_PSKY(ap)
+	if (sky == NULL)
+	    return
+	if (AP_SKYPIX(sky) != NULL)
+	    call mfree (AP_SKYPIX(sky), TY_REAL)
+	if (AP_INDEX(sky) != NULL)
+	    call mfree (AP_INDEX(sky), TY_INT)
+	if (AP_COORDS(sky) != NULL)
+	    call mfree (AP_COORDS(sky), TY_INT)
+	if (AP_SWGT(sky) != NULL)
+	    call mfree (AP_SWGT(sky), TY_REAL)
+	call mfree (AP_PSKY(ap), TY_STRUCT)
+end
diff --git a/noao/digiphot/apphot/fitsky/apsinit.x b/noao/digiphot/apphot/fitsky/apsinit.x
new file mode 100644
index 00000000..e180835f
--- /dev/null
+++ b/noao/digiphot/apphot/fitsky/apsinit.x
@@ -0,0 +1,116 @@
+include "../lib/apphotdef.h"
+include "../lib/fitskydef.h"
+include "../lib/fitsky.h"
+
+# APSINIT - Procedure to initialize the sky fitting structure and parameters.
+
+procedure apsinit (ap, function, annulus, dannulus, fwhmpsf, noise)
+
+pointer	ap		# pointer to the apphot structure
+int	function	# sky fitting algorithm
+real	annulus		# radius of sky annulus
+real	dannulus	# width of sky annulus
+real	fwhmpsf		# FWHM of the PSF
+int	noise		# noise function
+
+begin
+	# Set up the object parameters.
+	call malloc (ap, LEN_APSTRUCT, TY_STRUCT)
+
+	# Set up the global apphot package parameters.
+	call ap_defsetup (ap, fwhmpsf)
+
+	# Set up the noise model parameters.
+	call ap_noisesetup (ap, noise)
+
+	# Set up the display options.
+	call ap_dispsetup (ap)
+
+	# Initialize the sky fitting parameters.
+	call ap_skysetup (ap, function, annulus, dannulus)
+
+	# Unused structures are set to null.
+	AP_PCENTER(ap) = NULL
+	AP_PPHOT(ap) = NULL
+	AP_POLY(ap) = NULL
+	AP_PPSF(ap) = NULL
+	AP_RPROF(ap) = NULL
+end
+
+
+# AP_SKYSETUP -- Procedure to set up the sky fitting arrays and parameters.
+
+procedure ap_skysetup (ap, function, annulus, dannulus)
+
+pointer	ap		# pointer to apphot structure
+int	function	# sky fitting function
+real	annulus		# inner radius of sky annulus
+real	dannulus	# outer radius of sky annulus
+
+pointer	sky
+
+begin
+	call malloc (AP_PSKY(ap), LEN_SKYSTRUCT, TY_STRUCT)
+	sky = AP_PSKY(ap)
+	AP_SXCUR(sky) = INDEFR
+	AP_SYCUR(sky) = INDEFR
+
+	# Initialize the sky fitting parameters.
+	AP_SKYFUNCTION(sky) = function
+	switch (function) {
+	case AP_CONSTANT:
+	    call strcpy ("constant", AP_SSTRING(sky), SZ_FNAME)
+	case AP_MODE:
+	    call strcpy ("mode", AP_SSTRING(sky), SZ_FNAME)
+	case AP_CENTROID:
+	    call strcpy ("centroid", AP_SSTRING(sky), SZ_FNAME)
+	case AP_SKYFILE:
+	    call strcpy ("file", AP_SSTRING(sky), SZ_FNAME)
+	case AP_HISTPLOT:
+	    call strcpy ("histplot", AP_SSTRING(sky), SZ_FNAME)
+	case AP_RADPLOT:
+	    call strcpy ("radplot", AP_SSTRING(sky), SZ_FNAME)
+	case AP_MEDIAN:
+	    call strcpy ("median", AP_SSTRING(sky), SZ_FNAME)
+	case AP_GAUSS:
+	    call strcpy ("gauss", AP_SSTRING(sky), SZ_FNAME)
+	case AP_OFILT:
+	    call strcpy ("ofilt", AP_SSTRING(sky), SZ_FNAME)
+	case AP_CROSSCOR:
+	    call strcpy ("crosscor", AP_SSTRING(sky), SZ_FNAME)
+	case AP_MEAN:
+	    call strcpy ("mean", AP_SSTRING(sky), SZ_FNAME)
+	default:
+	    AP_SKYFUNCTION(sky) = DEF_SKYFUNCTION
+	    call strcpy ("mode", AP_SSTRING(sky), SZ_FNAME)
+	}
+
+	AP_SKYBACKGROUND(sky) = DEF_SKYVALUE
+	AP_ANNULUS(sky) = annulus
+	AP_DANNULUS(sky) = dannulus
+	AP_K1(sky) = DEF_K1
+	AP_BINSIZE(sky) = DEF_BINSIZE
+	AP_SMOOTH(sky) = DEF_SMOOTH
+	AP_SLOCLIP(sky) = DEF_SLOCLIP
+	AP_SHICLIP(sky) = DEF_SHICLIP
+	AP_SMAXITER(sky) = DEF_SMAXITER
+	AP_RGROW(sky) = DEF_RGROW
+	AP_SNREJECT(sky) = DEF_SNREJECT
+	AP_SLOREJECT(sky) = DEF_SLOREJECT
+	AP_SHIREJECT(sky) = DEF_SHIREJECT
+
+	# Initialize the sky pixel buffers.
+	AP_LENSKYBUF(sky) = 0
+	AP_NSKYPIX(sky) = 0
+	AP_SKYPIX(sky) = NULL
+	AP_INDEX(sky) = NULL
+	AP_COORDS(sky) = NULL
+	AP_SWGT(sky) = NULL
+
+	# Initialize results parameters.
+	AP_SKY_MODE(sky) = INDEFR
+	AP_SKY_SIG(sky) = INDEFR
+	AP_SKY_SKEW(sky) = INDEFR
+	AP_NSKY(sky) = 0
+	AP_NSKY_REJECT(sky) = 0
+end
diff --git a/noao/digiphot/apphot/fitsky/apsky.x b/noao/digiphot/apphot/fitsky/apsky.x
new file mode 100644
index 00000000..40faac3d
--- /dev/null
+++ b/noao/digiphot/apphot/fitsky/apsky.x
@@ -0,0 +1,346 @@
+include <ctype.h>
+include <gset.h>
+include <imhdr.h>
+include "../lib/apphot.h"
+include "../lib/display.h"
+include "../lib/fitsky.h"
+
+define	HELPFILE	"apphot$fitsky/fitsky.key"
+
+# APSKY -- Procedure to interactively determine sky values in an annular
+# region around a list of objects.
+
+int procedure apsky (ap, im, cl, sd, gd, mgd, id, out, stid, interactive, cache)
+
+pointer ap			# pointer to apphot structure
+pointer	im			# pointer to IRAF image
+int	cl			# starlist file descriptor
+int	sd			# the sky file descriptor
+pointer	gd			# pointer to graphcis descriptor
+pointer	mgd			# pointer to graphics metacode file
+pointer	id			# pointer to image display stream
+int	out			# output file descriptor
+int	stid			# output file sequence number
+int	interactive		# interactive mode
+int	cache			# cache the input image pixels
+
+real	wx, wy, xlist, ylist
+pointer sp, cmd
+int	wcs, key, colonkey, newimage, newobject, newsky, newlist, ier
+int	ip, ltid, oid, prev_num, req_num, buf_size, req_size, old_size
+int	memstat
+
+real	apstatr()
+int	clgcur(), apfitsky(), aprefitsky(), apgscur(), ctoi(), 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
+	newsky = 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
+		newsky = YES
+	    }
+
+	    # Loop over the colon commands.
+	    switch (key) {
+
+	    # Quit.
+	    case 'q':
+		if (interactive == YES) {
+		    if (apgqverify ("fitsky", ap, key) == YES) {
+			call sfree (sp)
+			return (apgtverify (key))
+		    }
+		} else {
+		    call sfree (sp)
+		    return (NO)
+		}
+
+	    # Print the error messages.
+	    case 'e':
+		if (interactive == YES)
+		    call ap_serrors (ap, ier)
+
+	    # Print the help page.
+	    case '?':
+		if ((id != NULL) && (id == gd))
+		    call gpagefile (id, HELPFILE, "")
+		else if (interactive == YES)
+		    call pagefile (HELPFILE, "[space=morehelp,q=quit,?=help]")
+
+	    # Rewind the coordinate file.
+	    case 'r':
+		if (cl != NULL) {
+		    call seek (cl, BOFL)
+		    ltid = 0
+		} else if (interactive == YES)
+		    call printf ("No coordinate list\n")
+
+	    # Verify the critical parameters.
+	    case 'v':
+		call ap_sconfirm (ap, out, stid)
+		newobject = YES
+		newsky = YES
+
+	    # Save the sky fitting parameters.
+	    case 'w':
+		call ap_pspars (ap)
+
+	    # Draw a centered radial profile plot.
+	    case 'd':
+	        if (interactive == YES) {
+		    call ap_qrad (ap, im, wx, wy, gd)
+		    newobject = YES
+		    newsky = YES
+	        }
+
+	    # Interactively set up sky fitting parameters.
+	    case 'i':
+		if (interactive == YES) {
+		    call ap_sradsetup (ap, im, wx, wy, gd, out, stid)
+		    newobject = YES
+		    newsky = YES
+		}
+
+	    # Process fitsky colon commands.
+	    case ':':
+		for (ip = 1; IS_WHITE(Memc[cmd+ip-1]); ip = ip + 1)
+		    ;
+		colonkey = Memc[cmd+ip-1]
+		switch (colonkey) {
+		case 'm', 'n':
+
+		    # Show/set fitsky commands.
+		    if (Memc[cmd+ip] != EOS && Memc[cmd+ip] != ' ') {
+		        call apskycolon (ap, im, cl, out, stid, ltid,
+			    Memc[cmd], newimage, newobject, newsky)
+			goto endswitch_
+		    }
+
+                    # No coordinate list.
+                    if (cl == NULL) {
+                        if (interactive == YES)
+                            call printf ("No coordinate list\n")
+                        goto endswitch_
+                    }
+
+		    # Get next object from the 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
+			newsky = YES
+                        goto endswitch_
+		    }
+
+		    # Measure the next object.	
+		    ier = apfitsky (ap, im, xlist, ylist, sd, gd)
+		    if (id != NULL) {
+		        call apmark (ap, id, NO, apstati (ap, MKSKY), NO)
+			if (id == gd)
+			    call gflush (id)
+			else
+			    call gframe (id)
+		    }
+		    call ap_splot (ap, stid, gd, apstati (ap, RADPLOTS))
+		    if (interactive == YES)
+		        call ap_qspsky (ap, ier)
+		    if (stid == 1)
+		        call ap_param (ap, out, "fitsky")
+		    call ap_splot (ap, stid, mgd, YES)
+		    call ap_pssky (ap, out, stid, ltid, ier)
+		    stid = stid + 1
+		    newobject = NO; newsky = NO
+
+		default:
+		    call apskycolon (ap, im, cl, out, stid, ltid, Memc[cmd],
+		        newimage, newobject, newsky)
+		}
+
+		if (newimage == YES) {
+		    if ((id != NULL) && (gd != id))
+		        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 sky and store the results.
+	    case 'f', ' ':
+		if (newobject == YES)
+		    ier = apfitsky (ap, im, wx, wy, sd, gd)
+		else if (newsky == YES)
+		    ier = aprefitsky (ap, im, gd)
+		if (id != NULL) {
+		    call apmark (ap, id, NO, apstati (ap, MKSKY), NO)
+		    if (id == gd)
+			call gflush (id)
+		    else
+		        call gframe (id)
+		}
+		call ap_splot (ap, stid, gd, apstati (ap, RADPLOTS))
+		if (interactive == YES)
+		    call ap_qspsky (ap, ier)
+		newobject = NO; newsky = NO
+
+		if (key == ' ') {
+		    if (stid == 1)
+		        call ap_param (ap, out, "fitsky")
+		    if (newlist == YES)
+		        call ap_pssky (ap, out, stid, ltid, ier)
+		    else
+		        call ap_pssky (ap, out, stid, 0, ier)
+		    call ap_splot (ap, stid, mgd, YES)
+		    stid = stid + 1
+		}
+
+	    # Get, fit the next object in the list.
+	    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 the next object.
+		newlist = YES
+		if (key == 'm') {
+		    newobject = YES
+		    newsky = YES
+                    goto endswitch_
+		}
+
+		# Measure the next object.
+		ier = apfitsky (ap, im, xlist, ylist, sd, gd)
+		if (id != NULL) {
+		    call apmark (ap, id, NO, apstati (ap, MKSKY), NO)
+		    if (id == gd)
+			call gflush (id)
+		    else
+			call gframe (id)
+		}
+		call ap_splot (ap, stid, gd, apstati (ap, RADPLOTS))
+		if (interactive == YES)
+		    call ap_qspsky (ap, ier)
+		if (stid == 1)
+		    call ap_param (ap, out, "fitsky")
+		call ap_pssky (ap, out, stid, ltid, ier)
+		call ap_splot (ap, stid, mgd, YES)
+		stid = stid + 1
+		newobject = NO
+		newsky = NO
+
+	    # Process the remainder of the list.
+	    case 'l':
+		if (cl != NULL) {
+		    ltid = ltid + 1
+		    oid = stid
+		    call apbsky (ap, im, cl, sd, out, stid, ltid, gd, mgd, id,
+		        YES)
+		    ltid = ltid + stid - oid + 1
+		    if (id != NULL) {
+		        if (id == gd)
+			    call gflush (id)
+		        else
+		            call gframe (id)
+		    }
+		} else if (interactive == YES)
+		    call printf ("No coordinate list\n")
+
+	    default:
+		# do nothing
+		call printf ("Unknown or ambiguous keystroke command\n")
+	    }
+
+endswitch_
+
+	    # Setup for the next object.
+	    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/fitsky/apskybuf.x b/noao/digiphot/apphot/fitsky/apskybuf.x
new file mode 100644
index 00000000..cacd4737
--- /dev/null
+++ b/noao/digiphot/apphot/fitsky/apskybuf.x
@@ -0,0 +1,254 @@
+include <imhdr.h>
+include <math.h>
+include <mach.h>
+include "../lib/apphotdef.h"
+include "../lib/fitskydef.h"
+include "../lib/fitsky.h"
+
+# APSKYBUF -- Procedure to fetch the sky pixels given the pointer to the
+# IRAF image, the coordinates of the center and the size of the apphot
+# sky annulus.
+
+int procedure apskybuf (ap, im, wx, wy)
+
+pointer	ap		# pointer to apphot structure
+pointer	im		# pointer to the IRAF image
+real	wx, wy		# center coordinates
+
+int	lenbuf
+pointer	sky
+real	annulus, dannulus, datamin, datamax
+int	ap_skypix(), ap_bskypix()
+
+begin
+	# Check for 0 radius annulus.
+	sky = AP_PSKY(ap)
+	annulus = AP_ANNULUS(sky) * AP_SCALE(ap)
+	dannulus = AP_DANNULUS(sky) * AP_SCALE(ap)
+	if (dannulus <= 0.0)
+	    return (AP_NOSKYAREA)
+
+	# Allocate space for sky pixels.
+	lenbuf = PI * (2.0 * annulus + dannulus + 1.0) * (dannulus + 0.5)
+
+	if (lenbuf != AP_LENSKYBUF(sky)) {
+	    if (AP_SKYPIX(sky) != NULL)
+		call mfree (AP_SKYPIX(sky), TY_REAL)
+	    call malloc (AP_SKYPIX(sky), lenbuf, TY_REAL)
+	    if (AP_COORDS(sky) != NULL)
+		call mfree (AP_COORDS(sky), TY_INT)
+	    call malloc (AP_COORDS(sky), lenbuf, TY_INT)
+	    if (AP_INDEX(sky) != NULL)
+		call mfree (AP_INDEX(sky), TY_INT)
+	    call malloc (AP_INDEX(sky), lenbuf, TY_INT)
+	    if (AP_SWGT(sky) != NULL)
+		call mfree (AP_SWGT(sky), TY_REAL)
+	    call malloc (AP_SWGT(sky), lenbuf, TY_REAL)
+	    AP_LENSKYBUF(sky) = lenbuf
+	}
+
+	# Fetch the sky pixels.
+	if (IS_INDEFR(AP_DATAMIN(ap)) && IS_INDEFR(AP_DATAMAX(ap))) {
+	    AP_NSKYPIX(sky) = ap_skypix (im, wx, wy, annulus, (annulus +
+	        dannulus), Memr[AP_SKYPIX(sky)], Memi[AP_COORDS(sky)],
+		AP_SXC(sky), AP_SYC(sky), AP_SNX(sky), AP_SNY(sky))
+	    AP_NBADSKYPIX(sky) = 0
+	} else {
+	    if (IS_INDEFR(AP_DATAMIN(ap)))
+		datamin = -MAX_REAL
+	    else
+		datamin = AP_DATAMIN(ap)
+	    if (IS_INDEFR(AP_DATAMAX(ap)))
+		datamax = MAX_REAL
+	    else
+		datamax = AP_DATAMAX(ap)
+	    AP_NSKYPIX(sky) = ap_bskypix (im, wx, wy, annulus, (annulus +
+	        dannulus), datamin, datamax, Memr[AP_SKYPIX(sky)],
+		Memi[AP_COORDS(sky)], AP_SXC(sky), AP_SYC(sky), AP_SNX(sky),
+		AP_SNY(sky), AP_NBADSKYPIX(sky))
+	}
+
+	if (AP_NSKYPIX(sky) <= 0) {
+	    if (AP_NBADSKYPIX(sky) <= 0)
+	        return (AP_SKY_OUTOFBOUNDS)
+	    else
+	        return (AP_NSKY_TOO_SMALL)
+	} else
+	    return (AP_OK)
+end
+
+
+# AP_SKYPIX -- Procedure to fetch the sky pixels from the image
+
+int procedure ap_skypix (im, wx, wy, rin, rout, skypix, coords, xc, yc,
+        nx, ny)
+
+pointer	im		# pointer to IRAF image
+real	wx, wy		# center of sky annulus
+real	rin, rout	# inner and outer radius of sky annulus
+real	skypix[ARB]	# skypixels
+int	coords[ARB]	# sky subraster coordinates [i + nx * (j - 1)]
+real	xc, yc		# center of sky subraster
+int	nx, ny		# max dimensions of sky subraster (output)
+
+int	i, j, ncols, nlines, c1, c2, l1, l2, nskypix
+pointer	buf
+real	xc1, xc2, xl1, xl2, rin2, rout2, rj2, r2
+pointer	imgs2r()
+
+#pointer	tbuf
+
+begin
+	if (rout <= rin)
+	    return (0)
+
+	# Test for out of bounds sky regions.
+	ncols = IM_LEN(im,1)
+	nlines = IM_LEN(im,2)
+	xc1 = wx - rout
+	xc2 = wx + rout
+	xl1 = wy - rout
+	xl2 = wy + rout
+	if (xc2 < 1.0 || xc1 > real (ncols) || xl2 < 1.0 || xl1 > real (nlines))
+	    return (0)
+
+	# Compute the column and line limits.
+	c1 = max (1.0, min (real (ncols), wx - rout)) + 0.5
+	c2 = min (real (ncols), max (1.0, wx + rout)) + 0.5
+	l1 = max (1.0, min (real (nlines), wy - rout)) + 0.5
+	l2 = min (real (nlines), max (1.0, wy + rout)) + 0.5
+	nx = c2 - c1 + 1
+	ny = l2 - l1 + 1
+	xc = wx - c1 + 1
+	yc = wy - l1 + 1
+
+	# Fetch the sky pixels.
+	rin2 = rin ** 2
+	rout2 = rout ** 2
+	nskypix = 0
+
+	do j = l1, l2 {
+	    buf = imgs2r (im, c1, c2, j, j)
+	    rj2 = (wy - j) ** 2
+	    do i = c1, c2 {
+	        r2 = (wx - i) ** 2 + rj2
+		if (r2 > rin2 && r2 <= rout2) {
+		    skypix[nskypix+1] = Memr[buf+i-c1]
+		    coords[nskypix+1] = (i - c1 + 1) + nx * (j - l1)
+		    nskypix = nskypix + 1
+		}
+	    }
+	}
+
+	#buf = imgs2r (im, c1, c2, l1, l2)
+	#tbuf = buf
+	#do j = l1,  l2 {
+	    #rj2 = (wy - j) ** 2
+	    #do i = c1, c2 {
+	        #r2 = (wx - i) ** 2 + rj2
+		#if (r2 > rin2 && r2 <= rout2) {
+		    #skypix[nskypix+1] = Memr[tbuf+i-c1]
+		    #coords[nskypix+1] = (i - c1 + 1) + nx * (j - l1)
+		    #nskypix = nskypix + 1
+		#}
+	    #}
+	    #tbuf = tbuf + nx
+	#}
+
+	return (nskypix)
+end
+
+
+# AP_BSKYPIX -- Procedure to fetch the sky pixels from the image
+
+int procedure ap_bskypix (im, wx, wy, rin, rout, datamin, datamax,
+	skypix, coords, xc, yc, nx, ny, nbad)
+
+pointer	im		# pointer to IRAF image
+real	wx, wy		# center of sky annulus
+real	rin, rout	# inner and outer radius of sky annulus
+real	datamin		# minimum good value
+real	datamax		# maximum good value
+real	skypix[ARB]	# skypixels
+int	coords[ARB]	# sky subraster coordinates [i + nx * (j - 1)]
+real	xc, yc		# center of sky subraster
+int	nx, ny		# max dimensions of sky subraster (output)
+int	nbad		# number of bad pixels
+
+int	i, j, ncols, nlines, c1, c2, l1, l2, nskypix
+pointer	buf
+real	xc1, xc2, xl1, xl2, rin2, rout2, rj2, r2, pixval
+pointer	imgs2r()
+
+#pointer	tbuf
+
+begin
+	if (rout <= rin)
+	    return (0)
+
+	# Test for out of bounds sky regions.
+	ncols = IM_LEN(im,1)
+	nlines = IM_LEN(im,2)
+	xc1 = wx - rout
+	xc2 = wx + rout
+	xl1 = wy - rout
+	xl2 = wy + rout
+	if (xc2 < 1.0 || xc1 > real (ncols) || xl2 < 1.0 || xl1 > real (nlines))
+	    return (0)
+
+	# Compute the column and line limits.
+	c1 = max (1.0, min (real (ncols), wx - rout)) + 0.5
+	c2 = min (real (ncols), max (1.0, wx + rout)) + 0.5
+	l1 = max (1.0, min (real (nlines), wy - rout)) + 0.5
+	l2 = min (real (nlines), max (1.0, wy + rout)) + 0.5
+	nx = c2 - c1 + 1
+	ny = l2 - l1 + 1
+	xc = wx - c1 + 1
+	yc = wy - l1 + 1
+
+	rin2 = rin ** 2
+	rout2 = rout ** 2
+	nskypix = 0
+	nbad = 0
+
+	# Fetch the sky pixels.
+	do j = l1, l2 {
+	    buf = imgs2r (im, c1, c2, j, j)
+	    rj2 = (wy - j) ** 2
+	    do i = c1, c2 {
+	        r2 = (wx - i) ** 2 + rj2
+		if (r2 > rin2 && r2 <= rout2) {
+		    pixval = Memr[buf+i-c1] 
+		    if (pixval < datamin || pixval > datamax)
+		        nbad = nbad + 1
+		    else {
+		        skypix[nskypix+1] = pixval
+		        coords[nskypix+1] = (i - c1 + 1) + nx * (j - l1)
+		        nskypix = nskypix + 1
+		    }
+		}
+	    }
+	}
+
+	#buf = imgs2r (im, c1, c2, l1, l2)
+	#tbuf = buf
+	#do j = l1, l2 {
+	    #rj2 = (wy - j) ** 2
+	    #do i = c1, c2 {
+	        #r2 = (wx - i) ** 2 + rj2
+		#if (r2 > rin2 && r2 <= rout2) {
+		    #pixval = Memr[tbuf+i-c1] 
+		    #if (pixval < datamin || pixval > datamax)
+		        #nbad = nbad + 1
+		    #else {
+		        #skypix[nskypix+1] = pixval
+		        #coords[nskypix+1] = (i - c1 + 1) + nx * (j - l1)
+		        #nskypix = nskypix + 1
+		    #}
+		#}
+	    #}
+	    #tbuf = tbuf + nx
+	#}
+
+	return (nskypix)
+end
diff --git a/noao/digiphot/apphot/fitsky/apskycolon.x b/noao/digiphot/apphot/fitsky/apskycolon.x
new file mode 100644
index 00000000..35b59aa9
--- /dev/null
+++ b/noao/digiphot/apphot/fitsky/apskycolon.x
@@ -0,0 +1,367 @@
+include "../lib/apphot.h"
+include "../lib/noise.h"
+include "../lib/fitsky.h"
+include "../lib/display.h"
+
+
+# APSKYCOLON -- Procedure to process the fitsky colon commands.
+
+procedure apskycolon (ap, im, cl, out, stid, ltid, cmdstr, newimage,
+	newskybuf, newsky)
+
+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				# coord list sequence number
+char	cmdstr[ARB]			# command string
+int	newimage			# new image ?
+int	newskybuf			# new sky buffer ?
+int	newsky				# new sky fit ?
+
+int	junk
+pointer	sp, incmd, outcmd
+int	strdic()
+
+begin
+	# Get the command.
+	call smark (sp)
+	call salloc (incmd, SZ_LINE, TY_CHAR)
+	call salloc (outcmd, SZ_LINE, TY_CHAR)
+	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, SCMDS) != 0)
+	    call apscolon (ap, out, stid, cmdstr, newskybuf, newsky)
+	else if (strdic (Memc[incmd], Memc[outcmd], SZ_LINE, APCMDS) != 0)
+	    call ap_apcolon (ap, im, cl, out, stid, ltid, cmdstr,
+	        newimage, junk, junk, newskybuf, newsky, junk, junk)
+	else if (strdic (Memc[incmd], Memc[outcmd], SZ_LINE, NCMDS) != 0)
+	    call apnscolon (ap, im, out, stid, cmdstr, junk, junk,
+	        newskybuf, newsky, junk, junk)
+	else
+	    call ap_simcolon (ap, cmdstr)
+
+	call sfree (sp)
+end
+
+
+# APSCOLON --  Procedure to examine and edit the sky fitting parameters.
+
+procedure apscolon (ap, out, stid, cmdstr, newbuf, newfit)
+
+pointer	ap		# pointer to the apphot structure
+int	out		# output file descriptor
+int	stid		# output file number
+char	cmdstr		# command string
+int	newbuf		# new sky buffer
+int	newfit		# new sky fit
+
+bool	bval
+int	ncmd, ival, stat
+pointer	sp, cmd
+real	rval
+
+bool	itob()
+int	nscan(), strdic(), btoi(), apstati()
+real	apstatr()
+
+
+begin
+	# Get the command
+	call smark (sp)
+	call salloc (cmd, SZ_LINE, TY_CHAR)
+	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, SCMDS)
+	switch (ncmd) {
+	case SCMD_ANNULUS:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("%s = %g %s\n")
+		    call pargstr (KY_ANNULUS)
+		    call pargr (apstatr (ap, ANNULUS))
+		    call pargstr (UN_SSCALEUNIT)
+	    } else {
+		call apsetr (ap, ANNULUS, rval)
+		if (stid > 1)
+		    call ap_rparam (out, KY_ANNULUS, rval, UN_SSCALEUNIT,
+			"inner radius of sky annulus")
+		newbuf = YES
+		newfit = YES
+	    }
+	case SCMD_DANNULUS:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("%s = %g %s\n")
+		    call pargstr (KY_DANNULUS)
+		    call pargr (apstatr (ap, DANNULUS))
+		    call pargstr (UN_SSCALEUNIT)
+	    } else {
+		call apsetr (ap, DANNULUS, rval)
+		if (stid > 1)
+		    call ap_rparam (out, KY_DANNULUS, rval, UN_SSCALEUNIT,
+			"width of the sky annulus")
+		newbuf = YES
+		newfit = YES
+	    }
+	case SCMD_SALGORITHM:
+	    call gargwrd (Memc[cmd], SZ_LINE)
+	    if (Memc[cmd] == EOS) {
+		call apstats (ap, SSTRING, Memc[cmd], SZ_FNAME)
+	        call printf ("%s = %s\n")
+		    call pargstr (KY_SSTRING)
+		    call pargstr (Memc[cmd])
+	    } else {
+		stat = strdic (Memc[cmd], Memc[cmd], SZ_LINE, SFUNCS)
+		if (stat > 0) {
+		    call apseti (ap, SKYFUNCTION, stat)
+		    call apsets (ap, SSTRING, Memc[cmd])
+		    if (stid > 1)
+		        call ap_sparam (out, KY_SSTRING, Memc[cmd],
+			    UN_SALGORITHM, "sky fitting algorithm")
+		    newfit = YES
+		}
+	    }
+	case SCMD_KHIST:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("%s = %g %s\n")
+		    call pargstr (KY_K1)
+		    call pargr (apstatr (ap, K1))
+		    call pargstr (UN_SSIGMA)
+	    } else {
+		call apsetr (ap, K1, rval)
+		if (stid > 1)
+		    call ap_rparam (out, KY_K1, rval, UN_SSIGMA,
+			"half width of sky histogram")
+		newfit = YES
+	    }
+	case SCMD_SLOREJECT:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("%s = %g %s\n")
+		    call pargstr (KY_SLOREJECT)
+		    call pargr (apstatr (ap, SLOREJECT))
+		    call pargstr (UN_SSIGMA)
+	    } else {
+		call apsetr (ap, SLOREJECT, rval)
+		if (stid > 1)
+		    call ap_rparam (out, KY_SLOREJECT, rval, UN_SSIGMA,
+			"lower k-sigma rejection criterion")
+		newfit = YES
+	    }
+	case SCMD_SHIREJECT:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("%s = %g %s\n")
+		    call pargstr (KY_SHIREJECT)
+		    call pargr (apstatr (ap, SHIREJECT))
+		    call pargstr (UN_SSIGMA)
+	    } else {
+		call apsetr (ap, SHIREJECT, rval)
+		if (stid > 1)
+		    call ap_rparam (out, KY_SHIREJECT, rval, UN_SSIGMA,
+			"upper k-sigma rejection criterion")
+		newfit = YES
+	    }
+	case SCMD_SLOCLIP:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("%s = %g %s\n")
+		    call pargstr (KY_SLOCLIP)
+		    call pargr (apstatr (ap, SLOCLIP))
+		    call pargstr (UN_SPERCENT)
+	    } else {
+		call apsetr (ap, SLOCLIP, rval)
+		if (stid > 1)
+		    call ap_rparam (out, KY_SLOCLIP, rval, UN_SPERCENT,
+			"lower k-sigma rejection criterion")
+		newfit = YES
+	    }
+	case SCMD_SHICLIP:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("%s = %g %s\n")
+		    call pargstr (KY_SHICLIP)
+		    call pargr (apstatr (ap, SHICLIP))
+		    call pargstr (UN_SPERCENT)
+	    } else {
+		call apsetr (ap, SHICLIP, rval)
+		if (stid > 1)
+		    call ap_rparam (out, KY_SHICLIP, rval, UN_SPERCENT,
+			"lower k-sigma rejection criterion")
+		newfit = YES
+	    }
+	case SCMD_SMAXITER:
+	    call gargi (ival)
+	    if (nscan () == 1) {
+		call printf ("%s = %d\n")
+		    call pargstr (KY_SMAXITER)
+		    call pargi (apstati (ap, SMAXITER))
+	    } else {
+		call apseti (ap, SMAXITER, ival)
+		if (stid > 1)
+		    call ap_iparam (out, KY_SMAXITER, ival, UN_SNUMBER,
+			"maximum number of iterations")
+		newfit = YES
+	    }
+	case SCMD_BINSIZE:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("%s = %g %s\n")
+		    call pargstr (KY_BINSIZE)
+		    call pargr (apstatr (ap, BINSIZE))
+		    call pargstr (UN_SSIGMA)
+	    } else {
+		call apsetr (ap, BINSIZE, rval)
+		if (stid > 1)
+		    call ap_rparam (out, KY_BINSIZE, rval, UN_SSIGMA,
+			"width of the sky histogram bin")
+		newfit = YES
+	    }
+	case SCMD_SMOOTH:
+	    call gargb (bval)
+	    if (nscan () == 1) {
+		call printf ("%s = %b\n")
+		    call pargstr (KY_SMOOTH)
+		    call pargb (itob (apstati (ap, SMOOTH)))
+	    } else {
+		call apseti (ap, SMOOTH, btoi (bval))
+		if (stid > 1)
+		    call ap_bparam (out, KY_SMOOTH, bval, UN_SSWITCH,
+			"Lucy smooth the histogram")
+		newfit = YES
+	    }
+	case SCMD_RGROW:
+	    call gargr (rval)
+	    if (nscan() == 1) {
+		call printf ("%s = %g %s\n")
+		    call pargstr (KY_RGROW)
+		    call pargr (apstatr (ap, RGROW))
+		    call pargstr (UN_SSCALEUNIT)
+	    } else {
+		call apsetr (ap, RGROW, rval)
+		if (stid > 1)
+		    call ap_rparam (out, KY_RGROW, rval, UN_SSCALEUNIT,
+			"region growing radius")
+		newfit = YES
+	    }
+	case SCMD_SNREJECT:
+	    call gargi (ival)
+	    if (nscan () == 1) {
+		call printf ("%s = %d\n")
+		    call pargstr (KY_SNREJECT)
+		    call pargi (apstati (ap, SNREJECT))
+	    } else {
+		call apseti (ap, SNREJECT, ival)
+		if (stid > 1)
+		    call ap_iparam (out, KY_SNREJECT, ival, UN_SNUMBER,
+			"maximum number of rejection cycles")
+		newfit = YES
+	    }
+	case SCMD_SKYVALUE:
+	    call gargr (rval)
+	    if (nscan () == 1) {
+		call printf ("%s = %g\n")
+		    call pargstr (KY_SKY_BACKGROUND)
+		    call pargr (apstatr (ap, SKY_BACKGROUND))
+		    call pargstr (UN_SCOUNTS)
+	    } else {
+		call apsetr (ap, SKY_BACKGROUND, rval)
+		if (stid > 1)
+		    call ap_rparam (out, KY_SKY_BACKGROUND, rval,
+		        UN_SCOUNTS, "user supplied sky value")
+		newfit = YES
+	    }
+	case SCMD_MKSKY:
+	    call gargb (bval)
+	    if (nscan () == 1) {
+		call printf ("%s = %b\n")
+		    call pargstr (KY_MKSKY)
+		    call pargb (itob (apstati (ap, MKSKY)))
+	    } else {
+		call apseti (ap, MKSKY, btoi (bval))
+	    }
+	default:
+	    # do nothing gracefully
+	    call printf ("Unrecognized command\7\n")
+	}
+
+	call sfree (sp)
+end
+
+
+# AP_SIMCOLON -- Procedure to process fitsky commands which alter parameters
+# other than the sky fitting parameters themselves.
+
+procedure ap_simcolon (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 commands.
+	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, SSHOWARGS)
+	    switch (ncmd) {
+	    case SCMD_DATA:
+		call printf ("\n")
+		call ap_nshow (ap)
+		call printf ("\n")
+	    case SCMD_SKY:
+		call printf ("\n")
+	        call ap_spshow (ap)
+		call printf ("\n")
+	    default:
+		call printf ("\n")
+	        call ap_sshow (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/fitsky/apsplot.x b/noao/digiphot/apphot/fitsky/apsplot.x
new file mode 100644
index 00000000..0a75b0c8
--- /dev/null
+++ b/noao/digiphot/apphot/fitsky/apsplot.x
@@ -0,0 +1,244 @@
+include <pkg/gtools.h>
+include <gset.h>
+include "../lib/apphotdef.h"
+include "../lib/apphot.h"
+include "../lib/noise.h"
+include "../lib/fitskydef.h"
+include "../lib/fitsky.h"
+
+# AP_SPLOT -- Procedure to compute radial profile plots for the sky fitting
+# routine.
+
+procedure ap_splot (ap, sid, gd, makeplot)
+
+pointer	ap		# pointer to the apphot structure
+int	sid		# id number of the star
+pointer	gd		# graphics stream
+int	makeplot	# make a plot
+
+int	nx, ny, nskypix
+pointer	sp, sky, str, r, gt
+real	xcenter, ycenter, rmin, rmax, imin, imax
+real	u1, u2, v1, v2, x1, x2, y1, y2
+int	apstati()
+pointer	ap_gtinit()
+real	apstatr()
+
+begin
+	# Initialize.
+	if (gd == NULL || makeplot == NO)
+	    return
+
+	# Check for defined center and non-constant algorithm.
+	xcenter = apstatr (ap, SXCUR)
+	ycenter = apstatr (ap, SYCUR)
+	if (IS_INDEFR(xcenter) || IS_INDEFR(ycenter))
+	    return
+	if (apstati (ap, SKYFUNCTION) == AP_CONSTANT)
+	    return
+
+	# Check that a buffer of sky pixels exists.
+	sky = AP_PSKY(ap)
+	nskypix = AP_NSKYPIX(sky)
+	nx = AP_SNX(sky)
+	ny = AP_SNY(sky)
+	if (nskypix <= 0 || nx <= 0 || ny <= 0)
+	    return
+
+	# Allocate working space
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+	call salloc (r, nskypix, TY_REAL)
+
+	# Compute the radii and the plot limits.
+	call ap_xytor (Memi[AP_COORDS(sky)], Memi[AP_INDEX(sky)],
+	    Memr[r], nskypix, AP_SXC(sky), AP_SYC(sky), nx)
+	call alimr (Memr[r], nskypix, rmin, rmax)
+	rmin = rmin - 1.0
+	rmax = rmax + 1.0
+	call alimr (Memr[AP_SKYPIX(sky)], nskypix, imin, imax)
+
+	# Reactivate the work station.
+	call greactivate (gd, 0)
+
+	# Store the 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,OSXCUR), apstatr(ap,OSYCUR))
+	
+	# Draw the plot.
+	call gclear (gd)
+	call ap_spset (gd, gt, ap, rmin, rmax, imin, imax)
+	call ap_spannotate (gd, ap, rmin, rmax, imin, imax)
+	call ap_plotrad (gd, gt, Memr[r], Memr[AP_SKYPIX(sky)], nskypix, "plus")
+
+	# 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_SPSET -- Procedure to set up the parameters for the fitsky radial profile
+# plot.
+
+procedure ap_spset (gd, gt, ap, xmin, xmax, ymin, ymax)
+
+pointer	gd		# graphics stream
+pointer	gt		# gtools pointer
+pointer	ap		# apphot pointer
+real	xmin, xmax	# minimum and maximum radial distance
+real	ymin, ymax	# minimum and maximum of the y axis
+
+int	fd
+pointer	sp, str, title
+real	aspect, scale, vx1, vx2, vy1, vy2
+int	apstati(), stropen()
+real	apstatr(), gstatr()
+
+begin
+	call smark (sp)
+	call salloc (str, 3 * SZ_LINE, TY_CHAR)
+	call salloc (title, SZ_LINE, TY_CHAR)
+
+	# Encode the parameter string.
+	fd = stropen (Memc[str], 3 * SZ_LINE, WRITE_ONLY)
+
+	call sysid (Memc[title], SZ_FNAME)
+	call fprintf (fd, "%s\n")
+	    call pargstr (Memc[title])
+
+	call fprintf (fd,
+	    "Sky: value=%0.2f sigma=%0.2f skew=%0.2f nsky=%d nreject=%d\n")
+	    call pargr (apstatr (ap, SKY_MODE))
+	    call pargr (apstatr (ap, SKY_SIGMA))
+	    call pargr (apstatr (ap, SKY_SKEW))
+	    call pargi (apstati (ap, NSKY))
+	    call pargi (apstati (ap, NSKY_REJECT))
+
+	call gt_gets (gt, GTTITLE, Memc[title], SZ_LINE)
+	call fprintf (fd, "%s\n")
+	    call pargstr (Memc[title])
+
+	call strclose (fd)
+
+	# Set the aspect ratio.
+	scale = apstatr (ap, SCALE)
+	aspect = gstatr (gd, G_ASPECT)
+	call gsetr (gd, G_ASPECT, 0.75)
+
+	# Set the labels and window.
+	call gseti (gd, G_XDRAWAXES, 2)
+	call gswind (gd, xmin / scale, xmax / scale, ymin, 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 gsview (gd, vx1, vx2, vy1, vy2)
+	call gswind (gd, xmin, xmax, ymin, ymax)
+	call glabax (gd, "",
+	    "Radial Distance (lower-pixels, upper-scale units)", "")
+
+	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 sfree (sp)
+end
+
+
+# AP_SPANNOTATE -- Procedure to annotate the radial plot in fitsky.
+
+procedure ap_spannotate (gd, ap, xmin, xmax, ymin, ymax)
+
+pointer	gd		# graphics stream
+pointer	ap		# apphot structure
+real	xmin, xmax	# min and max of x axis
+real	ymin, ymax	# min and max of y axis
+
+pointer	sp, str
+real	annulus, dannulus, sigma, skyval, skysigma
+real	apstatr ()
+
+begin
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+	call gseti (gd, G_PLTYPE, GL_DASHED)
+
+	# Mark the inner sky annulus.
+	annulus = apstatr (ap, SCALE) * apstatr (ap, ANNULUS)
+	if (annulus >= xmin && annulus <= xmax) {
+	    call gamove (gd, annulus, ymin)
+	    call gadraw (gd, annulus, ymax)
+	    call sprintf (Memc[str], SZ_LINE, "inner sky radius = %0.2f")
+	        call pargr (annulus)
+	    call gtext (gd, annulus, ymax, Memc[str], "q=h;u=180;v=t;p=r")
+	}
+ 
+	# Mark the outer sky annulus.
+	dannulus = annulus + apstatr (ap, SCALE) * apstatr (ap, DANNULUS)
+	if (dannulus >= xmin && dannulus <= xmax) {
+	    call gamove (gd, dannulus, ymin)
+	    call gadraw (gd, dannulus, ymax)
+	    call sprintf (Memc[str], SZ_LINE, "outer sky radius = %0.2f")
+	        call pargr (dannulus)
+	    call gtext (gd, dannulus, ymax, Memc[str], "q=h;u=180;v=t;p=r")
+	}
+
+	# Mark the sky sigma if defined.
+	sigma = apstatr (ap, SKY_SIGMA)
+	if (! IS_INDEFR(sigma)) {
+	    call gmark (gd, (xmin + xmax) / 2.0, (ymin + ymax) / 2.0, 
+		GM_VEBAR, -0.25, -sigma)
+	    call sprintf (Memc[str], SZ_LINE, "sigma = %g")
+		call pargr (sigma)
+	    call gtext (gd, (xmin + xmax) / 2.0, (ymin + ymax + sigma) / 2.0,
+		Memc[str], "q=h;h=c")
+	}
+
+	# Mark the sky value.
+	call gseti (gd, G_PLTYPE, GL_SOLID)
+	skyval = apstatr (ap, SKY_MODE)
+	if (skyval >= ymin && skyval <= ymax) {
+	    call gamove (gd, xmin, skyval)
+	    call gadraw (gd, xmax, skyval)
+	}
+
+	# Mark the upper sky sigma.
+	call gseti (gd, G_PLTYPE, GL_DASHED)
+	if (! IS_INDEFR(apstatr (ap, SKY_SIGMA)))
+	    skysigma = skyval + apstatr (ap, SHIREJECT) * apstatr (ap,
+	        SKY_SIGMA)
+	else
+	    skysigma = INDEFR
+	if (! IS_INDEFR(skysigma) && (skysigma >= ymin) && skysigma <= ymax) {
+	    call gamove (gd, xmin, skysigma)
+	    call gadraw (gd, xmax, skysigma)
+	    #call sprintf (Memc[str], SZ_LINE, "sky sigma= %g")
+	        #call pargr (skysigma)
+	}
+
+	# Mark the lower sky sigma
+	if (! IS_INDEFR(apstatr (ap, SKY_SIGMA)))
+	    skysigma = skyval - apstatr (ap, SLOREJECT) * apstatr (ap,
+	        SKY_SIGMA)
+	else
+	    skysigma = INDEFR
+	if (! IS_INDEFR(skysigma) && (skysigma >= ymin) && skysigma <= ymax) {
+	    call gamove (gd, xmin, skysigma)
+	    call gadraw (gd, xmax, skysigma)
+	}
+
+	call sfree (sp)
+end
diff --git a/noao/digiphot/apphot/fitsky/apsradsetup.x b/noao/digiphot/apphot/fitsky/apsradsetup.x
new file mode 100644
index 00000000..fd4039f8
--- /dev/null
+++ b/noao/digiphot/apphot/fitsky/apsradsetup.x
@@ -0,0 +1,97 @@
+include "../lib/display.h"
+include "../lib/fitsky.h"
+
+define	HELPFILE	"apphot$fitsky/ifitsky.key"
+
+# AP_SRADSETUP -- Procedure to set up the sky fitting interactively using
+# radial profile plot around the given coordinates.
+
+procedure ap_sradsetup (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	ier, key, wcs
+pointer	sp, cmd
+real	xcenter, ycenter, xc, yc, rmin, rmax, imin, imax, rval
+real	u1, u2, v1, v2, x1, x2, y1, y2
+
+int	apstati(), apfitsky(), clgcur(), ap_showplot()
+real	apstatr(), ap_cannulus(), ap_cdannulus(), ap_csigma()
+real	ap_cdatamin(), ap_cdatamax(), ap_crgrow()
+
+begin
+	if (gd == NULL)
+	    return
+	call greactivate (gd, 0)
+
+	# Store the old window and viewport coordinates.
+	call ggview (gd, u1, u2, v1, v2)
+	call ggwind (gd, x1, x2, y1, y2)
+
+	# Plot the profile.
+	if (ap_showplot (ap, im, wx, wy, gd, xcenter, ycenter, rmin, rmax,
+    	    imin, imax) == ERR) {
+	    call gdeactivate (gd, 0)
+	    return
+	}
+
+	# Allocate memory.
+	call smark (sp)
+	call salloc (cmd, SZ_LINE, TY_CHAR)
+
+	call printf (
+	    "Waiting for setup menu command (?=help, v=default setup, q=quit):")
+	while (clgcur ("gcommands", xc, yc, wcs, key, Memc[cmd],
+	    SZ_LINE) != EOF) {
+
+	switch (key) {
+
+	    case 'q':
+	        break
+	    case '?':
+		call gpagefile (gd, HELPFILE, "")
+	    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 'a':
+		rval = ap_cannulus (ap, gd, out, stid, rmin, rmax, imin, imax)
+	    case 'd':
+		rval = ap_cdannulus (ap, gd, out, stid, apstatr (ap, ANNULUS),
+		    rmin, rmax, imin, imax)
+	    case 'g':
+		rval = ap_crgrow (ap, gd, out, stid, rmin, rmax, imin, imax)
+	    case 'v':
+		rval = ap_cannulus (ap, gd, out, stid, rmin, rmax, imin, imax)
+		rval = ap_cdannulus (ap, gd, out, stid, apstatr (ap, ANNULUS),
+		    rmin, rmax, imin, imax)
+	        rval = ap_csigma (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):")
+	}
+	call printf (
+	    "Interactive setup is complete. Type w to store parameters.\n")
+
+	# Store the old window and viewport coordinates.
+	call gsview (gd, u1, u2, v1, v2)
+	call gswind (gd, x1, x2, y1, y2)
+
+	# Free the plotting space.
+	call sfree (sp)
+	call gdeactivate (gd, 0)
+
+	# Fit the new sky value and print it on the standard output.
+	ier = apfitsky (ap, im, xcenter, ycenter, NULL, gd)
+	call ap_splot (ap, 0, gd, apstati (ap, RADPLOTS))
+	call ap_qspsky (ap, ier)
+end
diff --git a/noao/digiphot/apphot/fitsky/apsshow.x b/noao/digiphot/apphot/fitsky/apsshow.x
new file mode 100644
index 00000000..da54bfbf
--- /dev/null
+++ b/noao/digiphot/apphot/fitsky/apsshow.x
@@ -0,0 +1,99 @@
+include "../lib/display.h"
+include "../lib/fitsky.h"
+
+
+# AP_SSHOW -- Procedure to print sky fitting parameters on the terminal.
+
+procedure ap_sshow (ap)
+
+pointer	ap	# pointer to the apphot strucuture
+
+bool	itob()
+int	apstati()
+
+begin
+	call ap_nshow (ap)
+	call printf ("\n")
+	call ap_spshow (ap)
+	call printf ("    %s = %b\n")
+	    call pargstr (KY_RADPLOTS)
+	    call pargb (itob (apstati (ap, RADPLOTS)))
+end
+
+
+# AP_SPSHOW -- Procedure to print sky fitting parameters on the terminal.
+
+procedure ap_spshow (ap)
+
+pointer	ap	# pointer to the apphot strucuture
+
+pointer	sp, str
+bool	itob()
+int	apstati()
+real	apstatr()
+
+begin
+	# Print the image characteristics
+	call smark (sp)
+	call salloc (str, SZ_LINE, TY_CHAR)
+
+	# Print the sky fitting parameters.
+	call printf ("Sky Fitting Parameters\n")
+	call apstats (ap, SSTRING, Memc[str], SZ_FNAME)
+	call printf ("    %s = %s %s    %s = %g %s\n")
+	    call pargstr (KY_SSTRING)
+	    call pargstr (Memc[str])
+	    call pargstr (UN_SALGORITHM)
+	    call pargstr (KY_SKY_BACKGROUND)
+	    call pargr (apstatr (ap, SKY_BACKGROUND))
+	    call pargstr (UN_SCOUNTS)
+
+	call printf ("    %s = %g %s    %s = %g %s\n")
+	    call pargstr (KY_ANNULUS)
+	    call pargr (apstatr (ap, ANNULUS))
+	    call pargstr (UN_SSCALEUNIT)
+	    call pargstr (KY_DANNULUS)
+	    call pargr (apstatr (ap, DANNULUS))
+	    call pargstr (UN_SSCALEUNIT)
+
+	call printf ("    %s = %g %s    %s = %g %s    %s = %b\n")
+	    call pargstr (KY_K1)
+	    call pargr (apstatr (ap, K1))
+	    call pargstr (UN_SSIGMA)
+	    call pargstr (KY_BINSIZE)
+	    call pargr (apstatr (ap, BINSIZE))
+	    call pargstr (UN_SSIGMA)
+	    call pargstr (KY_SMOOTH)
+	    call pargb (itob (apstati (ap, SMOOTH)))
+
+	call printf ("    %s = %g %s    %s = %g %s\n")
+	    call pargstr (KY_SLOCLIP)
+	    call pargr (apstatr (ap, SLOCLIP))
+	    call pargstr (UN_SPERCENT)
+	    call pargstr (KY_SHICLIP)
+	    call pargr (apstatr (ap, SHICLIP))
+	    call pargstr (UN_SPERCENT)
+
+	call printf ("    %s = %g %s    %s = %g %s    %s = %d\n")
+	    call pargstr (KY_SLOREJECT)
+	    call pargr (apstatr (ap, SLOREJECT))
+	    call pargstr (UN_SSIGMA)
+	    call pargstr (KY_SHIREJECT)
+	    call pargr (apstatr (ap, SHIREJECT))
+	    call pargstr (UN_SSIGMA)
+	    call pargstr (KY_SMAXITER)
+	    call pargi (apstati (ap, SMAXITER))
+
+	call printf ("    %s = %d    %s = %g %s\n")
+	    call pargstr (KY_SNREJECT)
+	    call pargi (apstati (ap, SNREJECT))
+	    call pargstr (KY_RGROW)
+	    call pargr (apstatr (ap, RGROW))
+	    call pargstr (UN_SSCALEUNIT)
+
+	call printf ("    %s = %b\n")
+	    call pargstr (KY_MKSKY)
+	    call pargb (itob (apstati (ap, MKSKY)))
+
+	call sfree (sp)
+end
diff --git a/noao/digiphot/apphot/fitsky/fitsky.key b/noao/digiphot/apphot/fitsky/fitsky.key
new file mode 100644
index 00000000..ee515bf1
--- /dev/null
+++ b/noao/digiphot/apphot/fitsky/fitsky.key
@@ -0,0 +1,81 @@
+	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 sky for current star
+spbar	Fit sky for current star, output results
+m	Move to next star in coordinate list
+m	Fit sky for next star in coordinate list, output results
+l	Fit sky for remaining stars in coordinate list, output results
+e	Print error messages
+r	Rewind the coordinate list
+q	Exit task
+
+
+        Colon commands
+
+:show	[data/sky]	List the parameters
+:m [n]	Move to the next [nth] star in coordinate list
+:n [n]	Fit sky to 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 PSF (scale units)
+:emission	[y/n]		Emission feature (y), absorption (n)
+:sigma		[value]		Standard deviation of sky (counts)
+:datamin	[value]		Minimum good pixel value (counts)
+:datamax	[value]		Maximum good pixel 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)
+
+# Sky fitting algorithm parameters
+
+:salgorithm	[string]	Sky fitting algorithm 
+:skyvalue	[value]		User supplied sky value (counts)
+:annulus 	[value]		Inner radius of sky annulus (scale units)
+:dannulus	[value]		Width of sky annulus (scale units)
+:khist		[value]		Sky histogram extent (+/- sky sigma)
+:smooth		[y/n]		Lucy smooth the sky histogram 
+:binsize	[value]		Resolution of sky histogram (sky sigma)
+:smaxiter	[value]		Maximum number of iterations
+:sloclip	[value]		Low side clipping factor (percent)
+:shiclip	[value]		High side clipping factor (percent)
+:snreject	[value]		Maximum number of rejection cycles
+:sloreject	[value]		Low side pixel rejection limits (sky sigma)
+:shireject	[value]		High side pixel rejection limits (sky sigma)
+:rgrow		[value]		Region growing radius (scale units)
+
+# Marking and plotting parameters
+
+:mksky		[y/n]		Mark sky annuli on the display
+:radplot	[y/n]		Plot radial profile of sky pixels
diff --git a/noao/digiphot/apphot/fitsky/ifitsky.key b/noao/digiphot/apphot/fitsky/ifitsky.key
new file mode 100644
index 00000000..bd651bc3
--- /dev/null
+++ b/noao/digiphot/apphot/fitsky/ifitsky.key
@@ -0,0 +1,11 @@
+                      Interactive Fitsky Setup Menu
+
+	v	Mark and verify the critical parameters (a,d,s)
+
+	s	Mark and verify the standard deviation of the sky
+	l	Mark and verify the minimum good data value
+	u	Mark and verify the maximum good data value
+
+	a	Mark and verify the inner radius of the sky annulus
+	d	Mark and verify the width of the sky annulus
+	g	Mark and verify the region growing radius
diff --git a/noao/digiphot/apphot/fitsky/mkpkg b/noao/digiphot/apphot/fitsky/mkpkg
new file mode 100644
index 00000000..a72edc57
--- /dev/null
+++ b/noao/digiphot/apphot/fitsky/mkpkg
@@ -0,0 +1,59 @@
+# FITSKY Task Tools
+
+$checkout libpkg.a ".."
+$update libpkg.a
+$checkin libpkg.a ".."
+$exit
+
+libpkg.a:
+	apavsky.x	../lib/display.h	../lib/fitsky.h
+	apbsky.x	<fset.h>		../lib/apphot.h         \
+			../lib/display.h
+	apcentroid.x	<mach.h>		../lib/fitsky.h
+	apcrosscor.x	<mach.h>		../lib/fitsky.h
+	apfitsky.x	../lib/apphotdef.h	../lib/noisedef.h	\
+			../lib/fitskydef.h	../lib/fitsky.h         \
+			../lib/apphot.h
+	apgauss.x	<mach.h>		../lib/fitsky.h         \
+			<math/nlfit.h>
+	apgrowhist.x
+	apgspars.x	../lib/display.h	../lib/fitsky.h         \
+			../lib/noise.h
+	aphgmsub.x
+	aphistplot.x	<gset.h>		<pkg/gtools.h>		\
+			../lib/fitsky.h
+	aplgsky.x	<mach.h>		../lib/fitsky.h
+	apmean.x	<mach.h>		../lib/fitsky.h
+	apmedian.x	<mach.h>		../lib/fitsky.h
+	apmode.x	<mach.h>		../lib/fitsky.h
+	appsky.x	../lib/fitsky.h		../lib/apphot.h
+	apradplot.x	../lib/fitsky.h
+	apreadsky.x
+	aprefitsky.x	../lib/apphotdef.h	../lib/noisedef.h	\
+			../lib/fitskydef.h	../lib/fitsky.h         \
+			../lib/apphot.h
+	aprgrow.x
+	apsconfirm.x	../lib/apphot.h		../lib/noise.h          \
+			../lib/fitsky.h
+	apserrors.x	../lib/fitsky.h
+	apsfree.x	../lib/apphotdef.h	../lib/fitskydef.h
+	apsinit.x	../lib/apphotdef.h	../lib/fitskydef.h      \
+			../lib/fitsky.h
+	apsky.x		<ctype.h>		<gset.h>		\
+			../lib/apphot.h		../lib/display.h        \
+			../lib/fitsky.h         <imhdr.h>
+	apskybuf.x	<imhdr.h>		<math.h>                \
+			../lib/apphotdef.h	../lib/fitskydef.h	\
+			../lib/fitsky.h         <mach.h>
+	apskycolon.x	../lib/apphot.h		../lib/fitsky.h         \
+			../lib/display.h	../lib/noise.h
+	apsradsetup.x	../lib/display.h	../lib/fitsky.h
+	appspars.x	../lib/display.h
+	apsplot.x	<pkg/gtools.h>		<gset.h>		\
+			../lib/apphot.h		../lib/noise.h          \
+			../lib/fitsky.h		../lib/fitskydef.h	\
+			../lib/apphotdef.h
+	apsshow.x	../lib/display.h	../lib/fitsky.h
+	t_fitsky.x	<fset.h>		<gset.h>		\
+			../lib/apphot.h		<imhdr.h>
+	;
diff --git a/noao/digiphot/apphot/fitsky/t_fitsky.x b/noao/digiphot/apphot/fitsky/t_fitsky.x
new file mode 100644
index 00000000..c9fb7db8
--- /dev/null
+++ b/noao/digiphot/apphot/fitsky/t_fitsky.x
@@ -0,0 +1,303 @@
+include <gset.h>
+include <fset.h>
+include <imhdr.h>
+include "../lib/apphot.h"
+
+# T_FITSKY -- Procedure to fit sky values in a an annular region around
+# a list of objects.
+
+procedure t_fitsky ()
+
+pointer image		# pointer to the name of the image
+pointer	output		# pointer to output file name
+pointer	coords		# coordinate file
+pointer	plotfile	# pointer to file of graphics metacode
+pointer graphics	# pointer to graphics display device
+pointer	display		# pointer to display device
+int	interactive	# mode of use
+int	cache		# cache the image pixels in memory
+int	verify		# verify critical parameters
+int	update		# update the critical parameter
+int	verbose		# verbose mode
+
+pointer	sp, outfname, cname, ap, im, mgd, gd, id, str
+int	sid, lid, limlist, lclist, lolist, out, cl, pfd, 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	strlen(), apsky(), imtopenp(), clpopnu(), open(), clgwrd(), ap_memstat()
+int	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 the standard output to flush on a 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, "Imcompatable 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"))
+	verify = btoi (clgetb ("verify"))
+	update = btoi (clgetb ("update"))
+	verbose = btoi (clgetb ("verbose"))
+
+	# Get the parameters.
+	call ap_sgpars (ap)
+	if (verify == YES && interactive == NO) {
+	    call ap_sconfirm (ap, NULL, 1)
+	    if (update == YES)
+		call ap_pspars (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 plot files.
+	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 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 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 the output text file, if output is "default", dir$default
+	    # or a directory specification then the extension "sky" is added
+	    # to the image name and a suitable version number is appended to
+	    # the output name. If the output string is null to 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], "sky",
+		        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 sky
+	    if (interactive == NO) {
+	        if (Memc[cname] != EOS)
+		    stat =  apsky (ap, im, cl, NULL, NULL, mgd, NULL, out,
+		        sid, NO, cache)
+	        else if (cl != NULL) {
+		    lid = 1
+	            call apbsky (ap, im, cl, NULL, out, sid, lid, gd, mgd, id,
+		        verbose)
+		    stat = NO
+		} else
+		    stat = NO
+	    } else
+		stat = apsky (ap, im, cl, NULL, gd, mgd, id, out, sid, YES,
+		    cache)
+
+	    # Cleanup.
+	    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.
+	if (cl != NULL && lclist == 1)
+	    call close (cl)
+
+	# If only one output file for a list of images close.
+	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 plot files.
+	if (id == gd && id != NULL) {
+	    call gclose (id)
+	} 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 sky fitting structure
+	call apsfree (ap)
+
+	# Close up lists.
+	call imtclose (imlist)
+	call clpcls (clist)
+	call clpcls (olist)
+
+	call sfree (sp)
+end