Repatch (from linux) of OSX IRAF

1 files changed, 154 insertions, 0 deletions

diff --git a/noao/digiphot/apphot/aputil/apmeds.x b/noao/digiphot/apphot/aputil/apmeds.x
new file mode 100644
index 00000000..0a206c8a
--- /dev/null
+++ b/noao/digiphot/apphot/aputil/apmeds.x
@@ -0,0 +1,154 @@
+# APSMED -- Compute the averaged median given a sorted array and an averaging
+# half width.
+
+real procedure apsmed (pix, index, npts, medcut)
+
+real	pix[ARB]	# array of sky pixels
+int	index[ARB]	# sorted index array
+int	npts		# number of pixels
+int	medcut		# averaging half width
+
+int	med, j, nmed, medlo, medhi
+real	sumed
+
+begin
+	med = (npts + 1) / 2
+	if (mod (med, 2) == 1) {
+	    medlo = max (1, med - medcut)
+	    medhi = min (npts, med + medcut)
+	} else {
+	    medlo = max (1, med - medcut)
+	    medhi = min (npts, med + medcut + 1)
+	}
+
+	sumed = 0.0
+	nmed = 0
+	do j = medlo, medhi {
+	    sumed = sumed + pix[index[j]]
+	    nmed = nmed + 1
+	}
+
+	return (sumed / nmed)
+end
+
+
+# APIMED -- Compute the index of new median value. Weight is an arbitrary
+# weight array which is assumed to be zero if the pixels has been rejected
+# and is positive otherwise.
+
+int procedure apimed (weight, index, lo, hi, nmed)
+
+real	weight[ARB]		# array of weights
+int	index[ARB]		# array of sorted indices
+int	lo, hi			# ranges of weights
+int	nmed			# number of good sky pixels
+
+int	npts, med
+
+begin
+	npts = 0
+	for (med = lo; med <= hi && npts < nmed; med = med + 1) {
+	    if (weight[index[med]] > 0.0)
+		npts = npts + 1
+	}
+
+	if (npts == 0)
+	    return (0)
+	else
+	    return (med)
+end
+
+
+# APWSMED -- Compute the new averaged median given a sorted input array,
+# an averaging half-width, and assuming that there has been pixel rejection.
+
+real procedure apwsmed (pix, index, weight, npix, med, medcut)
+
+real	pix[ARB]		# pixel values array
+int	index[ARB]		# sorted indices array
+real	weight[ARB]		# the weights array
+int	npix			# number of pixels
+int	med			# index of median value
+int	medcut			# of median cut
+
+int	j, nmed, maxmed
+real	sumed
+
+begin
+	sumed = pix[index[med]]
+	if (mod (med, 2) == 1)
+	    maxmed = 2 * medcut + 1 
+	else
+	    maxmed = 2 * medcut + 2 
+
+	nmed = 1
+	for (j = med - 1; j >= 1; j = j - 1) {
+	    if (nmed >=  medcut + 1)
+		break
+	    if (weight[index[j]] > 0.0) {
+		sumed = sumed + pix[index[j]]
+		nmed = nmed + 1
+	    }
+	}
+	for (j = med + 1; j <= npix; j = j + 1) {
+	    if (nmed >= maxmed)
+		break
+	    if (weight[index[j]] > 0.0) {
+		sumed = sumed + pix[index[j]]
+		nmed = nmed + 1
+	    }
+	}
+
+	return (sumed / nmed)
+end
+
+
+# APMEDR -- Vector median selection. The selection is carried out in a temporary
+# array, leaving the input vector unmodified.  Especially demanding applications
+# may wish to call the asok routine directory to avoid the call to the memory
+# allocator.
+
+real procedure apmedr (a, index, npix)
+
+real	a[ARB]				# input array of values
+int	index[ARB]			# sorted index array
+int	npix				# number of pixels
+
+int	i
+pointer	sp, aa
+real	median
+real	asokr()		# select the Kth smallest element from A
+
+begin
+	switch (npix) {
+	case 1, 2:
+	    return (a[1])
+
+	case 3:
+	    if (a[1] < a[2]) {
+		if (a[2] < a[3])
+		    return (a[2])
+		else if (a[1] < a[3])
+		    return (a[3])
+		else
+		    return (a[1])
+	    } else {
+		if (a[2] > a[3])
+		    return (a[2])
+		else if (a[1] < a[3])
+		    return (a[1])
+		else
+		    return (a[3])
+	    }
+
+	default:
+	    call smark (sp)
+	    call salloc (aa, npix, TY_REAL)
+	    do i = 1, npix
+		Memr[aa+i-1] = a[index[i]]
+	    median = asokr (Memr[aa], npix, (npix + 1) / 2)
+	    call sfree (sp)
+
+	    return (median)
+	}
+end