1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
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
|
