光学設計とその周辺、そしてたまに全く関係ないやつ

学んだことを書き留めていきたいと思いますが、ありふれたことを書いても人類の進歩に貢献しないので、専門的な事柄をメインにしたいと思います。なお私の専門とは光学設計とか画像処理とかです。

光学素子製造における補正(ターゲット)寸法の例

材料を加工をする際に補正(ターゲット)寸法という考え方があります. 厚み5±0.1mmの金属板を製造したい場合, 5.0mmジャストを狙った加工するのではなく, あえて狙いとしてやや厚め, 例えば5.02mmをターゲットして加工するようなことがあります. これはもし何かあったときに厚めなら追加工することで5.00mmに修正することが出来る一方, その逆は物理的に難しいからです.

このように指定された寸法に対して, 製造上の理由であえて別の寸法を目標にすることをChatGPT曰くどうやら補正寸法とかターゲット寸法と呼ぶらしいです. このように製造上の理由を考慮して加工することは光学素子の製造でもよくありますので, 今回簡単な例を紹介します.

レンズ, フィルター厚み

これは先ほどの金属板の例と同じです. あえて実際の厚みはやや厚めにすること(されていること)はあります.

曲率

曲率についてもあえて, 指定された設計値ではなく若干値の違うターゲット値を設定することがあります. これは製造したレンズを原器に押し付けて曲率誤差を見るときに, 完全に同じ曲率だと面が接し, 面に傷などがつく可能性があるからです. 以下の図ではこれを防ぐために, あえて本来の設計値とは別の曲率で球面を作り, 光束の通らないレンズの端のみが接するようにさせています. 製造レンズが凸面ならあえて原器より緩くなりますし, 凹面なら曲率が強くなります. この状態で検査をするとき, ニュートンリングが0ではなく, ある本数出ている状態が狙い通りの面になっていると考えます.

実際は素子そのものにターゲット寸法を設定するよりは, むしろ原器と検査用の原器を分けて, その検査用原器の曲率を指定値からずらし, レンズ自体は設計値どおりになるようにすることもあります.

熱膨張を考慮した金型設計

特にプラスチック素子のように射出成型で製造する場合, 高温の膨張した状態で金型に押し込まれるため, 冷めて常温に戻ると金型自体の寸法からずれます.
先ほどの曲率の件と同じく, 熱膨張が起きることを見越して, 金型の寸法自体を最終的に合うように設計することが多いはずですので, この手のことは加工屋と確認したほうが良いでしょう.

波動光学的にPSF, MTFをシミュレーションしてみる2(2/4) -光学系のモデル化とPSF, MTF-

前回の記事では自由空間の回折による光波の伝搬を扱いました. 2回目の今回は実際にレンズがあった場合にどのように像面の強度パターンが回折の影響を受けるか扱います.
波動光学的にPSF, MTFをシミュレーションしてみる1(1/4) -自由空間の回折伝搬- - 光学設計とその周辺、そしてたまに全く関係ないやつ

光学系のモデル化

カメラレンズのように何枚も光学素子があるような光学系の中の回折伝搬を考える際, もちろん一つ一つの素子で個別に伝搬を計算するのが発想としてはシンプルです. 一方で, モデル化をして, 光学系自体を単純化することができれば手間という意味ではシンプルで, 計算量も少なく済むかもしれません.
そのモデルとして用いられるのが, 以下のような入射・射出瞳です. そもそも瞳とは何かは他文献を見てもらうとして, 射出瞳を用いることで, 複雑な光学系も単一のレンズ系として扱うことができ, 最初に射出瞳位置やそのサイズを光線追跡で求める作業は必要ですが, 全ての回折伝搬計算は射出瞳位置での1回だけで行えばよくなります. また理論的にも見通しが良くなり, いわゆる瞳の相関がOTFだとかそういう重要な考察が得られるわけです. この後者の手法をCodeVの開発元は射出瞳回折モデルと呼んでいたような. 本記事では射出瞳モデルと呼びます.

さて2種類のアプローチを紹介しましたが, 世にある光学設計ソフトではどちらの機能も実装されていることが多いです. ただ, カメラレンズとか顕微鏡のような一般的なアプリケーションでPSFやMTFを計算する目的なら, 射出瞳モデルを使うことが多いはず. 設計ソフトにあるMTF解析機能も基本的にこれに沿っています. 逆に, 特殊な状況*1では各面で個別に回折伝搬を計算することもあり, 例えばCodeVのビーム伝搬解析とかZemax OSの物理光学伝搬などの機能がこちらに相当します. 本記事の目的はMTFを得ることなので, 射出瞳モデルで考えていきます.

射出瞳モデルを使った光の伝搬とレンズのフーリエ変換作用

射出瞳モデルを使うことで, 物体からの光が単一のレンズで集光されて像面に光が伝搬されますが, このフォーカシングによる光場への影響を見てみます. 改めて上の図を見ると, 物体からの発散した球面波が瞳に入射し, レンズが像点に新たな収束する球面波を作ります. 以降レンズのある座標を (\xi,\eta)で表します.
この物体-瞳の球面波は物体から瞳の距離をaとすると瞳直前での波面は以下のように表されます.


