Initial commit

1 files changed, 190 insertions, 0 deletions

diff --git a/noao/digiphot/apphot/wphot/apgmeasure.x b/noao/digiphot/apphot/wphot/apgmeasure.x
new file mode 100644
index 00000000..738630f2
--- /dev/null
+++ b/noao/digiphot/apphot/wphot/apgmeasure.x
@@ -0,0 +1,190 @@
+# AP_GMEASURE -- Procedure to measure the fluxes and effective areas of a set of
+# apertures assuming a Gaussian weighting function.
+
+procedure ap_gmeasure (im, wx, wy, c1, c2, l1, l2, aperts, sums, areas,
+        naperts, sigsq, gain, varsky)
+
+pointer	im			# pointer to image
+real	wx, wy			# center of subraster
+int	c1, c2			# column limits
+int	l1, l2			# line limits
+real	aperts[ARB]		# array of apertures
+double	sums[ARB]		# array of sums
+double	areas[ARB]		# aperture areas
+int	naperts			# number of apertures
+real	sigsq			# the profile widht squared
+real	gain			# the sky value
+real	varsky			# the sky variance
+
+int	i, j, k, nx, yindex
+double	fctn, weight, norm
+pointer	sp, buf, sump, sumpw
+real	xc, yc, apmaxsq, dy2, r2, r, prof, var
+pointer	imgs2r()
+
+begin
+	# Initialize.
+	call smark (sp)
+	call salloc (sump, naperts, TY_DOUBLE)
+	call salloc (sumpw, naperts, TY_DOUBLE)
+	call aclrd (Memd[sump], naperts)
+	call aclrd (Memd[sumpw], naperts)
+
+	# Get array boundary parameters.
+	nx = c2 - c1 + 1
+	xc = wx - c1 + 1
+	yc = wy - l1 + 1
+	apmaxsq = (aperts[naperts] + 0.5) ** 2
+
+	# Clear out the accumulaters
+	call aclrd (sums, naperts)
+	call aclrd (areas, naperts)
+
+	# Loop over the pixels.
+	do j = l1, l2 {
+	    buf = imgs2r (im, c1, c2, j, j)
+	    if (buf == EOF) {
+		call sfree (sp)
+		return
+	    }
+	    yindex = j - l1 + 1
+	    dy2 = (yindex - yc) ** 2
+	    do i = 1, nx {
+		r2 = (i - xc) ** 2 + dy2
+		if (r2 > apmaxsq)
+		    next
+		prof = max (exp (-r2 / sigsq), 0.0)
+		if (prof <= 0.0)
+		    next
+		var = max (0.0, Memr[buf+i-1])
+		var = var / gain + varsky
+		if (var <= 0.0)
+		    next
+		weight = prof / var
+		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] + weight * fctn * Memr[buf+i-1]
+		    areas[k] = areas[k] + weight * fctn
+		    Memd[sump+k-1] = Memd[sump+k-1] + prof
+		    Memd[sumpw+k-1] = Memd[sumpw+k-1] + weight * prof
+		}
+	    }
+	}
+
+	# Normalize.
+	do k = 1, naperts {
+	    if (Memd[sumpw+k-1] <= 0.0d0)
+	        norm = 0.0d0
+	    else
+	        norm = Memd[sump+k-1] / Memd[sumpw+k-1]
+	    sums[k] = sums[k] * norm
+	    areas[k] = areas[k] * norm
+	}
+
+	call sfree (sp)
+end
+
+
+# AP_BGMEASURE -- Procedure to measure the fluxes and effective areas of a set
+# of apertures assuming a Gaussian weighting function.
+
+procedure ap_bgmeasure (im, wx, wy, c1, c2, l1, l2, datamin, datamax,
+	aperts, sums, areas, naperts, minapert, sigsq, gain, varsky)
+
+pointer	im			# pointer to image
+real	wx, wy			# center of subraster
+int	c1, c2			# column limits
+int	l1, l2			# line limits
+real	datamin			# minimum good data
+real	datamax			# maximum good data
+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 aperture fo bad pixels
+real	sigsq			# the profile widht squared
+real	gain			# the image gain value
+real	varsky			# the sky variance
+
+int	i, j, k, nx, yindex, kindex
+double	fctn, weight, norm
+pointer	sp, buf, sump, sumpw
+real	xc, yc, apmaxsq, dy2, r2, r
+real	pixval, prof, var
+pointer	imgs2r()
+
+begin
+	# Initialize.
+	call smark (sp)
+	call salloc (sump, naperts, TY_DOUBLE)
+	call salloc (sumpw, naperts, TY_DOUBLE)
+	call aclrd (Memd[sump], naperts)
+	call aclrd (Memd[sumpw], naperts)
+	minapert = naperts + 1
+
+	# Get array boundary parameters.
+	nx = c2 - c1 + 1
+	xc = wx - c1 + 1
+	yc = wy - l1 + 1
+	apmaxsq = (aperts[naperts] + 0.5) ** 2
+
+	# Clear out the accumulaters
+	call aclrd (sums, naperts)
+	call aclrd (areas, naperts)
+
+	# Loop over the pixels.
+	do j = l1, l2 {
+	    buf = imgs2r (im, c1, c2, j, j)
+	    if (buf == EOF) {
+		call sfree (sp)
+		return
+	    }
+	    yindex = j - l1 + 1
+	    dy2 = (yindex - yc) ** 2
+	    do i = 1, nx {
+		r2 = (i - xc) ** 2 + dy2
+		if (r2 > apmaxsq)
+		    next
+		prof = max (exp (-r2 / sigsq), 0.0)
+		if (prof <= 0.0)
+		    next
+		pixval = Memr[buf+i-1]
+		var = max (0.0, pixval)
+		var = pixval / gain + varsky
+		if (var <= 0.0)
+		    next
+		weight = prof / var
+		r = sqrt (r2) - 0.5
+		kindex = naperts + 1
+		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] + weight * fctn * pixval
+		    areas[k] = areas[k] + weight * fctn
+		    Memd[sump+k-1] = Memd[sump+k-1] + prof
+		    Memd[sumpw+k-1] = Memd[sumpw+k-1] + weight * prof
+		}
+		if (kindex < minapert) {
+		    if (pixval < datamin || pixval > datamax)
+			minapert = kindex
+		}
+	    }
+	}
+
+	# Normalize.
+	do k = 1, naperts {
+	    if (Memd[sumpw+k-1] <= 0.0d0)
+	        norm = 0.0d0
+	    else
+	        norm = Memd[sump+k-1] / Memd[sumpw+k-1]
+	    sums[k] = sums[k] * norm
+	    areas[k] = areas[k] * norm
+	}
+
+	call sfree (sp)
+end