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
|
include <imhdr.h>
include "../../lib/astrom.h"
include "../../lib/aimpars.h"
include <pkg/skywcs.h>
# AT_PARWCS -- Compute a FITS WCS from an image using WCS definitions
# read from the AWCSPARS parameter file and stored in the astromery
# package descriptor. At the moment I am going to keep this routine simple
# by not worrying about the units of any quantities but the world coordinates
# of the reference point. This routine can be made more sophisticated later
# as time permits. The information is there ...
int procedure at_parwcs (im, at, update, verbose)
pointer im #I the input image descriptor
pointer at #I the astrometry package descriptor
bool update #I update rather than list the wcs
bool verbose #I verbose mode ?
double xref, yref, xmag, ymag, xrot, yrot, lngref, latref, dval
pointer sp, wfield, wtype, ctype, wcst, sym, coo, mw
int i, wkey, lngunits, latunits, coostat, stat
double at_imhms(), imgetd(), at_statd()
pointer at_statp(), stfind()
int at_wrdstr(), at_stati(), sk_decwcs()
bool streq()
errchk imgetd()
begin
# Return if the input is not 2D.
if (IM_NDIM(im) != 2)
return (ERR)
# Return if the wcs pointer is undefined.
if (at_statp (at, PWCS) == NULL)
return (ERR)
# Return if the keyword symbol table is undefined.
wcst = at_statp (at, WCST)
if (wcst == NULL)
return (ERR)
# Allocate working space.
call smark (sp)
call salloc (wfield, SZ_FNAME, TY_CHAR)
call salloc (wtype, SZ_FNAME, TY_CHAR)
call salloc (ctype, SZ_FNAME, TY_CHAR)
# Initialize.
xref = (1.0d0 + IM_LEN(im,1)) / 2.0d0
yref = (1.0d0 + IM_LEN(im,2)) / 2.0d0
xmag = INDEFD
ymag = INDEFD
xrot = 180.0d0
yrot= 0.0d0
lngref = INDEFD
latref = INDEFD
lngunits = 0
latunits = 0
call strcpy ("tan", Memc[wtype], SZ_FNAME)
call strcpy ("J2000", Memc[ctype], SZ_FNAME)
do i = 1, AT_NWFIELDS {
# Which keyword have we got ?
wkey = at_wrdstr (i, Memc[wfield], SZ_FNAME, AT_WFIELDS)
switch (wkey) {
# Get the x reference point in pixels.
case WCS_WXREF:
sym = stfind (wcst, Memc[wfield])
if (sym != NULL) {
if (streq (AT_WCSTKVAL(sym), "INDEF"))
dval = at_statd (at, WXREF)
else iferr (dval = imgetd (im, AT_WCSTKVAL(sym)))
dval = at_statd (at, WXREF)
} else
dval = at_statd (at, WXREF)
if (! IS_INDEFD(dval))
xref = dval
case WCS_WYREF:
sym = stfind (wcst, Memc[wfield])
if (sym != NULL) {
if (streq (AT_WCSTKVAL(sym), "INDEF"))
dval = at_statd (at, WYREF)
else iferr (dval = imgetd (im, AT_WCSTKVAL(sym)))
dval = at_statd (at, WYREF)
} else
dval = at_statd (at, WYREF)
if (! IS_INDEFD(dval))
yref = dval
case WCS_WXMAG:
sym = stfind (wcst, Memc[wfield])
if (sym != NULL) {
if (streq (AT_WCSTKVAL(sym), "INDEF"))
dval = at_statd (at, WXMAG)
else iferr (dval = imgetd (im, AT_WCSTKVAL(sym)))
dval = at_statd (at, WXMAG)
} else
dval = at_statd (at, WXMAG)
if (! IS_INDEFD(dval))
xmag = dval
case WCS_WYMAG:
sym = stfind (wcst, Memc[wfield])
if (sym != NULL) {
if (streq (AT_WCSTKVAL(sym), "INDEF"))
dval = at_statd (at, WYMAG)
else iferr (dval = imgetd (im, AT_WCSTKVAL(sym)))
dval = at_statd (at, WYMAG)
} else
dval = at_statd (at, WYMAG)
if (! IS_INDEFD(dval))
ymag = dval
case WCS_WXROT:
sym = stfind (wcst, Memc[wfield])
if (sym != NULL) {
if (streq (AT_WCSTKVAL(sym), "INDEF"))
dval = at_statd (at, WXROT)
else iferr (dval = imgetd (im, AT_WCSTKVAL(sym)))
dval = at_statd (at, WXROT)
} else
dval = at_statd (at, WXROT)
if (! IS_INDEFD(dval))
xrot = dval
case WCS_WYROT:
sym = stfind (wcst, Memc[wfield])
if (sym != NULL) {
if (streq (AT_WCSTKVAL(sym), "INDEF"))
dval = at_statd (at, WYROT)
else iferr (dval = imgetd (im, AT_WCSTKVAL(sym)))
dval = at_statd (at, WYROT)
} else
dval = at_statd (at, WYROT)
if (! IS_INDEFD(dval))
yrot = dval
case WCS_WRAREF:
sym = stfind (wcst, Memc[wfield])
if (sym != NULL) {
if (streq (AT_WCSTKVAL(sym), "INDEF"))
dval = at_statd (at, WRAREF)
else {
dval = at_imhms (im, AT_WCSTKVAL(sym))
if (IS_INDEFD(dval)) {
iferr (dval = imgetd (im, AT_WCSTKVAL(sym)))
dval = at_statd (at, WRAREF)
}
}
} else
dval = at_statd (at, WRAREF)
if (! IS_INDEFD(dval))
lngref = dval
case WCS_WDECREF:
sym = stfind (wcst, Memc[wfield])
if (sym != NULL) {
if (streq (AT_WCSTKVAL(sym), "INDEF"))
dval = at_statd (at, WDECREF)
else {
dval = at_imhms (im, AT_WCSTKVAL(sym))
if (IS_INDEFD(dval)) {
iferr (dval = imgetd (im, AT_WCSTKVAL(sym)))
dval = at_statd (at, WDECREF)
}
}
} else
dval = at_statd (at, WDECREF)
if (! IS_INDEFD(dval))
latref = dval
case WCS_WRAUNITS:
lngunits = at_stati (at, WRAUNITS)
case WCS_WDECUNITS:
latunits = at_stati (at, WDECUNITS)
case WCS_WPROJ:
sym = stfind (wcst, Memc[wfield])
if (sym != NULL) {
if (streq (AT_WCSTKVAL(sym), "INDEF"))
call at_stats (at, WPROJ, Memc[wtype], SZ_FNAME)
else iferr (call imgstr (im, AT_WCSTKVAL(sym), Memc[wtype],
SZ_FNAME))
call at_stats (at, WPROJ, Memc[wtype], SZ_FNAME)
} else
call at_stats (at, WPROJ, Memc[wtype], SZ_FNAME)
if (streq (Memc[wtype], "INDEF"))
call strcpy ("tan", Memc[wtype], SZ_FNAME)
case WCS_WSYSTEM:
sym = stfind (wcst, Memc[wfield])
if (sym != NULL) {
if (streq (AT_WCSTKVAL(sym), "INDEF"))
call at_stats (at, WSYSTEM, Memc[ctype], SZ_FNAME)
else iferr (call imgstr (im, AT_WCSTKVAL(sym), Memc[ctype],
SZ_FNAME))
call at_stats (at, WSYSTEM, Memc[ctype], SZ_FNAME)
} else
call at_stats (at, WSYSTEM, Memc[ctype], SZ_FNAME)
if (streq (Memc[ctype], "INDEF"))
call strcpy ("J2000", Memc[ctype], SZ_FNAME)
default:
;
}
}
# Update the header.
if (IS_INDEFD(xmag) || IS_INDEFD(ymag) || IS_INDEFD(lngref) ||
IS_INDEFD(latref)) {
stat = ERR
} else {
# Open coordinate system struct
coostat = sk_decwcs (Memc[ctype], mw, coo, NULL)
if (coostat == ERR || mw != NULL) {
if (mw != NULL)
call mw_close (mw)
stat = ERR
} else {
if (verbose)
call printf (
" Writing FITS wcs using default parameters\n")
if (lngunits > 0)
call sk_seti (coo, S_NLNGUNITS, lngunits)
if (latunits > 0)
call sk_seti (coo, S_NLATUNITS, latunits)
call at_uwcs (im, coo, Memc[wtype], lngref, latref, xref,
yref, xmag, ymag, xrot, yrot, false, update)
stat = OK
}
# Close the coordinate structure
if (coo != NULL)
call sk_close (coo)
}
call sfree (sp)
return (stat)
end
|