\displaystyle U_1(\xi,\eta) = A \exp{ \left[ i \frac{\pi}{\lambda a} (\xi^2+\eta^2) \right] } \tag{1-1} \label{1-1}

また瞳-像を伝搬する光はこの距離をbとすると瞳直後の波面は以下のように同じ形の式となります*2.

\displaystyle U_2(\xi,\eta)= A' \exp{ \left[ -i \frac{\pi}{\lambda a} (\xi^2+\eta^2) \right] } \tag{1-2} \label{1-2}

A->A'と振幅が変化したのはレンズによる透過率などを考慮したことによります.
この時, レンズによる波面の作用を t(x,y)とすると, U_2(\xi,\eta)=t(\xi,\eta) U_1(\xi,\eta)となるため, t(x,y)は係数比を T_0とまとめると

\displaystyle t(\xi,\eta)= \frac{U_2(\xi,\eta)}{U_1(\xi,\eta)}=T_0 \exp{\left[ -i \frac{\pi}{\lambda f} (\xi^2+\eta^2) \right] } \tag{1-3} \label{1-3}

ここで結像の式 1/a+1/b=1/fを使いました. ここからわかる通り, レンズによる位相への影響は物体像配置にかかわらず, そのレンズの固有パラメーターである焦点距離のみで表すことが出来ます.

レンズの有限の開口を表す関数を瞳関数と呼びますがこれを p(\xi,\eta)とすると, この関数も式(1-3)にそのまま加えることが出来ます. また収差の影響もこの t(\xi,\eta)に含めることができ, 以下の図のように収差というのは射出瞳位置における理想球面の波面と実際の収差がある波面の光路差となるため, これを W(\xi,\eta)とすると, 先の瞳関数も合わせて, 以下式(1-4)のようにまとめられます*3.


\displaystyle t(\xi,\eta)= \frac{U_2(\xi,\eta)}{U_1(\xi,\eta)}=T_0 p(\xi,\eta) \exp{ \left[ -i \frac{\pi}{\lambda f} (\xi^2+\eta^2) \right] } \exp{ \left[ -i W(\xi,\eta)\right] }
\displaystyle = T_0 P(\xi,\eta) \exp{ \left[ -i \frac{\pi}{\lambda f} (\xi^2+\eta^2) \right] } \tag{1-4} \label{1-4}

となります.

点像分布関数 PSF(Point Spread Function)と周波数応答関数

以上の内容をベースに点像分布関数PSF, つまり点光源を結像させた際の像面での広がり具合, を見ていきます.
まず結像にはコヒーレント結像とインコヒーレント結像の2種類があります. 通常PSFとかMTFとかいうときはインコヒーレント結像を想定することが多いのですが, 進行の都合上コヒーレント結像をまず考えます. コヒーレント結像とは物体からの光波の振幅と位相がそのまま光学系を通って像面で干渉・回折し, 振幅の重ね合わせで像を形成する状況です. つまり像面のある点における振幅, 強度は物体面全体からの影響を受けるということです.
以下の図のような状況を考えていきますが, 物体面(座標 xy_1)の位置( x_o,y_o)にある点 P_1から距離 z_1の位置に焦点距離fの単レンズがあり, この点 P_1から伝搬する光波の振幅と位相がそのまま光学系を通ってレンズから距離 z_2を通って像面(座標 xy_2)に干渉・回折し位置( x_i,y_i)にある点 P_2に像が形成されています.

重ね合わせ積分を使うと, 像面における振幅 U_2(x_i,y_i)は, 物体面の振幅 U_1(x_o,y_o)を使って


\displaystyle U_2(x_i,y_i)= \iint_{\infty} U_1(x_o,y_o) h(x_i,y_i;x_o,y_o) dx_o dy_o \tag{2-1} \label{2-1}

となります. ここでhはインパルス応答であり, 点像つまりPSFそのものです. ここからの目的はこのPSFであるhをこれまでの回折計算結果を使って理論的に求めていきます.

物体からの1点が作る振幅分布は U_1=\delta(x_o, y_o) とデルタ関数を使うと, 像面での振幅分布 U_2自体が点像分布関数PSFとなります.

前回の記事含めこれまでの内容を振り返りP1からP2の伝搬をステップ毎にまとめると,
1)物体面のデルタ関数が距離z1の分フレネル伝搬し, レンズまで到達(U1->Ui)

  • >デルタ関数と前回記事の式(2-3)のインパルス応答を使い式(1-3)の畳み込み積分で計算

2)レンズ位置で式(1-4)で表される位相関数 t(\xi,\eta)が付加(Ui->Ui')

  • >1)の振幅と式(1-4)のレンズ作用tの掛け算を計算

3)その光が再度距離z2の分像面まで到達(Ui'->U2)

  • >2)の振幅を入力として再度1)の計算

