Skip to content

Commit

Permalink
switch cos/sin for the sign of m in spherical harmonic function and d…
Browse files Browse the repository at this point in the history 
…erivatives, update testing scripts
  • Loading branch information
Mark Baum committed Dec 15, 2020
1 parent 8a9ff46 commit e26c458
Show file tree
Hide file tree
Showing 9 changed files with 140 additions and 75 deletions.
7 changes: 6 additions & 1 deletion docs/_static/pygments.css
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,10 @@
pre { line-height: 125%; }
td.linenos pre { color: #000000; background-color: #f0f0f0; padding-left: 5px; padding-right: 5px; }
span.linenos { color: #000000; background-color: #f0f0f0; padding-left: 5px; padding-right: 5px; }
td.linenos pre.special { color: #000000; background-color: #ffffc0; padding-left: 5px; padding-right: 5px; }
span.linenos.special { color: #000000; background-color: #ffffc0; padding-left: 5px; padding-right: 5px; }
.highlight .hll { background-color: #ffffcc }
.highlight { background: #f8f8f8; }
.highlight { background: #f8f8f8; }
.highlight .c { color: #8f5902; font-style: italic } /* Comment */
.highlight .err { color: #a40000; border: 1px solid #ef2929 } /* Error */
.highlight .g { color: #000000 } /* Generic */
Expand Down
10 changes: 5 additions & 5 deletions docs/spherical_harmonic.html
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,9 +44,9 @@
\begin{equation}
Y_n^m(\theta,\phi) =
\begin{cases}
\sqrt{2} \cos(m\phi) P_n^{|m|}(\cos\theta) & m < 0 \\
\sqrt{2} \sin(m\phi) P_n^{|m|}(\cos\theta) & m < 0 \\
P_n^{m}(\cos\theta) & m = 0 \\
\sqrt{2} \sin(m\phi) P_n^{m}(\cos\theta) & m > 0
\sqrt{2} \cos(m\phi) P_n^{m}(\cos\theta) & m > 0
\end{cases}
\end{equation}</div><p>where <span class="math notranslate nohighlight">\(P_n^m\)</span> is an associated Legendre polynomial (<code class="xref py py-func docutils literal notranslate"><span class="pre">legen_theta</span></code>). The harmonics here are orthonormal.</p>
<p>This module contains functions for evaluating the <a class="reference internal" href="#orthopoly.spherical_harmonic.sph_har" title="orthopoly.spherical_harmonic.sph_har"><code class="xref py py-func docutils literal notranslate"><span class="pre">spherical</span> <span class="pre">harmonics</span></code></a>, their <a class="reference internal" href="#orthopoly.spherical_harmonic.grad_sph_har" title="orthopoly.spherical_harmonic.grad_sph_har"><code class="xref py py-func docutils literal notranslate"><span class="pre">gradients</span></code></a>, and their <a class="reference internal" href="#orthopoly.spherical_harmonic.lap_sph_har" title="orthopoly.spherical_harmonic.lap_sph_har"><code class="xref py py-func docutils literal notranslate"><span class="pre">laplacians</span></code></a>. It doesn’t contain functions to perform transforms from values on the sphere to spherical harmonic expansion coefficients. The module also contains some functions for generating grids in spherical space and for generating random spherical harmonic expansions with specific power density properties (<a class="reference internal" href="#orthopoly.spherical_harmonic.noise" title="orthopoly.spherical_harmonic.noise"><code class="xref py py-func docutils literal notranslate"><span class="pre">noise</span></code></a>). The <a class="reference internal" href="#orthopoly.spherical_harmonic.Expansion" title="orthopoly.spherical_harmonic.Expansion"><code class="xref py py-class docutils literal notranslate"><span class="pre">Expansion</span></code></a> class stores harmonic coefficients, evaluates them, and may be multiplied/divided by scalars. For evaluating only a few expansions at a few sets of points, using the <a class="reference internal" href="#orthopoly.spherical_harmonic.sph_har" title="orthopoly.spherical_harmonic.sph_har"><code class="xref py py-func docutils literal notranslate"><span class="pre">sph_har</span></code></a> function or <a class="reference internal" href="#orthopoly.spherical_harmonic.Expansion" title="orthopoly.spherical_harmonic.Expansion"><code class="xref py py-class docutils literal notranslate"><span class="pre">Expansion</span></code></a> class should be fine. For evaluating many expansions (of the same size) at the same set of points, consider using the <a class="reference internal" href="#orthopoly.spherical_harmonic.sph_har_matrix" title="orthopoly.spherical_harmonic.sph_har_matrix"><code class="xref py py-func docutils literal notranslate"><span class="pre">sph_har_matrix</span></code></a> function.</p>
Expand Down Expand Up @@ -228,9 +228,9 @@
\begin{equation}
Y_n^m(\theta,\phi) =
\begin{cases}
\sqrt{2} \cos(m\phi) P_n^{|m|}(\cos\theta) &amp; m &lt; 0 \\
\sqrt{2} \sin(m\phi) P_n^{|m|}(\cos\theta) &amp; m &lt; 0 \\
P_n^{m}(\cos\theta) &amp; m = 0 \\
\sqrt{2} \sin(m\phi) P_n^{m}(\cos\theta) &amp; m &gt; 0
\sqrt{2} \cos(m\phi) P_n^{m}(\cos\theta) &amp; m &gt; 0
\end{cases}
\end{equation}</div><dl class="simple">
<dt>where <span class="math notranslate nohighlight">\(P_n^m\)</span> is the normalized associated Legendre function implemented in this module as <code class="xref py py-func docutils literal notranslate"><span class="pre">legen_hat</span></code> or <code class="xref py py-func docutils literal notranslate"><span class="pre">legen_theta</span></code>. The normalization ensures that the functions are orthonormal. More info can be found in the references below, among many other places</dt><dd><ul class="simple">
Expand Down Expand Up @@ -624,4 +624,4 @@ <h3 id="searchlabel">Quick search</h3>


</body>
</html>
</html>
135 changes: 80 additions & 55 deletions orthopoly/spherical_harmonic.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,9 +11,9 @@
\begin{equation}
Y_n^m(\theta,\phi) =
\begin{cases}
\sqrt{2} \cos(m\phi) P_n^{|m|}(\cos\theta) & m < 0 \\
\sqrt{2} \sin(m\phi) P_n^{|m|}(\cos\theta) & m < 0 \\
P_n^{m}(\cos\theta) & m = 0 \\
\sqrt{2} \sin(m\phi) P_n^{m}(\cos\theta) & m > 0
\sqrt{2} \cos(m\phi) P_n^{m}(\cos\theta) & m > 0
\end{cases}
\end{equation}
Expand Down Expand Up @@ -246,9 +246,9 @@ def sph_har(t, p, n, m):
\begin{equation}
Y_n^m(\theta,\phi) =
\begin{cases}
\sqrt{2} \cos(m\phi) P_n^{|m|}(\cos\theta) & m < 0 \\
\sqrt{2} \sin(m\phi) P_n^{|m|}(\cos\theta) & m < 0 \\
P_n^{m}(\cos\theta) & m = 0 \\
\sqrt{2} \sin(m\phi) P_n^{m}(\cos\theta) & m > 0
\sqrt{2} \cos(m\phi) P_n^{m}(\cos\theta) & m > 0
\end{cases}
\end{equation}
Expand All @@ -268,9 +268,9 @@ def sph_har(t, p, n, m):
_check_sph_har_args(t, p, n, m)

#handle different cases for m
if m < 0: Y = sqrt(2.0)*cos(m*p)
if m < 0: Y = sqrt(2.0)*sin(m*p)
elif m == 0: Y = 1.0 + 0.0*p #is one by default, but the same shape as p
else: Y = sqrt(2.0)*sin(m*p)
else: Y = sqrt(2.0)*cos(m*p)
#apply the associated Legendre function
Y *= legen_theta(t, n, abs(m))

Expand Down Expand Up @@ -316,14 +316,14 @@ def grad_sph_har(t, p, n, m, R=1):

#handle different cases for m
if m < 0:
dt *= sqrt(2.0)*cos(m*p)
dp *= -sqrt(2.0)*m*sin(m*p)
dt *= sqrt(2.0)*sin(m*p)
dp *= sqrt(2.0)*m*cos(m*p)
elif m == 0:
dt *= 1.0 + 0.0*p #forece the same shape as p
dt *= 1.0 + 0.0*p #force the same shape as p
dp *= 0.0*p #force the same shape as p
else:
dt *= sqrt(2.0)*sin(m*p)
dp *= sqrt(2.0)*m*cos(m*p)
dt *= sqrt(2.0)*cos(m*p)
dp *= -sqrt(2.0)*m*sin(m*p)

#apply associated Legendre functions
dt *= dlegen_theta(t, n, abs(m))
Expand Down Expand Up @@ -352,14 +352,14 @@ def lap_sph_har(t, p, n, m, R=1):

#handle different cases for m
if m < 0:
ddt *= cos(m*p)*sqrt(2.0)*(sin(t)*ddlegen_theta(t,n,abs(m)) + cos(t)*dlegen_theta(t,n,abs(m)))
ddp *= -sqrt(2.0)*legen_theta(t,n,abs(m))*(m**2)*cos(m*p)
ddt *= sin(m*p)*sqrt(2.0)*(sin(t)*ddlegen_theta(t,n,abs(m)) + cos(t)*dlegen_theta(t,n,abs(m)))
ddp *= -sqrt(2.0)*legen_theta(t,n,abs(m))*(m**2)*sin(m*p)
elif m == 0:
ddt *= sin(t)*ddlegen_theta(t,n,m) + cos(t)*dlegen_theta(t,n,m)
ddp *= 0.0*t #force the same shape as t
else:
ddt *= sin(m*p)*sqrt(2.0)*(sin(t)*ddlegen_theta(t,n,m) + cos(t)*dlegen_theta(t,n,m))
ddp *= -sqrt(2.0)*legen_theta(t,n,m)*(m**2)*sin(m*p)
ddt *= cos(m*p)*sqrt(2.0)*(sin(t)*ddlegen_theta(t,n,m) + cos(t)*dlegen_theta(t,n,m))
ddp *= -sqrt(2.0)*legen_theta(t,n,m)*(m**2)*cos(m*p)

#sum the components
lap = ddt + ddp
Expand Down Expand Up @@ -515,12 +515,11 @@ def _find_el(x, y, z):
el - list of length 3 index arrays indicating vertices of each el"""

try:
tri = SphericalVoronoi(np.stack((x, y, z)).T)._tri
sv = SphericalVoronoi(np.stack((x, y, z)).T)
except MemoryError:
raise MemoryError('There are too many vertices to perform Delaunay triangulation with available memory. You may need to use a bigger computer.')
else:
el = tri.simplices.astype(np.uint32)
return(el)
return(sv._simplices.astype(np.uint32))

def _subdivide_elements(t, p, el):
"""Subdivide each triangular element in the mesh into four new triangles
Expand Down Expand Up @@ -552,46 +551,72 @@ def _subdivide_elements(t, p, el):
return(t, p, el)

def _icosahedron_base(R=1):
"""Generate the triangles forming an icosahedron, returned in cartesian
"""Generates the triangles forming an icosahedron, returned in cartesian
or spherical coordinates
optional args:
R - radius of sphere to put vertices on
returns:
x - x coordinates of icosahedron vertices
y - y coordinates of icosahedron vertices
z - z coordinates of icosahedron vertices
el - list of triangular element indices"""
#golden ratio
PHI = (1.0 + sqrt(5.0))/2.0
:param float R: radius of sphere to put vertices on
#icosahedron vertices (12 of them)
i, j, k = zip(*product([-1, 1], [-PHI, PHI], [0]))
pts = list(zip(i,j,k)) + list(zip(j,k,i)) + list(zip(k,i,j))
x, y, z = zip(*pts)
x, y, z = array(x), array(y), array(z)

#assemble edges (30 of them)
D = zeros((12,12))
for i in range(12):
for j in range(i+1,12):
D[i,j] = sqrt((x[i] - x[j])**2 + (y[i] - y[j])**2 + (z[i] - z[j])**2)
edg = [array(i, dtype=int_) for i in zip(*np.nonzero(isclose(D,2)))]

#assemble elements (20 of them)
el = set()
for i in range(30):
for j in range(i+1,30):
for k in range(j+1,30):
s = set(edg[i]) | set(edg[j]) | set(edg[k])
if len(s) == 3:
el.add(tuple(s))
el = [array(e, dtype=int_) for e in el]

#normalize vertices for the proper radius
x *= R/sqrt(1 + PHI**2)
y *= R/sqrt(1 + PHI**2)
z *= R/sqrt(1 + PHI**2)
:return: tuple containing
- array of x coordinates of icosahedron vertices
- array of y coordinates of icosahedron vertices
- array of z coordinates of icosahedron vertices
- list of length 3 arrays of triangular element indices
- list of length 3 arrays of element neighbor indices"""

#golden ratio
Φ = (1.0 + sqrt(5.0))/2.0

x = array([-1, -1, 1, 1, -Φ, Φ, -Φ, Φ, 0, 0, 0, 0])
y = array([-Φ, Φ, -Φ, Φ, 0, 0, 0, 0, -1, -1, 1, 1])
z = array([0, 0, 0, 0, -1, -1, 1, 1, -Φ, Φ, -Φ, Φ])
el = array([
[0, 2, 9],
[1, 4, 10],
[0, 2, 8],
[2, 7, 9],
[7, 9, 11],
[3, 5, 7],
[5, 8, 10],
[2, 5, 8],
[3, 5, 10],
[6, 9, 11],
[1, 6, 11],
[2, 5, 7],
[3, 7, 11],
[1, 3, 11],
[0, 4, 8],
[4, 8, 10],
[1, 4, 6],
[0, 4, 6],
[0, 6, 9],
[1, 3, 10]])
ne = array([
[ 2, 3, 18],
[15, 16, 19],
[ 0, 7, 14],
[ 0, 4, 11],
[ 3, 9, 12],
[ 8, 11, 12],
[ 7, 8, 15],
[ 2, 6, 11],
[ 5, 6, 19],
[ 4, 10, 18],
[ 9, 13, 16],
[ 3, 5, 7],
[ 4, 5, 13],
[10, 12, 19],
[ 2, 15, 17],
[ 1, 6, 14],
[ 1, 10, 17],
[14, 16, 18],
[ 0, 9, 17],
[ 1, 8, 13]])

#normalize for proper radius
x *= R/sqrt(1 + Φ**2)
y *= R/sqrt(1 + Φ**2)
z *= R/sqrt(1 + Φ**2)

return(x, y, z, el)

Expand Down
4 changes: 3 additions & 1 deletion test.sh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,8 +8,10 @@ cd tests/chebyshev
python test_bvp.py
python test_interpolation.py

cd ../spherical-harmonic
cd ../legendre
python test_associated_legendre.py

cd ../spherical-harmonic
python test_spherical_harmonic.py
python test_noise.py

Expand Down
4 changes: 4 additions & 0 deletions tests/gegenbauer/plot_gegenbauer.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,10 @@
sys.path.insert(0, join('..', '..'))
from orthopoly.gegenbauer import *

plt.style.use('dark_background')
#plt.rc('text', usetex=True)
#plt.rc('font', family='serif')

xhat = linspace(-1, 1, 1000)
ylims = [(-3,4), (-8,12),(-20,25)]

Expand Down
4 changes: 2 additions & 2 deletions tests/legendre/test_associated_legendre.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,7 +52,7 @@ def test_dlegen_theta(nmax, N=100, h=1e-3):
err[n,m] = np.abs(dP - dP_fd).max()
#print('(%d,%d) %g' % (n,m,err[n,m]))
fig, ax = plt.subplots(1,1)
r = ax.pcolormesh(range(nmax+1), range(nmax+1), np.log10(err))
r = ax.pcolormesh(range(nmax+2), range(nmax+2), np.log10(err))
cb = plt.colorbar(r, ax=ax)
cb.set_label('log$_{10}$(Maximum Absolute Error)')
ax.set_ylabel("$n$")
Expand Down Expand Up @@ -86,7 +86,7 @@ def test_ddlegen_theta(nmax, N=100, h=1e-2):
err[n,m] = np.abs(ddP - ddP_fd).max()
#print('(%d,%d) %g' % (n,m,err[n,m]))
fig, ax = plt.subplots(1,1)
r = ax.pcolormesh(range(nmax+1), range(nmax+1), np.log10(err))
r = ax.pcolormesh(range(nmax+2), range(nmax+2), np.log10(err))
cb = plt.colorbar(r, ax=ax)
cb.set_label('log$_{10}$(Maximum Absolute Error)')
ax.set_ylabel("$n$")
Expand Down
27 changes: 22 additions & 5 deletions tests/spherical-harmonic/test_analysis.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,11 +7,15 @@
sys.path.insert(0, join('..', '..'))
from orthopoly.spherical_harmonic import *

#plt.style.use('dark_background')
#plt.rc('text', usetex=True)
#plt.rc('font', family='serif')

#-------------------------------------------------------------------------------
# INPUT

#degree of truncation (should be an even number)
T = 16
T = 20

#number of grid points for plotting
nphi = 500
Expand Down Expand Up @@ -42,25 +46,38 @@
ex = Expansion(a, yn, ym)

#make a grid for sampling
ntheta = nphi//2
pg = linspace(0, 2*pi, nphi)
tg = linspace(0, pi, nphi//2)
tg = linspace(0, pi, ntheta)
Pg, Tg = meshgrid(pg, tg)

#plot results
fig, (axa, axb, axc) = plt.subplots(3,1)

r = axa.pcolormesh(Pg, Tg, f_z(Tg, Pg))
r = axa.pcolormesh(
linspace(0, 2*pi, nphi+1),
linspace(0, pi, ntheta+1),
f_z(Tg, Pg))
plt.colorbar(r, ax=axa)
axa.scatter(p, t, c='k', s=1)
axa.set_title('Exact Function')

r = axb.pcolormesh(Pg, Tg, ex(Tg, Pg))
r = axb.pcolormesh(
linspace(0, 2*pi, nphi+1),
linspace(0, pi, ntheta+1),
ex(Tg, Pg))
plt.colorbar(r, ax=axb)
axb.scatter(p, t, c='k', s=1)
axb.set_title('Reproduced Function')

err = ex(Tg, Pg) - f_z(Tg, Pg)
r = axc.pcolormesh(Pg, Tg, err, vmin=-abs(err).max(), vmax=abs(err).max(), cmap='RdBu')
r = axc.pcolormesh(
linspace(0, 2*pi, nphi+1),
linspace(0, pi, ntheta+1),
err,
vmin=-abs(err).max(),
vmax=abs(err).max(),
cmap='RdBu')
plt.colorbar(r, ax=axc)
axc.scatter(p, t, c='k', s=1)
axc.set_title('Error with Dense Sampling')
Expand Down
18 changes: 15 additions & 3 deletions tests/spherical-harmonic/test_noise.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,13 @@ def test_noise(N, ntheta=201, nphi=401, cmap='RdBu'):
ex = noise(N, p)
Y = M.dot(ex.a).reshape(ntheta, nphi)
v = abs(Y).max()
r = axa.pcolormesh(P - pi, -T + pi/2, Y, cmap='RdBu', vmin=-v, vmax=v)
r = axa.pcolormesh(#P - pi, -T + pi/2, Y,
linspace(-pi, pi, nphi+1),
linspace(-pi/2, pi/2, ntheta+1),
Y,
cmap='RdBu',
vmin=-v,
vmax=v)
axa.grid(False)
axa.set_xticks([])
axa.set_yticks([])
Expand All @@ -57,7 +63,7 @@ def test_noise(N, ntheta=201, nphi=401, cmap='RdBu'):
axb.legend()
figb.tight_layout()

def plot_noise(N, p='red', ntheta=151, nphi=301, ni=7, nj=8, cmap='RdBu'):
def plot_noise(N, p='red', ntheta=151, nphi=301, ni=4, nj=5, cmap='RdBu'):
"""Generate noise, from red to blue, and plot the spectra along with
the expansions
args:
Expand All @@ -82,7 +88,13 @@ def plot_noise(N, p='red', ntheta=151, nphi=301, ni=7, nj=8, cmap='RdBu'):
ex = noise(N, p)
Y = M.dot(ex.a).reshape(ntheta, nphi)
v = abs(Y).max()
ax.pcolormesh(P - pi, T - pi/2, Y, cmap=cmap, vmin=-v, vmax=v)
ax.pcolormesh(#P - pi, T - pi/2, Y,
linspace(-pi, pi, nphi+1),
linspace(-pi/2, pi/2, ntheta+1),
Y,
cmap=cmap,
vmin=-v,
vmax=v)
ax.set_xticks([])
ax.set_yticks([])

Expand Down
6 changes: 3 additions & 3 deletions tests/spherical-harmonic/test_spherical_harmonic.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,9 +8,9 @@
sys.path.insert(0, join('..', '..'))
from orthopoly.spherical_harmonic import *

#plt.style.use('dark_background')
plt.rc('text', usetex=True)
plt.rc('font', family='serif')
plt.style.use('dark_background')
#plt.rc('text', usetex=True)
#plt.rc('font', family='serif')

#-------------------------------------------------------------------------------
#FUNCTIONS
Expand Down

0 comments on commit e26c458

Please sign in to comment.