光学設計の文脈でよく出てくる点像分布関数(PSF)と、MTF曲線の具体的な数値計算の方法について最近調べていたので、それらについてまとめてみました。私は光学設計の専門ではなく、そもそも物理学科出身でもないので、間違っているところなどあるかと思いますが、お気づきの方はご指摘いただけると幸いです。
球面波の場合、波源からの距離を$r$として $$ u(r, t) = \frac{A}{r}\exp[i(kr - \omega t + \delta)] $$ と表されます。ここで$k$は波数です。
以下の計算では$\omega t$の項と$\delta$の項は影響してこないので無視されます。
点像分布関数(PSF)について
点像分布関数(Point Spread Function - PSF)とは点光源が像面上で結像したときの、像面上における光の強度分布を表す関数です。 |
| 図1: 点像分布関数 |
図1はPSFの例です。赤色に近いほど光の強度が高いことを表します。理想的な結像では点光源は1点に集光するので、PSFはデルタ関数のような1点のみにピークを持つ関数になります。しかし、実際には収差や回折の影響で一点には集光せず、図1のようなある程度広がりを持った像が得られます。
PSFには回折の影響を考慮しない幾何光学的PSFと、回折の影響を考慮した波動光学的PSFの2種類があります。一般的にはPSFと言うと波動光学的PSFの方を指すようです。図1は回折を考慮した波動光学的PSFで、中心の周りにピークがもう一つあり、波打っているように見えます。
縦軸がコントラストを表し、横軸は像高と言って像面中心からの距離を表します。赤線・緑線は特定の空間周波数に対応し、実線・破線はサジタル(S)方向とメリジオナル(M)方向に対応します。ここで言うサジタルとメリジオナル方向の意味がよく分かっていないのですが、おそらく像面を極座標系$(r, \theta)$で表したときの動経方向$r$と方位角方向$\theta$に対応するのではないかと思われます。
要約すると、MTF曲線は動経方向と方位角方向の特定の空間周波数におけるコントラストを像高を変えながら表したものということになります。
MTFもPSF同様、光の回折の影響を考慮するかしないかで、幾何光学的MTFと波動光学的MTFがあります。
要約すると、MTF曲線は動経方向と方位角方向の特定の空間周波数におけるコントラストを像高を変えながら表したものということになります。
MTFもPSF同様、光の回折の影響を考慮するかしないかで、幾何光学的MTFと波動光学的MTFがあります。
PSFの数値計算
PSFには幾何光学的PSFと波動光学的PSFの2種類がありました。計算方法が変わってくるので分けて説明します。
幾何光学的PSFの数値計算
幾何光学的PSFは点光源が像面上でどのように結像するかを、回折の影響を考慮しないで表したものでした。これは光線追跡(レイトレーシング)を用いて比較的簡単に計算することができます。
点光源を適当な位置(一般的には無限遠?)に配置し、そこからレンズ前面に向けて多数の光線を発射します。これらの光線がレンズに当たるたびに屈折させ次の方向を計算、さらに次のレンズとの衝突位置の計算を繰り返します。最後に光線が像面と交わる位置を計算し、その位置に光が持つ放射輝度を加算します。
このようにして多数の光線に対しレイトレーシングを繰り返し、サンプル点の間は適当に補間することで、像面上における光の強度分布を得ることができます。こうして幾何光学的PSFを得ることができます。
なお、光線が像面と交わった位置を表すグラフはスポットダイアグラム(Spot Diagram)と呼ばれます。図2はその例です。
 |
| 図2: スポットダイアグラムの例 |
波動光学的PSFの数値計算
波動光学的PSFを計算するためには光を波として表現し、レンズを通過したときの回折を考慮し、像面上の波を足し合わせたものを求めることで像面上の強度を計算できます。これを愚直にやろうとするとかなり大変そうです。そこで各種の近似を用いることで計算を実行可能な形に落とし込みます。
まずは基本となる回折の計算方法について説明していきます。
光の波としての表現方法
光を波として考え、数式で表すには複素数を用いた表現が利用されます。
平面波の場合、位置$\boldsymbol{r}$, 時刻$t$における複素振幅$u(\boldsymbol{r}, t)$は $$ u(\boldsymbol{r}, t) = A\exp[i(\boldsymbol{k}\cdot\boldsymbol{r} - \omega t + \delta)] $$
と表されます。ここで$A$は振幅, $i$は虚数単位, $\boldsymbol{k}$は波数ベクトル, $\omega$は角周波数, $\delta$は初期位相です。平面波の場合、位置$\boldsymbol{r}$, 時刻$t$における複素振幅$u(\boldsymbol{r}, t)$は $$ u(\boldsymbol{r}, t) = A\exp[i(\boldsymbol{k}\cdot\boldsymbol{r} - \omega t + \delta)] $$
球面波の場合、波源からの距離を$r$として $$ u(r, t) = \frac{A}{r}\exp[i(kr - \omega t + \delta)] $$ と表されます。ここで$k$は波数です。
以下の計算では$\omega t$の項と$\delta$の項は影響してこないので無視されます。
回折積分
回折の計算を説明するにあたって、以下のような座標系で考えます。
 |
