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 <imhdr.h>
include <error.h>
# T_EXTINCTION -- CL task for applying extinction corrections to images.
#
# The image headers must contain the parameters DISPAXIS, CRVALn,
# CRPIXn, and CDELTn to define the wavelength coordinates and
# either AIRMASS, ZD, or information needed to compute the zenith
# distance (HA, LATITUDE, RA, DEC, ST).
#
# The extinction table contains wavelengths and extinctions in
# magnitudes such that the multiplicative extinction correction
# is given by:
#
# correction = 10 ** (0.4 * airmass * extinction value)
#
# The extinction table need not be sorted.
procedure t_extinction()
int list1 # List of images to be corrected
int list2 # List of extinction corrected images
char table[SZ_FNAME] # Extinction table filename
bool extcor
char image1[SZ_FNAME], image2[SZ_FNAME]
int fd, nalloc, len_table
real wavelen, ext
pointer im1, im2, w, e
int clpopnu(), fscan(), nscan(), open(), clgfil()
bool imgetb(), streq()
pointer immap()
errchk ext_cor()
begin
# Get the list of images and the extinction table.
list1 = clpopnu ("input")
list2 = clpopnu ("output")
call clgstr ("extinction", table, SZ_FNAME)
# Read the extinction table. Dynamically allocate memory for the
# table.
fd = open (table, READ_ONLY, TEXT_FILE)
nalloc = 100
call malloc (w, nalloc, TY_REAL)
call malloc (e, nalloc, TY_REAL)
len_table = 0
while (fscan (fd) != EOF) {
call gargr (wavelen)
call gargr (ext)
if (nscan() < 2)
next
if (len_table == nalloc) {
nalloc = nalloc + 100
call realloc (w, nalloc, TY_REAL)
call realloc (e, nalloc, TY_REAL)
}
Memr[w + len_table] = wavelen
Memr[e + len_table] = ext
len_table = len_table + 1
}
call close (fd)
# If there are no extinction values in the table then return an error.
# Sort the extinction values by wavelength.
if (len_table > 0) {
call realloc (w, len_table, TY_REAL)
call realloc (e, len_table, TY_REAL)
call xt_sort2 (Memr[w], Memr[e], len_table)
} else {
call mfree (w, TY_REAL)
call mfree (e, TY_REAL)
call error (0, "No extinction values extinction table")
}
# Loop through each pair of input and output images. Check if
# the input image has been corrected previously. If TRUE then
# print message and go on to the next input image. If FALSE
# print message and apply extinction corrections.
# Missing information in the image header will return an error
# which will warn the user and go on to the next image.
while (clgfil (list1, image1, SZ_FNAME) != EOF) {
if (clgfil (list2, image2, SZ_FNAME) == EOF) {
call eprintf ("No output image for %s.\n")
call pargstr (image1)
next
}
if (streq (image1, image2)) {
im1 = immap (image1, READ_WRITE, 0)
im2 = im1
} else {
im1 = immap (image1, READ_ONLY, 0)
im2 = immap (image2, NEW_COPY, im1)
}
iferr (extcor = imgetb (im1, "extcor"))
extcor = false
if (extcor) {
call printf ("Image %s is extinction corrected.\n")
call pargstr (image1)
} else {
call printf ("Extinction correction: %s -> %s.\n")
call pargstr (image1)
call pargstr (image2)
call flush (STDOUT)
iferr (call do_extinct(im1, im2, Memr[w], Memr[e], len_table)) {
call printf ("!!No extinction correction for %s!!\n")
call pargstr (image1)
call flush (STDOUT)
call erract (EA_WARN)
}
}
if (im2 != im1)
call imunmap (im2)
call imunmap (im1)
}
# Finish up.
call mfree (w, TY_REAL)
call mfree (e, TY_REAL)
call clpcls (list1)
call clpcls (list2)
end
# DO_EXTINCT -- Apply extinction correction.
define SZ_FIELD 8 # Size of field string
procedure do_extinct (im1, im2, w, e, len_table)
pointer im1 # Input IMIO pointer
pointer im2 # Output IMIO pointer
real w[len_table] # Wavelengths
real e[len_table] # Extinction values
int len_table # Length of extinction table
char field[SZ_FIELD]
int laxis, paxis, npix, i, flag, dcflag
real crval, cdelt, crpix, airmass, wavelen, extval
long v1[IM_MAXDIM], v2[IM_MAXDIM]
pointer sp, ext, pix1, pix2
int imgeti(), imgnlr(), impnlr()
real imgetr(), img_airmass()
errchk get_daxis, imgeti, imgetr, img_airmass
begin
# Determine the dispersion axis and linear coordinates.
call get_daxis (im1, laxis, paxis)
call sprintf (field, SZ_FIELD, "crval%d")
call pargi (laxis)
crval = imgetr (im1, field)
call sprintf (field, SZ_FIELD, "crpix%d")
call pargi (laxis)
crpix = imgetr (im1, field)
call sprintf (field, SZ_FIELD, "cdelt%d")
call pargi (laxis)
iferr (cdelt = imgetr (im1, field)) {
call sprintf (field, SZ_FIELD, "cd%d_%d")
call pargi (laxis)
call pargi (laxis)
cdelt = imgetr (im1, field)
}
dcflag = imgeti (im1, "dc-flag")
# Determine the airmass.
airmass = img_airmass (im1)
# Determine the extinction values at each pixel.
npix = IM_LEN (im1, laxis)
call smark (sp)
call salloc (ext, npix, TY_REAL)
do i = 1, npix {
wavelen = crval + (i - crpix) * cdelt
if (dcflag == 1)
wavelen = 10. ** wavelen
call intrp (1, w, e, len_table, wavelen, extval, flag)
Memr[ext+i-1] = 10. ** (0.4 * airmass * extval)
}
# Loop through the image applying the extinction correction to each
# pixel.
call amovkl (long (1), v1, IM_MAXDIM)
call amovkl (long (1), v2, IM_MAXDIM)
while ((imgnlr(im1, pix1, v1) != EOF) &&
(impnlr(im2, pix2, v2) != EOF)) {
switch (laxis) {
case 1:
call amulr (Memr[pix1], Memr[ext], Memr[pix2], IM_LEN (im1, 1))
default:
extval = Memr[ext+v1[laxis]-2]
call amulkr (Memr[pix1], extval, Memr[pix2], IM_LEN (im1, 1))
}
}
call sfree (sp)
# Add the extinction correction flag, history, and return.
# The parameter ex-flag is added for compatibility with onedspec.
call imaddb (im2, "extcor", true)
call imaddi (im2, "ex-flag", 0)
call xt_phistory (im2, "Extinction correction applied.")
end
|
