Mathematicaの回折計算用のちょい便利な関数

MathematicaでFraunhoffer回折を計算するのは簡単で
1.瞳関数となる2次元matrixを作る
2.Fourier[]で瞳関数をFourier変換する
3.例えば強度分布（絶対値の2乗）を表示する
とすればよい。
FFTを行うライブラリを呼ぶプログラムを書くことを考えればずいぶん楽である。

Mathematicaでは数値的に2次元Fourier変換した後は[[1,1]]の位置が原点（DC成分）となり、これを表示すると回折像が領域の四隅に切り取られてしまう。これを簡単に表示する関数を作った。

matrixTakeFourier2D[m_?MatrixQ, halfSpan_] :=
        Take[
            RotateRight[
                Map[Take[RotateRight[#, halfSpan], 2halfSpan + 1] &, m],
            halfSpan],
        2halfSpan + 1]

表示したいmatrixと、原点の周りの何ポイント分を表示するかを指定する。
例えば

In[40]:= pupil = With[{r = 32},
    Table[If[x^2 + y^2 <= r^2, 1., 0.], {y, -4r, 4r}, {x, -4r, 4r}]];
In[41]:= ld[Abs[pupil]];

で無収差瞳関数。これをfourier変換すると

In[42]:= spot = Fourier[pupil];
In[43]:= ld[Abs[spot]^2];

となって、スポットが四隅に泣き別れ。

In[44]:= ld[matrixTakeFourier2D[Abs[spot]^2, 10]];

となって見やすくなる、というもの。
上のld[]は単なるDensityPlotのオプションを決め打ちして見やすくする関数

ld[m_,opt___]:=ListDensityPlot[m,Mesh->False,AspectRatio->Automatic,opt]

ついでにvectorやmatrixを補間する関数

vectorReduce[v_, nReduct_:2] :=
    Map[Apply[Plus, #]/nReduct &, Partition[v, nReduct]]

matrixReduce[m_?MatrixQ, nReduct_:2] :=
    vectorReduce[Map[vectorReduce[#, nReduct] &, m], nReduct]

一部を取り出す関数

matrixTake[m_?MatrixQ, {{xs_, ys_}, {xe_, ye_}}] :=
    Module[{xsi, ysi, xei, yei},
        {{xsi, ysi}, {xei, yei}} = Ceiling[{{xs, ys}, {xe, ye}}];
        Take[Map[Take[#, {xsi, xei}] &, m], {ysi, yei}]
    ]

など。これはAddOnにあるか。

コメント 2

こんにちは、
pigpioで検索して来ました。Mathematicaによる偏光にも興味があり、ファイルをダウンロードさせて頂きたいのですが、404 not foundになってしまいます。確認頂けないでしょうか？
by ぶなのもり (2018-08-16 11:39) 

コメントありがとうござます。

実は、このブログにあるリンクのほどんとは切れてしまっています。すみません。このファイルは改めて復活させますので少しお待ちください。

by decafish (2018-08-16 19:32)