初めに
初めまして．今年の4月まで波動光学を専攻していた新社会人です．
コロナ嫌ですね．マスクが熱いです．熱いといえば最近はリアルタイムレイトレーシングが熱いですね．「写実的な絵をリアルタイムで！」．なんとも良い時代に生まれました．
しかし，幾何光学で実現できる物理表現にも限界があります．例えば回折ですね．遮蔽物で遮られた光がその遮蔽物の裏側に回り込むことです．このような現象には数値的な計算手法が確立されており，DOE(回折光学素子)の設計時のシミュレーションや計算機合成ホログラム(CGH)の計算に使用されています．
というわけで，この記事では波動光学の伝搬計算をしてみたいという方向けに，有名な計算手法の説明を記そうと思います．
本記事での演算結果は
CPUであれば
こちら
GPUであれば
こちら
にあるプログラムにて確認できます．
初投稿どうぞお手柔らかに．
波動光学とは
「光は波である」という視点に立って光の理論を突き詰め，応用する学問です．反射・回折・波長分散・屈折など様々な現象を扱います．
特に回折は計算や理論が少し難しい部分がありますが，使いこなせば強力な武器にもなります．
早速，伝搬計算手法について説明していきます．
伝搬計算手法
最初にこの記事を読んで計算できるであろう平行平面間の伝搬計算の結果をgif画像で示しておきます．
伝搬計算手法には主に以下の3通りが挙げられます．
- フレネル回折
- フラウンホーファー回折
- 角スペクトル法
最初の二つは近似を使用して導き出されたものであり，伝搬距離によって重大な計算誤差を生じます．ありがたいことに
こちら
で近似について記されているので，ぜひ確認してみてください．
角スペクトル法の概要
例えば，フレネル回折の式は，開口面(遮蔽物の面)から離れた位置で観測される回折波の複素振幅分布を与えます．離れた位置での計算はできますが，開口に近接する領域での結果はレイリーゾンマーフェルトの回折積分式のようなごちゃごちゃした積分を直接法で計算する必要があります．
直接法は時間がかかるので，近傍領域でも高速に計算できる手法はないのでしょうか?
それが角スペクトル法です．
この手法はヘルムホルツ方程式に基づいた回折伝搬計算手法です．というかレイリーゾンマーフェルトの回折積分式による厳密解を与えます．従って近似などを使用せず，大変精度の良い結果が得られる手法です．しかも高速!なんという!
早速この手法に迫っていきましょう!!!!!
角スペクトル法の説明
この部分も先ほど紹介した記事に載っているのですが，折角なので，角スペクトル法が実はそのまま用いていては完ぺきな手法とはならない場合があることも踏まえて，最初から説明していきます．以下にまず，使用する座標系を示します．
位置$z = 0$にある平面$(x,y;z = 0)$に光波分布$g(x,y;z = 0)$を考え，そこから距離$z$離れた位置の平面$(x,y;z = z)$における回折波の光波分布$g(x,y;z = z)$を計算することを考えます．いきなりですがここで，$g(x,y;z = 0)$のフーリエ変換$G(u,v;z=0)$を示します．
G(u,v;z=0) = \iint g(x,y;z = 0)\exp[-i2\pi(ux+vy)]dxdy
更に，逆フーリエ変換により次が成り立ちます．
g(x,y;z = 0) = \iint G(u,v;z=0)\exp[i2\pi(ux+vy)]dudv
ここで，$u$と$v$は空間周波数と呼ばれるもので，それぞれ$x$軸$y$軸に沿った平面波の波動ベクトルの各成分です．平面波は波動ベクトルの方向に進行する波であり，
A(x,y)\exp[i2\pi(ux+vy)]
の形式で表される波です．平面波を知ったところで，逆フーリエ変換の式を見てください．平面波が被積分関数になっていますね．
つまりどういうことかといいますと，元の光波はそのスペクトルを振幅として，様々な方向に回折する平面波$G(u,v;z=0)\exp[i2\pi(ux+vy)]$を無限に足し合わせたものであると解釈できるのです．図に謎の黄色い何かが記されていましたが，実はいろいろな方向に進む平面波を示していました!!
ちなみに，$\mathbf{k}=(k_x,k_y,k_z)$は波動ベクトルを示しており，$(k_x,k_y,k_z)=2\pi(u,v,w)$と$|\mathbf{k}|=2\pi/\lambda$を満足します．
逆に，フーリエ変換の式は，光波が平面波の群に展開されていると解釈することもできます．
となれば，伝搬計算というのは，この展開されて得た平面波群の各平面波を観測面までそれぞれ伝搬して，再度足し合わせればよいということになります．これが角スペクトル法の原理です．
平面波$G(u,v;z=0)\exp[i2\pi(ux+vy)]$を距離$z$進めるには$\exp[i2\pi\sqrt{(1/\lambda)^2 - u^2 -v^2} \cdot z]$をかければよいです．これは，$g(x,y;z = z)=\iint G(u,v;z)\exp[i2\pi(ux+vy)]dudv$がヘルムホルツ方程式を満足するという条件から導出されます．
以上より，観測面$(x,y;z=z)$における光波のフーリエ変換結果は
G(u,v;z=z) = G(u,v;z=0)\exp[i2\pi\sqrt{(1/\lambda)^2 - u^2 -v^2} \cdot z]
となりますから，この逆フーリエ変換
\iint G(u,v;z=0)\exp[i2\pi\sqrt{(1/\lambda)^2 - u^2 -v^2} \cdot z]\exp[i2\pi(ux+vy)]dudv
が，求めたい光波分布$g(x,y;z = z)$になります．
ここで，$\exp[i2\pi\sqrt{(1/\lambda)^2 - u^2 -v^2} \cdot z]$は周波数応答関数と呼ばれ，周波数空間で光波を輸送する働きをします．ここでは
H(u,v;z=z) = \exp[i2\pi w(u,v)z]\\
w(u,v) = \left\{
\begin{array}{ll}
\sqrt{(1/\lambda)^2 - u^2 -v^2} & (u^2 +v^2 \leq (1/\lambda)^2) \\
0 & (\text{otherwize})
\end{array}
\right.
です．以上をまとめると，フーリエ変換の記号を$\mathcal{F}[\cdot]$として，
g(x,y;z = z) = \mathcal{F}^{-1}[\mathcal{F}[g(x,y;z = 0)]H(u,v;z=z)]
で計算が実行されます．フーリエ変換(離散フーリエ変換)はFFTにより高速に処理されますから，相当の速度で伝搬計算が実行できます．うれしいですね!!!
角スペクトル法の計算結果
ここで，いったん数値計算をしてみます．私のプログラムを使用します．
こちら
こんな開口を用意します．縦横512ピクセル，標本間隔1[μm]，波長は630[nm]で設定します．
1[mm]伝搬してみると
となりました．もう少し長距離いっちゃいましょうか．10[mm]
おや?かなり汚いですね．なぜでしょうか．何が光波に起きているのでしょうか．プログラムのミス??はたまたバグ??
周波数応答関数とエイリアシング
先に原因を述べると，角スペクトル法による長距離伝搬では，周波数応答関数(先ほどの$H(u,v;z=z) = \exp[i2\pi w(u,v)z]$)が，エイリアシングしてしまうからです．
この問題を解決した手法として，Band-Limited Angular Spectrum Method(原論文)があります．
ここでは論文の内容を踏まえて簡単に角スペクトル法のエイリアシングとBand-Limited Angular Spectrum Methodについて，説明していきます．
エイリアシングはなぜ起きる?
角スペクトル法は，原理的には完全な伝搬計算手法です．しかし，数値シミュレーションを行う以上は，我々は離散化された信号にFFTを用いて，畳み込みを実行することになります．ここで，周波数応答関数の式をじっくり眺めていただくと，伝搬距離に対して周波数が増加する関数であることがわかります(これをアップチャープな関数といいます)．
このことから，長距離伝搬を数値計算する場合，周波数応答関数を適当な周波数で帯域制限しなければなりません．周波数が上がりすぎると，標本間隔を変えなかった場合，ナイキスト定理を満たさなくなってしまうからです．
Band-Limited Angular Spectrum Method
では，いったいどのような周波数帯域で制限すればよいのでしょうか?
$H(u,v;z=z)$の局所的な信号の周波数はそれぞれ$u,v$方向に対して$ \left|\frac{\partial (w(u,0)z)}{\partial u}\right|$と$ \left|\frac{\partial (w(0,v)z)}{\partial v}\right|$となります．計算すると，$ \left|\frac{uz}{w(u,0)}\right|$と$ \left|\frac{vz}{w(0,v)}\right|$となりますね．
この周波数がナイキスト定理を満たすとき，周波数$u,v$の標本間隔をそれぞれ$\Delta u, \Delta v$とすれば，標本周波数は${\Delta u}^{-1}, {\Delta v}^{-1}$ですから(今は周波数空間を基準にしているのでこちらが「周波数」扱いです．)，各軸方向に対して
2\left|\frac{uz}{w(u,0)}\right|\leq {\Delta u}^{-1}\; \text {and}\; 2\left|\frac{vz}{w(0,v)}\right|\leq {\Delta v}^{-1}
を満足する必要があります．変形すると次のようになります．
|u|\leq \frac{1}{[(2\Delta u z)^2+1]^{1/2}\lambda}\; \text {and}\; |v|\leq \frac{1}{[(2\Delta v z)^2+1]^{1/2}\lambda}
従って，これが求める制限帯域となります．これを導き出し，実践したものが先ほどの論文であり，この手法はBand-Limited Angular Spectrum Methodと呼ばれます．
早速，導入してみましょう．主要な部分を抜粋しました．書き方などは目をつぶってください(震え声)．
void WaveFront::generateFRF(double distance)//周波数応答関数の生成
{
double k = 1 / w_lambda / w_lambda;
int i, j;
for (i = 0; i < w_nx; i++)
{
for (j = 0; j < w_ny; j++)
{
double kx = itox(i);
double ky = jtoy(j);
double w = k - kx * kx - ky * ky;
if (w > 0)
{
double phase = 2 * PI * sqrt(w) * distance;
SetPixel(i, j, complex<double>(cos(phase), sin(phase)));
}
else//evanescent element
{
SetPixel(i, j, complex<double>(1.0, 0.0));
}
}
}
}
void WaveFront::bandlimit(double uband, double vband)//指定された周波数帯域で制限
{
int i, j;
for (j = 0; j < w_ny; ++j)
{
for (i = 0; i < w_nx; ++i)
{
if (abs(itox(i)) <= uband / 2 && abs(jtoy(j)) <= vband / 2)
{
continue;
}
else
{
SetPixel(i, j, complex<double>(0.0, 0.0));
}
}
}
}
void WaveFront::AsmProp(const double R)//角スペクトル法で指定した距離分伝搬計算
{
fft2D(-1);
AsmPropInFourierSpace(R);
fft2D(1);
*this /= GetN();
}
void WaveFront::AsmPropInFourierSpace(const double R)//周波数空間での伝搬
{
WaveFront h(w_nx, w_ny, w_px, w_py, w_lambda);//光波オブジェクトの生成
h.generateFRF(R);//周波数応答関数を設定
double uband = 2 / (sqrt((2 * h.w_px * R) * (2 * h.w_px * R) + 1) * h.w_lambda);//u方向の制限帯域の計算
double vband = 2 / (sqrt((2 * h.w_py * R) * (2 * h.w_py * R) + 1) * h.w_lambda);//v方向の制限帯域の計算
h.bandlimit(uband, vband);// 帯域制限
int x, y;
for (y = 0; y < w_ny; ++y)
{
for (x = 0; x < w_nx; ++x)
{
SetPixel(x, y, GetPixel(x, y) * h.GetPixel(x, y));
}
}
SetOrigin(GetOrigin() + GetNormal() * R);
}
いい感じです．因みに，制限帯域は少し複雑な形状をした領域になるのですが，矩形の領域で近似してもあまり影響ありません．このプログラムでは矩形領域で帯域制限しています．早速計算結果のおかしかった10[mm]の伝搬計算を再計算しましょう．
[結果]
[先ほどの結果]
ふつくしい．．．大変理想的な結果が得られました．
と書いたのはいいものの，これもより遠距離の伝搬になると結果が望まぬものになります．その場合は，伝搬元の標本窓を縦横倍ずつにゼロパディングで拡張した後に伝搬計算を行い，その後元のサイズに縮小することで理想の結果を獲得できます．是非やってみてください．私のプログラムでどんな感じなのか見るのもよいかもしれません．
最後にもう一度徐々に伝搬した場合の変化をgif画像で示しておきます．
おわりに
ここまで読んでいただいた方，稚拙な文章に長々とお付き合いいただきありがとうございました．
これから波動光学の道に進まれる方や単純に興味のある方に少しでも役立てる記事となっていることを願います．
「光あれ!!!!!!」
ではこれにて失礼します．お疲れさまでした．
次