をすることになります.

まず1)の計算をすると, 物体面上の点P1が座標 (x_o,y_o)にあるとすると, 到達面のレンズ座標 (\xi,\eta)上の振幅位相 U_i(\xi,\eta)は, 物体強度分布を表す入射光波はデルタ関数を使って U_1(x_1,y_1)=A \delta(x_1-x_o, y_1-y_o)となるため,

\displaystyle U_i(\xi,\eta)= A \delta(x_1-x_o, y1-y_o) \ast \frac{\exp{(ikz_1)}}{i \lambda z_1} \exp{ \left[ i \frac{\pi}{\lambda z_1} (\xi^2+\eta^2) \right] }
\displaystyle = \frac{A \exp{(ikz_1)}}{i\lambda z_1} \iint_{\Sigma} \delta(x_1-x_o, y_1-y_o) \exp{ \left[ i \frac{\pi}{\lambda z_1} (\xi-x_1)^2+(\eta-y_1)^2 \right] } d\xi d\eta
\displaystyle = \frac{A \exp{(ikz_1)}}{i\lambda z_1} \exp{ \left[ i \frac{\pi}{\lambda z_1} (\xi-x_o)^2+(\eta-y_o)^2 \right] } \tag{2-2} \label{2-2}

と計算できます.

次に2)の計算ですが, これは単に式(1-4)を掛け算をするだけです. ということでレンズ直後の振幅位相分布 U_i'(\xi,\eta)は


\displaystyle U_i'(\xi,\eta)=t(\xi,\eta) U_i(\xi,\eta)
\displaystyle =T_0 \frac{A \exp{(ikz_1)}}{i\lambda z_1} \exp{ \left[ i \frac{\pi}{\lambda z_1} (\xi-x_o)^2+(\eta-y_o)^2 \right]} P(\xi,\eta)\exp{ \left[ -i \frac{\pi}{\lambda f} (\xi^2+\eta^2) \right] } \tag{2-3} \label{2-3}

です.

最後に3)の U_2(x_i,y_i)計算ですが, 1)の計算と同様に,

\displaystyle U_2(x_i,y_i)=U_i'(\xi,\eta) \ast \frac{\exp{(ikz_2)}}{i \lambda z_2} \exp{ \left[ i \frac{\pi}{\lambda z_2} (x_i^2+y_i^2) \right] }
\displaystyle = T_0 \frac{A \exp{ \{ik (z_1+z_2) \} } }{(i\lambda)^2 z_1 z_2} \iint_{\infty} P(\xi,\eta) \exp{ \left[ -i \frac{\pi}{\lambda f} (\xi^2+\eta^2) \right] }
\displaystyle \exp{ \left[ i \frac{\pi}{\lambda z_1} (\xi-x_o)^2+(\eta-y_o)^2 \right] } \exp{ \left[ i \frac{\pi}{\lambda z_2} (x_i-\xi)^2+(y_i-\eta)^2 \right] } d\xi d\eta
\displaystyle = T_0 \frac{A \exp{ \{ik (z_1+z_2) \} } }{(i\lambda)^2 z_1 z_2} \exp{ \left[ i \frac{\pi}{\lambda z_1} (x_o^2+y_o^2) \right] } \exp{ \left[ i \frac{\pi}{\lambda z_2} (x_i^2+y_i^2) \right] }
\displaystyle \iint_{\Sigma} P(\xi,\eta) \exp{ \left[ i \frac{\pi}{\lambda} \left( \frac{1}{z_1}+\frac{1}{z_2}-\frac{1}{f} \right) (\xi^2+\eta^2) \right] } \exp{ \left[-i \frac{2\pi}{\lambda} \left( \left( \frac{x_o}{z_1}+\frac{x_i}{z_2} \right) \xi+\left( \frac{y_o}{z_1}+\frac{y_i}{z_2} \right) \eta \right) \right] } d\xi d\eta \tag{2-4} \label{2-4}

ここで結像の式 1/z_1+1/z_2=1/fを使えば, 瞳座標の2乗の項はキャンセルされるため, 式(2-4)はフーリエ変換の形になる. これがレンズのフーリエ変換と呼ばれる所以です.
また通常最終的には振幅ではなくその絶対値2乗である強度が重要なため, 積分の前にある \exp{ \{ik (z_1+z_2) \} }や \exp{ \left[ (x_i^2+y_i^2) \right] } は無視できます. 注意が必要なのは物体座標に関係する項 \exp{ \left[ (x_o^2+y_o^2) \right] } です. というのは式(2-1)の畳み込み積分の積分変数と同じ変数となっているため, 単純に無視することはできません. 良く説明されるのがMTFのよい光学系なら物体点と像点がほぼ1対1となるため, 事実上無視できるということです. そんなに単純な話では実際はないですが, インパルス応答の絶対値2乗が伝達関数となるインコヒーレント結像ではどうせキャンセルされるため, これ以上は深堀せず以降は無視して考えていきます*4.
以上の議論をまとめ, この光学系の倍率mが m=-z_2/z_1となることも利用すると式(2-4)を記述しなおすと,

