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
|
# Copyright(c) 1986 Association of Universities for Research in Astronomy Inc.
include <gio.h>
# GCTRAN -- Transform a point in world coordinate system WCS_A to world
# coordinate system WCS_B. The transformation is performed by transforming
# from WCS_A to NDC, followed by a transformation from NDC to WCS_B. The
# transformation parameters are cached for efficiency when transforming
# multiple points in the same pairs of coordinate systems. Three types of
# transformations are supported: linear, log, and "elog". The latter is
# a logarithmic function defined for all X, i.e., for negative as well as
# positive X.
procedure gctran (gp, x1,y1, x2,y2, wcs_a, wcs_b)
pointer gp # graphics descriptor
real x1, y1 # coords of point in WCS_A (input)
real x2, y2 # coords of point in WCS_B (output)
int wcs_a # input WCS
int wcs_b # output WCS
int w, a
int wcsord, tran[2,2], wcs[2]
real morigin[2,2], worigin[2,2], scale[2,2], ds
real w1[2,2], w2[2,2], s1[2,2], s2[2,2], p1[2], p2[2]
pointer wp
bool fp_nondegenr()
real elogr(), aelogr()
begin
# Verify that the WCS has not changed since we were last called.
# WCSORD is a unique integer (ordinal) assigned by GIO each time a
# WCS is fixed.
if (GP_WCSSTATE(gp) != FIXED || GP_WCSORD(gp) != wcsord)
wcs[1] = -1
# Verify that cached transformation parameters are up to date, and if
# not, recompute them.
if (wcs[1] != wcs_a || wcs[2] != wcs_b) {
wcsord = GP_WCSORD(gp)
wcs[1] = wcs_a
wcs[2] = wcs_b
# Copy the WCS parameters into 2-dim arrays so that we can use the
# same code for both axes.
do w = 1, 2 {
wp = GP_WCSPTR (gp, wcs[w])
tran[1,w] = WCS_XTRAN(wp)
tran[2,w] = WCS_YTRAN(wp)
# If the window is degenerate enlarge the window until there
# is enough range to make a plot.
if (fp_nondegenr (WCS_WX1(wp), WCS_WX2(wp)))
GP_WCSSTATE(gp) = MODIFIED
if (fp_nondegenr (WCS_WY1(wp), WCS_WY2(wp)))
GP_WCSSTATE(gp) = MODIFIED
w1[1,w] = WCS_WX1(wp)
w2[1,w] = WCS_WX2(wp)
w1[2,w] = WCS_WY1(wp)
w2[2,w] = WCS_WY2(wp)
s1[1,w] = WCS_SX1(wp)
s2[1,w] = WCS_SX2(wp)
s1[2,w] = WCS_SY1(wp)
s2[2,w] = WCS_SY2(wp)
}
# Compute the transformation parameters for both axes and both
# world coordinate systems.
do w = 1, 2 # w = wcs index
do a = 1, 2 { # a = axis
morigin[a,w] = s1[a,w]
ds = s2[a,w] - s1[a,w]
if (tran[a,w] == LINEAR) {
worigin[a,w] = w1[a,w]
scale[a,w] = ds / (w2[a,w] - w1[a,w])
} else if (tran[a,w] == LOG && w1[a,w] > 0 && w2[a,w] > 0) {
worigin[a,w] = log10 (w1[a,w])
scale[a,w] = ds / (log10(w2[a,w]) - worigin[a,w])
} else {
worigin[a,w] = elogr (w1[a,w])
scale[a,w] = ds / (elogr(w2[a,w]) - worigin[a,w])
}
}
}
p1[1] = x1
p1[2] = y1
# Forward Transformation. Transform P1 (point-A) in wcs_a to NDC
# coordinates, if the input WCS is not number zero (the NDC coordinate
# system).
if (wcs_a != 0)
do a = 1, 2 {
if (tran[a,1] != LINEAR)
if (tran[a,1] == LOG) {
if (p1[a] <= 0)
p1[a] = INDEF
else
p1[a] = log10 (p1[a])
} else
p1[a] = elogr (p1[a])
p1[a] = ((p1[a] - worigin[a,1]) * scale[a,1]) + morigin[a,1]
}
# Inverse Transformation. Transform point P1, now in NDC coordinates,
# to WCS-B. If WCS-B is zero (NDC), we need only copy the points.
if (wcs_b == 0) {
p2[1] = p1[1]
p2[2] = p1[2]
} else {
do a = 1, 2 {
if (IS_INDEF (p1[a]))
p2[a] = INDEF
else {
p2[a] = (p1[a] - morigin[a,2]) / scale[a,2] + worigin[a,2]
if (tran[a,2] != LINEAR)
if (tran[a,2] == LOG)
p2[a] = 10.0 ** p2[a]
else
p2[a] = aelogr (p2[a])
}
}
}
x2 = p2[1]
y2 = p2[2]
end
|
