球面調和関数の回転 - 驚異のアニヲタ社会復帰の予備

MikuHatsune2016-12-07

球面調和関数をRでやったりPythonでやったりしていたのだが、球面調和関数展開して係数をもっているときに、物体が回転してしまうとルジャンドル関数が回転してしまうので、係数も回転に応じて適当に変換しないといけない。
これはたいてい、回転後のパラメータがカンマつきで表すと
$f(\theta^{'},\phi{'})=\sum_{l=0}^\infty \sum_{m{'}=-l}^l c_l^{m^{'}} Y_l^{{'}}(\theta{'},\phi{'})$
という、回転後の係数 $c_l^{m^{'}}$ で考えないといけない。ちなみに回転は、オイラー角 $(\alpha, \beta, \gamma)$ で、zyz 軸がどの教科書や論文をみてもそうかいてある。
回転自体はこれがわかりやすい。
　
回転しても、もとの係数 $c_l^m$ から計算できるように偉い人が考えていて、
$c_l^{m^{'}}=\sum_{m=-l}^l D_{l}^{m^{'},m} c_l^m$
$D_l^{m^{'},m}=e^{im\gamma}d_l^{m^{'},m}(\beta)e^{im\alpha}$
$d_l^{m^{'},m}(\beta)=\sum_{t=max(0,m-m^{'})}^{min(l+m,l-m^{'})} (-1)^t \frac{\sqrt{(l+m)!(l-m)!(l+m^{'})!(l-m^{'})!}}{(l+m-t)!(l-m^{'}-t)!(t+m^{'}-m)!t!}(\cos\frac{\beta}{2})^{2l+m-m^{'}-2t}(\sin\frac{\beta}{2})^{2t+m^{'}-m}$
とすごいことになっている。
R でがんばってみたが、計算は遅いし結果もおかしいしで、Python に頼る。pyshtools はもともとfortran でshtools というものがあるっぽいのだが、それのPython 用モジュールである。例を参考にやってみる。

import pyshtools as p
import numpy as np
import math
import itertools

pi = math.pi

# 球面調和関数係数をおさめる箱を適当につくる
x = p.SHCoeffs.from_random([10, 1, 0])
x.get_coeffs()

lmax = 2
x = p.SHCoeffs.from_zeros(lmax)
ls = [0, 1, 1, 1, 2, 2, 2, 2, 2]     # l の番地
ms = [0, -1, 0, 1, -2, -1, 0, 1, 2]  # m の番地
value = [10, 1, 2, 0, 1, 5, 3, 4, 0] # (l, m) の値

# sphericalharmonics coefficients
for (l, m, v) in zip(ls, ms, value):
  x.set_coeffs(v, l, m)

x.get_coeffs() でnumpy array float 型のリストを取得できる。 $l=2$ のデータを作ってみると、長さ2、中身はl+1 * l+1 次元の行列になっている。

array([[[ 10.,   0.,   0.],
        [  2.,   0.,   0.],
        [  3.,   4.,   0.]],

       [[  0.,   0.,   0.],
        [  0.,   1.,   0.],
        [  0.,   5.,   1.]]])

リストの最初は $m\geq 0$ の項で、2つめは $m<0$ の項である。ひとつめのリストには、 $l=0,1,2$ の要素がはいっていて、各 $l$ のリストには $m=0,1,2$ の要素が入っている。
$m<0$ のほうは少し変わっていて、各 $l$ のなかに $m=0,-1,-2$ の順番で要素がはいっている。要素をいじるにはset_coeffs() で値、 $l$ , $m$ の番地を指定すれば任意にいじることができる。
係数を取り出したかったら

def shc_extract(coef, lmax):
  c = [coef.get_coeffs()[0][0][0]]
  for i in range(1, lmax+1):
    for k in [1, 0]:
      if k == 0:
        for j in range(k, i+1):
          c += [coef.get_coeffs()[k][i][j]]
      else:
        for j in reversed(range(k, i+1)):
          c += [coef.get_coeffs()[k][i][j]]
  return c

shc_extract(x_rot, lmax)

　
球面調和関数は球面上のポテンシャルの(数学的には厳密な)級数展開なので、球面上のなにかを比較したかったら、係数を比較してもいいことになっている。もとの形を球面調和関数で記述してあるとして、適当にぐるぐる回転してみて、もっとも近いもしくは遠い回転を探してみる。
たぶんこれも、何かいい方法があるのだろうが、愚直に回転させまくってみる。10度ずつの回転をzyz に対して総当りする。総当りにはitertools が使い勝手がよく、めっちゃ速い。ちなみに上でzip を使って一括iterate しているが、これもパクっている。
　

div = 10                  # 何度ずつ回転するか
degs = range(1, 360, div) # 回転角度
rot_res = []
for i, j, k in itertools.product(degs, degs, degs):
  irad, jrad, krad = i*pi/180, j*pi/180, k*pi/180 # 回転角度をradian で記録しておく
  x_rot = x.rotate(irad, jrad, krad, degrees=False)
  rot_res += [[irad, jrad, krad, shc_extract(x_rot, lmax)]]


# R じゃないといろいろできない情弱なのでファイルに書き出す
w0 = open("rotate_shc.txt", "w")
for z in rot_res:
  txt = " ".join(map(str, z[:3])) + " " + " ".join(map(str, z[3])) + "\n"
  w0.write(txt)

w0.close()

　
実際は単位球面上に模様のように色が付いているが、その値に応じて半径を決めると、球面調和関数で画像をググったときにでてくるようなチョウチョみたいなアレができる。
差が最小の時の係数の回転である。ほとんど重なっているように見えるが、これでも(321, 1, 41)度回転しているらしい。

　
差が最大の時の係数の回転である。このときは(51, 161, 131)度回転している。

　
回転自体は、どれくらい細かくグリッドサーチするかと、 $d_l^{m,m^{'}}$ が $O(l^3)$ のオーダーで計算量がかかるらしいので、 $d_l^{m,m^{'}}$ を再帰的に求める方法もあるらしい。
Journal of Computational Physics 228(16):5621-5627 September 2009
　
$d_l^{m,m^{'}}$ を計算したかったらdjpi2 という関数が使える。これを使うと、最終的な回転はrotate 関数でできるが、

alpha = pi/8
beta = pi/9
gamma = pi/10
ang = [pi/1.2, pi/1.8, pi/3]
dj = p.shtools.djpi2(lmax)

x_rot = p.shtools.SHRotateRealCoef(x.get_coeffs(), np.array(ang), dj)
# 同じ
x_rot = x.rotate(alpha, beta, gamma, degrees=False)

とやっても同様な結果が得られる。