\displaystyle U_2(x_i,y_i)=h(x_i,y_i;x_o,y_o)
\displaystyle = \frac{A'}{\lambda^2 z_1 z_2} \iint_{\infty} P(\xi,\eta) \exp{ \left[-i \frac{2\pi}{\lambda z_2} \left[ (x_i-mx_o) \xi+ (y_i-my_o) \eta \right] \right] } d\xi d\eta \tag{2-5} \label{2-5}

つまり, インパルス応答は像座標 (mx_o, my_o)を中心とした開口関数pのフラウンホーファー回折像となります.

PSFをもとめるだけならこの式(2-5)で終わりなのですが, OTFを求めるには式(2-1)が畳み込みの形で表す必要がありますので, もうちょい式(2-5)を変えていきましょう.
まず式(2-5)を理想的な結像状態で幾何光学的に表すとどうなるかを見てみます. \xi' = \xi/\lambda z_2, \eta'=\eta/\lambda z_2とし式を変形すると

\displaystyle h(x_i,y_i;x_o,y_o)= A' m \iint_{\infty} P(\lambda z_2 \xi', \lambda z_2 \eta') \exp{ \left[-i 2\pi \left[ (x_i-mx_o) \xi'+ (y_i-my_o) \eta' \right] \right] } d\xi' d\eta' \tag{2-6} \label{2-6}

幾何光学かつ理想的な結像となると瞳関数は瞳座標すべてで1, \lambda \to 0となるため, この時式(2-6)はデルタ関数の定義そのものとなるため,

\displaystyle h(x_i,y_i;x_o,y_o) \to A' m \delta( x_i-mx_o, y_i-my_o) \tag{2-7} \label{2-7}

となります.
これを式(2-1)に代入すると,

\displaystyle U_{2g}(x_i,y_i) = \frac{1}{m}U_1(x_i/m,y_i/m) \tag{2-8} \label{2-8}

です. 回折も考慮した振幅U_2とこの幾何光学的な振幅を区別するためこの式(2-8)はU_2gとしています.
mx_o=x_o', my_o=y_o', h'=h/mとし式を変形すると式(2-8)も使うと,

\displaystyle U_2(x_i,y_i)= A' \iint_{\infty} h'(x_i-x_o',y_i-y_o') \frac{1}{m} U_1\left( \frac{x_o'}{m},\frac{y_o'}{m} \right) dx_o' dy_o' = h'(x_o',y_o') \ast U_{2g}(x_o',y_o') \tag{2-9} \label{2-9}

強度の場合は

\displaystyle I_2(x_i,y_i) = | U_2(x_i,y_i)|^2= |h'(x_o',y_o') \ast U_{2g}(x_o',y_o') |^2 \tag{2-10} \label{2-10}

です.
また改めてh'は

\displaystyle h'(x_i-x_o',y_i-y_o')= A' \iint_{\infty} P(\lambda z_2 \xi', \lambda z_2 \eta') \exp{ \left[-i 2\pi \left[ (x_i-x_o') \xi'+ (y_i-y_o') \eta' \right] \right] } d\xi' d\eta' \tag{2-11} \label{2-11}

です.

以上の内容をまとめると, 式(2-9)のように回折限界の光学系がつくる理想像は幾何光学に理想な像とインパルス応答の畳み込みとなる, そしてそのインパルス応答は式(2-6)のように回折の効果は瞳関数のフラウンホーファー回折像となる.

ここまではコヒーレント結像の場合でした. ではインコヒーレント結像ではどうなるでしょうか. インコヒーレント結像では物体各点からの振幅はインコヒーレントゆえに像面で干渉しません. その代わり, 各点からの強度が像面でそのまま足し算として最終的な強度を作ります.
つまり非常にざっくりと言えば式(2-9)の振幅の関係がそのまま強度の関係となります. 式(2-6)のインパルス応答自体は振幅に対する作用なので, 強度への作用はは式(2-6)を二乗だとすると,

\displaystyle I_2(x_i,y_i)= | A'|^2 \iint_{\infty} |h'(x_i-x_o',y_i-y_o')|^2 I_{2g}(x_o',y_o' dx_o' dy_o' = |h'(x_o',y_o') |^2 \ast | I_{2g}(x_o',y_o')|^2 \tag{2-12} \label{2-12}

となります.

光学伝達関数OTFとMTF

改めて前節の内容をまとめると, コヒーレント結像の場合は式(2-9)のように振幅の入出力が畳み込み積分で表される, つまり線形時不変システムだということ. インコヒーレントの場合は式(2-11)のように強度となります.
それぞれの式をフーリエ変換すると

\displaystyle G_{U2}(f_X,f_Y)=H(f_X,f_Y)G_{Ug}(f_X,f_Y) \tag{3-1} \label{3-1}
\displaystyle G_{I2}(f_X,f_Y)=H'(f_X,f_Y)G_{Ig}(f_X,f_Y) \tag{3-2} \label{3-2}

となります.

ここでHとH'はそれぞれの伝達関数となり, 光学ではそのまま光学伝達関数(OTF, Optical transfer function)と呼びます. コヒーレント結像の方のHを特に振幅伝達関数(ATF)と呼ぶこともあります.
インコヒーレント結像の場合さらにOTFの絶対値を取ったものを変調伝達関数(MTF, Modulation transfer function)と呼びます.
ここでさらにインコヒーレント結像の伝達関数は正規化することが普通なので, 式(3-2)のH'は

\displaystyle H'(f_X,f_Y)=\frac{\mathcal{F} [ |h|^2 ] }{\iint_{\infty} |h'(x_i,y_I)|^2 dx_i dy_i} \tag{3-3} \label{3-3}

となります. OTFを絶対値と位相部分で分解すると, 以下の式(3-4)が得られます. PTFが位相部分となりあまり光学で議論されることはないですが位相伝達関数と呼びます(Phase Transfer Function).

\displaystyle H'(f_X,f_Y)=OTF(f_X,f_Y)=MTF(f_X,f_Y) \exp{ \left[ i PTF(f_X,f_Y) \right] } \tag{3-4} \label{3-4}

数値計算例

せっかくなので数値計算例を示したいと思います. いったん波面収差はないとして(W=0), 焦点距離100mm, 瞳径5mm, F値20の軸上のPSFとMTFを計算するPythonコードが以下です.

import numpy as np
from matplotlib import pyplot as plt

def circ(t, width):
    #矩形開口を作成.
    return np.where(np.abs(t) <= width / 2, 1.0, 0.0)

M=1024
L=1e-3
dx=L/M

x=np.arange(-L/2,L/2-dx,dx)
y=x

wl=0.55e-6
k=2*np.pi/wl
Dxp=5e-3
wxp=Dxp/2
zxp=100e-3
lz=wl*zxp

fx=np.arange(-1/(2*dx),1/(2*dx),1/L)
Fx,Fy=np.meshgrid(fx,fx)

H=circ(np.sqrt(Fx*Fx+Fy*Fy),2*wxp/lz)

h2=np.abs(np.fft.ifftshift(np.fft.ifft2(np.fft.fftshift(H))))**2

OTF=np.fft.fft2(np.fft.fftshift(h2))
MTF=np.abs(OTF)
MTF=np.fft.ifftshift(MTF)

PSF_cross_section= h2[int(M/2), :]
MTF_cross_section= MTF[int(M/2), :]

mask = fx >= 0
fx_filtered = fx[mask]
MTF_Normalized_temp = MTF_cross_section[mask]
MTF_Normalized = MTF_Normalized_temp/MTF_Normalized_temp[0]

# 2D PSF plot
plt.imshow(h2,extent=[x.min(), x.max(), y.min(), y.max()])
plt.colorbar(label="Value")
plt.title("2D PSF")
plt.show()

# PSF(cross section) plot
plt.scatter(x, PSF_cross_section)
plt.title("PSF")
plt.xlabel("cyc/m")
plt.ylabel("Intensity")
plt.show()

# MTF plot
plt.scatter(fx_filtered, MTF_Normalized)
plt.title("MTF")
plt.xlabel("cyc/m")
plt.ylabel("MTF")
plt.show()

まずPSFは以下の結果. 無収差なのでこのスケールだとグラフにする意味もあまりないですけど.

以下はMTFです. 瞳径を変えれば結果も変わりますので試してみください.

今回はここまで. 次回は収差があるときを詳しく扱う予定です.

*1:例えば絞りが特殊な形状しているとか中間像があるときなど.

*2:この時プラスマイナスが違うことに注意.

*3:式(1-4)の内,瞳関数と波面収差の項だけを一般化瞳関数と呼びます.

*4:逆にいえばコヒーレント結像では要注意. ここらへんはGoodmanの本に説明があります.

波動光学的にPSF, MTFをシミュレーションしてみる1(1/4) -自由空間の回折伝搬-

光学設計でとっても重要なPSF(Point Spread Function)やMTF(Modulation Transfer Function)を扱います. そもそもPSFとかMTFって何というところは本ブログ以上にうまく説明されている情報源はたくさんあるため割愛します. 理論的なところも詳しい式の導出は本格的な教科書に譲るとして, 一方設計ソフトに完全にお任せの実際のコードとかを示しながら説明していければなと思っています. 最終的なゴールは実際の収差も考慮してPSFとかMTFの数値計算を自分でもできるようになることです. 本記事では自由空間を光が伝搬するときの回折の影響を扱います. 数値計算も交えながら, 回折シミュレーションの雰囲気を感じ, 次の記事ではレンズなど光学素子があるときの特性を扱います.

重要なポイント

今回扱うのはスカラー近似に基づいた波動光学の方です. 波動光学ということは回折が考慮されるということですが, その回折の扱いも近似のレベルから何通りかありますが, それぞれ個別に扱います. 一通り説明し終えた後, 基本すぎて逆に取り上げられることの少ない幾何光学MTFについても考えたいと思います.

回折と各種伝搬の概要

最初に回折による光波の伝搬とその代表的モデルであるフレネル伝搬, フラウンホーファー伝搬をまとめます.
回折という現象を図にすると以下の通りになります. ホイヘンスの2次波の原理に基づき, 各点で波面が2次波を作り, そこが観測点で重なり合うことにより新たな光場が得られます.

この現象を数式で表してみましょう. 座標は以下の図のような取り方, つまり左側から平行光が \Sigmaの開口をもつ遮蔽物に入射したときに回折により観測面O'-xyの点P(x,y)でどう光場を作るかを表しています.

これを式にすると観測点P'の振幅 U_{2}(x,y)は,


\displaystyle U_{2}(x,y) = \frac{z}{i λ} \iint_{\Sigma} U_{1}(\eta,\xi) \frac{ \exp(ikr_{12})}{r_{12}^2} d\xi d\eta \tag{1-1} \label{1-1}
\displaystyle r_{12}^2=\sqrt{z^2+(x-\xi)^2+(y-\eta)^2} \tag{1-2} \label{1-2}

となります. \sigma面で球面波である2次波が最終的にP'に寄与すると考えると式の内容がわかりやすいかなと思います.

ここで線形時不変システムだと考えると


\displaystyle U_{2}(x,y) = \iint_{\Sigma} U_{1}(\eta,\xi) h(x-\xi,y-\eta) d\xi d\eta \tag{1-3} \label{1-3}
\displaystyle h(x,y) = \frac{z}{i λ} \frac{\exp(ikr)}{r^2} \tag{1-4} \label{1-4}
\displaystyle r^2=\sqrt{z^2+x^2+y^2} \tag{1-5} \label{1-5}

とインパルス応答hを使って表すことができます.
よって最終的にU2を計算するときも, フーリエ変換の関係を利用して

\displaystyle U_{2}(x,y)=\mathcal{F}^{-1}[\mathcal{F}\{U_{1}(x,y)\} \mathcal{F}\{h(x,y)\}] =\mathcal{F}^{-1}[\mathcal{F}\{U_{1}(x,y)\} H(f_{X},f_{y})] \tag{1-6} \label{1-6}

と考えても良いです. ここでHはインパルス応答hをフーリエ変換した伝達関数で,

\displaystyle H(f_{X},f_{y})= \mathcal{F}\{h(x,y) \}= \exp{ \left( jkz \sqrt{1-(\lambda f_X)^2-(\lambda f_Y)^2 } \right) } \tag{1-7} \label{1-7}

となります( f_X=x/(\lambda z)). ここまでの厳密な解析積分を使った展開をRayleigh-Sommerfeld解とよばれ, この定式化に基づく光の伝搬の扱いをホイヘンス伝搬とか呼ばれることがあります*1.

次に近似を使った別の伝搬の定式化を考えていきます. 観測面は十分遠くにある場合, 式(1-2)内のzは x-\xiや y-\etaより大きくなり, 式(1-2)は


\displaystyle r_{12}^2=\sqrt{z^2+(x-\xi)^2+(y-\eta)^2} \approx z \left[ 1+\frac{1}{2z^2} (x-\xi)^2+\frac{1}{2z^2} (y-\eta)^2 \right] \tag{2-1} \label{2-1}

と近似できます.

これを式(1-1)に代入すると, 指数の部分についてはこの近似をそのまま使い, 分母の r_{12}^2はそのままzと置き換えても比較的影響は小さいため*2, 式は


\displaystyle U_{2}(x,y) = \frac{ \exp(ikz)}{iλz} \iint_{\Sigma} U_{1}(\eta,\xi) \exp{ \left[ i\frac{k}{2z}[ (x-\xi)^2+(y-\eta)^2 ] \right] } d\xi d\eta \tag{2-2} \label{2-2}


となります.
ここでインパルス応答と伝達関数を以下のように定義すると, 式(1-3)と同じ畳み込みで表されます.

\displaystyle h(x,y) = \frac{ \exp{(ikz)}}{iλz} \exp{ \left[ \frac{ik}{2z} (x^2+y^2 ) \right] } \tag{2-3} \label{2-3}
\displaystyle H(f_X,f_Y) = \exp{(jkz)} \exp{ \left[ -i\pi \lambda z (f_X^2+f_Y^2) \right] } \tag{2-4} \label{2-4}

ここで式(2-1)をフレネル近似, そしてこの近似に基づく伝搬式(2-2)をフレネル伝搬と呼びます.

最後にもう1段階近似した伝搬式を紹介します.
式(2-2)の2次式を展開すると以下の式(2-5)となります.


\displaystyle U_{2}(x,y) = \frac{ \exp(ikz)}{iλz} \exp{ \left[ i \frac{k}{2z}(x^2+y^2) \right] } \iint_{\Sigma} U_{1}(\eta,\xi) \exp{ \left[ i \frac{k}{2z}(\xi^2+\eta^2) \right] } \exp{ \left[ - i\frac{2 \pi}{\lambda z}[ x\xi+y\eta ] \right] } d\xi d\eta \tag{2-5} \label{2-5}

ここで以下の式(3-1)の近似をすると, これはフレネル近似よりさらに

\displaystyle z>>k(\xi^2+\eta^2)/2 \tag{3-1} \label{3-1}

式(2-4)は以下の式(3-2)となります.

\displaystyle U_{2}(x,y) = \frac{ \exp(ikz)}{iλz} \exp{ \left[ i \frac{k}{2z}(x^2+y^2) \right]} \iint_{\Sigma} U_{1}(\eta,\xi) \exp{ \left[ - i\frac{2 \pi}{\lambda z}[ x\xi+y\eta ] \right] } d\xi d\eta \tag{3-2} \label{3-2}

この近似(3-1)と基づく光の伝搬をフラウンホーファー近似, フラウンホーファー伝搬と呼ばれます.

式(3-2)は入射場のフーリエ変換そのものです. ホイヘンス伝搬やフレネル伝搬違い, フラウンホーファー伝搬の場合は伝達関数は存在しません. 式(3-2)のフーリエ変換演算そのものは線形時不変ではないので.

理論は以上として, 数値計算をしてみます. いろいろ紹介しましたが, ここでは
A) 式(1-3)のRS積分を直接計算したホイヘンス伝搬
B) 式(1-6)のフーリエ変換の関係を使ったホイヘンス伝搬
C) 同じくフーリエ変換を使ったフレネル伝搬計算
D) そして式(3-2)のフラウンホーファー伝搬の数値計算
の計4つを扱います.
以下のPythonコードが最初の3つを実装した例です. RSPropagationDirectがホイヘンス伝搬の直接積分, RSPropagationがホイヘンス伝搬のフーリエ変換を使った計算, FresnelPropagationが同じくフーリエ変換を使ったフレネル伝搬で, 行76で矩形である開口面上の入射場を定義しており, 80-82行で回折伝搬を計算しています.

