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
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
|
include <error.h>
include <imhdr.h>
include <imio.h>
include <mwset.h>
# APMW_OPEN -- Open APMW structure.
# APMW_CLOSE -- Close APMW structure.
# APMW_SETAP -- Set aperture values in APMW structure.
# APMW_SAVEIM -- Set WCS in image header.
# APMW_WCSFIX -- Fix up WCS
# Output formats
define ONEDSPEC 1 # Individual 1D spectra
define MULTISPEC 2 # Multiple spectra
define ECHELLE 3 # Echelle spectra
define STRIP 4 # Strip spectra
define NORM 5 # Normalized spectra
define FLAT 6 # Flat spectra
define RATIO 7 # Ratio of data to model
define DIFF 8 # Difference of data and model
define FIT 9 # Model
define NOISE 10 # Noise calculation
# Data structure for the apertures. This version assumes the coordinates
# are the same for all the apertures.
define APMW_LEN (8 + $1 * 6) # Structure length
define APMW_LABEL Memi[$1] # WCS label
define APMW_UNITS Memi[$1+1] # WCS units
define APMW_DTYPE Memi[$1+2] # Dispersion type
define APMW_NW Memi[$1+3] # Number of pixels
define APMW_W1 Memd[P2D($1+4)] # Starting coordinate
define APMW_DW Memd[P2D($1+6)] # Coordinate per pixel
define APMW_AP Memi[$1+6*($2-1)+8] # Aperture
define APMW_BEAM Memi[$1+6*($2-1)+9] # Beam
define APMW_APLOW Memd[P2D($1+6*($2-1)+10)] # Aperture low
define APMW_APHIGH Memd[P2D($1+6*($2-1)+12)] # Aperture high
# APMW_OPEN -- Open APMW structure.
pointer procedure apmw_open (in, out, dispaxis, naps, nw)
pointer in #I Input IMIO pointer
pointer out #I Output IMIO pointer
int dispaxis #I Input dispersion axis
int naps #I Number of apertures
int nw #I Number of dispersion pixels
pointer apmw #O Returned APMW pointer
int imgeti()
double mw_c1trand()
pointer mw, ct, mw_openim(), mw_sctran()
errchk mw_openim, mw_sctran, mw_c1trand, apmw_wcsfix
begin
# Allocate data structure.
call malloc (apmw, APMW_LEN(naps), TY_STRUCT)
call malloc (APMW_LABEL(apmw), SZ_LINE, TY_CHAR)
call malloc (APMW_UNITS(apmw), SZ_LINE, TY_CHAR)
# Set defaults.
call strcpy ("Pixel", Memc[APMW_LABEL(apmw)], SZ_LINE)
Memc[APMW_UNITS(apmw)] = EOS
APMW_DTYPE(apmw,i) = -1
APMW_NW(apmw,i) = nw
APMW_W1(apmw,i) = 1.
APMW_DW(apmw,i) = 1.
# Get WCS info from input image.
iferr {
mw = mw_openim (in)
iferr (APMW_DTYPE(apmw) = imgeti (in, "DC-FLAG"))
APMW_DTYPE(apmw) = -1
iferr (call mw_gwattrs (mw, dispaxis, "label",
Memc[APMW_LABEL(apmw)], SZ_LINE)) {
if (APMW_DTYPE(apmw) == -1)
call strcpy ("Pixel", Memc[APMW_LABEL(apmw)], SZ_LINE)
else
call strcpy ("Wavelength", Memc[APMW_LABEL(apmw)], SZ_LINE)
}
iferr (call mw_gwattrs (mw, dispaxis, "units",
Memc[APMW_UNITS(apmw)], SZ_LINE)) {
if (APMW_DTYPE(apmw) == -1)
Memc[APMW_UNITS(apmw)] = EOS
else
call strcpy ("Angstroms", Memc[APMW_UNITS(apmw)], SZ_LINE)
}
call apmw_wcsfix (in, mw)
iferr (ct = mw_sctran (mw, "logical", "world", dispaxis))
call error (1,
"Coordinate system ignored (rotated?). Using pixel coordinates.")
APMW_W1(apmw) = mw_c1trand (ct, 1D0)
APMW_DW(apmw) = mw_c1trand (ct, double (nw))
APMW_DW(apmw) = (APMW_DW(apmw)-APMW_W1(apmw))/(nw-1)
} then
call erract (EA_WARN)
call mw_close (mw)
return (apmw)
end
# APMW_CLOSE -- Close APMW structure.
procedure apmw_close (apmw)
pointer apmw # APMW pointer
begin
call mfree (APMW_LABEL(apmw), TY_CHAR)
call mfree (APMW_UNITS(apmw), TY_CHAR)
call mfree (apmw, TY_STRUCT)
end
# APMW_SETAP -- Set aperture values in APMW structure.
procedure apmw_setap (apmw, line, ap, beam, aplow, aphigh)
pointer apmw # APMW pointer
int line # Image line
int ap # Aperture
int beam # Beam
real aplow # Aperture lower limit
real aphigh # Aperture upper limit
begin
APMW_AP(apmw,line) = ap
APMW_BEAM(apmw,line) = beam
APMW_APLOW(apmw,line) = aplow
APMW_APHIGH(apmw,line) = aphigh
end
# APMW_SAVEIM -- Save WCS in image header.
procedure apmw_saveim (apmw, im, fmt)
pointer apmw #I APMW pointer
pointer im #I IMIO pointer
int fmt #I Output format
int i, naps, wcsdim, axes[3], imaccf()
double r[3], w[3], cd[9]
bool strne()
pointer sp, key, str, mw, list, mw_open(), imofnlu(), imgnfn()
errchk imdelf
data axes/1,2,3/
begin
call smark (sp)
call salloc (key, SZ_FNAME, TY_CHAR)
call salloc (str, SZ_LINE, TY_CHAR)
if (fmt == STRIP)
naps = 1
else
naps = IM_LEN(im, 2)
# Workaround for truncation of header during image header copy.
IM_HDRLEN(im) = IM_LENHDRMEM(im)
# Delete keywords.
list = imofnlu (im, "SLFIB[0-9]*")
while (imgnfn (list, Memc[key], SZ_FNAME) != EOF)
call imdelf (im, Memc[key])
call imcfnl (list)
# Add aperture parameters to image header.
do i = 1, naps {
call sprintf (Memc[key], SZ_FNAME, "APNUM%d")
call pargi (i)
call sprintf (Memc[str], SZ_LINE, "%d %d %.2f %.2f")
call pargi (APMW_AP(apmw,i))
call pargi (APMW_BEAM(apmw,i))
call pargd (APMW_APLOW(apmw,i))
call pargd (APMW_APHIGH(apmw,i))
call imastr (im, Memc[key], Memc[str])
if (naps == 1) {
call sprintf (Memc[key], SZ_FNAME, "APID%d")
call pargi (i)
ifnoerr (call imgstr (im, Memc[key], Memc[str], SZ_LINE)) {
if (strne (Memc[str], IM_TITLE(im))) {
call imastr (im, "MSTITLE", IM_TITLE(im))
call strcpy (Memc[str], IM_TITLE(im), SZ_IMTITLE)
}
call imdelf (im, Memc[key])
}
}
}
# Add dispersion parameters to image header.
if (APMW_DTYPE(apmw) != -1)
call imaddi (im, "DC-FLAG", APMW_DTYPE(apmw))
else if (imaccf (im, "DC-FLAG") == YES)
call imdelf (im, "DC-FLAG")
if (APMW_NW(apmw) < IM_LEN(im,1))
call imaddi (im, "NP2", APMW_NW(apmw))
else if (imaccf (im, "NP2") == YES)
call imdelf (im, "NP2")
iferr (call imdelf (im, "dispaxis"))
;
if (fmt == STRIP)
call imaddi (im, "dispaxis", 1)
# Set WCS in image header.
wcsdim = IM_NPHYSDIM(im)
mw = mw_open (NULL, wcsdim)
if (fmt == STRIP)
call mw_newsystem (mw, "linear", wcsdim)
else
call mw_newsystem (mw, "equispec", wcsdim)
call mw_swtype (mw, axes, wcsdim, "linear", "")
if (Memc[APMW_LABEL(apmw)] != EOS)
call mw_swattrs (mw, 1, "label", Memc[APMW_LABEL(apmw)])
if (Memc[APMW_UNITS(apmw)] != EOS)
call mw_swattrs (mw, 1, "units", Memc[APMW_UNITS(apmw)])
call aclrd (r, 3)
call aclrd (w, 3)
call aclrd (cd, 9)
r[1] = 1.
w[1] = APMW_W1(apmw)
cd[1] = APMW_DW(apmw)
if (wcsdim == 2)
cd[4] = 1.
if (wcsdim == 3) {
cd[5] = 1.
cd[9] = 1.
}
call mw_swtermd (mw, r, w, cd, wcsdim)
call mw_saveim (mw, im)
call mw_close (mw)
call sfree (sp)
end
# APMW_WCSFIX -- Fix up WCS to avoid CDELT=0 which occurs if there are WCS
# keywords in the header but no CDELT.
procedure apmw_wcsfix (im, mw)
pointer im # IMIO pointer
pointer mw # MWCS pointer
int i, ndim, mw_stati()
double val
pointer sp, r, w, cd
errchk mw_gwtermd, mw_swtermd
begin
call mw_seti (mw, MW_USEAXMAP, NO)
ndim = mw_stati (mw, MW_NDIM)
call smark (sp)
call salloc (r, ndim, TY_DOUBLE)
call salloc (w, ndim, TY_DOUBLE)
call salloc (cd, ndim*ndim, TY_DOUBLE)
# Check cd terms. Assume no rotation.
call mw_gwtermd (mw, Memd[r], Memd[w], Memd[cd], ndim)
do i = 0, ndim-1 {
val = Memd[cd+i*(ndim+1)]
if (val == 0D0) {
Memd[w+i] = 1D0
Memd[cd+i*(ndim+1)] = 1D0
}
}
call mw_swtermd (mw, Memd[r], Memd[w], Memd[cd], ndim)
call sfree (sp)
end
|
