Skip to content
GitLab
Menu
Projects
Groups
Snippets
/
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
Simon Perkins
ducc
Commits
e0c87ef7
Commit
e0c87ef7
authored
Feb 14, 2020
by
Martin Reinecke
Browse files
enhance demo
parent
6a5dd0cf
Changes
1
Hide whitespace changes
Inline
Sidebyside
pysharp/demo.py
View file @
e0c87ef7
...
...
@@ 7,9 +7,10 @@
import
pysharp
import
numpy
as
np
from
numpy.testing
import
assert_allclose
from
time
import
time
# set maximum multipole moment
lmax
=
4095
lmax
=
2047
# maximum m. For SHTOOLS this is alway equal to lmax, if I understand correctly.
mmax
=
lmax
...
...
@@ 18,7 +19,7 @@ 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.
nlon
=
8192
nlon
=
4096
# create an object which will do the SHT work
job
=
pysharp
.
sharpjob_d
()
...
...
@@ 34,13 +35,10 @@ job = pysharp.sharpjob_d()
nalm
=
((
mmax
+
1
)
*
(
mmax
+
2
))
//
2
+
(
mmax
+
1
)
*
(
lmax

mmax
)
# number of realvalued random numbers to draw
nalm_r
=
nalm
*
2

lmax

1
# get random numbers
alm_r
=
np
.
random
.
uniform
(

1.
,
1.
,
nalm_r
)
# create the complexvalued a_lm array
alm
=
np
.
empty
(
nalm
,
dtype
=
np
.
complex128
)
alm
[
0
:
lmax
+
1
]
=
alm_r
[
0
:
lmax
+
1
]
alm
[
lmax
+
1
:]
=
np
.
sqrt
(
0.5
)
*
(
alm_r
[
lmax
+
1
::
2
]
+
1j
*
alm_r
[
lmax
+
2
::
2
])
# get random a_lm
alm
=
np
.
random
.
uniform
(

1.
,
1.
,
nalm
)
+
1j
*
np
.
random
.
uniform
(

1.
,
1.
,
nalm
)
# make a_lm with m==0 realvalued
alm
[
0
:
lmax
+
1
].
imag
=
0.
# describe the a_lm array to the job
job
.
set_triangular_alm_info
(
lmax
,
mmax
)
...
...
@@ 48,8 +46,20 @@ job.set_triangular_alm_info(lmax, mmax)
# describe the GaussLegendre geometry to the job
job
.
set_Gauss_geometry
(
nlat
,
nlon
)
# go from a_lm to map and back
alm2
=
job
.
map2alm
(
job
.
alm2map
(
alm
))
# 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 realvalued 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 isolatitude 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
assert_allclose
(
alm
,
alm2
)
Write
Preview
Supports
Markdown
0%
Try again
or
attach a new file
.
Attach a file
Cancel
You are about to add
0
people
to the discussion. Proceed with caution.
Finish editing this message first!
Cancel
Please
register
or
sign in
to comment