照度分布計算その14 - デバグの続き

回折による照度分布の計算の続き。とりあえずコンパイルは通って一見動いているようになった。ちゃんとしておかなければいけない仕上げのデバグ。
とりあえず、結果が解ってる条件で計算する。まずピンホールの手前までの回折計算。

円形の瞳内で強度分布が均一の場合、その焦点位置での振幅分布A(r)は、当然回転対称になって解析的に計算できて

となる。J₁(x)は第1種のベッセル関数の1次のやつ。これをプロットすると

この形を「ベッセル・シンク（Bessel-sinc）」と言う人もいる。

rは極座標の径方向の変数でもちろん

rの単位はちょっと気をつける必要があるけど回折計算ではよく出てくる

を単位にはかる。

普通実際に観測できるのは振幅ではなくて強度なのでざっくりA(r)²になって、これがどんな形をしているかと言うと

みたいなの。この同心円状の強度を図中の矢印の位置で一旦強度が0になる（J₁(x)の最初の零点）。この内側には全強度の83.8%が集中するのでこれを「エアリーのディスク」と呼ぶ人がいる。矢印の位置を"えいや"で「ディスク半径」、というよりは「スポット半径」と言う。これはだいたい0.61λ/N.A.になる。

さて、こないだ書いた数値計算のピンホールの手前でこれにあっているかどうかを見てみる。瞳内の一様な振幅分布のためには専用の振幅関数を作った方がいいけど、とりあえずFWHMを大きな値にしてしまえば計算上は問題ない。波長を1μmとして、絞り込みレンズのN.A.が0.122になるように例えば焦点距離10mmのアパチャ径2.44mmにすればスポット半径は5μm、直径で10μmになるはず。これを試してみよう。

Mathematicaの上から

In[86]:= params = {fftSize -> 512, fwhmh -> 1000, fwhmv -> 1000, lambda -> 1, 
   pupilAreaSize -> 40, focalLengthOfConversion -> 10, 
   apertureDiameter -> 1.22*2, pinholeDiameter -> 10, 
   outputLaserPlane -> 0, outputAperturePlane -> 1, 
   outputPrePinholePlane -> 1, outputPinholePlane -> 1, 
   outputPostPinholePlane -> 0, outputIntensityPlane -> 0, 
   outputDirectIllumination -> 0};
In[87]:= t = RunThrough["./a.out", params];

としてみる。帰ってくるのが遅い。どうも、計算そのものよりも出力が大きすぎてMathematicaに送っているpipeの方が律速しているみたい。転送サイズを設定するパラメータを作った方がいいな。

まず、アパチャでの振幅分布を表示すると

In[88]:= With[{r = pupilAreaSize*dsize/fftSize /. params}, 
 ListPlot3D[
  Chop[matrixTakeFourier2D[outputAperturePlane /. t, dsize]], 
  PlotRange -> All, DataRange -> {{-r, r}, {-r, r}}]]

半径1.22mm。よくわからんけどだいたいあってそうか。

ピンホール前の強度分布を表示させる。これが式-43のA(r)²で、エアリーの半径が5μmになっていると正しい。

With[{r = focalLengthOfConversion/pupilAreaSize*dsize /. params}, 
 ListDensityPlot[
  Chop[matrixTakeFourier2D[outputPrePinholePlane /. t, dsize]]^2, 
  PlotRange -> Automatic, DataRange -> {{-r, r}, {-r, r}}]]

まあ、こんなもんか。

わかりにくいので、プロファイルの表示にしてみる。

In[91]:= With[{r = focalLengthOfConversion/pupilAreaSize*dsize /. params}, 
 ListPlot[Chop[
     matrixTakeFourier2D[outputPrePinholePlane /. t, 
      dsize]][[dsize]]^2, DataRange -> {-r, r}, Joined -> True, 
  PlotStyle -> Thick]]

とすると

最初の零点は5のところになる。OKじゃないでしょうか。

座標は、Mathematicaで与えているので同時に間違っていたら意味が無い。径を10μmに設定してあったピンホールの振幅分布を表示させてみると

In[92]:= With[{r = focalLengthOfConversion/pupilAreaSize*dsize /. params}, 
 ListDensityPlot[
  Chop[matrixTakeFourier2D[outputPinholePlane /. t, dsize]], 
  PlotRange -> All, DataRange -> {{-r, r}, {-r, r}}]]

あってるね。たぶん、だいたい。

ここまでであっていれば、次のピンホールの透過率分布を合成したり、その後の伝播も同じクラスを使っているので問題ないはず。残りもさくっとやってみよう。