diff options
Diffstat (limited to 'noao/digiphot/apphot/wphot/apgmeasure.x')
| -rw-r--r-- | noao/digiphot/apphot/wphot/apgmeasure.x | 190 |
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 |