| 図3: 回折計算の座標系 |
$\Sigma$は開口部分を表し、$P(x_p, y_p, 0)$は$\Sigma$上の点です。$Q(x_q, y_q, z_q)$は像面上の点です。$\theta$は線分$PQ$と$z$軸のなす角を表します。
開口部分$\Sigma$を通過する平面波が、回折しながら点$Q$に到達する状況を考えます。フレネル・ホイヘンスの原理を用いると、点$Q$の複素振幅は$\Sigma$上で発生する二次波を重ね合わせたものとして表すことができます。
まずは、点$P$から点$Q$に到達する二次波の複素振幅を考えると
$$
\frac{1}{i\lambda}\frac{e^{ikr}}{r}U(x_p, y_p, 0)\cos{\theta}
$$
となります。ここで$\lambda$は光の波長, $k$は波数, $U(x_p, y_p, 0)$は点$P$における複素振幅です。(この部分の導出は専門書などを参照ください。ここでは結果だけ使わせてもらいました。)
点$Q$の複素振幅$U(x_q, y_q, z_q)$はこれらの二次波を$\Sigma$上全体に渡って重ね合わせたものなので $$ \begin{align} U(x_q, y_q, z_q) &= \frac{1}{i\lambda}\int\int_{\Sigma} U(x_p, y_p, 0)\frac{e^{ikr}}{r}\cos{\theta}dx_p dy_p \tag{1} \\ &= \frac{z_q}{i\lambda}\int\int_{\Sigma} U(x_p, y_p, 0)\frac{e^{ikr}}{r^2}dx_p dy_p \end{align} $$ となります。式変形では$\cos{\theta} = \frac{z_q}{r}$を用いました。
この積分はフレネル・キルヒホッフの回折積分と呼ばれます。これを数値積分すれば点$Q$における複素振幅が計算できるので、その絶対値の2乗を計算することで像面上の強度を得ることができます。しかし、この計算方法は開口面と像面の分割数を$N$とすると$O(N^4)$の計算量になるため、実用的ではありません。
点$Q$の複素振幅$U(x_q, y_q, z_q)$はこれらの二次波を$\Sigma$上全体に渡って重ね合わせたものなので $$ \begin{align} U(x_q, y_q, z_q) &= \frac{1}{i\lambda}\int\int_{\Sigma} U(x_p, y_p, 0)\frac{e^{ikr}}{r}\cos{\theta}dx_p dy_p \tag{1} \\ &= \frac{z_q}{i\lambda}\int\int_{\Sigma} U(x_p, y_p, 0)\frac{e^{ikr}}{r^2}dx_p dy_p \end{align} $$ となります。式変形では$\cos{\theta} = \frac{z_q}{r}$を用いました。
この積分はフレネル・キルヒホッフの回折積分と呼ばれます。これを数値積分すれば点$Q$における複素振幅が計算できるので、その絶対値の2乗を計算することで像面上の強度を得ることができます。しかし、この計算方法は開口面と像面の分割数を$N$とすると$O(N^4)$の計算量になるため、実用的ではありません。
フレネル近似
$\theta$が小さく、$\cos{\theta} \approx 1$となる状況(近軸近似)を考えます。
式1に注目すると、$\cos{\theta}$の部分は1になり、$\frac{1}{r}$は$\frac{1}{z}$と近似することができます。一方で$e^{ikr}$の項は$k = \frac{2\pi}{\lambda}$の値が大きいため、$r$を$z$で置き換えてしまうと近似の精度が良くありません。
そこで$r$を級数展開して第2項まで使うことを考えます。$r = \sqrt{(x_q - x_p)^2 + (y_q - y_p)^2 + z_q^2}$を第2項まで級数展開すると
$$
r = z_q + \frac{(x_q - x_p)^2 + (y_q - y_p)^2}{2z_q} + \cdots
$$
これを回折積分の式に代入することで
$$
U(x_q, y_q, z_q) = \frac{1}{i\lambda z_q}e^{ikz_q} \int\int_{\Sigma} U(x_p, y_p, 0)e^{ik\frac{(x_q - x_p)^2 + (y_q - y_p)^2}{2z_q}}dx_pdx_q \tag{2}
$$
となります。この近似をフレネル近似といい、式2をフレネル回折積分といいます。
式2の$e$の肩の中身を展開すると $$ U(x_q, y_q, z_q) = \frac{1}{i\lambda z_q}e^{ikz_q}e^{ik\frac{x_q^2 + y_q^2}{2z_q}} \int\int_{\Sigma} U(x_p, y_p, 0)e^{-ik\frac{x_px_q + y_py_q}{z_q}}e^{ik\frac{x_p^2 + y_p^2}{2z_q}}dx_p dx_q \tag{3} $$ となります。式3において$f_x = \frac{x_q}{\lambda z_q}$, $f_y = \frac{y_q}{\lambda z_q}$とおくと $$ \begin{align} U(x_q, y_q, z_q) &= \frac{1}{i\lambda z_q}e^{ikz_q}e^{ik\frac{x_q^2 + y_q^2}{2z_q}} \int\int_{\Sigma} U(x_p, y_p, 0)e^{-ik\frac{x_px_q + y_py_q}{z_q}}e^{ik\frac{x_p^2 + y_p^2}{2z_q}}dx_p dx_q \\ &= \frac{1}{i\lambda z_q}e^{ikz_q}e^{ik\frac{x_q^2 + y_q^2}{2z_q}} \int\int_{\Sigma} U(x_p, y_p, 0)e^{ik\frac{x_p^2 + y_p^2}{2z_q}}e^{-2\pi i(f_x x_p + f_y y_p)}dx_p dx_q \\ &= \frac{1}{i\lambda z_q}e^{ikz_q}e^{ik\frac{x_q^2 + y_q^2}{2z_q}} \mathcal{F}\{U(x_p, y_p, 0)e^{ik\frac{x_p^2 + y_p^2}{2z_q}}\}_{f_x = \frac{x_q}{\lambda z_q}, f_y = \frac{y_q}{\lambda z_q}} \end{align} $$ となり、$U(x_p, y_p, 0)e^{ik\frac{x_p^2 + y_p^2}{2z_q}}$の2次元フーリエ変換で計算できることが分かります。
フレネル近似が成り立つためには、$r$を級数展開した第3項と$k$をかけたものが $$ \frac{k}{8z_q^3}[(x_q - x_p)^2 + (y_q - y_p)^2]^2 \ll 1 $$ を満たす必要があり、これを$z_q$に関して書き直すと $$ z_q^3 \gg \frac{\pi}{4\lambda}[(x_q - x_p)^2 + (y_q - y_p)^2]^2 $$ となります。
式2の$e$の肩の中身を展開すると $$ U(x_q, y_q, z_q) = \frac{1}{i\lambda z_q}e^{ikz_q}e^{ik\frac{x_q^2 + y_q^2}{2z_q}} \int\int_{\Sigma} U(x_p, y_p, 0)e^{-ik\frac{x_px_q + y_py_q}{z_q}}e^{ik\frac{x_p^2 + y_p^2}{2z_q}}dx_p dx_q \tag{3} $$ となります。式3において$f_x = \frac{x_q}{\lambda z_q}$, $f_y = \frac{y_q}{\lambda z_q}$とおくと $$ \begin{align} U(x_q, y_q, z_q) &= \frac{1}{i\lambda z_q}e^{ikz_q}e^{ik\frac{x_q^2 + y_q^2}{2z_q}} \int\int_{\Sigma} U(x_p, y_p, 0)e^{-ik\frac{x_px_q + y_py_q}{z_q}}e^{ik\frac{x_p^2 + y_p^2}{2z_q}}dx_p dx_q \\ &= \frac{1}{i\lambda z_q}e^{ikz_q}e^{ik\frac{x_q^2 + y_q^2}{2z_q}} \int\int_{\Sigma} U(x_p, y_p, 0)e^{ik\frac{x_p^2 + y_p^2}{2z_q}}e^{-2\pi i(f_x x_p + f_y y_p)}dx_p dx_q \\ &= \frac{1}{i\lambda z_q}e^{ikz_q}e^{ik\frac{x_q^2 + y_q^2}{2z_q}} \mathcal{F}\{U(x_p, y_p, 0)e^{ik\frac{x_p^2 + y_p^2}{2z_q}}\}_{f_x = \frac{x_q}{\lambda z_q}, f_y = \frac{y_q}{\lambda z_q}} \end{align} $$ となり、$U(x_p, y_p, 0)e^{ik\frac{x_p^2 + y_p^2}{2z_q}}$の2次元フーリエ変換で計算できることが分かります。
フレネル近似が成り立つためには、$r$を級数展開した第3項と$k$をかけたものが $$ \frac{k}{8z_q^3}[(x_q - x_p)^2 + (y_q - y_p)^2]^2 \ll 1 $$ を満たす必要があり、これを$z_q$に関して書き直すと $$ z_q^3 \gg \frac{\pi}{4\lambda}[(x_q - x_p)^2 + (y_q - y_p)^2]^2 $$ となります。
フラウンホーファー近似
式3において$\frac{x_p^2 + y_p^2}{2z_q}$が十分に小さく、$e^{ik\frac{x_p^2 + y_p^2}{2z_q}} \approx 1$と近似すると
$$
U(x_q, y_q, z_q) = \frac{1}{i\lambda z_q}e^{ikz_q}e^{ik\frac{x_q^2 + y_q^2}{2z_q}} \int\int_{\Sigma} U(x_p, y_p, 0)e^{-ik\frac{x_px_q + y_py_q}{z_q}}dx_p dx_q \tag{4}
$$
となります。この近似をフラウンホーファー近似といい、式4をフラウンホーファー回折積分といいます。
式4において$f_x = \frac{x_q}{\lambda z_q}$, $f_y = \frac{y_q}{\lambda z_q}$とすると $$ \begin{align} U(x_q, y_q, z_q) &= \frac{1}{i\lambda z_q}e^{ikz_q}e^{ik\frac{x_q^2 + y_q^2}{2z_q}} \int\int_{\Sigma} U(x_p, y_p, 0)e^{-2\pi i(f_x x_p + f_y y_p)}dx_p dy_p \\ &= \frac{1}{i\lambda z_q}e^{ikz_q}e^{ik\frac{x_q^2 + y_q^2}{2z_q}} \mathcal{F}\{U(x_p, y_p, 0)\}_{f_x = \frac{x_q}{\lambda z_q}, f_y = \frac{y_q}{\lambda z_q}} \end{align} $$ となり、$U(x_p, y_p, 0)$の2次元フーリエ変換で計算できることが分かります。 フラウンホーファー近似が成り立つためには $$ \frac{k}{2z_q}(x_p^2 + y_p^2) \ll 1 $$ を満たす必要があり、これを$z_q$に関して書き直すと $$ z_q \gg \frac{\pi}{\lambda}(x_p^2 + y_p^2) $$ となります。
式4において$f_x = \frac{x_q}{\lambda z_q}$, $f_y = \frac{y_q}{\lambda z_q}$とすると $$ \begin{align} U(x_q, y_q, z_q) &= \frac{1}{i\lambda z_q}e^{ikz_q}e^{ik\frac{x_q^2 + y_q^2}{2z_q}} \int\int_{\Sigma} U(x_p, y_p, 0)e^{-2\pi i(f_x x_p + f_y y_p)}dx_p dy_p \\ &= \frac{1}{i\lambda z_q}e^{ikz_q}e^{ik\frac{x_q^2 + y_q^2}{2z_q}} \mathcal{F}\{U(x_p, y_p, 0)\}_{f_x = \frac{x_q}{\lambda z_q}, f_y = \frac{y_q}{\lambda z_q}} \end{align} $$ となり、$U(x_p, y_p, 0)$の2次元フーリエ変換で計算できることが分かります。 フラウンホーファー近似が成り立つためには $$ \frac{k}{2z_q}(x_p^2 + y_p^2) \ll 1 $$ を満たす必要があり、これを$z_q$に関して書き直すと $$ z_q \gg \frac{\pi}{\lambda}(x_p^2 + y_p^2) $$ となります。
回折の数値計算
開口$\Sigma$を半径$r = 1~\mathrm{[mm]}$の円形開口、$\lambda = 550~\mathrm{[nm]}$、$z_q = 5~\mathrm{[m]}$として数値計算してみます。
まずはフラウンホーファー近似で計算してみます。近似の条件について確認してみると、$x_p$と$y_p$は最大で$r$になるので
$$
z_q \gg \frac{\pi}{\lambda}r^2 \approx 6.28
$$
となります。仮定の条件を満たしていないですが、このまま計算を進めることにします。
Comments
Post a Comment