Skip to content
GitLab
Menu
Projects
Groups
Snippets
Loading...
Help
Help
Support
Community forum
Keyboard shortcuts
?
Submit feedback
Contribute to GitLab
Sign in
Toggle navigation
Menu
Open sidebar
Martin Reinecke
ducc
Commits
e05a563f
Commit
e05a563f
authored
Feb 29, 2020
by
Martin Reinecke
Browse files
improve demo
parent
2396eb78
Pipeline
#69985
passed with stages
in 8 minutes and 57 seconds
Changes
1
Pipelines
1
Hide whitespace changes
Inline
Side-by-side
interpol_ng/demo.py
View file @
e05a563f
...
@@ -10,25 +10,27 @@ import pysharp
...
@@ -10,25 +10,27 @@ import pysharp
import
time
import
time
np
.
random
.
seed
(
42
)
np
.
random
.
seed
(
42
)
lmax
=
999
lmax
=
2047
mmax
=
lmax
mmax
=
lmax
kmax
=
18
kmax
=
2
def
nalm
(
lmax
,
mmax
):
return
((
mmax
+
1
)
*
(
mmax
+
2
))
//
2
+
(
mmax
+
1
)
*
(
lmax
-
mmax
)
def
idx
(
l
,
m
):
def
idx
(
l
,
m
):
return
(
m
*
((
2
*
lmax
+
1
)
-
m
))
//
2
+
l
return
(
m
*
((
2
*
lmax
+
1
)
-
m
))
//
2
+
l
def
deltabeam
(
lmax
,
kmax
):
def
deltabeam
(
lmax
,
kmax
):
nalm_beam
=
((
kmax
+
1
)
*
(
kmax
+
2
))
//
2
+
(
kmax
+
1
)
*
(
lmax
-
kmax
)
beam
=
np
.
zeros
(
nalm
(
lmax
,
kmax
))
+
0j
beam
=
np
.
zeros
(
nalm_beam
)
+
0j
for
l
in
range
(
lmax
+
1
):
for
l
in
range
(
lmax
+
1
):
beam
[
l
]
=
np
.
sqrt
((
2
*
l
+
1.
)
/
(
4
*
np
.
pi
))
beam
[
l
]
=
np
.
sqrt
((
2
*
l
+
1.
)
/
(
4
*
np
.
pi
))
return
beam
return
beam
# number of required a_lm coefficients
nalm
=
((
mmax
+
1
)
*
(
mmax
+
2
))
//
2
+
(
mmax
+
1
)
*
(
lmax
-
mmax
)
# get random a_lm
# get random a_lm
slmT
=
np
.
random
.
uniform
(
-
1.
,
1.
,
nalm
)
+
1j
*
np
.
random
.
uniform
(
-
1.
,
1.
,
nalm
)
slmT
=
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
# make a_lm with m==0 real-valued
slmT
[
0
:
lmax
+
1
].
imag
=
0.
slmT
[
0
:
lmax
+
1
].
imag
=
0.
...
@@ -39,9 +41,9 @@ t0=time.time()
...
@@ -39,9 +41,9 @@ t0=time.time()
foo
=
interpol_ng
.
PyInterpolator
(
slmT
,
blmT
,
lmax
,
kmax
,
1e-5
)
foo
=
interpol_ng
.
PyInterpolator
(
slmT
,
blmT
,
lmax
,
kmax
,
1e-5
)
print
(
"setup: "
,
time
.
time
()
-
t0
)
print
(
"setup: "
,
time
.
time
()
-
t0
)
# evaluate total convolution on a sufficiently resolved
Clenshaw-Curtis
grid
# evaluate total convolution on a sufficiently resolved
Gauss-Legendre
grid
nth
=
1000
nth
=
lmax
+
1
nph
=
2
000
nph
=
2
*
mmax
+
1
ptg
=
np
.
zeros
((
nth
,
nph
,
3
))
ptg
=
np
.
zeros
((
nth
,
nph
,
3
))
th
,
wgt
=
np
.
polynomial
.
legendre
.
leggauss
(
nth
)
th
,
wgt
=
np
.
polynomial
.
legendre
.
leggauss
(
nth
)
th
=
np
.
arccos
(
th
)
th
=
np
.
arccos
(
th
)
...
...
Write
Preview
Markdown
is supported
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