healpix_test.py 3.46 KB
Newer Older
Martin Reinecke's avatar
Martin Reinecke committed
1
import ducc_0_1.healpix as ph
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
import numpy as np
import math

def random_ptg(vlen):
  res = np.empty((vlen, 2), dtype=np.float64)
  res[:,0] = np.arccos((np.random.random_sample(vlen)-0.5)*2)
  res[:,1] = np.random.random_sample(vlen)*2*math.pi
  return res

def check_pixangpix(vlen,ntry,nside,isnest):
  base = ph.Healpix_Base (nside, "NEST" if isnest else "RING")
  cnt = 0
  while cnt < ntry:
    cnt += 1
    inp = np.random.randint(low=0, high=12*nside*nside-1, size=vlen)
    out = base.ang2pix(base.pix2ang(inp))
    if not np.array_equal(inp, out):
      raise ValueError("Test failed")

def check_vecpixvec(vlen, ntry, nside, isnest):
  base = ph.Healpix_Base (nside, "NEST" if isnest else "RING")
  cnt = 0
  while cnt < ntry:
    cnt += 1
    inp = ph.ang2vec(random_ptg(vlen))
    out = base.pix2vec(base.vec2pix(inp))
    if np.any(ph.v_angle(inp,out) > base.max_pixrad()):
      raise ValueError("Test failed")

def check_pixangvecpix(vlen, ntry, nside, isnest):
  base = ph.Healpix_Base (nside, "NEST" if isnest else "RING")
  cnt = 0
  while cnt < ntry:
    cnt += 1
    inp = np.random.randint(low=0, high=12*nside*nside-1, size=vlen)
    out = base.vec2pix(ph.ang2vec(base.pix2ang(inp)))
    if not np.array_equal(inp,out):
      raise ValueError("Test failed")

def check_pixvecangpix(vlen, ntry, nside, isnest):
  base = ph.Healpix_Base (nside, "NEST" if isnest else "RING")
  cnt = 0
  while cnt < ntry:
    cnt += 1
    inp = np.random.randint(low=0, high=12*nside*nside-1, size=vlen)
    out = base.ang2pix(ph.vec2ang(base.pix2vec(inp)))
    if not np.array_equal(inp,out):
      raise ValueError("Test failed")

def check_pixvecpix(vlen,ntry,nside,isnest):
  base=ph.Healpix_Base (nside,"NEST" if isnest else "RING")
  cnt=0
  while (cnt<ntry):
    cnt+=1
    inp=np.random.randint(low=0,high=12*nside*nside-1,size=vlen)
    out=base.vec2pix(base.pix2vec(inp))
    if (np.array_equal(inp,out)==False):
      raise ValueError("Test failed")

def check_ringnestring(vlen,ntry,nside):
  base=ph.Healpix_Base (nside,"NEST")
  cnt=0
  while (cnt<ntry):
    cnt+=1
    inp=np.random.randint(low=0,high=12*nside*nside-1,size=vlen)
    out=base.nest2ring(base.ring2nest(inp))
    if (np.array_equal(inp,out)==False):
      raise ValueError("Test failed")

def check_pixxyfpix(vlen,ntry,nside,isnest):
  base=ph.Healpix_Base (nside,"NEST" if isnest else "RING")
  cnt=0
  while (cnt<ntry):
    cnt+=1
    inp=np.random.randint(low=0,high=12*nside*nside-1,size=vlen)
    out=base.xyf2pix(base.pix2xyf(inp))
    if (np.array_equal(inp,out)==False):
      raise ValueError("Test failed")

def check_vecangvec(vlen,ntry):
  cnt=0
  while (cnt<ntry):
    cnt+=1
    inp=random_ptg(vlen)
    out=ph.vec2ang(ph.ang2vec(inp))
    if (np.any(np.greater(np.abs(inp-out),1e-10))):
      raise ValueError("Test failed")

check_vecangvec(1000,1000)

for nside in (1,32,512,8192,32768*8):
  check_ringnestring(1000,1000,nside)
  for isnest in (False,True):
    check_vecpixvec(1000,1000,nside,isnest)
    check_pixangpix(1000,1000,nside,isnest)
    check_pixvecpix(1000,1000,nside,isnest)
    check_pixxyfpix(1000,1000,nside,isnest)
    check_pixangvecpix(1000,1000,nside,isnest)
    check_pixvecangpix(1000,1000,nside,isnest)

isnest=False
for nside in (3,7,514,8167,32768*8+7):
  check_vecpixvec(1000,1000,nside,isnest)
  check_pixangpix(1000,1000,nside,isnest)
  check_pixvecpix(1000,1000,nside,isnest)
  check_pixxyfpix(1000,1000,nside,isnest)
  check_pixangvecpix(1000,1000,nside,isnest)
  check_pixvecangpix(1000,1000,nside,isnest)