sht_demo.py 3.21 KB
Newer Older
Martin Reinecke's avatar
Martin Reinecke committed
1
# Elementary demo for the ducc_0_1.sht interface using a Gauss-Legendre grid
Martin Reinecke's avatar
Martin Reinecke committed
2
3
4
5
6
# I'm not sure I have a perfect equivalent for the DH grid(s) at the moment,
# since they apparently do not include the South Pole. The Clenshaw-Curtis
# and Fejer quadrature rules are very similar (see the documentation in
# sharp_geomhelpers.h). An exact analogon to DH can be added easily, I expect.

Martin Reinecke's avatar
Martin Reinecke committed
7
import ducc_0_1.sht as sht
Martin Reinecke's avatar
Martin Reinecke committed
8
import numpy as np
Martin Reinecke's avatar
Martin Reinecke committed
9
from time import time
Martin Reinecke's avatar
Martin Reinecke committed
10

Martin Reinecke's avatar
stage 2  
Martin Reinecke committed
11
12
13
def _l2error(a, b):
    return np.sqrt(np.sum(np.abs(a-b)**2)/np.sum(np.abs(a)**2))

Martin Reinecke's avatar
Martin Reinecke committed
14
# set maximum multipole moment
Martin Reinecke's avatar
Martin Reinecke committed
15
lmax = 2047
Martin Reinecke's avatar
Martin Reinecke committed
16
17
18
19
20
# maximum m. For SHTOOLS this is alway equal to lmax, if I understand correctly.
mmax = lmax

# Number of pixels per ring. Must be >=2*lmax+1, but I'm choosing a larger
# number for which the FFT is faster.
Martin Reinecke's avatar
Martin Reinecke committed
21
nlon = 4096
Martin Reinecke's avatar
Martin Reinecke committed
22
23

# create an object which will do the SHT work
Martin Reinecke's avatar
Martin Reinecke committed
24
job = sht.sharpjob_d()
Martin Reinecke's avatar
Martin Reinecke committed
25
26
27
28
29
30

# create a set of spherical harmonic coefficients to transform
# Libsharp works exclusively on real-valued maps. The corresponding harmonic
# coefficients are termed a_lm; they are complex numbers with 0<=m<=lmax and
# m<=l<=lmax.
# Symmetry: a_l,-m = (-1)**m*conj(a_l,m).
Martin Reinecke's avatar
Martin Reinecke committed
31
# The symmetry implies that all coefficients with m==0 are purely real-valued.
Martin Reinecke's avatar
Martin Reinecke committed
32
33
# The a_lm are stored in a 1D complex-valued array, in the following order:
# a_(0,0), a(1,0), ..., a_(lmax,0), a(1,1), a(2,1), ... a(lmax,1), ..., a(lmax, mmax)
Martin Reinecke's avatar
Martin Reinecke committed
34
35
36

# number of required a_lm coefficients
nalm = ((mmax+1)*(mmax+2))//2 + (mmax+1)*(lmax-mmax)
Martin Reinecke's avatar
Martin Reinecke committed
37
# get random a_lm
38
39
rng = np.random.default_rng(42)
alm = rng.uniform(-1., 1., nalm) + 1j*rng.uniform(-1., 1., nalm)
Martin Reinecke's avatar
Martin Reinecke committed
40
41
# make a_lm with m==0 real-valued
alm[0:lmax+1].imag = 0.
Martin Reinecke's avatar
Martin Reinecke committed
42
43
44
45

# describe the a_lm array to the job
job.set_triangular_alm_info(lmax, mmax)

46
47
48
49
50
51

print("testing Gauss-Legendre grid")

# Number of iso-latitude rings required for Gauss-Legendre grid
nlat = lmax+1

Martin Reinecke's avatar
Martin Reinecke committed
52
# describe the Gauss-Legendre geometry to the job
Martin Reinecke's avatar
stage 2  
Martin Reinecke committed
53
job.set_gauss_geometry(nlat, nlon)
Martin Reinecke's avatar
Martin Reinecke committed
54

Martin Reinecke's avatar
Martin Reinecke committed
55
56
57
58
59
60
61
62
63
64
65
66
67
68
# go from a_lm to map
t0=time()
map = job.alm2map(alm)
print("time for map synthesis: {}s".format(time()-t0))

# map is a 1D real-valued array with (nlat*nlon) entries. It can be reshaped
# to (nlat, nlon) for plotting.
# Libsharp woks on "1D" maps because it apso supports pixelizations that varying
# number of pixels on each iso-latitude ring, which cannot be represented by 2D
# arrays (e.g. Healpix)

t0=time()
alm2 = job.map2alm(map)
print("time for map analysis: {}s".format(time()-t0))
Martin Reinecke's avatar
Martin Reinecke committed
69
70

# make sure input was recovered accurately
Martin Reinecke's avatar
stage 2  
Martin Reinecke committed
71
print("L2 error: ", _l2error(alm,alm2))
72
73
74
75


print("testing Driscoll-Healy grid")

Martin Reinecke's avatar
stage 2  
Martin Reinecke committed
76
77
# Number of iso-latitude rings required for Driscoll-Healy grid
nlat = 2*lmax+2
78
79

# describe the Gauss-Legendre geometry to the job
Martin Reinecke's avatar
stage 2  
Martin Reinecke committed
80
job.set_dh_geometry(nlat, nlon)
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97

# go from a_lm to map
t0=time()
map = job.alm2map(alm)
print("time for map synthesis: {}s".format(time()-t0))

# map is a 1D real-valued array with (nlat*nlon) entries. It can be reshaped
# to (nlat, nlon) for plotting.
# Libsharp woks on "1D" maps because it apso supports pixelizations that varying
# number of pixels on each iso-latitude ring, which cannot be represented by 2D
# arrays (e.g. Healpix)

t0=time()
alm2 = job.map2alm(map)
print("time for map analysis: {}s".format(time()-t0))

# make sure input was recovered accurately
Martin Reinecke's avatar
stage 2  
Martin Reinecke committed
98
print("L2 error: ", _l2error(alm,alm2))