import numpy as np
from matplotlib import pyplot as plt

def rectpuls(t, width):
    #矩形開口を作成
    return np.where(np.abs(t) <= width / 2, 1.0, 0.0)

def RSPropagationDirect(u1, L, wl, z):
    #直接Rayleigh-Sommefeldの畳み込みを計算
    M, N = u1.shape
    dx = L / M
    k = 2 * np.pi / wl

    x=-L/2+dx*np.arange(M)
    y=-L/2+dx*np.arange(M)
    X1, Y1 = np.meshgrid(x, y)

    u2 = np.zeros_like(u1, dtype=complex)

    for i in range(M):
        for j in range(N):
            x2 = x[i]
            y2 = y[j]

            R = np.sqrt((X1 - x2)**2 + (Y1 - y2)**2 + z**2)
            h = (1 / (1j * wl)) * (np.exp(1j * k * R) / R) * (z / R)
            u2[j, i] = np.sum(u1 * h) * dx * dx

    return u2

def RSPropagation(u1,L,wl,z):
    #フーリエ変換を使ってRayleigh-Sommefeldの畳み込みを計算
    M,N=u1.shape
    dx=L/M
    k=2*np.pi/wl
    
    fx=-1/(2*dx)+(1/L)*np.arange(M)
    FX,FY=np.meshgrid(fx,fx)

    H=np.exp(1j*k*z*np.sqrt(1-wl*wl*FX*FX-wl*wl*FY*FY))
    H=np.fft.fftshift(H)
    U1=np.fft.fft2(np.fft.fftshift(u1))
    U2=H*U1
    u2=np.fft.ifftshift(np.fft.ifft2(U2))
    return u2

