球面上のポテンシャルなど極座標表示できるものを関数と表すとする。離散データは容量がかさばるのと特徴をつかみにくいという問題点があるので何らかの方法で連続化して表せると良い。そこで用いられる方法の一つが球面調和関数による展開である(他には球面ウェーブレット変換などもあるらしい)。
以下は球面調和関数の概要とその可視化のメモ。
1. 球面調和関数とは
関数をという級数を用いて表示することを考える。が一次独立なら基底として扱えて、さらに内積$<F_{i},F_{j}> = \int F_{i}F_{j} dx = \delta_{ij}$なら正規直交基底である。フーリエ展開や球面調和関数展開は例の一つ。は$<f, F_{i}>$として求められる。
ルジャンドル陪関数
ルジャンドル陪関数は $$ \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} $$
球面調和関数
球面調和関数は、ルジャンドル関数の入力の代わりにを使ったものを用いて表される。はz軸からの角度、は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} $$ 可視化すると
(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)
を可視化した。左上はそのまま、右上は正規化したもの、下は半径としてそれぞれ可視化した。
- 実数版
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
左上が、右上が正規化したもの、左下が符号付き半径、右下が符号なし半径をそれぞれ可視化した。
- 適当に係数を与えて足しあわせてみると何らかの形を得ることができる。
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)
発展的な内容
球面調和関数は立体に効率よく影を落としたりするのに用いられているらしい。
- https://basesandframes.wordpress.com/2016/05/11/spherical-harmonic-lighting-the-gritty-details/: 球面調和関数の全てが分かる上に実装の参考になる
- OpenGLで球面調和関数を描画したきゅーぶろぐ(仮) | きゅーぶろぐ(仮): 動くデモがある
- ☆PROJECT ASURA☆ [XNA GS] 『球面調和関数』
- ☆PROJECT ASURA☆ [OpenGL] 『球面調和関数(1)』
- ☆PROJECT ASURA☆ [OpenGL] 『球面調和関数(2)』
球面調和関数で与えられた磁気ポテンシャルのデータから極の位置を計算できる