sht_demo.py 3.19 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
38
39
40
# get random a_lm
alm = np.random.uniform(-1., 1., nalm) + 1j*np.random.uniform(-1., 1., nalm)
# make a_lm with m==0 real-valued
alm[0:lmax+1].imag = 0.
Martin Reinecke's avatar
Martin Reinecke committed
41
42
43
44

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

45
46
47
48
49
50

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
51
# describe the Gauss-Legendre geometry to the job
Martin Reinecke's avatar
stage 2    
Martin Reinecke committed
52
job.set_gauss_geometry(nlat, nlon)
Martin Reinecke's avatar
Martin Reinecke committed
53

Martin Reinecke's avatar
Martin Reinecke committed
54
55
56
57
58
59
60
61
62
63
64
65
66
67
# 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
68
69

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


print("testing Driscoll-Healy grid")

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

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

# 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
97
print("L2 error: ", _l2error(alm,alm2))