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
|
include <mwset.h>
# RG_RXYL -- Compute the grid of logical coordinates.
procedure rg_rxyl (xl, yl, nx, ny, x1, x2, y1, y2)
double xl[ARB] #O array of output x coordinates
double yl[ARB] #O array of output y coordinates
int nx #I the size of the grid in x
int ny #I the size of the grid in y
double x1 #I the lower limit of the grid in x
double x2 #I the upper limit of the grid in x
double y1 #I the lower limit of the grid in y
double y2 #I the upper limit of the grid in y
double xstep, ystep, x, y
int i, j, npts
begin
if (nx == 1)
xstep = 0.0d0
else
xstep = (x2 - x1) / (nx - 1)
if (ny == 1)
ystep = 0.0d0
else
ystep = (y2 - y1) / (ny - 1)
npts = 0
y = y1
do j = 1, ny {
x = x1
do i = 1, nx {
npts = npts + 1
xl[npts] = x
yl[npts] = y
x = x + xstep
}
y = y + ystep
}
end
# RG_XYTOXY -- Compute the world coordinate list give the wcs descriptor
# and the logical coordinates.
pointer procedure rg_xytoxy (mw, xl, yl, xw, yw, npts, inwcs, outwcs, ax1, ax2)
pointer mw #I the wcs descriptor
double xl[ARB] #I the input logical x coordinate
double yl[ARB] #I the input logical y coordinate
double xw[ARB] #O the output world x coordinate
double yw[ARB] #O the output world y coordinate
int npts #I the number of coordinates.
char inwcs[ARB] #I the input wcs
char outwcs[ARB] #I the output wcs
int ax1 #I the logical x axis
int ax2 #I the logical y axis
int i, axbits
pointer ct
double mw_c1trand()
int mw_stati()
pointer mw_sctran()
errchk mw_sctran()
begin
# Compile the transformation.
if (mw == NULL) {
ct = NULL
} else if (mw_stati (mw, MW_NDIM) >= 2) {
axbits = 2 ** (ax1 - 1) + 2 ** (ax2 - 1)
iferr (ct = mw_sctran (mw, inwcs, outwcs, axbits))
ct = NULL
} else {
axbits = 2 ** (ax1 - 1)
iferr (ct = mw_sctran (mw, inwcs, outwcs, axbits))
ct = NULL
}
# Compute the world coordinates.
if (ct == NULL) {
call amovd (xl, xw, npts)
call amovd (yl, yw, npts)
} else if (mw_stati (mw, MW_NDIM) == 2) {
do i = 1, npts
call mw_c2trand (ct, xl[i], yl[i], xw[i], yw[i])
} else {
do i = 1, npts {
xw[i] = mw_c1trand (ct, xl[i])
yw[i] = yl[i]
}
}
return (ct)
end
|
