Maximaで球面調和関数展開 - 次世代3Dコンテントパイプライン開発室

Maximaという数式処理ソフトを知っていますか？これ、フリーながらかなりの高機能で、個人的にかなり重宝しています。そんなMaximaの威力を、CG屋さんに馴染み深い例で紹介しましょう。SH(Spherical Harmonics)の関数展開です。

それではまずMaximaを立ち上げて、ためしにY(l=3,m=3)を求めてみましょう。
SHの場合、プリセットで定義されているので、

INPUT:spherical_harmonic(3,3,theta,phi);

と入力するだけで

$\frac{5\,\sqrt{7}\,{e}^{3\,i\,\phi}\,{sin\left( \theta) }^{3}}{8\,\sqrt{5}\,\sqrt{\pi }}$

と出てきます。複素数がそのまま出てきてしまっているので、ここでは実数（m<0なら虚数）を抽出しておきましょう。ついでにsqrt(2)もかけておきます。

INPUT:realpart(%) * sqrt(2);
$\frac{5\,\sqrt{2}\,\sqrt{7}\,cos\left( 3\,\phi) \,{sin\left( \theta) }^{3}}{8\,\sqrt{5}\,\sqrt{\pi }}$

でました。これだけでも十分に実用的な計算式なのですが、極座標系から直交座標系への変換にも挑戦してみましょう。Maximaではまず以下のおまじないが必要です。

INPUT:%,trigexpand=true;
$\frac{5\,\sqrt{2}\,\sqrt{7}\,\left( {cos\left( \phi) }^{3}-3\,cos\left( \phi) \,{sin\left( \phi) }^{2}) \,{sin\left( \theta) }^{3}}{8\,\sqrt{5}\,\sqrt{\pi }}$

これは三角関数の展開を意味しているのですが、Maximaではこの手続きを踏まないと以下の変換が何故かうまくいきませんでした。唯一のハマりポイントと言えるでしょう。それではtheta,piの順に関数合成します。

INPUT:ratsubst(acos(z),theta,%);
$-\frac{\left( 15\,\sqrt{2}\,\sqrt{7}\,cos\left( \phi) \,{sin\left( \phi) }^{2}-5\,\sqrt{2}\,\sqrt{7}\,{cos\left( \phi) }^{3}) \,{\left( 1-{z}^{2}) }^{\frac{3}{2}}}{8\,\sqrt{5}\,\sqrt{\pi }}$
INPUT:ratsubst(atan2(y,x),phi,%);
$-\frac{\left( 15\,\sqrt{2}\,\sqrt{7}\,x\,{y}^{2}-5\,\sqrt{2}\,\sqrt{7}\,{x}^{3}) \,{\left( 1-{z}^{2}) }^{\frac{3}{2}}}{8\,\sqrt{5}\,\sqrt{\pi }\,{\left( {y}^{2}+{x}^{2}) }^{\frac{3}{2}}}$