Zernike多項式 - その7 Mathematica関数を使う

ちょっと時間があいたけど前回の続き。今日は使用例。

使用例

いつものようにZernikePolynomials.mを$PathのどこかにOpticsフォルダを作ってそこにコピーする。そのあと

Needs["Optics`ZernikePolynomials`"]; 

として読み込む。

Packageで定義されたシンボルは

で確認できる。usage messageも書いてあるのでシンボルをクリックすると簡単なヘルプが表示される。 とりあえず例えばfringe orderで最初の16項を極座標とデカルト座標で書いてみると

あってるかどうかは、まあ僕が20年近く使っているのでおっけーでしょう、たぶん。

トリビアルな関数の展開例

実際に単位円の中で定義される簡単な関数をZernike多項式で展開してみよう。 例のひとつ目は、

1. 原点を中心として半径1/2の内部は1
2. その外は0

というような円筒型の関数を考える。

これは回転対称なので、Zernike多項式のうち回転対称でないものは使えない（係数を計算しても0になるだけ）なので回転対称な式だけ抜き出す。たとえばfringe orderで最初の80のうち回転対称な式は

In[68]:= sphs=Select[Range[0, 80], sphericalQ]
Out[68]= {0, 3, 8, 15, 24, 35, 48, 63, 80}

の9つだけである。これを使って展開すると、

In[69]:= coef=Map[Integrate[
     zernikePolarPolynomial[#][r, t] r, {r, 0, 1/2}, {t, 0, 
      2 Pi}]/(Pi/rmsCoefficient[#]^2) &, sphs]
Out[69]= {1/4, -(9/16), 15/32, -(21/256), -(135/512), 627/2048, -(273/
  4096), -(13005/65536), 32181/131072}

となる。

この係数を使ってZernike多項式を足し合わせる

In[70]:= Plot3D[Evaluate[
  Map[zernikeCartesianPolynomial[#][x, y] &, sphs].coef*circle[x, y]],
    {x, -1.1, 1.1}, {y, -1.1, 1.1}, PlotRange -> All]

と、

のようになる。ちょうどFourier級数で矩形波を展開したときのGibbsの現象と同じようなことが起こるのはまあ、あたりまえ。

表示のところで使っているcircle[x,y]という関数は半径1の内部で1を返してその外では0を返すユーティリティ関数で、表示を見やすくするためだけのもの。

上の例では単位円内の積分は極座標の方が簡単で、表示するときはデカルト座標の方が簡単なので使い分けている。数値的にではなく代数式として等しいので置き換えても問題ない。こういうところはMathematicaを使う醍醐味。

つぎはもうひとつ別の例。