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
|
include <error.h>
include <syserr.h>
define EXTN_LOOKUP 10 # Interp index for de-extinction
define DEXTN_LOOKUP 11 # Interp index for differential extn table
define TEMP_SPACE 100 # Amount of temporary space to allocate
# GET_EXTN -- Get extinction from calibration file and
# any update as indicated from the SENSITIVITY
# computation
procedure get_extn (wave_tbl, extn_tbl, nwaves)
pointer wave_tbl, extn_tbl
int nwaves
pointer waves, extns
begin
# Get standard extinction values
call ext_load (waves, extns, nwaves)
# Copy values to external pointer.
# Use of salloc is incorrect but this is a hack on old code. FV
call salloc (extn_tbl, nwaves, TY_REAL)
call salloc (wave_tbl, nwaves, TY_REAL)
call amovr (Memr[waves], Memr[wave_tbl], nwaves)
call amovr (Memr[extns], Memr[extn_tbl], nwaves)
call mfree (waves, TY_REAL)
call mfree (extns, TY_REAL)
end
# DE_EXT_SPEC -- Apply extinction correction to a spectrum
procedure de_ext_spec (spec, airm, w0, wpc, wave_tbl, extn_tbl, nwaves, len)
real spec[ARB], wave_tbl[ARB], extn_tbl[ARB]
real airm, w0, wpc
int nwaves, len
int i, ierr
real wave, ext
bool lin_log
begin
# Assume linear dispersion, but possibly in LOG10
if (w0 < 5.0 && wpc < 0.05)
lin_log = true
else
lin_log = false
# Initialize interpolator
call intrp0 (EXTN_LOOKUP)
do i = 1, len {
wave = w0 + (i-1) * wpc
if (lin_log)
wave = 10.0 ** wave
# Table must be in wavelength, not log[]
call intrp (EXTN_LOOKUP, wave_tbl, extn_tbl, nwaves,
wave, ext, ierr)
spec[i] = spec[i] * 10.0 ** (0.4 * airm * ext)
}
end
# SUM_SPEC -- Add up counts within a specified region of a spectrum
# denoted by a wavelength range.
# The summation is active only over those pixels which
# are completely within the range specification.
# Data referenced outside the spectrum is ignored.
procedure sum_spec (spec, w1, w2, w0, wpc, counts, len)
real spec[ARB], w1, w2, w0, wpc, counts
int len
int i, pix1, pix2
real pix_index()
begin
# Compute pixel numbers from w1 to w2
pix1 = max (int (pix_index (w0, wpc, w1) + 1.0), 1)
pix2 = max (int (pix_index (w0, wpc, w2) ), pix1)
pix2 = min (pix2, len)
counts = 0.0
do i = pix1, pix2
counts = counts + spec[i]
# Guarantee that there are no negative counts
if (counts < 0.0)
counts = 0.0
end
# PIX_INDEX -- Returns the pixel index at wavelength for linearly
# dispersion corrected spectra
#
# The "Guess" is made that if the start wavelength for the
# spectrum is less than 5.0 and the dispersion is less than
# 0.05, the spectrum has been linearized in LOG10 space.
#
# Note that in IRAF, a pixel index effectively refers to the center of a pixel.
# So a spectrum must actually extend from w0-0.5*wpc to w0+(len+0.5)*wpc
real procedure pix_index (w0, wpc, w)
real w0, wpc, w
real xw
begin
# Check for LOG10 dispersion
if (w0 < 5.0 && wpc < 0.05)
xw = log10 (w)
else
xw = w
pix_index = (xw - w0) / wpc + 1.0
end
define NALLOC 128 # Allocation block size
# EXT_LOAD -- Read extinction data from database directory.
procedure ext_load (waves, extns, nwaves)
pointer waves, extns
int nwaves
real wave, extn
int fd, nalloc
pointer sp, file
int open(), fscan(), nscan(), errcode()
begin
call smark (sp)
call salloc (file, SZ_FNAME, TY_CHAR)
# Get the extinction file and open it.
call clgstr ("extinction", Memc[file], SZ_FNAME)
iferr (fd = open (Memc[file], READ_ONLY, TEXT_FILE)) {
switch (errcode()) {
case SYS_FNOFNAME:
nwaves = 2
call malloc (waves, nwaves, TY_REAL)
call malloc (extns, nwaves, TY_REAL)
Memr[waves] = 1000.
Memr[extns] = 0.
Memr[waves+1] = 10000.
Memr[extns+1] = 0.
call eprintf ("No extinction correction applied\n")
return
default:
call erract (EA_ERROR)
}
}
# Read the extinction data.
nalloc = 0
nwaves = 0
while (fscan (fd) != EOF) {
call gargr (wave)
call gargr (extn)
if (nscan() != 2)
next
if (nalloc == 0) {
nalloc = nalloc + NALLOC
call malloc (waves, nalloc, TY_REAL)
call malloc (extns, nalloc, TY_REAL)
} else if (nwaves == nalloc) {
nalloc = nalloc + NALLOC
call realloc (waves, nalloc, TY_REAL)
call realloc (extns, nalloc, TY_REAL)
}
Memr[waves+nwaves] = wave
Memr[extns+nwaves] = extn
nwaves = nwaves + 1
}
call close (fd)
if (nwaves == 0)
call error (1, "No extinction data found")
call realloc (waves, nwaves, TY_REAL)
call realloc (extns, nwaves, TY_REAL)
call sfree (sp)
end
# EXT_FREE -- Free extinction data arrays.
procedure ext_free (waves, extns)
pointer waves, extns
begin
call mfree (waves, TY_REAL)
call mfree (extns, TY_REAL)
end
|
