usage.py 2.1 KB
Newer Older
Martin Reinecke's avatar
comment    
Martin Reinecke committed
1
2
# Short usage demo for the interpol_ng module

Martin Reinecke's avatar
Martin Reinecke committed
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
import pyinterpol_ng
import numpy as np


# establish a random number generator
rng = np.random.default_rng(np.random.SeedSequence(42))


def nalm(lmax, mmax):
    return ((mmax+1)*(mmax+2))//2 + (mmax+1)*(lmax-mmax)


def random_alm(lmax, mmax, ncomp):
    res = rng.uniform(-1., 1., (nalm(lmax, mmax), ncomp)) \
     + 1j*rng.uniform(-1., 1., (nalm(lmax, mmax), ncomp))
    # make a_lm with m==0 real-valued
    res[0:lmax+1,:].imag = 0.
    return res


# simulation parameters
lmax=1024 # highest l and m moment for sky, highest l moment for beam
kmax=13 # highest m moment for beam
ncomp=3 # T, E and B

# get random sky a_lm
# the a_lm arrays follow the same conventions as those in healpy
slm = random_alm(lmax, lmax, ncomp)

# build beam a_lm
blm = random_alm(lmax, kmax, ncomp)

# build random pointings
npnt = 1000000 # play with this to measure performance
ptg = rng.uniform(0., 1., (npnt,3))
ptg[:,0]*=np.pi   # theta
ptg[:,1]*=2*np.pi # phi
ptg[:,2]*=2*np.pi # psi

# build the interpolator
# For a "classic" CMB experiment we want to combine the T, E and B signals,
# so we set `separate` to `False`

print("classic interpolator setup...")
inter_classic = pyinterpol_ng.PyInterpolator(
    slm,blm,separate=False,lmax=lmax, kmax=kmax, epsilon=1e-4, nthreads=2)
print("...done")

# get interpolated values
print("interpolating...")
res = inter_classic.interpol(ptg)
print("...done")

# res is an array of shape (nptg, 1).
# If we had specified `separate=True`, it would be of shape(nptg, 3).
print(res.shape)

Martin Reinecke's avatar
Martin Reinecke committed
60
61
# Since the interpolator object holds large data structures, it should be
# deleted once it is no longer needed
Martin Reinecke's avatar
Martin Reinecke committed
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
del inter_classic

# Now the same thing for an experiment with HWP. In this case we need the
# interpolated T, E and B signals separate.
separate = True

print("HWP interpolator setup...")
inter_hwp = pyinterpol_ng.PyInterpolator(
    slm,blm,separate=True,lmax=lmax, kmax=kmax, epsilon=1e-4, nthreads=2)
print("...done")

# get interpolated values
print("interpolating...")
res = inter_hwp.interpol(ptg)
print("...done")

# now res has shape(nptg,3)
print(res.shape)