Googleが出しているC++向けのSpherical Harmonics(SPH)ライブラリの簡単な使い方です．(Spherical HarmonicsをSPHと訳すの，粒子法のSPHと混ざってややこしいですね…)
球面調和関数の各項を計算するだけでなく，3次元極座標関数の球面調和変換(球面調和展開, Spherical Harmonics Transform=SHT)も行うことができます．
ダウンロード
以下のページからダウンロードできます．zipでDLするもよし，git cloneするもよしです．
導入
dllやlibに依存していない.cpp(.cc)と.hのみのライブラリなので手軽に導入できます．
自分のプロジェクトのinclude用ディレクトリにshディレクトリを移動させるだけです．楽！
ただしこのときEigenも必要になるのでこちらもDLして導入しておきます．
(Eigen自体もDLしてきてEigenディレクトリをinclude用ディレクトリに移動させるだけ)
使い方
自分のコードに
#include "sh/spherical_harmonics.h"
を記入して開始．
いちいち名前空間を書きたくなければ
using namespace sh;
と書いておきましょう．(以下はusingした状態でのコードとなります．)
球面調和関数(Y_l^m(Φ,θ))を求める
超簡単です。以下の関数にl, m, phi, thetaを与えて実行するだけです．
double EvalSH(int l, int m, double phi, double theta);
球面調和変換
テイラー展開の球面調和関数版です．
まずは変換(展開)したい関数をつくって関数ポインタの形で渡せるようにします．
そのあとProjectFunction関数に渡すことによって，展開した際の係数のベクトルを得ることができます．
// 楕円体出力用の関数
double a = 1.0, b = 1.0, c = 0.5; // 楕円体のx, y, z 方向の長さの半分
SphericalFunction ellipsoidFn = [a, b, c](double phi, double theta) {
return ellipsoid(theta, phi, a, b, c);
};
// 球面座標変換(SPT)
auto result = ProjectFunction(8, ellipsoidFn, kDefaultSampleCount);
for (int i = 0; i < result->size(); i++)
{
std::cout << "C(" << i << "): " << result->at(i) << std::endl;
}
ここではあらかじめ楕円体関数(theta, phi, a, b, cを与えるとrが帰ってくるような関数)を他に作成しておいて，ラムダ式の中にそれを突っ込んでいます．
ProjectFunction(int order, const SphericalFunction &func, int sample_count);
関数の引数についてですが，orderは何次元まで計算させるかの数字を入れます．理論的にはここに∞を入れれば元の関数と一致します．よって大きくすればするほど近似度合いは高くなりますが，計算量も増えていきます．なお出てくる係数の数は当然order^2個となります．funcは計算させたい極座標関数の関数ポインタを与えます．sample_countは計算で使用するモンテカルロ積分のサンプル数です．kDefaultSampleCountにデフォルト値が設定されていますが，使用用途によってはこれが大きすぎることもあるようです．