def FresnelPropagation(u1,L,wl,z):
    #フーリエ変換を使ってフレネル伝搬を計算
    M,N=u1.shape
    dx=L/M
    k=2*np.pi/wl
    
    #fx=array_with_interval(-1/(2*dx),1/L,1/(2*dx)-2)
    fx=-1/(2*dx)+(1/L)*np.arange(M)
    FX,FY=np.meshgrid(fx,fx)

    H=np.exp(-1j*np.pi*wl*z*(FX*FX+FY*FY))
    H=np.fft.fftshift(H)
    U1=np.fft.fft2(np.fft.fftshift(u1))
    U2=H*U1
    u2=np.fft.ifftshift(np.fft.ifft2(U2))
    return u2

L1=0.5
M=250
dx1=L1/M
x1=-L1/2+dx1*np.arange(M)
y1=-L1/2+dx1*np.arange(M)
wl=0.5*pow(10,-6)
k=2*np.pi/wl

w=0.051
z=2,000

X1,Y1=np.meshgrid(x1,y1)
u1=rectpuls(X1,2*w)*rectpuls(Y1,2*w)  #開口面の振幅
I1=np.abs(u1)**2  #開口面の強度

# 伝搬計算
#u2=RSPropagationDirect(u1,L1,wl,z) #RS伝搬の直接計算
#u2=RSPropagation(u1,L1,wl,z) #RS伝搬のFFTを使った計算
u2=FresnelPropagation(u1,L1,wl,z) #フレネル伝搬のFFTを使った計算

