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
|
include <math.h>
define DIRS "|ang2row|row2ang|"
define ANG2PIX 1
define PIX2ANG 2
define CUNITS "|hourdegree|degrees|radians|"
define H 1
define D 2
define R 3
define MTYPES "|nest|ring|"
define NEST 1
define RING 2
# T_HPCTRAN -- Convert between HEALPix rows and spherical coordinates.
#
# It is up to the user to know the coordinate and map type; e.g.
# galactic/nested, equatorial/ring. However, the use can use
# whatever units for the coordinate type; e.g. hours/degrees, radians.
#
# The HEALPix row is 1 indexed to be consistent with IRAF conventions.
# This row can be used to access the map data with TTOOLS tasks.
procedure t_hpctran ()
int dir # Direction (ang2row|row2ang)
int row # HEALpix map row (1 indexed)
double lng # RA/longitude
double lat # DEC/latitude
int nside # Resolution parameter
int cunits # Coordinate units
int mtype # HEALpix map type
char str[10]
int clgeti(), clgwrd()
double clgetd()
errchk ang2row, row2ang
begin
# Get parameters.
dir = clgwrd ("direction", str, 10, DIRS)
nside = clgeti ("nside")
cunits = clgwrd ("cunits", str, 10, CUNITS)
mtype = clgwrd ("maptype", str, 10, MTYPES)
switch (dir) {
case ANG2PIX:
lng = clgetd ("lng")
lat = clgetd ("lat")
switch (cunits) {
case 0:
call error (1, "Unknown coordinate units")
case H:
lng = lng * 15D0
case R:
lng = RADTODEG(lng)
lat = RADTODEG(lat)
}
call ang2row (row, lng, lat, mtype, nside)
call clputi ("row", row)
case PIX2ANG:
row = clgeti ("row")
call row2ang (row, lng, lat, mtype, nside)
switch (cunits) {
case 0:
call error (1, "Unknown coordinate units")
case H:
lng = lng / 15D0
case R:
lng = DEGTORAD(lng)
lat = DEGTORAD(lat)
}
call clputd ("lng", lng)
call clputd ("lat", lat)
}
# Output the map row.
call printf ("%d %g %g\n")
call pargi (row)
call pargd (lng)
call pargd (lat)
end
# TEST_HEALPIX2 -- Test routine as in the HEALPix distribution.
procedure test_healpix2 ()
double theta, phi
int nside
int ipix, npix, dpix, ip1
begin
call printf("Starting C Healpix pixel routines test\n")
nside = 1024
dpix = 23
# Find the number of pixels in the full map
npix = 12*nside*nside
call printf("Number of pixels in full map: %d\n")
call pargi (npix)
call printf("dpix: %d\n")
call pargi (dpix)
call printf("Nest -> ang -> Ring -> ang -> Nest\n")
# call printf("Nest -> ang -> Nest\n")
# call printf("Ring -> ang -> Ring\n")
for (ipix = 0; ipix < npix; ipix = ipix + dpix) {
call pix2ang_nest(nside, ipix, theta, phi)
call ang2pix_ring(nside, theta, phi, ip1)
call pix2ang_ring(nside, ip1, theta, phi)
call ang2pix_nest(nside, theta, phi, ip1)
# call pix2ang_ring(nside, ipix, theta, phi)
# call ang2pix_ring(nside, theta, phi, ip1)
if (ip1 != ipix) {
call printf("Error: %d %d %d\n")
call pargi (nside)
call pargi (ipix)
call pargi (ip1)
}
}
call printf("%d\n")
call pargi (nside)
call printf("test completed\n\n")
end
|
