はじめに
スペックル(Speckle)をPython使って計算してみるシリーズの第3回です。
前回は実際の試料表面の高さ場を求めて、スペックルの計算を行って高さが1nm程度のばらつきでもスペックルが発生することを確認しました。
今回は試料の反射率特性を取り込んで計算してみたいと思います。
前回記事
取り込む要素について
計算条件は下記とします。
・波長532nm
・試料はSi(n=4.1+0.03i)
・76度(ブリュスター角)からP偏光の斜方照明
その他の条件と計算方法を最初にまとめます。
フレネル係数
界面で反射した場合、光の電場の大きさだけでなく位相も変化します。この2つを合わせて、複素数であらわしたものをフレネル係数と呼びます。
界面の物質の屈折率が分かれば、入射角の情報があれば計算可能です。
今回はP偏光照明で計算してみます。照明光がP偏光であれば、下式で書けます。
$$r_p = \frac{n_2 \cos \theta_i – n_1 \cos \theta_t}{n_2 \cos \theta_i + n_1 \cos \theta_t}$$
屈折率nも$\theta$も一般に複素数です。
$\theta_t$は屈折角なのでスネルの法則から求める必要があります。
反射角について
試料の表面は完全に平坦でないため、各点で傾きを持つことになります。下式で各高さ場における法線を求められます。
$$\vec{n}_{unit} = \frac{1}{\sqrt{(\frac{\partial h}{\partial x})^2 + (\frac{\partial h}{\partial y})^2 + 1}} \left( -\frac{\partial h}{\partial x}, -\frac{\partial h}{\partial y}, 1 \right)$$
照明光の規格化波数ベクトルと上記で内積をとってなす角を求めて、フレネル係数の値を高さ場の各点に割り当てていきます。
斜方照明
前回記事の結果から正反射光の影響が大きいことが分かりました。正反射光を切り離すために試料を斜方から照明してみます。斜方照明を表現するにはどうしたらいいでしょうか？
試料表面の電場分布において、斜方照明方向に付加する位相が増えていくようにしていけばよいです。
詳細はコードを確認ください。
計算コード
まず高さ場を求めます。前回記事から流用です。
import numpy as np
import matplotlib.pyplot as plt
# 1. パラメータ設定
N = 1024
L = 100e-6
wavelength=532e-9
dx = L / N
x = np.linspace(-L/2,L/2,N)
y = np.copy(x)
X,Y = np.meshgrid(x,y)
# 周波数空間のPSDを (N, N//2 + 1) の形状で定義
# kx_r は np.fft.rfftfreq(N, d=dx) を使用
kx_r = np.fft.rfftfreq(N, d=dx)
ky_r = np.fft.fftfreq(N, d=dx)
KXr, KYr = np.meshgrid(kx_r, ky_r) # shape: (N, N//2+1)
# PSDと複素振幅の生成
def Gaussian_PSD(KX, KY, corr=5e-6):
return np.exp(-(KX**2 + KY**2)*corr**2)
PSD_r = Gaussian_PSD(KXr, KYr)
phase = 2 * np.pi * np.random.rand(N, N // 2 + 1)
F = np.sqrt(PSD_r) * np.exp(1j * phase)
# 逆フーリエ変換
h = np.fft.irfft2(F, s=(N, N))
#高さスケールを1nmに調整
h = (h-h.mean())/h.std() * 1e-9次にフレネル係数を返す関数を用意します。
def fresnel_reflection(theta_i, n1=1.0, n2=4.1+0.03j):
"""
theta_i: 入射角 [rad]
n1: 入射側の屈折率 (空気なら1.0)
n2: 試料の複素屈折率 (4.1+0.03j)
"""
# スネルの法則で屈折角を計算
# sin(theta_t) = (n1/n2) * sin(theta_i)
sin_theta_t = (n1 / n2) * np.sin(theta_i)
# 全反射の可能性を考慮してcosを計算
cos_theta_i = np.cos(theta_i)
cos_theta_t = np.sqrt(1 - sin_theta_t**2 + 0j) # 複素数対応
# s偏光とp偏光の反射係数
r_s = (n1 * cos_theta_i - n2 * cos_theta_t) / (n1 * cos_theta_i + n2 * cos_theta_t)
r_p = (n2 * cos_theta_i - n1 * cos_theta_t) / (n2 * cos_theta_i + n1 * cos_theta_t)
return r_p次に高さ場の傾きから照明光との角度を求めていきます。
np.gradientは便利ですね。自前でやろうとするとセル数が変わるとか意外と面倒。
# 勾配の計算
grad_y, grad_x = np.gradient(h,dx,dx)
# 正規化された法線ベクトルの算出
norm = np.sqrt(grad_x**2 + grad_y**2 + 1)
nx, ny, nz = -grad_x / norm, -grad_y / norm, 1 / norm
# 入射光ベクトル k_in (76度入射)
deg = np.pi/180
k_in = np.array([np.sin(76*deg), 0, -np.cos(76*deg)])
# 内積をとって粗面とのなす角theta_surface を算出
cos_theta = -(nx * k_in[0] + ny * k_in[1] + nz * k_in[2])
theta_surface = np.arccos(np.clip(cos_theta, -1, 1))
# フレネル係数マップの作成
r_map = fresnel_reflection(theta_surface) 1nm程度の高さばらつきの粗面でどの程度の角度がつくのか、興味あるので照明光と粗面の間の反射角のずれ(コード中”theta_surface”)を図示してみます。
約±0.2度ほど照明角の76度からばらつきがあるようです。また、相関長さによってばらつきスケールも変わるようです。
フレネル係数の角度依存性が大きい条件だと特性が敏感に変わるかもしれません。
フレネル係数マップが求まったので、位相スクリーン近似(Phase Screen Approximation)などを行い手計算します。
# 入射角に応じた空間周波数シフト量
k_shift = np.sin(76*deg) / wavelength
# 空間周波数平面 KX, KY に位相傾斜を加える
illumination_phase = np.exp(1j * 2 * np.pi * k_shift * X)
# 位相スクリーン近似、高さばらつき分
phi = 4* np.pi / wavelength * h * np.cos(76*deg)
# 試料面の複素電場
field = np.exp(1j * phi) * illumination_phase * r_map
# フーリエ変換（遠方回折 / フラウンフォーファー回折）
# FFTのシフトを行い、強度（絶対値の2乗）を計算
fft_field = np.fft.fft2(field)
speckle = np.abs(fft_field)**2
speckle_shifted = np.fft.fftshift(speckle)
# スペックルコントラストの計算 (標準偏差 / 平均値)
contrast = np.std(speckle) / np.mean(speckle)
#波数を計算
k = 1 / wavelength#532nmの波数[cycle/m]
kx_s = np.fft.fftshift(ky_r)
kx_s /= k
ky_s = np.copy(kx_s)
KX_s,KY_s = np.meshgrid(kx_s,ky_s)
#k<1でマスク
mask = KX_s**2 +KY_s**2 <1
#グラフ化
#上位99.9%の強度を求めて、その値でクリップ
vmax_val = np.percentile(speckle_shifted*mask, 99.9)
clipped_intensity = np.clip(speckle_shifted * mask, 0, vmax_val)
fig, ax = plt.subplots(figsize=(6,6))
c = ax.contourf(KX_s,KY_s, clipped_intensity ,levels=100, cmap="gray",vmax=vmax_val)
ax.set_aspect("equal")
ax.set_xlabel("KX")
ax.set_ylabel("KY")
fig.colorbar(c, ax=ax, label="Intensity",shrink=0.7)
plt.title(f"Speckle Pattern\nContrast: {contrast:.2f}", fontsize=14)
plt.tight_layout()
plt.show() 結果は下図です。正反射光が強すぎるため、上位99.9%の強度でクリップして線形スケールで示しました。
正反射光を表す輝点が右側に来ているので、斜方照明の効果を正しく計算出来ています。
また斜方照明としたことで、正反射光の輝点を右側に配置できたので、より広い空間周波数までみることができます。
正反射光付近にスペックルが発生しています。おそらくガウシアンPSDだと高周波成分がそぎ落とされているためで、高周波成分を許すともっと広範囲にスペックルが広がります。
KY=0の位置に輝線が入っていますが、これは原因がよく分かりませんでした。PSDの低周波成分が多いとこうなるのでしょうか？
まとめ
PSDから高さ場を形成し、表面の傾きによるフレネル係数への影響をとりこみ、ブリュスター角での照明条件でスペックルを計算しました。
斜方照明とすることで、より広い範囲の空間周波数成分まで見ることができます。
状況に合わせて以下の効果を取り入れると面白そうです。
・高さ場のばらつきを変える
・PSDを非対称なものを使って、スペックルへの影響を計算する
・試料(屈折率)や波長を変えてみる
・吸収の強い物質で計算する
・照明光の形状の効果を取り入れてみる(今は100μm角が均一に照明)
次回記事
シミュレーション方法はここで一区切りとします。
次回からはスペックルの統計的性質について別シリーズでまとめていこうと思います。
参考文献
スペックルに関係する論文の引用元や海外コードの解説で必ず登場するグッドマンの不朽の名著。英語の専門用語や理論的なニュアンスをそのまま正確に理解したい、本格志向の技術者におすすめの一冊です。
散乱光学について体系的に学びたい方には、John C. Stover 著 Optical Scattering: Measurement and Analysis がおすすめです。実際の測定を意識した解説が多く、表面粗さや散乱理論、BRDFなど、研究や開発で役立つ内容が幅広くまとまっています。
数式だけでなく物理的な解釈にも触れられているため、散乱現象をより深く理解したい方に適した一冊です。
回折計算（FFT）をコードに落とし込む際の物理的・数学的な根拠を最も丁寧に解説してくれる不朽の名著です。光学シミュレーションで迷ったときの道標として、手元に置きたい一冊。
コメント