x2=-L1/2+dx1*np.arange(M)
y2=-L1/2+dx1*np.arange(M)
I2=np.abs(u2)**2 

Cross_section= I2[int(M/2), :] #断面取得

fig, axs = plt.subplots(1, 3, figsize=(15, 6))  

# 入力面の強度
im1 = axs[0].imshow(I1, extent=[x1.min(), x1.max(), y1.min(), y1.max()],
origin='lower', aspect='equal', cmap='gray')
axs[0].set_title('Input Intensity')
axs[0].set_xlabel('x (m)')
axs[0].set_ylabel('y (m)')
plt.colorbar(im1, ax=axs[0], label='Intensity')

# スクリーン面の強度
im2 = axs[1].imshow(I2, extent=[x1.min(), x1.max(), y1.min(), y1.max()],
origin='lower', aspect='equal', cmap='gray')
axs[1].set_title('Output Intensity (after propagation)')
axs[1].set_xlabel('x (m)')
axs[1].set_ylabel('y (m)')
plt.colorbar(im2, ax=axs[1], label='Intensity')

# 強度断面図
axs[2].plot(x1, Cross_section, color='black')
axs[2].set_title('1D Profile at y = 0 (center slice of Output)')
axs[2].set_xlabel('x (m)')
axs[2].set_ylabel('Intensity')
axs[2].grid(True)

