demo.py 1.52 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(20909)
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
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
20
21
22
23
24
25
26
27
28
29
30
31
def compress_alm(alm,lmax):
    res = np.empty(2*len(alm)-lmax-1, dtype=np.float64)
    res[0:lmax+1] = alm[0:lmax+1].real
    res[lmax+1::2] = np.sqrt(2)*alm[lmax+1:].real
    res[lmax+2::2] = np.sqrt(2)*alm[lmax+1:].imag
    return res


def myalmdot(a1,a2,lmax,mmax,spin):
    return np.vdot(compress_alm(a1,lmax),compress_alm(a2,lmax))
def mydot(a1,a2,spin):
    return np.vdot(theta_extend(a1,spin),theta_extend(a2,spin))
Martin Reinecke's avatar
Martin Reinecke committed
32

33
lmax=512
Martin Reinecke's avatar
Martin Reinecke committed
34
mmax=lmax
35
kmax=12
Martin Reinecke's avatar
Martin Reinecke committed
36
37
38
39


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

Martin Reinecke's avatar
Martin Reinecke committed
42
# build beam a_lm
Martin Reinecke's avatar
Martin Reinecke committed
43
blmT = random_alm(lmax, kmax)
44
ptg = np.array([[1.2,0.01,1.],[1.4,0.7,2.]])
Martin Reinecke's avatar
Martin Reinecke committed
45
foo = interpol_ng.PyInterpolator(slmT,blmT,lmax, kmax, epsilon=1e-6, nthreads=1)
46
47
48
bar=foo.interpol(ptg)
print(foo.Nphi(),foo.Nphi0())
fake = np.random.uniform(0.,1., ptg.shape[0])
Martin Reinecke's avatar
Martin Reinecke committed
49
50
51
foo2 = interpol_ng.PyInterpolator(lmax, kmax, epsilon=1e-6, nthreads=2)
foo2.deinterpol(ptg.reshape((-1,3)), fake)
bla=foo2.getSlm(blmT)
52
print(myalmdot(slmT, np.conj(bla), lmax, lmax, 0))
Martin Reinecke's avatar
Martin Reinecke committed
53
print(np.vdot(fake,bar))
54
print(myalmdot(slmT, np.conj(bla), lmax, lmax, 0)/np.vdot(fake,bar))