導入

散乱光分布を計算していて、ふとこれレンズで結像するとどう見えるんだろう？と思ったので計算してみたシリーズの第二回。

「瞳面上に電場を並べてFFTすればいいだけでは？」

と思ってやっていました。

しかし実際にやってみると、いくつか問題が出てきます。

(θ,φ)で等間隔にサンプリングしたデータは、そのままでは瞳面で等間隔にならない
そのため補間が必要になる
さらに、電場の強度に補正が必要になる

つまり、

「方向空間での分布」を「空間周波数空間（瞳面）」に正しく変換する必要がある

ということです。

このあたりが、実際に手を動かしてみて一番ハマったポイントでした。

前回の記事

計算の全体フロー

計算全体の流れはこんな感じです。

散乱光のダミーデータを準備(方向空間)　
瞳面座標に変換(k空間)
等間隔グリッドに再配置(FFTのため)
強度補正
FFTして結像

FFTかけるだけでしょ、と思っていましたが意外とやることが多いです。
1個ずつ説明するので安心してください。

1. 散乱光のダミーデータを準備(方向空間)　

散乱光の計算結果のダミーデータを作ります。
散乱光計算は通常、角度に対して等間隔で散乱光強度を計算します。

import numpy as np
from scipy.interpolate import griddata #補間用
import matplotlib.pyplot as plt
plt.style.use("ggplot")

#度から弧度法への変換係数
deg = np.pi / 180

#dummyデータ生成
theta = np.linspace(0,90*deg,90)
phi = np.linspace(0,360*deg,360)
theta,phi = np.meshgrid(theta,phi)
E = np.cos(theta)

#コンター図作成
theta_cont = (90*deg-theta)/deg #散乱角から仰角に。コンター図用に度に直す
phi_cont = phi / deg#コンター図用に度に直す
plt.figure(figsize=(6,5))
plt.contourf(phi_cont,theta_cont,E,levels=50)
plt.xlabel("azimuth angle(deg)")
plt.ylabel("elevation angle(deg)")
plt.show()

上方への散乱光強度が強いです。
thetaは散乱角なので直上方向が基準ですが、コンター図にする際は90度から引いて
仰角に変換すればメルカトル図法っぽくなるので見やすいと思います。

2. 瞳面座標に変換(k空間)

theta,phiからk空間に座標を変換します。

#theta,phiからk空間に変換
kx = np.sin(theta)*np.cos(phi)
ky = np.sin(theta)*np.sin(phi)

plt.figure(figsize=(6,5))
plt.contourf(kx,ky,E,levels=50)
plt.xlabel("kx")
plt.ylabel("ky")
plt.xlim(-1,1)
plt.ylim(-1,1)
plt.gca().set_aspect('equal')
plt.show()

電場は何も変えていないのですが、
座標をk空間に変換したことで輝点を真上から見たようなコンター図になります。

ここで重要なのは、

👉 k空間の座標がそのまま「瞳面の位置」に対応するという点です。

つまり、光の進行方向（θ,φ）は、
レンズによって「どの位置を通るか」に変換されます。

この対応付けが、結像計算の出発点になります。

3. 等間隔グリッドに再配置(FFTのため)

(θ,φ)で等間隔にサンプリングしたデータをそのままk空間に変換すると、
瞳面上では点の密度が不均一になります。

しかしFFTは等間隔サンプリングを前提としているため、

一度、等間隔グリッドに再配置する必要があります。

# FFTするデータ数、2の累乗にする。
N=512 
#kを等間隔で作成
kx_im = np.linspace(-1,1,N)
ky_im = np.linspace(-1,1,N)
KX_im,KY_im =np.meshgrid(kx_im,ky_im)

#電場強度を補間
points = np.column_stack([kx.ravel(),ky.ravel()])
values = E.ravel()

#後でtheta使うので補間
values_theta=theta.ravel()

E_interpolate = griddata(points,values,(KX_im,KY_im),method='linear')
theta_interpolate = griddata(points,values_theta,(KX_im,KY_im),method='linear')

#nanが含まれるので0にする
E_interpolate[np.isnan(E_interpolate)]=0
theta_interpolate[np.isnan(theta_interpolate)]=0

E_interpolateが512x512の等間隔グリッドに補間したものです。
見た目は同じなので図示省略します。
thetaも次で使うので同じように等間隔グリッドに補間しています。

4. 強度補正

ここが今回の中で一番重要なポイントです。
少し直感と違う処理になるので、順番に説明します。

レンズで拾える範囲に電場を制限する必要があります。
ここでは仮にNA = 0.9以下が検出可能、それ以上検出できないので0を掛けます。
(θ,φ)からk空間に変換したのでヤコビアンである$\sin\theta$を掛けます
直感的にはθが大きい所の光が瞳面に変換した際に密度が高くなるのでそれを補正します。
レンズが適切に収差補正されている(アプラナート条件)場合は、$\sqrt{\cos\theta}$を掛ける必要があります。
　検出NAが小さいときは省略してもあまり影響がありません。

これらの補正は、「方向空間でのサンプリングの歪み」を補うためのものです。

NA=0.9
mask = (np.sqrt(KX_im**2+KY_im**2) < NA).astype(int)
weight =np.sin(theta_interpolate)*np.cos(theta_interpolate)**0.5
E_interpolate *= mask*weight

図示してみるとこんな感じです。
ここ確認する必要はあまりないんですが、逐一みるのが好きです。

5. FFTして結像

後はFFTするだけです。
瞳面の電場分布と像面の電場はフーリエ変換の関係にあります。

今回は「瞳面 → 像面」への変換なので、逆フーリエ変換を用います。

結像の強度のみを考えたいのであれば、FFTと逆FFTは区別しなくてよいですが位相も考慮したい時は区別する必要があります。

E_img = np.fft.ifft2(E_interpolate)
#FFT結果を並べなおす
E_img = np.fft.fftshift(E_img)

nx=ny= N
x = np.fft.fftfreq(nx, d = 1/N)
y = np.fft.fftfreq(ny, d = 1/N)
x = np.fft.fftshift(x)
y = np.fft.fftshift(y)
X, Y = np.meshgrid(y, x)

#図示用のコード
plt.figure(figsize=(6,5))
plt.contourf(X,Y,np.abs(E_img)**2,levels=50)
plt.xlabel("X")
plt.ylabel("Y")
plt.gca().set_aspect('equal')
plt.show()

ちょっと見づらいですが、中央に小さい点が見えます。

散乱光はある一点から発生していて、それを結像したので小さいな1点に結像されています。

ちょっと見づらいので、logをとってプロットしてみるのも良いです。

ここまでずっと電場を計算していましたが、実際に観測できるのはパワーです。
そのため電場の絶対値の2乗をプロットします。
電場は一般に複素数なのでnp.absを付けるのが安全です。

ここで疑問

これで目標達成！と思ったのですが、図を見た時にいろんな疑問がわいてきました。

結像の横軸の値は意味は？どれくらいの大きさに結像されているの？
レンズの特性やフォーカスずれの影響はどう取り込めばいいの？
位相とか偏光の影響はどう取り込めばいいの？

今回は結像計算の最小セット、ここからが本番です。

次回は結像画像の大きさをどう考えればいいかについて解説したいと思います。

雑談

瞳面をうまく説明できません。

教科書の説明は正しいですが直感的に理解しづらく、
初心者に向けて説明する時はいつも困ってしまいます。

こう説明するといいよっていうのがあれば、教えてください。