demo.py 1.97 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

7
np.random.seed(48)
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

Martin Reinecke's avatar
Martin Reinecke committed
13
14
15
16
17
18
19
20
def random_alm(lmax, mmax):
    res = np.random.uniform(-1., 1., nalm(lmax, mmax)) \
     + 1j*np.random.uniform(-1., 1., nalm(lmax, mmax))
    # make a_lm with m==0 real-valued
    res[0:lmax+1].imag = 0.
    return res


Martin Reinecke's avatar
Martin Reinecke committed
21
def deltabeam(lmax,kmax):
Martin Reinecke's avatar
Martin Reinecke committed
22
    beam=np.zeros(nalm(lmax, kmax))+0j
Martin Reinecke's avatar
Martin Reinecke committed
23
24
25
26
    for l in range(lmax+1):
        beam[l] = np.sqrt((2*l+1.)/(4*np.pi))
    return beam

27
28
29
30
31
32
33
def convolve(alm1, alm2, lmax):
    job = pysharp.sharpjob_d()
    job.set_triangular_alm_info(lmax, lmax)
    job.set_gauss_geometry(lmax+1, 2*lmax+1)
    map = job.alm2map(alm1)*job.alm2map(alm2)
    job.set_triangular_alm_info(0,0)
    return job.map2alm(map)[0]*np.sqrt(4*np.pi)
Martin Reinecke's avatar
Martin Reinecke committed
34

35
lmax=10
Martin Reinecke's avatar
Martin Reinecke committed
36
mmax=lmax
37
kmax=lmax
Martin Reinecke's avatar
Martin Reinecke committed
38
39
40
41


# get random sky a_lm
# the a_lm arrays follow the same conventions as those in healpy
Martin Reinecke's avatar
Martin Reinecke committed
42
slmT = random_alm(lmax, mmax)
Martin Reinecke's avatar
Martin Reinecke committed
43
44
45

# build beam a_lm (pencil beam for now)
blmT = deltabeam(lmax,kmax)
46
blmT = random_alm(lmax, mmax)
Martin Reinecke's avatar
Martin Reinecke committed
47
48

t0=time.time()
Martin Reinecke's avatar
Martin Reinecke committed
49
# build interpolator object for slmT and blmT
Martin Reinecke's avatar
Martin Reinecke committed
50
51
foo = interpol_ng.PyInterpolator(slmT,blmT,lmax, kmax, epsilon=1e-6, nthreads=2)
print("setup time: ",time.time()-t0)
52
nth = 2*lmax+1
Martin Reinecke's avatar
Martin Reinecke committed
53
nph = 2*mmax+1
54

Martin Reinecke's avatar
Martin Reinecke committed
55
ptg = np.zeros((nth,nph,3))
56
ptg[:,:,0] = (np.pi*np.arange(nth)/(nth-1)).reshape((-1,1))
Martin Reinecke's avatar
Martin Reinecke committed
57
ptg[:,:,1] = (2*np.pi*np.arange(nph)/nph).reshape((1,-1))
Martin Reinecke's avatar
Martin Reinecke committed
58
ptg[:,:,2] = np.pi*0.2
Martin Reinecke's avatar
Martin Reinecke committed
59
t0=time.time()
Martin Reinecke's avatar
Martin Reinecke committed
60
# do the actual interpolation
Martin Reinecke's avatar
Martin Reinecke committed
61
bar=foo.interpol(ptg.reshape((-1,3))).reshape((nth,nph))
Martin Reinecke's avatar
Martin Reinecke committed
62
print("interpolation time: ", time.time()-t0)
63
64
65
66
67
68
69
70
71
72
73
plt.subplot(2,2,1)
plt.imshow(bar.reshape((nth,nph)))
bar2 = np.zeros((nth,nph))
for ith in range(nth):
    for iph in range(nph):
        rbeam=interpol_ng.rotate_alm(blmT, lmax, ptg[ith,iph,2],ptg[ith,iph,0],ptg[ith,iph,1])
        bar2[ith,iph] = convolve(slmT, rbeam, lmax)
plt.subplot(2,2,2)
plt.imshow(bar2)
plt.subplot(2,2,3)
plt.imshow((bar2-bar.reshape((nth,nph))))
Martin Reinecke's avatar
Martin Reinecke committed
74
plt.show()