1 files changed, 96 insertions, 0 deletions
diff --git a/noao/digiphot/apphot/radprof/aprmmeasure.x b/noao/digiphot/apphot/radprof/aprmmeasure.x
new file mode 100644
index 00000000..fe455008
--- /dev/null
+++ b/noao/digiphot/apphot/radprof/aprmmeasure.x
@@ -0,0 +1,96 @@
+include <math/curfit.h>
+
+# AP_RMMEASURE -- Procedure to compute the sums and areas inside the apertures.
+
+procedure ap_rmmeasure (pixels, nx, ny, wx, wy, aperts, sums, areas, naperts)
+
+real	pixels[nx,ARB]		# subraster pixel values
+int	nx, ny			# dimensions of the subraster
+real	wx, wy			# center of subraster
+real	aperts[ARB]		# array of apertures
+double	sums[ARB]		# array of sums
+double	areas[ARB]		# aperture areas
+int	naperts			# number of apertures
+
+int	i, j, k
+double	fctn
+real	apmaxsq, dy2, r2, r
+
+begin
+	# Initialize.
+	apmaxsq = (aperts[naperts] + 0.5) ** 2
+	call aclrd (sums, naperts)
+	call aclrd (areas, naperts)
+
+	# Loop over the pixels.
+	do j = 1, ny {
+	    dy2 = (j - wy) ** 2
+	    do i = 1, nx {
+		r2 = (i - wx) ** 2 + dy2
+		if (r2 > apmaxsq)
+		    next
+		r = sqrt (r2) - 0.5
+		do k = 1, naperts {
+		    if (r > aperts[k])
+			next
+		    fctn = max (0.0, min (1.0, aperts[k] - r))
+		    sums[k] = sums[k] + fctn * pixels[i,j]
+		    areas[k] = areas[k] + fctn
+		}
+	    }
+	}
+end
+
+
+# AP_BRMMEASURE -- Procedure to compute the sums and areas inside the apertures.
+
+procedure ap_brmmeasure (pixels, nx, ny, wx, wy, datamin, datamax, aperts,
+	sums, areas, naperts, minapert)
+
+real	pixels[nx,ARB]		# subraster pixel values
+int	nx, ny			# dimensions of the subraster
+real	wx, wy			# center of subraster
+real	datamin			# minimum good data value
+real	datamax			# maximum good data value
+real	aperts[ARB]		# array of apertures
+double	sums[ARB]		# array of sums
+double	areas[ARB]		# aperture areas
+int	naperts			# number of apertures
+int	minapert		# minimum number of apertures
+
+int	i, j, k, kindex
+double	fctn
+real	apmaxsq, dy2, r2, r, pixval
+
+begin
+	# Initialize.
+	apmaxsq = (aperts[naperts] + 0.5) ** 2
+	call aclrd (sums, naperts)
+	call aclrd (areas, naperts)
+	minapert = naperts + 1
+
+	# Loop over the pixels.
+	do j = 1, ny {
+	    dy2 = (j - wy) ** 2
+	    do i = 1, nx {
+		r2 = (i - wx) ** 2 + dy2
+		if (r2 > apmaxsq)
+		    next
+		r = sqrt (r2) - 0.5
+		kindex = naperts + 1
+		pixval = pixels[i,j]
+		do k = 1, naperts {
+		    if (r > aperts[k])
+			next
+		    kindex = min (k, kindex)
+		    fctn = max (0.0, min (1.0, aperts[k] - r))
+		    sums[k] = sums[k] + fctn * pixval
+		    areas[k] = areas[k] + fctn
+		}
+		if (kindex < minapert) {
+		    if (pixval < datamin || pixval > datamax)
+			minapert = kindex
+		}
+	    }
+	}
+end