aboutsummaryrefslogtreecommitdiff
path: root/noao/onedspec/sensfunc/sfextinct.x
diff options
context:
space:
mode:
Diffstat (limited to 'noao/onedspec/sensfunc/sfextinct.x')
-rw-r--r--noao/onedspec/sensfunc/sfextinct.x226
1 files changed, 226 insertions, 0 deletions
diff --git a/noao/onedspec/sensfunc/sfextinct.x b/noao/onedspec/sensfunc/sfextinct.x
new file mode 100644
index 00000000..f7b95326
--- /dev/null
+++ b/noao/onedspec/sensfunc/sfextinct.x
@@ -0,0 +1,226 @@
+include <pkg/gtools.h>
+include "sensfunc.h"
+
+define RANGE_AIRMASS 0.1 # Minimum airmass range
+define SIGMA_AIRMASS 0.05 # Minimum sigma in airmass
+
+# SF_EXINCT -- Determine a residual extinction curve. At each wavelength
+# for which there are multiple observations or neighboring wavelengths
+# such that there is a sufficient airmass range determine the slope
+# of the sensitivity vs airmass. Residual sensitivity is used to minimize
+# wavelength scatter when multiple wavelengths are needed because of
+# nonoverlapping standard star data. Each such slope is a measure of the
+# residual extinction at that wavelength. To make the residual extinction
+# curve fit the extinction vs. wavelength using ICFIT.
+
+procedure sf_extinct (gp, stds, nstds, cv, ecv, function, order)
+
+pointer gp # Graphics structure
+pointer stds[nstds] # Standard star data
+int nstds # Number of standard stars
+pointer cv # Sensitivity function curve
+pointer ecv # Residual extinction curve
+char function[ARB] # Fitting function
+int order # Function order
+
+bool ans
+int i, j, n, nwaves, sum, npts, scan()
+real a, amin, amax, rms, rms1, r2, sig, cveval()
+double x, y, z, sumx, sumy, sumz, sumx2, sumxy
+pointer sp, waves, sens, airm, xp, yp, fp, wp, ic
+pointer gt, gt_init()
+errchk salloc, xt_sort3, icg_fit, ic_open
+
+define cancel_ 99
+
+begin
+ # Cancel previous extinction if defined.
+ if (ecv != NULL)
+ goto cancel_
+
+ # Check for minimum airmass range and determine number of points.
+ # Ignore added objects and composite data.
+ amin = 100.
+ amax = 0.
+ nwaves = 0
+ do i = 1, nstds - 2 {
+ if (STD_FLAG(stds[i]) != SF_INCLUDE)
+ next
+ nwaves = nwaves + STD_NWAVES(stds[i])
+ a = STD_AIRMASS(stds[i])
+ amin = min (amin, a)
+ amax = max (amax, a)
+ }
+ if (amax - amin < RANGE_AIRMASS) {
+ call printf (
+ "Insufficient airmass coverage for extinction determination")
+ return
+ }
+
+ # Extract data to be used and sort by wavelength.
+ # The data is wavelength, airmass, and residual sensitivity.
+ call smark (sp)
+ call salloc (waves, nwaves, TY_REAL)
+ call salloc (sens, nwaves, TY_REAL)
+ call salloc (airm, nwaves, TY_REAL)
+
+ nwaves = 0
+ do i = 1, nstds-2 {
+ if (STD_FLAG(stds[i]) != SF_INCLUDE)
+ next
+ n = STD_NWAVES(stds[i])
+ a = STD_AIRMASS(stds[i])
+ xp = STD_WAVES(stds[i])
+ yp = STD_SENS(stds[i])
+ fp = STD_FIT(stds[i])
+ wp = STD_WTS(stds[i])
+ do j = 1, n {
+ if (Memr[wp] != 0.) {
+ Memr[airm+nwaves] = a
+ Memr[waves+nwaves] = Memr[xp]
+ Memr[sens+nwaves] = Memr[yp] - Memr[fp]
+ nwaves = nwaves + 1
+ }
+ xp = xp + 1
+ yp = yp + 1
+ fp = fp + 1
+ wp = wp + 1
+ }
+ }
+
+ call xt_sort3 (Memr[waves], Memr[sens], Memr[airm], nwaves)
+
+ # Bin points with common wavelengths or at least two points.
+ sum = 0
+ sumx = 0.
+ sumy = 0.
+ sumz = 0.
+ sumx2 = 0.
+ sumxy = 0.
+ n = 0
+ do i = 0, nwaves-2 {
+ x = Memr[airm+i]
+ y = Memr[sens+i]
+ z = Memr[waves+i]
+ sum = sum + 1
+ sumx = sumx + x
+ sumy = sumy + y
+ sumx2 = sumx2 + x * x
+ sumxy = sumxy + x * y
+ sumz = sumz + z
+
+ if ((z == Memr[waves+i+1]) || (sum < 2))
+ next
+
+ x = sumx2 - sumx * sumx / sum
+ if (x > SIGMA_AIRMASS) {
+ Memr[waves+n] = sumz / sum
+ Memr[sens+n] = (sumx * sumy / sum - sumxy) / x
+ Memr[airm+n] = 1.
+ n = n + 1
+ sum = 0
+ sumx = 0.
+ sumy = 0.
+ sumz = 0.
+ sumx2 = 0.
+ sumxy = 0.
+ }
+ }
+ if (sum > 1) {
+ x = sumx2 - sumx * sumx / sum
+ if (x > SIGMA_AIRMASS) {
+ Memr[waves+n] = sumz / sum
+ Memr[sens+n] = (sumx * sumy / sum - sumxy) / x
+ Memr[airm+n] = 1.
+ n = n + 1
+ }
+ }
+
+ if (n < 2) {
+ call printf ("Cannot determine residual extinction")
+ call sfree (sp)
+ return
+ }
+
+ # Fit residual extinction curve using ICFIT.
+ gt = gt_init()
+ call gt_sets (gt, GTTYPE, "mark")
+ call gt_seti (gt, GTCOLOR, GP_PLCOLOR(gp))
+ call ic_open (ic)
+ call ic_putr (ic, "xmin", min (STD_WSTART(stds[1]), STD_WEND(stds[1])))
+ call ic_putr (ic, "xmax", max (STD_WSTART(stds[1]), STD_WEND(stds[1])))
+ call ic_pstr (ic, "function", "chebyshev")
+ call ic_puti (ic, "order", 1)
+ call ic_pstr (ic, "xlabel", "wavelength")
+ call ic_pstr (ic, "ylabel", "residual extinction")
+ call ic_pstr (ic, "yunits", "mag")
+ call icg_fit (ic, GP_GIO(gp), "cursor", gt, ecv, Memr[waves],
+ Memr[sens], Memr[airm], n)
+ call ic_closer (ic)
+ call gt_free (gt)
+
+ # Determine significance of the fit.
+ call sf_fit (stds, nstds, cv, function, order,
+ min (GP_WSTART(gp), GP_WEND(gp)), max (GP_WSTART(gp), GP_WEND(gp)))
+ call sf_rms (stds, nstds, rms1, npts)
+ do i = 1, nstds - 2 {
+ if (STD_FLAG(stds[i]) != SF_INCLUDE)
+ next
+ n = STD_NWAVES(stds[i])
+ a = STD_AIRMASS(stds[i])
+ xp = STD_WAVES(stds[i])
+ yp = STD_SENS(stds[i])
+ call cvvector (ecv, Memr[xp], Memr[sens], n)
+ call amulkr (Memr[sens], a, Memr[sens], n)
+ call aaddr (Memr[yp], Memr[sens], Memr[yp], n)
+ }
+ call sf_fit (stds, nstds, cv, function, order,
+ min (GP_WSTART(gp), GP_WEND(gp)), max (GP_WSTART(gp), GP_WEND(gp)))
+ call sf_rms (stds, nstds, rms, npts)
+ do i = 1, SF_NGRAPHS
+ if (GP_SHDR(gp,i) != NULL)
+ call shdr_close (GP_SHDR(gp,i))
+
+ r2 = 1 - rms ** 2 / rms1 ** 2
+ sig = r2 * (nwaves - 2) / max (0.01, 1. - r2)
+ if (sig <= 0.0)
+ sig = 0.
+ else
+ sig = sqrt (sig)
+
+ # Apply to data if desired.
+ call printf (
+ "Significance = %4.1f sigma: Apply residual extinction correction? ")
+ call pargr (sig)
+ call flush (STDOUT)
+
+ ans = false
+ if (scan() != EOF)
+ call gargb (ans)
+
+ # Undo last fit if not applying correction.
+ if (!ans)
+ goto cancel_
+
+ call printf ("Residual extinction correction applied")
+ call sfree (sp)
+ return
+
+cancel_
+ do i = 1, nstds - 2 {
+ if (STD_FLAG(stds[i]) == SF_EXCLUDE)
+ next
+ n = STD_NWAVES(stds[i])
+ a = STD_AIRMASS(stds[i])
+ xp = STD_WAVES(stds[i])
+ yp = STD_SENS(stds[i])
+ do j = 1, n {
+ Memr[yp] = Memr[yp] - a * cveval (ecv, Memr[xp])
+ xp = xp + 1
+ yp = yp + 1
+ }
+ }
+ call cvfree (ecv)
+ call printf ("Residual extinction correction canceled")
+ call sfree (sp)
+end
generated by cgit (git 2.34.1) at 2025-09-20 07:19:37 -0400