plt.tight_layout()
plt.show()

まず開口面の強度は以下の通り(I2).

そしてz=2,000の時の計算結果は以下の通り(上からホイヘンス直接, ホイヘンスFFT, フレネルFFT). 3つの手法で特に違いはありませんが, 回折の影響で開口面の矩形分布が干渉縞のような変化をしています.

次にz=20,000の時は以下の通り. 先ほどと違い, ホイヘンス直接とそれ以外は明らかに違いがありますが, FFTを使った2つの方法の結果は不自然です.

今回の一連の趣旨と外れるためあまり深くは扱いませんが, これはサンプリングの条件がサンプリング定理を満たしていないからです.

いろいろめんどくさいことを考えたくない場合はホイヘンス伝搬の直接計算を使えばよいですが, こちらは圧倒的に遅いです. そこで対策としてFFTのほうはいろんな手法が研究されています. そちらは光学設計とは別の話となるためこれ以上は扱いませんし, 私は個人的には本当に正しい設定をしているか自信がない場合は畳み込み計算を使うのもよいと思います.

なおFFTを使った場合, ホイヘンス伝搬とフレネル伝搬はほとんど差はありません. フレネル近似の制約がある分, ホイヘンス伝搬のほうが正確かと思いきや, 意外とその差はなく, 次回の記事で扱う光学系のフーリエ変換作用を考えやすいため, 基本的にはフレネル伝搬で考えます.

さて最後にフラウンホーファー伝搬を扱います.
この場合の計算は

def FraunhoferPropagation(u1,L,wl,z):
    #フラウンホーファー伝搬を計算
    M,N=u1.shape
    dx=L/M
    k=2*np.pi/wl
    
    #x2=array_with_interval(-L/2,dx,L/2-dx)
    x2=-L/2+dx*np.arange(M)
    X2,Y2=np.meshgrid(x2,x2)

    c=1/(1j*wl*z)*np.exp(1j*k/(2*z)*(X2**2+Y2**2))

    u2=c*np.fft.ifftshift(np.fft.fft2(np.fft.fftshift(u1)))*(dx**2) 

    return u2

結果を見やすくするためにw=0.011, z=20,000としたとき, 以下の通りu2を計算します.
u2=FraunhoferPropagation(u1,L1,wl,z) #フラウンホーファー伝搬の計算

結果は次の通り. 強度分布のコントラストでわかりづらいですが, sinc関数の2乗のようになっていることが確認できると思います.

フラウンホーファー伝搬は伝搬距離の制約が非常に強いですが, 実は光学系のフーリエ変換作用により, PSFやMTFを計算するときはフラウンホーファー伝搬を使うことが可能であり, 先述したややこしいサンプリングの問題はあまり気にしなくて良いです. そういえばいわゆる角スペクトル法というのもありますが, これは本質的にはホイヘンス伝搬と同じことをしています.

内容は以上. 今回は自由空間上での回折伝搬でしたが, 次回は光学系がある場合はどう伝搬するかを扱います.

参考文献
Computational Fourier Optics: A MATLAB Tutorial

*1:後述するフレネル伝搬とかフラウンホーファー伝搬ほど一般的ではない.

*2:指数の部分については少しの違いでも指数の値は大きく変わる.