Kumpei Shiraishi's blog

物理、プログラミング、クライミングに関する雑記

Pythonで球面調和関数を描画

目標

Pythonを使って、下図のようなものを描く。

Spherical Harmonics

結論

好きな自然数{l}{m}{m=-l,-l+1,\dots,l-1,1})を以下のコード冒頭で指定して、実行する。

Draw spherical harmonics with your preferred quant …

そして以下のような図を得る。

20170530_sph_harm_3_1

能書き

私は、この図を一枚描くのにかなりの時間を取られたので、バカな自分向けに、極めて初歩的な所から備忘録的にメモをしていく。

モジュールのimport

  • numpy:自明
  • scipy.special.sph_harm:球面調和関数。後述。
  • matplotlib.pyplot:自明
  • mpl_toolkits.mplot3d.Axes3D:3次元の図を描くのに必要のようだ。

scipysph_harmについて

公式のreferenceを見ると、引数は以下の4つを取る。これがWikipedia等で用いられているnotation(球面調和関数球面座標系と異なることが、混乱の原因その1だった。

  1. m:普通の{m}と一緒。即ち、{m=-l,\dots,l}となる整数。
  2. n:Wikipedia等では{l}とされている自然数。
  3. theta:定義域は[0,2*pi]。つまり、Wikipedia等では{\phi}とされている角度成分。
  4. phi:定義域は[0,pi]。つまり、Wikipedia等では{\theta}とされている角度成分。

Tipとしては、用意したデータセットの範囲をよく確認して、どっちの引数に食わせるか考えると良い。後に参考文献で挙げるページでは、再定義しているが、私にはこれは余計混乱するだけだった。

データセットの構成

サンプルコードでは、極座標の角度のデータ数を100とした。このことは、絵的に100個の点を用意したことにはならない。10,000個の点を用意しているのである。極座標変換の式

{\displaystyle
x = r\sin\theta\cos\phi \\
y = r\sin\theta\sin\phi \\
z = r\cos\theta
}

を見れば分かるように、{x}{y}{\theta}{\phi}の積で得られる。よって、絵的には10,000点を用意しているのである。{z}はそれに合わせて10,000点用意する。

これをコードにすると、直積を取るnp.outerを用いることになる。これにより、100行1列のベクトル{\sin\theta}と1行100列のベクトル{\cos\phi}の積である100行100列の行列が生成される({x}の場合)。

np.meshgridで鮮やかに書けるが、とりあえず1ステップずつ、自分が何をやっているのか分かる方が良いと思ってnp.outerを使った。

これに対応する形で球面調和関数も作らないといけない。まずthetaphiの直積も作り、それらをsph_harmに放り込んで、用意した。もちろん最初に指定した量子数を使った。

球面調和関数が何を表しているか

目的の「例の図」で、球面調和関数は何に対応しているのか?それは、中心からの距離である。半径rの3次元球に、生成した球面調和関数を掛けると、「例の図」の1パーツが出来上がる。

注意したいのは、一般に球面調和関数は複素数値を取るということだ。実部と虚部の物理的な意味はよく把握していないが、「例の図」を描くには実部を選べばよい。また、中心からの距離なので、実部を選んだあとで、絶対値を取る(順番を間違えないように)。

3次元描画

3次元描画の方法は、色々ある。まずfigureのインスタンスを生成、サブプロットでprojection='3d'とすれば、3次元で描ける。

サンプルで生成した図では、wireframeにした。他の事例は参考文献参照。

それを保存すれば、めでたしめでたし。

参考文献