demo.py 2.3 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(20)
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

Martin Reinecke's avatar
Martin Reinecke committed
35
lmax=40
Martin Reinecke's avatar
Martin Reinecke committed
36
mmax=lmax
Martin Reinecke's avatar
Martin Reinecke committed
37
kmax=10
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

Martin Reinecke's avatar
Martin Reinecke committed
44
# build beam a_lm
Martin Reinecke's avatar
Martin Reinecke committed
45
blmT = random_alm(lmax, kmax)
Martin Reinecke's avatar
Martin Reinecke committed
46
47

t0=time.time()
Martin Reinecke's avatar
Martin Reinecke committed
48
# build interpolator object for slmT and blmT
Martin Reinecke's avatar
Martin Reinecke committed
49
foo = interpol_ng.PyInterpolator(slmT,blmT,lmax, kmax, epsilon=1e-6, nthreads=1)
Martin Reinecke's avatar
Martin Reinecke committed
50
print("setup time: ",time.time()-t0)
51
nth = 2*lmax+1
Martin Reinecke's avatar
Martin Reinecke committed
52
nph = 2*mmax+1
Martin Reinecke's avatar
Martin Reinecke committed
53
#exit()
Martin Reinecke's avatar
Martin Reinecke committed
54
ptg = np.zeros((nth,nph,3))
55
ptg[:,:,0] = (np.pi*np.arange(nth)/(nth-1)).reshape((-1,1))
Martin Reinecke's avatar
Martin Reinecke committed
56
ptg[:,:,1] = (2*np.pi*np.arange(nph)/nph).reshape((1,-1))
Martin Reinecke's avatar
Martin Reinecke committed
57
ptg[:,:,2] = np.pi*0.2
Martin Reinecke's avatar
Martin Reinecke committed
58
t0=time.time()
Martin Reinecke's avatar
Martin Reinecke committed
59
# do the actual interpolation
Martin Reinecke's avatar
Martin Reinecke committed
60
bar=foo.interpol(ptg.reshape((-1,3))).reshape((nth,nph))
Martin Reinecke's avatar
Martin Reinecke committed
61
print("interpolation time: ", time.time()-t0)
62
63
64
plt.subplot(2,2,1)
plt.imshow(bar.reshape((nth,nph)))
bar2 = np.zeros((nth,nph))
Martin Reinecke's avatar
Martin Reinecke committed
65
66
blmTfull = np.zeros(slmT.size)+0j
blmTfull[0:blmT.size] = blmT
67
68
for ith in range(nth):
    for iph in range(nph):
Martin Reinecke's avatar
Martin Reinecke committed
69
        rbeam=interpol_ng.rotate_alm(blmTfull, lmax, ptg[ith,iph,2],ptg[ith,iph,0],ptg[ith,iph,1])
Martin Reinecke's avatar
Martin Reinecke committed
70
        bar2[ith,iph] = convolve(slmT, rbeam, lmax).real
71
72
73
74
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
75
plt.show()
Martin Reinecke's avatar
Martin Reinecke committed
76
77
78
79
80
81
82
83
84
85

fake = np.random.uniform(-1.,1., bar.size)
foo2 = interpol_ng.PyInterpolator(lmax, kmax, epsilon=1e-6, nthreads=2)
foo2.deinterpol(ptg.reshape((-1,3)), fake)
bla=foo2.getSlm(blmT)
print(np.vdot(slmT,bla))
slmT[lmax+1:]*=np.sqrt(2)
bla[lmax+1:]*=np.sqrt(2)
print(np.vdot(slmT,bla))
print(np.vdot(fake,bar))