demo.py 2.32 KB
Newer Older
Martin Reinecke's avatar
Martin Reinecke committed
1
2
3
4
5
6
7
8
9
# Elementary demo for pysharp interface using a Gauss-Legendre grid
# 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.

import pysharp
import numpy as np
from numpy.testing import assert_allclose
Martin Reinecke's avatar
Martin Reinecke committed
10
from time import time
Martin Reinecke's avatar
Martin Reinecke committed
11
12

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

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

# 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
22
nlon = 4096
Martin Reinecke's avatar
Martin Reinecke committed
23
24
25
26
27
28
29
30
31

# create an object which will do the SHT work
job = pysharp.sharpjob_d()

# 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
32
# The symmetry implies that all coefficients with m==0 are purely real-valued.
Martin Reinecke's avatar
Martin Reinecke committed
33
34
# 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
35
36
37

# number of required a_lm coefficients
nalm = ((mmax+1)*(mmax+2))//2 + (mmax+1)*(lmax-mmax)
Martin Reinecke's avatar
Martin Reinecke committed
38
39
40
41
# 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
42
43
44
45
46
47
48

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

# describe the Gauss-Legendre geometry to the job
job.set_Gauss_geometry(nlat, nlon)

Martin Reinecke's avatar
Martin Reinecke committed
49
50
51
52
53
54
55
56
57
58
59
60
61
62
# 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
63
64
65

# make sure input was recovered accurately
assert_allclose(alm, alm2)