aboutsummaryrefslogtreecommitdiff
path: root/noao/twodspec/apextract/apvariance.x
diff options
context:
space:
mode:
Diffstat (limited to 'noao/twodspec/apextract/apvariance.x')
-rw-r--r--noao/twodspec/apextract/apvariance.x420
1 files changed, 420 insertions, 0 deletions
diff --git a/noao/twodspec/apextract/apvariance.x b/noao/twodspec/apextract/apvariance.x
new file mode 100644
index 00000000..015eed74
--- /dev/null
+++ b/noao/twodspec/apextract/apvariance.x
@@ -0,0 +1,420 @@
+include <gset.h>
+include "apertures.h"
+
+
+# AP_VARIANCE -- Variance weighted extraction based on profile and CCD noise.
+# If desired reject deviant pixels. In addition to the variance weighted
+# spectrum, the unweighted and uncleaned "raw" spectrum is extracted and
+# a sigma spectrum is returned. Wavelengths with saturated pixels are
+# flagged with 0 value and negative sigma if cleaning.
+
+procedure ap_variance (im, ap, dbuf, nc, nl, c1, l1, sbuf, svar, profile,
+ nx, ny, xs, ys, spec, raw, specsig, nsubaps, asi)
+
+pointer im # IMIO pointer
+pointer ap # Aperture structure
+pointer dbuf # Data buffer
+int nc, nl # Size of data buffer
+int c1, l1 # Origin of data buffer
+pointer sbuf # Sky values (NULL if none)
+pointer svar # Sky variance
+real profile[ny,nx] # Profile (returned)
+int nx, ny # Size of profile array
+int xs[ny], ys # Origin of profile array
+real spec[ny,nsubaps] # Spectrum
+real raw[ny,nsubaps] # Raw spectrum
+real specsig[ny,nsubaps] # Sky variance in, spectrum sigma out
+int nsubaps # Number of subapertures
+pointer asi # Image interpolator for edge pixel weighting
+
+real rdnoise # Readout noise in RMS data numbers.
+real gain # Gain in photons per data number.
+real saturation # Maximum value for an unsaturated pixel.
+bool clean # Clean cosmic rays?
+real nclean # Number of pixels to clean
+real lsigma, usigma # Rejection sigmas.
+
+bool sat
+int fd, iterate, niterate, nrej, irej, nreject
+int i, ix, iy, ix1, ix2
+real low, high, step, shift, x1, x2, wt1, wt2, s, w, dat, sk, var, var0
+real sum, wsum, wvsum, sum1, sum2, total1, total2
+real vmin, resid, rrej
+pointer cv, gp
+pointer sp, str, work, wt, xplot, yplot, eplot, fplot, data, sky, data1
+
+real apgetr(), apgimr(), ap_cveval()
+bool apgetb()
+errchk apgimr, ap_asifit
+
+begin
+ # Get task parameters.
+ gain = apgimr ("gain", im)
+ rdnoise = apgimr ("readnoise", im) ** 2
+ saturation = apgetr ("saturation")
+ if (!IS_INDEF(saturation))
+ saturation = saturation * gain
+ clean = apgetb ("clean")
+ lsigma = apgetr ("lsigma")
+ usigma = apgetr ("usigma")
+ call ap_popen (gp, fd, "clean")
+
+ # Allocate memory and one index.
+ call smark (sp)
+ call salloc (str, SZ_LINE, TY_CHAR)
+ call salloc (work, 6*nx, TY_REAL)
+ wt = work - 1
+ xplot = wt + nx
+ yplot = xplot + nx
+ eplot = yplot + nx
+ fplot = eplot + nx
+ data1 = fplot + nx
+ if (sbuf == NULL) {
+ call salloc (sky, nx, TY_REAL)
+ call aclrr (Memr[sky], nx)
+ sky = sky - 1
+ var0 = rdnoise
+ }
+
+ # Initialize
+ if (rdnoise == 0.)
+ vmin = 1.
+ else
+ vmin = rdnoise
+ if (clean) {
+ nclean = apgetr ("nclean")
+ if (nclean < 1.)
+ niterate = max (1., nclean * nx)
+ else
+ niterate = max (1., min (real (nx), nclean))
+ } else
+ niterate = 0
+
+ call aclrr (spec, ny * nsubaps)
+ call aclrr (raw, ny * nsubaps)
+ call amovkr (-1., specsig, ny * nsubaps)
+
+ i = AP_AXIS(ap)
+ low = AP_CEN(ap,i) + AP_LOW(ap,i)
+ high = AP_CEN(ap,i) + AP_HIGH(ap,i)
+ step = (high - low) / nsubaps
+ cv = AP_CV(ap)
+
+ # For each line compute the weighted spectrum and then iterate
+ # to reject deviant pixels. Rejected pixels are flagged by negative
+ # variance.
+
+ nreject = 0
+ total1 = 0.
+ total2 = 0.
+ do iy = 1, ny {
+ shift = ap_cveval (cv, real (iy + ys - 1)) - c1 + 1
+ x1 = max (0.5, low + shift) + c1 - xs[iy]
+ x2 = min (nc + 0.49, high + shift) + c1 - xs[iy]
+ if (x1 >= x2)
+ next
+ ix1 = nint (x1)
+ ix2 = nint (x2)
+
+ call ap_asifit (dbuf+(iy+ys-1-l1)*nc, nc, xs[iy]-c1+1,
+ low+shift, high+shift, data, asi)
+# data = dbuf + (iy + ys - 1 - l1) * nc + xs[iy] - c1 - 1
+ if (sbuf != NULL) {
+ sky = sbuf + (iy - 1) * nx - 1
+ var0 = rdnoise + Memr[svar+iy-1]
+ }
+
+ # Set pixel weights for summing.
+# if (asi != NULL)
+# call asifit (asi, Memr[data], nc-xs[iy]+c1)
+ call ap_edge (asi, x1+1, x2+1, wt1, wt2)
+
+ # First estimate spectrum by summing across the aperture.
+ # Accumulate the raw spectrum and set up various arrays for
+ # plotting and later access.
+
+ sat = false
+ sum = 0.
+ wsum = 0.
+ wvsum = 0.
+ do ix = ix1, ix2 {
+ if (ix1 == ix2)
+ w = wt1
+ else if (ix == ix1)
+ w = wt1
+ else if (ix == ix2)
+ w = wt2
+ else
+ w = 1.
+ dat = Memr[data+ix]
+ if (!IS_INDEF(saturation))
+ if (dat > saturation)
+ sat = true
+ sk = Memr[sky+ix]
+ raw[iy,1] = raw[iy,1] + w * (dat - sk)
+
+ Memr[xplot+ix] = ix + xs[iy] - 1
+ Memr[yplot+ix] = dat - sk
+ Memr[data1+ix] = dat - sk
+ Memr[wt+ix] = w
+ var = max (vmin, var0 + max (0., dat))
+ w = profile[iy,ix] / var
+ var = sqrt (var)
+ Memr[eplot+ix] = var
+ sum = sum + w * (dat - sk)
+ wsum = wsum + w * profile[iy,ix]
+ wvsum = wvsum + (w * var) ** 2
+ }
+ if (wsum > 0.) {
+ spec[iy,1] = sum / wsum
+ specsig[iy,1] = sqrt (wvsum) / abs (wsum)
+ } else {
+ spec[iy,1] = 0.
+ specsig[iy,1] = -1.
+ }
+
+ sum = 0.
+ wsum = 0.
+ wvsum = 0.
+ sum1 = 0.
+ sum2 = 0.
+ do ix = ix1, ix2 {
+ sum1 = sum1 + Memr[wt+ix] * Memr[data1+ix]
+ if (Memr[eplot+ix] <= 0.)
+ next
+ sk = Memr[sky+ix]
+ s = max (0., spec[iy,1]) * profile[iy,ix]
+ var = max (vmin, var0 + (s + sk))
+ w = profile[iy,ix] / var
+ var = sqrt (var)
+ Memr[eplot+ix] = var
+ Memr[fplot+ix] = s
+ sum = sum + w * Memr[data1+ix]
+ wsum = wsum + w * profile[iy,ix]
+ wvsum = wvsum + (w * var) ** 2
+ }
+ if (wsum > 0.) {
+ spec[iy,1] = sum / wsum
+ specsig[iy,1] = sqrt (wvsum) / abs (wsum)
+ sum2 = sum2 + spec[iy,1]
+ } else {
+ spec[iy,1] = 0.
+ specsig[iy,1] = -1.
+ sum1 = 0.
+ sum2 = 0.
+ }
+
+ # Reject cosmic rays one at a time.
+ nrej = 0
+ do iterate = 1, niterate {
+ irej = 0
+ rrej = 0.
+
+ # Compute revised variance estimate using profile model
+ # skip rejected pixels, find worst pixel.
+
+ do ix = ix1, ix2 {
+ if (Memr[eplot+ix] <= 0.)
+ next
+ s = max (0., spec[iy,1]) * profile[iy,ix]
+ sk = Memr[sky+ix]
+ var = sqrt (max (vmin, var0 + max (0., s + sk)))
+ Memr[fplot+ix] = s
+ Memr[eplot+ix] = var
+
+ resid = (Memr[data1+ix] - Memr[fplot+ix]) / var
+ if (abs (resid) > abs (rrej)) {
+ rrej = resid
+ irej = ix
+ }
+ }
+
+ # Reject worst outlier.
+
+ if (rrej <= -lsigma || rrej >= usigma) {
+ Memr[eplot+irej] = -Memr[eplot+irej]
+ Memr[data1+irej] = Memr[fplot+irej]
+ nrej = nrej + 1
+ } else
+ break
+
+ # Update spectrum estimate excluding rejected pixels.
+ sum = 0.
+ wsum = 0.
+ wvsum = 0.
+ sum1 = 0.
+ sum2 = 0.
+ do ix = ix1, ix2 {
+ sum1 = sum1 + Memr[wt+ix] * Memr[data1+ix]
+ if (Memr[eplot+ix] <= 0.)
+ next
+ w = profile[iy,ix] / Memr[eplot+ix]**2
+ sum = sum + w * Memr[data1+ix]
+ wsum = wsum + w * profile[iy,ix]
+ wvsum = wvsum + (w * Memr[eplot+ix]) ** 2
+ }
+
+ if (wsum > 0.) {
+ spec[iy,1] = sum / wsum
+ specsig[iy,1] = sqrt (wvsum) / abs (wsum)
+ sum2 = sum2 + spec[iy,1]
+ } else {
+ spec[iy,1] = 0.
+ specsig[iy,1] = -1.
+ sum1 = 0.
+ sum2 = 0.
+ }
+ }
+
+ nreject = nreject + nrej
+ total1 = total1 + sum1
+ total2 = total2 + sum2
+
+ # Calculate subapertures if desired.
+ if (nsubaps > 1) {
+ do i = 1, nsubaps {
+ x1 = max (0.5, low + (i - 1) * step + shift) + c1 - xs[iy]
+ x2 = min (nc + 0.49, low + i * step + shift) + c1 - xs[iy]
+ if (x1 >= x2) {
+ spec[iy,i] = 0.
+ raw[iy,i] = 0.
+ specsig[iy,i] = -1.
+ next
+ }
+ ix1 = nint (x1)
+ ix2 = nint (x2)
+ call ap_edge (asi, x1+1, x2+1, wt1, wt2)
+
+ sum = 0.
+ wvsum = 0.
+ raw[iy,i] = 0.
+ do ix = ix1, ix2 {
+ if (ix1 == ix2)
+ w = wt1
+ else if (ix == ix1)
+ w = wt1
+ else if (ix == ix2)
+ w = wt2
+ else
+ w = 1.
+ raw[iy,i] = raw[iy,i] + w * Memr[yplot+ix]
+ if (Memr[eplot+ix] <= 0.)
+ next
+ w = profile[iy,ix] / Memr[eplot+ix]**2
+ sum = sum + w * Memr[data1+ix]
+ wvsum = wvsum + (w * Memr[eplot+ix]) ** 2
+ }
+
+ if (wsum > 0.) {
+ spec[iy,i] = sum / wsum
+ specsig[iy,i] = sqrt (wvsum) / abs (wsum)
+ } else {
+ spec[iy,i] = 0.
+ specsig[iy,i] = -1.
+ }
+ }
+ }
+
+ # Flag points with saturated pixels.
+ if (sat)
+ do i = 1, nsubaps
+ specsig[iy,i] = -specsig[iy,i]
+
+ # Plot profile with cosmic rays if desired.
+ if (gp != NULL && nrej > 0 && spec[iy,1] > 0.) {
+ call sprintf (Memc[str], SZ_LINE, "Profile %4d")
+ call pargi (iy)
+ s = Memr[yplot+ix1] - abs (Memr[eplot+ix1])
+ w = Memr[yplot+ix1] + abs (Memr[eplot+ix1])
+ do ix = ix1+1, ix2 {
+ s = min (s, Memr[yplot+ix] - abs (Memr[eplot+ix]))
+ w = max (w, Memr[yplot+ix] + abs (Memr[eplot+ix]))
+ }
+ sum = w - s
+ x1 = ix1 + xs[iy] - 2
+ x2 = ix2 + xs[iy]
+ s = s - 0.1 * sum
+ w = w + 0.1 * sum
+ call gclear (gp)
+ call gswind (gp, x1, x2, s, w)
+ call glabax (gp, Memc[str], "", "")
+
+ do ix = ix1, ix2 {
+ if (Memr[eplot+ix] > 0.) {
+ call gmark (gp, Memr[xplot+ix], Memr[yplot+ix],
+ GM_PLUS, 2., 2.)
+ call gmark (gp, Memr[xplot+ix], Memr[yplot+ix],
+ GM_VEBAR, 2., -6.*Memr[eplot+ix])
+ } else {
+ call gmark (gp, Memr[xplot+ix], Memr[yplot+ix],
+ GM_CROSS, 2., 2.)
+ call gmark (gp, Memr[xplot+ix], Memr[yplot+ix],
+ GM_VEBAR, 1., 6.*Memr[eplot+ix])
+ }
+ }
+ call gpline (gp, Memr[xplot+ix1], Memr[fplot+ix1], ix2-ix1+1)
+ }
+ }
+
+ # To avoid any bias, scale weighted extraction to same total flux
+ # as raw spectrum (with rejected pixels replaced by fit).
+
+ if (total1 * total2 <= 0.) {
+ call sprintf (Memc[str], SZ_LINE,
+ "EXTRACT: WARNING - Aperture %d:")
+ call pargi (AP_ID(ap))
+ call ap_log (Memc[str], YES, NO, YES)
+ call sprintf (Memc[str], SZ_LINE,
+ " Total variance weighted spectrum flux is %g")
+ call pargr (total2)
+ call ap_log (Memc[str], YES, NO, YES)
+ call sprintf (Memc[str], SZ_LINE,
+ " Total unweighted spectrum flux is %g")
+ call pargr (total1)
+ call ap_log (Memc[str], YES, NO, YES)
+ call sprintf (Memc[str], SZ_LINE,
+ " Variance spectrum bias factor ignored")
+ call ap_log (Memc[str], YES, NO, YES)
+ } else {
+ sum = total1 / total2
+ call amulkr (spec, sum, spec, ny * nsubaps)
+ call amulkr (specsig, sum, specsig, ny * nsubaps)
+ if (sum < .5 || sum > 2) {
+ call sprintf (Memc[str], SZ_LINE,
+ "EXTRACT: WARNING - Aperture %d:")
+ call pargi (AP_ID(ap))
+ call ap_log (Memc[str], YES, NO, YES)
+ call sprintf (Memc[str], SZ_LINE,
+ " Total variance weighted spectrum flux is %g")
+ call pargr (total2)
+ call ap_log (Memc[str], YES, NO, YES)
+ call sprintf (Memc[str], SZ_LINE,
+ " Total unweighted spectrum flux is %g")
+ call pargr (total1)
+ call ap_log (Memc[str], YES, NO, YES)
+ call sprintf (Memc[str], SZ_LINE,
+ "EXTRACT: Aperture %d variance spectrum bias factor is %g")
+ call pargi (AP_ID(ap))
+ call pargr (total1 / total2)
+ call ap_log (Memc[str], YES, NO, YES)
+ } else {
+ call sprintf (Memc[str], SZ_LINE,
+ "EXTRACT: Aperture %d variance spectrum bias factor is %g")
+ call pargi (AP_ID(ap))
+ call pargr (total1 / total2)
+ call ap_log (Memc[str], YES, NO, NO)
+ }
+ }
+
+ # Log the number of rejected pixels.
+ if (clean) {
+ call sprintf (Memc[str], SZ_LINE,
+ "EXTRACT: %d pixels rejected from aperture %d")
+ call pargi (nreject)
+ call pargi (AP_ID(ap))
+ call ap_log (Memc[str], YES, NO, NO)
+ }
+
+ call ap_pclose (gp, fd)
+ call sfree (sp)
+end
generated by cgit (git 2.34.1) at 2025-09-20 08:57:04 -0400