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
|
include <mach.h>
include <math.h>
include <math/gsurfit.h>
include <pkg/skywcs.h>
# CC_INIT_STD -- Get the parameter values relevant to the transformation from
# the cl or the database file.
#
procedure cc_init_std (dt, record, geometry, lngunits, latunits, sx1,
sy1, sx2, sy2, mw, coo)
pointer dt #I pointer to database file produced by geomap
char record[ARB] #I the name of the database record
int geometry #I the type of geometry to be computed
int lngunits #I the input ra / longitude units
int latunits #I the input dec / latitude units
pointer sx1, sy1 #O pointers to the linear x and y surfaces
pointer sx2, sy2 #O pointers to the x and y distortion surfaces
pointer mw #O pointer to the mwcs structure
pointer coo #O pointer to the coordinate structure
double lngref, latref
int recstat, proj
pointer sp, projstr, projpars
int cc_dtrecord(), strdic()
pointer cc_celwcs()
begin
call smark (sp)
call salloc (projstr, SZ_FNAME, TY_CHAR)
call salloc (projpars, SZ_LINE, TY_CHAR)
if (dt == NULL) {
call cc_rinit (lngunits, latunits, sx1, sy1, mw, coo)
sx2 = NULL
sy2 = NULL
} else {
recstat = cc_dtrecord (dt, record, geometry, coo, Memc[projpars],
lngref, latref, sx1, sy1, sx2, sy2)
if (recstat == ERR) {
coo = NULL
sx1 = NULL
sy1 = NULL
sx2 = NULL
sy2 = NULL
mw = NULL
} else {
call sscan (Memc[projpars])
call gargwrd (Memc[projstr], SZ_FNAME)
proj = strdic (Memc[projstr], Memc[projstr], SZ_FNAME,
WTYPE_LIST)
if (proj <= 0 || proj == WTYPE_LIN)
Memc[projpars] = EOS
mw = cc_celwcs (coo, Memc[projpars], lngref, latref)
}
}
call sfree (sp)
end
# CC_FREE_STD -- Free the previously defined transformation.
procedure cc_free_std (sx1, sy1, sx2, sy2, mw, coo)
pointer sx1, sy1 #U pointers to the linear x and y surfaces
pointer sx2, sy2 #U pointers to the x and y distortion surfaces
pointer mw #U pointer to the mwcs structure
pointer coo #U pointer to the celestial coordinate structure
begin
if (sx1 != NULL)
call dgsfree (sx1)
if (sy1 != NULL)
call dgsfree (sy1)
if (sx2 != NULL)
call dgsfree (sx2)
if (sy2 != NULL)
call dgsfree (sy2)
if (mw != NULL)
call mw_close (mw)
if (coo != NULL)
call sk_close (coo)
end
# CC_RINIT -- Compute the required wcs structure from the input parameters.
procedure cc_rinit (lngunits, latunits, sx1, sy1, mw, coo)
int lngunits #I the input ra / longitude units
int latunits #I the input dec / latitude units
pointer sx1 #O pointer to the linear x coordinate surface
pointer sy1 #O pointer to the linear y coordinate surface
pointer mw #O pointer to the mwcs structure
pointer coo #O pointer to the celestial coordinate structure
double xref, yref, xscale, yscale, xrot, yrot, lngref, latref
int coostat, proj, tlngunits, tlatunits, pfd
pointer sp, projstr
double clgetd()
double dgseval()
int sk_decwcs(), sk_stati(), strdic(), open()
pointer cc_celwcs(), cc_rdproj()
errchk open()
begin
# Allocate some workin space.
call smark (sp)
call salloc (projstr, SZ_LINE, TY_CHAR)
# Get the reference point pixel coordinates.
xref = clgetd ("xref")
if (IS_INDEFD(xref))
xref = 0.0d0
yref = clgetd ("yref")
if (IS_INDEFD(yref))
yref = 0.0d0
# Get the scale factors.
xscale = clgetd ("xmag")
if (IS_INDEFD(xscale))
xscale = 1.0d0
yscale = clgetd ("ymag")
if (IS_INDEFD(yscale))
yscale = 1.0d0
# Get the rotation angles.
xrot = clgetd ("xrotation")
if (IS_INDEFD(xrot))
xrot = 0.0d0
xrot = -DEGTORAD(xrot)
yrot = clgetd ("yrotation")
if (IS_INDEFD(yrot))
yrot = 0.0d0
yrot = -DEGTORAD(yrot)
# Initialize the linear part of the solution.
call dgsinit (sx1, GS_POLYNOMIAL, 2, 2, NO, double (-MAX_REAL),
double (MAX_REAL), double (-MAX_REAL), double (MAX_REAL))
call dgsinit (sy1, GS_POLYNOMIAL, 2, 2, NO, double (-MAX_REAL),
double (MAX_REAL), double (-MAX_REAL), double (MAX_REAL))
call geo_rotmagd (sx1, sy1, xscale, yscale, xrot, yrot)
call geo_xyshiftd (sx1, sy1, -dgseval (sx1, xref, yref),
-dgseval (sy1, xref, yref))
lngref = clgetd ("lngref")
if (IS_INDEFD(lngref))
lngref = 0.0d0
latref = clgetd ("latref")
if (IS_INDEFD(latref))
latref = 0.0d0
coostat = sk_decwcs ("j2000", mw, coo, NULL)
if (coostat == ERR || mw != NULL) {
if (mw != NULL)
call mw_close (mw)
}
if (lngunits <= 0)
tlngunits = sk_stati (coo, S_NLNGUNITS)
else
tlngunits = lngunits
call sk_seti (coo, S_NLNGUNITS, tlngunits)
if (latunits <= 0)
tlatunits = sk_stati (coo, S_NLATUNITS)
else
tlatunits = latunits
call sk_seti (coo, S_NLATUNITS, tlatunits)
call clgstr ("projection", Memc[projstr], SZ_LINE)
iferr {
pfd = open (Memc[projstr], READ_ONLY, TEXT_FILE)
} then {
proj = strdic (Memc[projstr], Memc[projstr], SZ_LINE, WTYPE_LIST)
if (proj <= 0 || proj == WTYPE_LIN)
Memc[projstr] = EOS
} else {
proj = cc_rdproj (pfd, Memc[projstr], SZ_LINE)
call close (pfd)
}
mw = cc_celwcs (coo, Memc[projstr], lngref, latref)
call sfree (sp)
end
define MAX_NITER 20
# CC_DO_STD -- Transform the coordinates using the full transformation
# computed by CCMAP.
procedure cc_do_std (x, y, xt, yt, sx1, sy1, sx2, sy2, forward)
double x, y #I initial positions
double xt, yt #O transformed positions
pointer sx1, sy1 #I pointer to linear surfaces
pointer sx2, sy2 #I pointer to distortion surfaces
bool forward #I forward transform
double f, fx, fy, g, gx, gy, denom, dx, dy
int niter
pointer newsx, newsy
double dgseval()
begin
if (forward) {
xt = dgseval (sx1, x, y)
if (sx2 != NULL)
xt = xt + dgseval (sx2, x, y)
yt = dgseval (sy1, x, y)
if (sy2 != NULL)
yt = yt + dgseval (sy2, x, y)
} else {
xt = x / 1.0
yt = y / 1.0
call dgsadd (sx1, sx2, newsx)
call dgsadd (sy1, sy2, newsy)
niter = 0
repeat {
f = dgseval (newsx, xt, yt) - x
call dgsder (newsx, xt, yt, fx, 1, 1, 0)
call dgsder (newsx, xt, yt, fy, 1, 0, 1)
g = dgseval (newsy, xt, yt) - y
call dgsder (newsy, xt, yt, gx, 1, 1, 0)
call dgsder (newsy, xt, yt, gy, 1, 0, 1)
denom = fx * gy - fy * gx
dx = (-f * gy + g * fy) / denom
dy = (-g * fx + f * gx) / denom
xt = xt + dx
yt = yt + dy
if (max (abs (dx), abs (dy), abs(f), abs(g)) < 1.0e-5)
break
niter = niter + 1
} until (niter >= MAX_NITER)
call dgsfree (newsx)
call dgsfree (newsy)
}
end
|
