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
|
# The code here is taken from t_transform.x in the longslit package. The
# changes are to use a sum instead of an average when multiple surfaces
# are given and not to use the xgs interface. Also the convergence
# tolerance is user specified since in this application the units might
# not be pixels.
define MAX_ITERATE 20
define ERROR 0.05
define FUDGE 0.5
# TR_INVERT -- Given user coordinate surfaces U(X,Y) and V(X,Y)
# (if none use one-to-one mapping and if more than one sum)
# corresponding to a given U and V and also the various partial
# derivatives. This is done using a gradient following interative
# method based on evaluating the partial derivative at each point
# and solving the linear Taylor expansions simultaneously. The last
# point sampled is used as the starting point. Thus, if the
# input U and V progress smoothly then the number of iterations
# can be small. The output is returned in x and y and in the derivative array
# DER. A point outside of the surfaces is returned as the nearest
# point at the edge of the surfaces in the DER array.
procedure tr_invert (usf, nusf, vsf, nvsf, u, v, x, y, der,
xmin, xmax, ymin, ymax, tol)
pointer usf[ARB], vsf[ARB] # User coordinate surfaces U(X,Y) and V(X,Y)
int nusf, nvsf # Number of surfaces for each coordinate
double u, v # Input U and V to determine X and Y
double x, y # Output X and Y
double der[8] # Last result as input, new result as output
# 1=X, 2=Y, 3=U, 4=DUDX, 5=DUDY, 6=V,
# 7=DVDX, 8=DVDY
double xmin, xmax, ymin, ymax # Limits of coordinate surfaces.
double tol # Tolerance
int i, j, nedge
double fudge, du, dv, dx, dy, tmp[3]
begin
# Use the last result as the starting point for the next position.
# If this is near the desired value then the interation will converge
# quickly. Allow a iteration to go off the surface twice.
# Quit when DX and DY are within tol.
nedge = 0
do i = 1, MAX_ITERATE {
du = u - der[3]
dv = v - der[6]
dx = (der[8] * du - der[5] * dv) /
(der[8] * der[4] - der[5] * der[7])
dy = (dv - der[7] * dx) / der[8]
fudge = 1 - FUDGE / i
x = der[1] + fudge * dx
y = der[2] + fudge * dy
der[1] = max (xmin, min (xmax, x))
der[2] = max (ymin, min (ymax, y))
if ((abs (dx) < tol) && (abs (dy) < tol))
break
if (nusf == 0)
der[3] = der[1]
else if (nusf == 1) {
call dgsder (usf[1], der[1], der[2], der[3], 1, 0, 0)
call dgsder (usf[1], der[1], der[2], der[4], 1, 1, 0)
call dgsder (usf[1], der[1], der[2], der[5], 1, 0, 1)
} else {
call dgsder (usf[1], der[1], der[2], der[3], 1, 0, 0)
call dgsder (usf[1], der[1], der[2], der[4], 1, 1, 0)
call dgsder (usf[1], der[1], der[2], der[5], 1, 0, 1)
do j = 2, nusf {
call dgsder (usf[j], der[1], der[2], tmp[1], 1, 0, 0)
call dgsder (usf[j], der[1], der[2], tmp[2], 1, 1, 0)
call dgsder (usf[j], der[1], der[2], tmp[3], 1, 0, 1)
der[3] = der[3] + tmp[1]
der[4] = der[4] + tmp[2]
der[5] = der[5] + tmp[3]
}
}
if (nvsf == 0)
der[6] = der[2]
else if (nvsf == 1) {
call dgsder (vsf[1], der[1], der[2], der[6], 1, 0, 0)
call dgsder (vsf[1], der[1], der[2], der[7], 1, 1, 0)
call dgsder (vsf[1], der[1], der[2], der[8], 1, 0, 1)
} else {
call dgsder (vsf[1], der[1], der[2], der[6], 1, 0, 0)
call dgsder (vsf[1], der[1], der[2], der[7], 1, 1, 0)
call dgsder (vsf[1], der[1], der[2], der[8], 1, 0, 1)
do j = 2, nvsf {
call dgsder (vsf[j], der[1], der[2], tmp[1], 1, 0, 0)
call dgsder (vsf[j], der[1], der[2], tmp[2], 1, 1, 0)
call dgsder (vsf[j], der[1], der[2], tmp[3], 1, 0, 1)
der[6] = der[6] + tmp[1]
der[7] = der[7] + tmp[2]
der[8] = der[8] + tmp[3]
}
}
}
end
# TR_INIT -- Since the inversion iteration always begins from the last
# point we need to initialize before the first call to TR_INVERT.
procedure tr_init (usf, nusf, vsf, nvsf, x, y, der)
pointer usf[ARB], vsf[ARB] # User coordinate surfaces
int nusf, nvsf # Number of surfaces for each coordinate
double x, y # Starting X and Y
double der[8] # Inversion data
int j
double tmp[3]
begin
der[1] = x
der[2] = y
if (nusf == 0) {
der[3] = der[1]
der[4] = 1.
der[5] = 0.
} else if (nusf == 1) {
call dgsder (usf[1], der[1], der[2], der[3], 1, 0, 0)
call dgsder (usf[1], der[1], der[2], der[4], 1, 1, 0)
call dgsder (usf[1], der[1], der[2], der[5], 1, 0, 1)
} else {
call dgsder (usf[1], der[1], der[2], der[3], 1, 0, 0)
call dgsder (usf[1], der[1], der[2], der[4], 1, 1, 0)
call dgsder (usf[1], der[1], der[2], der[5], 1, 0, 1)
do j = 2, nusf {
call dgsder (usf[j], der[1], der[2], tmp[1], 1, 0, 0)
call dgsder (usf[j], der[1], der[2], tmp[2], 1, 1, 0)
call dgsder (usf[j], der[1], der[2], tmp[3], 1, 0, 1)
der[3] = der[3] + tmp[1]
der[4] = der[4] + tmp[2]
der[5] = der[5] + tmp[3]
}
}
if (nvsf == 0) {
der[6] = der[2]
der[7] = 0.
der[8] = 1.
} else if (nvsf == 1) {
call dgsder (vsf[1], der[1], der[2], der[6], 1, 0, 0)
call dgsder (vsf[1], der[1], der[2], der[7], 1, 1, 0)
call dgsder (vsf[1], der[1], der[2], der[8], 1, 0, 1)
} else {
call dgsder (vsf[1], der[1], der[2], der[6], 1, 0, 0)
call dgsder (vsf[1], der[1], der[2], der[7], 1, 1, 0)
call dgsder (vsf[1], der[1], der[2], der[8], 1, 0, 1)
do j = 2, nvsf {
call dgsder (vsf[j], der[1], der[2], tmp[1], 1, 0, 0)
call dgsder (vsf[j], der[1], der[2], tmp[2], 1, 1, 0)
call dgsder (vsf[j], der[1], der[2], tmp[3], 1, 0, 1)
der[6] = der[6] + tmp[1]
der[7] = der[7] + tmp[2]
der[8] = der[8] + tmp[3]
}
}
end
|
