demo.py 1.88 KB
Newer Older
Martin Reinecke's avatar
Martin Reinecke committed
1
2
3
4
import interpol_ng
import numpy as np
import pysharp
import time
Martin Reinecke's avatar
Martin Reinecke committed
5
import matplotlib.pyplot as plt
Martin Reinecke's avatar
Martin Reinecke committed
6

Martin Reinecke's avatar
Martin Reinecke committed
7
np.random.seed(42)
Martin Reinecke's avatar
Martin Reinecke committed
8
9
10
11

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

Martin Reinecke's avatar
Martin Reinecke committed
12
13

def deltabeam(lmax,kmax):
Martin Reinecke's avatar
Martin Reinecke committed
14
    beam=np.zeros(nalm(lmax, kmax))+0j
Martin Reinecke's avatar
Martin Reinecke committed
15
16
17
18
19
    for l in range(lmax+1):
        beam[l] = np.sqrt((2*l+1.)/(4*np.pi))
    return beam


Martin Reinecke's avatar
Martin Reinecke committed
20
21
22
23
24
25
26
27
28
lmax=2047
mmax=lmax
kmax=2  # doesn't make any sense for the beam we are using, but just for demonstration purposes ...


# get random sky a_lm
# the a_lm arrays follow the same conventions as those in healpy
slmT =      np.random.uniform(-1., 1., nalm(lmax, mmax)) \
       + 1j*np.random.uniform(-1., 1., nalm(lmax, mmax))
Martin Reinecke's avatar
Martin Reinecke committed
29
30
31
32
33
34
35
# make a_lm with m==0 real-valued
slmT[0:lmax+1].imag = 0.

# build beam a_lm (pencil beam for now)
blmT = deltabeam(lmax,kmax)

t0=time.time()
Martin Reinecke's avatar
Martin Reinecke committed
36
# build interpolator object for slmT and blmT
Martin Reinecke's avatar
Martin Reinecke committed
37
38
foo = interpol_ng.PyInterpolator(slmT,blmT,lmax, kmax, epsilon=1e-6, nthreads=2)
print("setup time: ",time.time()-t0)
Martin Reinecke's avatar
Martin Reinecke committed
39

Martin Reinecke's avatar
Martin Reinecke committed
40
41
42
43
# assemble a set of (theta, phi, psi) pointings, at which we want to evaluate
# the interpolated signal.
# For testig purposes, we choose theta and phi such that they form a
# Gauss-Legendre grid of sufficient resolution to recover the input a_lm
Martin Reinecke's avatar
Martin Reinecke committed
44
45
nth = lmax+1
nph = 2*mmax+1
Martin Reinecke's avatar
Martin Reinecke committed
46
ptg = np.zeros((nth,nph,3))
Martin Reinecke's avatar
Martin Reinecke committed
47
48
49
50
51
th, _ = np.polynomial.legendre.leggauss(nth)
th = np.arccos(-th)
ptg[:,:,0] = th.reshape((-1,1))
ptg[:,:,1] = (2*np.pi*np.arange(nph)/nph).reshape((1,-1))
ptg[:,:,2] = 0
Martin Reinecke's avatar
Martin Reinecke committed
52
t0=time.time()
Martin Reinecke's avatar
Martin Reinecke committed
53
# do the actual interpolation
Martin Reinecke's avatar
Martin Reinecke committed
54
bar=foo.interpol(ptg.reshape((-1,3))).reshape((nth,nph))
Martin Reinecke's avatar
Martin Reinecke committed
55
print("interpolation time: ", time.time()-t0)
Martin Reinecke's avatar
Martin Reinecke committed
56

Martin Reinecke's avatar
Martin Reinecke committed
57
# get a_lm back from interpolation results
Martin Reinecke's avatar
Martin Reinecke committed
58
59
60
61
62
job = pysharp.sharpjob_d()
job.set_triangular_alm_info(lmax, mmax)
job.set_gauss_geometry(nth, nph)
alm2 = job.map2alm(bar.reshape((-1,)))

Martin Reinecke's avatar
Martin Reinecke committed
63
64
65
#compare with original a_lm
plt.plot(np.abs(alm2-slmT))
plt.suptitle("Deviations between original and reconstructed a_lm")
Martin Reinecke's avatar
Martin Reinecke committed
66
plt.show()