球面上のポテンシャルなど極座標表示できるものを関数 $f(\theta, \phi)$ と表すとする。離散データは容量がかさばるのと特徴をつかみにくいという問題点があるので何らかの方法で連続化して表せると良い。そこで用いられる方法の一つが球面調和関数による展開である(他には球面ウェーブレット変換などもあるらしい)。

以下は球面調和関数の概要とその可視化のメモ。

1. 球面調和関数とは

関数を $f(x) = \sum_{i}^{\infty} c_{i}F_{i}(x)$ という級数を用いて表示することを考える。 $F_{i}(x)$ が一次独立なら基底として扱えて、さらに内積$<F_{i},F_{j}> = \int F_{i}F_{j} dx = \delta_{ij}$なら正規直交基底である。フーリエ展開や球面調和関数展開は例の一つ。 $c_{i}$ は$<f, F_{i}>$として求められる。

ルジャンドル陪関数

ルジャンドル陪関数 $P_{l}^{m}$ は $$ \begin{align} &P_{l}^{m}(x) &= &(-1)^{l} (1 - x^{2})^{m/2} \frac{d^{m}}{dx^{m}} P_{l}(x) \\ &P_{l}(x) &= &\frac{1}{2^{n}n!} \frac{d^{l}}{dx^{l}} (x^{2}-1)^{l} \end{align} $$ と表される。

内積は $$ \begin{align} <P_{l}^{m}, P_{l'}^{m}> &= &\int_{-1}^{1} P_{l}^{m} P_{l'}^{m}dx \\ &= &\frac{2}{2n+1}\frac{l+m}{l-m}\delta_{ll'} \end{align} $$

球面調和関数

球面調和関数 $Y_{l}^{m}$ は、ルジャンドル関数の入力 $x$ の代わりに $\cos\theta$ を使ったものを用いて表される。 $\theta$ はz軸からの角度、 $\phi$ はx軸からの角度を表す。 $$ Y_{l}^{m}(\theta, \phi) = (-1)^{(m + |m|)/2} \sqrt{\frac{2l+1}{4\pi} \frac{(l-|m|)!}{(l+|m|)!}} P_{l}^{|m|} (\cos\theta) e^{im\phi} $$ 右辺のルジャンドル陪関数以前の部分は正規直交基底にするための規格化定数なのであまり気にしなくて良さそう。

内積は、 $$ \begin{align} <Y_{l}^{m}, Y_{l'}^{m'}> &= &\int_{0}^{2\pi}d\phi \int_{0}^{\pi}d\theta \sin\theta Y_{l}^{m} Y_{l'}^{m'} \\ &= &\int_{0}^{2\pi}d\phi \int_{-1}^{1}d\cos\theta Y_{l}^{m} Y_{l'}^{m'} \\ &= &\delta_{ll'}\delta_{mm'} \end{align} $$

実数版は $$ Y_{l}^{m} = \begin{cases} \sqrt{2}(-1)^{m}Re(Y_{l}^{m}) &\text{if} &m>0\\ Y_{l}^{0} &\text{if} &m=0\\ \sqrt{2}(-1)^{m}Im(Y_{l}^{-m}) &\text{if} &m<0 \end{cases} $$ 可視化すると

f:id:kuyata:20160217125723p:plain f:id:kuyata:20160217130414p:plain

(wikipediaより改変したもの)

行がl(0,1,2...)、列がm(...,-1,0,1,...)を表している。実数版か複素数版どちらを可視化したのかは不明。

参考文献
- t-pot『テイラー、フーリエ、球面調和関数』: 分かりやすい解説
- Spherical Harmonic -- from Wolfram MathWorld
- 予報ができない人の気象メモ — 球面調和関数の添字: 添字の名前について

具体例

基底となる球面調和関数をプロットする。

複素数版

以下で用いるscipyのsph_harmはwikipedia等と表記が異なり、使いにくいので、

from scipy.special import sph_harm as _sph_harm
def sph_harm(l, m, theta, phi):
    return _sph_harm(m, l, phi, theta)

としておく。

import numpy as np
from scipy.special import sph_harm as _sph_harm
def sph_harm(l, m, theta, phi):
    return _sph_harm(m, l, phi, theta)
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import matplotlib.cm as cm

# R^2(u=θ, v=Φ)
u = np.linspace(0, np.pi, 80)
v = np.linspace(0, 2 * np.pi, 80)

# R^3
x = np.outer(np.sin(u), np.cos(v))
y = np.outer(np.sin(u), np.sin(v))
z = np.outer(np.cos(u), np.ones(v.size))

fig = plt.figure()
Lmax = 3
for l in range(Lmax + 1):
    for m in range(-l, l+1):
        # 球面調和関数の値を色で表す
        sph = np.abs(sph_harm(l, m, np.array(u, ndmin=2).T, np.array(v)))
# 正規化
#        if sph.max() == sph.min():
#            sph.fill(0.5)
#        else:
#            sph = (sph - sph.min()) / (sph.max() - sph.min())
        colors = cm.jet(sph) #cm.bwr
        ax = fig.add_subplot(Lmax + 1, 2*Lmax+1, l*(2*Lmax+1)+m+Lmax+1, projection="3d")
        ax.set_axis_off()
        surf = ax.plot_surface(x, y, z, facecolors=colors,
                rstride=2, cstride=2, linewidth=10, shade=False)

plt.savefig("sh_c.png", transparent=False)

$|Y_{l}^{m}|$ を可視化した。左上はそのまま、右上は正規化したもの、下は半径としてそれぞれ可視化した。

実数版

def sph_harm(l, m, theta, phi):
    if m > 0:
        return np.sqrt(2) * (-1)**m * _sph_harm(m, l, phi, theta).real
    if m == 0:
        return _sph_harm(m, l, phi, theta).real
    else:
        return np.sqrt(2) * (-1)**m * _sph_harm(-m, l, phi, theta).imag

左上が $Y_{l}^{m}$ 、右上が正規化したもの、左下が符号付き半径、右下が符号なし半径をそれぞれ可視化した。

適当に係数を与えて足しあわせてみると何らかの形を得ることができる。

import numpy as np
from scipy.special import sph_harm as _sph_harm
def sph_harm(l, m, theta, phi):
    if m >= 0:
        return np.sqrt(2) * _sph_harm(m, l, phi, theta).real
    else:
        return np.sqrt(2) * _sph_harm(-m, l, phi, theta).imag
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import matplotlib.cm as cm

# R^2
u = np.linspace(0, np.pi, 80)
v = np.linspace(0, 2 * np.pi, 80)

# lの最大値
Lmax = 2
# rの初期値
r = np.zeros([u.size, v.size])
# 係数
c = np.concatenate([[10], [1, 2, 0], [1, 5, 3, 4, 0]])

for l in range(Lmax + 1):
    for m in range(-l, l + 1):
        # 球面調和関数を足し合わせる
        r += c[l**2 + l + m] * sph_harm(l, m, np.array(u, ndmin=2).T, v).real

# R^3
x = r * np.outer(np.sin(u), np.cos(v))
y = r * np.outer(np.sin(u), np.sin(v))
z = r * np.outer(np.cos(u), np.ones(v.size))

fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")
ax.set_axis_off()
surf = ax.plot_surface(x, y, z, color='b',
        rstride=2, cstride=2, linewidth=1, shade=False)
plt.savefig("shape.png", transparent=False)