導入
散乱光分布を計算していて、ふとこれレンズで結像するとどう見えるんだろう？と思ったので計算してみたシリーズの第6回。最終回です。
今回は偏光の影響を取り込み、これまでの総まとめとします。
過去記事
偏光の計算方法
前回pyGDMで散乱光の電場を求めました。
この時x,y,zの3成分が得られました。
この座標系は局所座標でなく、グローバル座標であり、x,y,zは直行しています。
FFTは線形な演算なので、x,y,z成分に対してそれぞれ結像計算を行い、最後に強度を合成すればよいのです。
計算コード
散乱光のデータは"E"というリスト形式で得られました。
詳細は前回記事参照ください。
https://qiita.com/opticnomad/items/d8fe078e0b276eea5008
Eは所与として散乱計算は省略します。
前回記事までのコードを活用しています。
import numpy as np
import matplotlib.pyplot as plt
plt.style.use("ggplot")
from scipy.interpolate import griddata
#x,y,z成分の電場取得、整形
Ex =E[2][:,0].reshape(E[0].shape)
Ey =E[2][:,1].reshape(E[0].shape)
Ez= E[2][:,2].reshape(E[0].shape)
#後の計算用
theta =np.copy(E[0])
kx = np.sin(E[0])*np.cos(E[1])
ky = np.sin(E[0])*np.sin(E[1])
#電場の大きさをグラフ化
fig, axes = plt.subplots(1, 3, figsize=(15,5))
components = [Ex, Ey, Ez]
titles = ["Ex", "Ey", "Ez"]
for ax, comp, title in zip(axes, components, titles):
im = ax.contourf(
kx,
ky,
np.abs(comp),
levels=50)
ax.set_title(title)
ax.set_xlabel("kx")
ax.set_ylabel("ky")
ax.set_xlim(-1,1)
ax.set_ylim(-1,1)
ax.set_aspect('equal')
fig.colorbar(im, ax=ax)
plt.tight_layout()
plt.show()
各方向における電場の大きさはこんな感じです。
検出NA=0.9、収差は無いものとします。
前回言及したように瞳面では、512x512のうち瞳の半径は512の1/4くらいにしないとうまく結像計算ができません。FFTは周期境界条件を仮定するため、瞳を配列いっぱいに広げると折り返しの影響がでます。
そのため、NAのほかにmask_ratioという変数を導入します。
N=512
NA=0.9
mask_ratio =0.5
#kを等間隔で作成
kx_im = np.linspace(-1/mask_ratio*NA,1/mask_ratio*NA,N)
ky_im = np.linspace(-1/mask_ratio*NA,1/mask_ratio*NA,N)
KX_im,KY_im =np.meshgrid(kx_im,ky_im)
#マスク
mask = (np.sqrt(KX_im**2+KY_im**2) < NA).astype(int)
#電場強度を補間
points = np.column_stack([kx.ravel(),ky.ravel()])
values_x = Ex.ravel()
values_y = Ey.ravel()
values_z = Ez.ravel()
values_theta=theta.ravel()
#電場とthetaを補間
Ex_interpolate = griddata(points,values_x,(KX_im,KY_im),method='linear')
Ey_interpolate = griddata(points,values_y,(KX_im,KY_im),method='linear')
Ez_interpolate = griddata(points,values_z,(KX_im,KY_im),method='linear')
theta_interpolate = griddata(points,values_theta,(KX_im,KY_im),method='linear')
#nanが含まれるので0にする
Ex_interpolate[np.isnan(Ex_interpolate)]=0
Ey_interpolate[np.isnan(Ey_interpolate)]=0
Ez_interpolate[np.isnan(Ez_interpolate)]=0
theta_interpolate[np.isnan(theta_interpolate)]=0
#瞳関数とアプラナート条件を適用
weight =np.cos(theta_interpolate)**0.5
Ex_interpolate *= mask*weight
Ey_interpolate *= mask*weight
Ez_interpolate *= mask*weight
#FFTを実行
Ex_img = np.fft.fftshift(
np.fft.ifft2(
np.fft.ifftshift(Ex_interpolate)
)
)
Ey_img = np.fft.fftshift(
np.fft.ifft2(
np.fft.ifftshift(Ey_interpolate)
)
)
Ez_img = np.fft.fftshift(
np.fft.ifft2(
np.fft.ifftshift(Ez_interpolate)
)
)
FFTまで終わりました。
結像面のスケーリングを行います。レンズ焦点距離は100mmとします。
#像面のスケーリング
f = 100 #レンズ焦点距離
lam = 0.532 *10**-3#波長
#瞳面での1ピクセルの物理的大きさを求める
dx = f * NA / (N/2* mask_ratio)
k = np.fft.fftfreq(N, d=dx)
x_real = lam * f * k
x_real = np.fft.fftshift(x_real)
y_real = np.copy(x_real)
X_real,Y_real=np.meshgrid(x_real,y_real)
#グラフ化
plt.figure(figsize=(6,5))
pix =20
plt.contourf(X_real[-pix+N//2:pix+N//2,-pix+N//2:pix+N//2],Y_real[-pix+N//2:pix+N//2,-pix+N//2:pix+N//2],np.abs(Ex_img[-pix+N//2:pix+N//2,-pix+N//2:pix+N//2])**2+np.abs(Ey_img[-pix+N//2:pix+N//2,-pix+N//2:pix+N//2])**2+np.abs(Ez_img[-pix+N//2:pix+N//2,-pix+N//2:pix+N//2])**2,levels=50)
plt.xlabel("X(mm)")
plt.ylabel("Y(mm)")
plt.gca().set_aspect('equal')
plt.show()
集光スポットはこんな感じです。
像がボケて左右方向に広がっていますが、
元の散乱光の強度分布を反映しています。
例えば電場を一様分布(等方的)としてみると像はもっとシャープに見えます。
散乱光強度分布によるボケの影響もあることがわかります。
#これで電場置き換え
Ex =np.ones_like(E[0])
Ey =np.ones_like(E[0])
Ez =np.ones_like(E[0])
偏光の取り扱い
最初に思い違いをしていました。
瞳面ではレンズで曲げられてz方向に進みます。
光の波数ベクトルと偏光振動面は直行するため、瞳面ではx,y方向に電場振動面がくるため、
光を3次元行列で回転させて瞳面上での電場を出した後にFFTと思っていました。
コーディングがとても大変そうです。
実際はx,y,zはグローバル座標で直行して独立、更にpyGDMが3成分全部計算してくれるため、それぞれ計算すれば十分でした。
最初は理解できなかったのですが、FFT結像はフラウンフォーファー回折と同じ形で書けて、無限遠での回折現象を扱っていると思うと腹落ちしました。
最後に
FFT結像のシリーズはここまでです。
収差の入れ方を工夫する、散乱光の計算に意味を持たせる、瞳面で偏光フィルターを入れてみる、位相操作をするなどすれば色々面白い現象を計算出来そうです。
結像計算はフラウンフォーファー回折によるFFT計算が必要ですが、
レーザ光の伝播を計算する時は、角スペクトル法という別の方法が必要です。以下のページで解説しています。