はじめに
以前の記事において、\( q\)パラメータを用いた光線行列計算と幾何光学計算について紹介した。しかし、これらの計算は前者においては基本ガウシアンビームしかカバーされず、後者にしては有限のビーム幅を持つレーザーにおいては計算が不可能だった。
そこで、今回はこれらの問題を解決するために、幾何光学計算とモンテカルロ法を用いて\(q\)パラメータを用いた光線行列計算を再現し、さらに任意のレーザーでシミュレーションを行えるようにした。
高次ガウシアンモードにおける\(q\)パラメータの限界
\( q\)パラメータを用いたガウシアンビームの計算には高次ガウシアンモードへの拡張といった意味では限界がある。この理由としては、\(q\)パラメータは基本ガウシアンビームを想定した計算であること。また、高次ガウシアンモードの\(q\)パラメータへの拡張を考察した記事おいて[平等1998]、\(q\)パラメータは次のように表される;
\begin{equation}
\frac{1}{q(z)} = \frac{1}{R(z)} – i \frac{M^2 \lambda}{\pi w_{high}(0)^2}
\end{equation}
この式を見ると\(M^2\)を大きくすると計算におけるビーム径\(w_{high}(0)\)は小さくなり、逆もまたしかりである。このことから同じビーム径で異なる\(M^2\)を持つビームの計算が不可能であることがわかる。
そこで本記事ではモンテカルロ法を用いて力業で解こうといった考えである。
モンテカルロ法
モンテカルロ法(Monte Carlo method)[wikipedia]はフォン・ノイマンらの核分裂における中性子の拡散現象のシミュレーションへの応用に端を発するといわれている。特に光学との親和性は高く、zemaxでのノンシーケンシャルモードのように光線を多く飛ばす計算ではモンテカルロ法が用いられている(と思います)[参考]。特に今回の検討にあたってはシミュレーション光学といった教科書を参考にした。[牛山2003]
モンテカルロ法とは多くの場合、疑似乱数を用いてランダムに事象を解析し、解析的に解けない問題に使われていることが多いです。そこで今回は高次ガウシアンモードの解析にむけてプログラムを書いてみました。
特に照明光学系に利用する場合は以下の手順で考えていく[牛山2003]。
①疑似乱数の発生分布を決定する(今回はガウシアン分布)
②光源面上に、疑似乱数による座標を決定する
③同様に疑似乱数によって出射角を決定する
④発生した光線を被照明面まで光線追跡する
⑤被照明面上での微小区分に到達する光線本数をカウントする
⑥ ②から⑤までを任意回数行う
⑦放射照度を計算する
特に今回はレーザーであり、ビーム径だけが知りたかったので①-④⑥のみを行うと方針立てて計算を行っていきます。
また、Juliaでの実装は可読性を第一に考え、プログラムの速さは一切考慮していません。
for文によるABCD行列の計算はfor文内で解かずに、あらかじめ解析的に解きmap関数で代入処理したほうが早いと思います。
モンテカルロ法を用いた光線追跡
光源の用意
今回モンテカルロ法による計算を行うための乱数発生部分は光源部分にある。一般的な証明光学系と異なり、レーザー光源を考えるために様々な因子を簡略化/具体化して計算が行えるのが最大の特徴です。光線追跡を行う光線行列計算の式を思い出す[過去記事]と
\begin{equation} \left[ \begin{array}{c} r_o\\ \theta_i \end{array} \right] = \left[ \begin{array}{c c} A & B\\ C & D \end{array} \right] \left[ \begin{array}{c} r_i\\ \theta_i \end{array} \right] \end{equation}
と表すことができる。\(r_i, \theta_i \)が光源に因子となる。これらの値は光源観点からNFP(Near Field Pattern)とFFP(Far Field Pattern)によって決定される。 NFPの分布を\(D_{NFP}\)、DFPの分布を\(D_{FFP}\)としたときに次式のように光源の因子が決められる。
\begin{equation}
r_i = \sqrt{x_i^2+y_i^2} \sim D_{NFP}\\
\theta_i \sim D_{FFP}
\end{equation}
より具体的に、今回はガウシアンビームを仮定して話を進めていきます。まず、Juliaを含めた一般的な標準ガウス分布は次のように表せます。
\begin{equation}
f(x) = \frac{1}{\sqrt{2\pi\sigma}}\exp \left\{ – \frac{(x-\mu)^2}{2\sigma^2} \right\}
\end{equation}
また、ビームの強度分布が次式で表せます[Yariv2010]。
\begin{equation}
I(x) \propto \exp \left\{ – \frac{(x-\mu)^2}{w^2} \right\}
\end{equation}
この式からビーム半径\( w \)を用いて\(r_i\)の従うガウス分布の標準偏差\(\sigma_{NFP}\)を表すと
\begin{equation}
\sigma_{NFP} = \sqrt{2}w
\end{equation}
同様に広がり角(全角)は以下の様に表せ[平等1998]
\begin{equation}
\Theta = 2\theta = M^2 \frac{4\lambda}{2 \pi w}
\end{equation}
\(\theta_i\)の従うガウス分布の標準偏差\(\sigma_{FFP}\)は
\begin{equation}
\sigma_{FFP} = \sqrt{2}\theta = \frac{\Theta}{\sqrt{2}}
\end{equation}
となる。これをJuliaで実装すると
using Distributions
# パラメータ設定
rin = 8e-3 # 初期ビームウエスト半径 (mm) : w0
λ = 940e-6 # 波長 (mm)
M² = 1.0
#calculate spread angle
θ = 4λ*M² / (2π * rin) # ビームの発散角度 (rad)
sigma_d = rin * sqrt(2) # 1/e²幅=rinのときの標準偏差
N = 2000 # サンプル数
dist = Normal(0, sigma_d)
ri_distibution = rand(dist, N)
sigma_theta = θ / sqrt(2) # ビームの発散角度の標準偏差
dist2 = Normal(0,sigma_theta) # ビームの発散角度の分布
theta_distrbution = rand(dist2, N)基本ガウシアンビームにおける光線追跡
今回もコリメータ光学系をいったん考える。基本は過去記事と同じの場合を考える。
\begin{equation} \left[ \begin{array}{c} r_o\\ \theta_i \end{array} \right] = \left[ \begin{array}{c c} 1 & z\\ 0 & 1 \end{array} \right] \left[ \begin{array}{c c} 1 & 0\\ -\frac{1}{f} & 1 \end{array} \right] \left[ \begin{array}{c c} 1 & f\\ 0 & 1 \end{array} \right] \left[ \begin{array}{c} r_i\\ \theta_i \end{array} \right] \end{equation}
以上のような立式が可能となり、\(r_i, \theta_i \)は上記で定めたパラメータを代入する。そのコードを以下に示す。また、基本ガウシアンビーム(\(M^2 = 1\) )のときは\(q\)パラメータを用いた計算が理論的に正しいと確認されているので、ベンチマークとして\(q\)パラメータを用いた計算も並行して行っている。
using PyPlot
using Distributions
using PyCall
plt = PyPlot
PyCall.PyDict(matplotlib["rcParams"])["font.family"] = "Meiryo";
# パラメータ設定
f = 0.9 # レンズの焦点距離 (mm)
rin = 8e-3 # 初期ビームウエスト半径 (mm) : w0
λ = 940e-6 # 波長 (mm)
M² = 1.0
#calculate spread angle
θ = 4λ*M² / (2π * rin) # ビームの発散角度 (rad)
############# Monte Carlo法によるビーム半径の計算 ##############
sigma_d = rin * sqrt(2) # 1/e²幅=rinのときの標準偏差
N = 2000 # サンプル数
dist = Normal(0, sigma_d)
ri_distibution = rand(dist, N)
sigma_theta = θ / sqrt(2) # ビームの発散角度の標準偏差
dist2 = Normal(0,sigma_theta) # ビームの発散角度の分布
theta_distrbution = rand(dist2, N)
z_list_MC = range(0.0,10.0,step = 0.1) # distance from input ray in mm
r_results = zeros(length(ri_distibution), length(z_list_MC))
theta_results = zeros(length(ri_distibution), length(z_list_MC))
for (i,( rin_sample, theta_sample )) in enumerate(zip(ri_distibution, theta_distrbution))
input_matrix = [rin_sample ; theta_sample] # input ray matrix
lens_matrix = [1 0; -1/f 1] # lens matrix
for (j, z) in enumerate(z_list_MC)
if z <= 0.9 # before lens
M = [1 z; 0 1]
else # after lens
before_lens = [1 0.9; 0 1]
z_matrix = [1 z-0.9; 0 1]
M = z_matrix * lens_matrix * before_lens
end
output_matrix = M * input_matrix
r_out = output_matrix[1]
theta_out = output_matrix[2]
r_results[i, j] = r_out
theta_results[i, j] = rad2deg(theta_out)
end
end
# calculation of 1/e^2 radius
mean_rms = sqrt.(mean(r_results.^2, dims=1))[:]
mean_r_1e2 = mean_rms / sqrt(2)
#############ここまでがMC法##############
#############ここからがqパラメータを使ったビーム半径の計算##############
# 初期qパラメータの計算
q_in = 1 / (1im * lam / (π * rin^2))
# ABCD行列の定義
input_matrix = [q_in ; 1] # input ray matrix (q, 1)
lens_matrix = [1 0; -1 / f 1] # lens matrix
z_list = range(0.0,10.0,step = 0.01)
q_result_radius = []
# レンズの位置
z_lens = 0.9 # mm
for z_val in z_list
if z_val <= z_lens # レンズの手前
# 自由空間伝播
M_prop = [1 z_val; 0 1]
M_total = M_prop
else # レンズの後
# レンズまでの伝播
M_before_lens = [1 z_lens; 0 1]
# レンズ通過後の残りの距離
z_after_lens = z_val - z_lens
M_after_lens_prop = [1 z_after_lens; 0 1]
# システム全体のABCD行列
M_total = M_after_lens_prop * lens_matrix * M_before_lens
end
A, B, C, D = M_total[1,1], M_total[1,2], M_total[2,1], M_total[2,2]
q_out = (A * q_in + B) / (C * q_in + D)
q_out = 1 / q_out
#ここまでq_outの計算終わり
# qパラメータを使ったビーム半径の計算
q_out_im = imag(q_out)
beam_radius = √(lam / (π * abs(q_out_im)))
push!(q_result_radius,beam_radius)
end
#############ここまでがqパラメータを使ったビーム半径の計算##############
plt.figure(figsize=(10, 6))
plt.plot(z_list, q_result_radius, label="q-parameter (理論値)")
plt.plot(z_list_MC, mean_r_1e2, label="Monte Carlo", linestyle="--")
plt.xlabel("z (mm)")
plt.ylabel("ビーム半径 (mm)")
plt.legend()
plt.grid()
display(gcf())
この実行結果が以上のようなグラフとなる。理論値とモンテカルロ法による計算値が非常に一致していることがわかる。また、モンテカルロ法の計算結果については
\begin{equation}
RMS = \sigma = \sqrt{2}w
\end{equation}
の関係を用いた。
また、この結果はSample数(N)を増やせばより厳密に一致していく。今回のN数(=2000)では600msほどの実行時間だった。
高次ガウスモードでの光線追跡
上記の計算に対して高次ガウスモードでは\(M^2\)の値を変更する。今回は2を代入して計算を行った。
using PyPlot
using Distributions
using PyCall
plt = PyPlot
PyCall.PyDict(matplotlib["rcParams"])["font.family"] = "Meiryo";
# パラメータ設定
f = 0.9 # レンズの焦点距離 (mm)
rin = 8e-3 # 初期ビームウエスト半径 (mm) : w0
λ = 940e-6 # 波長 (mm)
M² = 2.0
#calculate spread angle
θ = 4λ*M² / (2π * rin) # ビームの発散角度 (rad)
############# Monte Carlo法によるビーム半径の計算 ##############
sigma_d = rin * sqrt(2) # 1/e²幅=rinのときの標準偏差
N = 2000 # サンプル数
dist = Normal(0, sigma_d)
ri_distibution = rand(dist, N)
sigma_theta = θ / sqrt(2) # ビームの発散角度の標準偏差
dist2 = Normal(0,sigma_theta) # ビームの発散角度の分布
theta_distrbution = rand(dist2, N)
z_list_MC = range(0.0,10.0,step = 0.1) # distance from input ray in mm
r_results = zeros(length(ri_distibution), length(z_list_MC))
theta_results = zeros(length(ri_distibution), length(z_list_MC))
for (i,( rin_sample, theta_sample )) in enumerate(zip(ri_distibution, theta_distrbution))
input_matrix = [rin_sample ; theta_sample] # input ray matrix
lens_matrix = [1 0; -1/f 1] # lens matrix
for (j, z) in enumerate(z_list_MC)
if z <= 0.9 # before lens
M = [1 z; 0 1]
else # after lens
before_lens = [1 0.9; 0 1]
z_matrix = [1 z-0.9; 0 1]
M = z_matrix * lens_matrix * before_lens
end
output_matrix = M * input_matrix
r_out = output_matrix[1]
theta_out = output_matrix[2]
r_results[i, j] = r_out
theta_results[i, j] = rad2deg(theta_out)
end
end
# calculation of 1/e^2 radius
mean_rms = sqrt.(mean(r_results.^2, dims=1))[:]
mean_r_1e2 = mean_rms / sqrt(2)
#############ここまでがMC法##############
#############ここからがqパラメータを使ったビーム半径の計算##############
# 初期qパラメータの計算
q_in = 1 / (1im * lam / (π * rin^2))
# ABCD行列の定義
input_matrix = [q_in ; 1] # input ray matrix (q, 1)
lens_matrix = [1 0; -1 / f 1] # lens matrix
z_list = range(0.0,10.0,step = 0.01)
q_result_radius = []
# レンズの位置
z_lens = 0.9 # mm
for z_val in z_list
if z_val <= z_lens # レンズの手前
# 自由空間伝播
M_prop = [1 z_val; 0 1]
M_total = M_prop
else # レンズの後
# レンズまでの伝播
M_before_lens = [1 z_lens; 0 1]
# レンズ通過後の残りの距離
z_after_lens = z_val - z_lens
M_after_lens_prop = [1 z_after_lens; 0 1]
# システム全体のABCD行列
M_total = M_after_lens_prop * lens_matrix * M_before_lens
end
A, B, C, D = M_total[1,1], M_total[1,2], M_total[2,1], M_total[2,2]
q_out = (A * q_in + B) / (C * q_in + D)
q_out = 1 / q_out
#ここまでq_outの計算終わり
# qパラメータを使ったビーム半径の計算
q_out_im = imag(q_out)
beam_radius = √(lam / (π * abs(q_out_im)))
push!(q_result_radius,beam_radius)
end
#############ここまでがqパラメータを使ったビーム半径の計算##############
plt.figure(figsize=(10, 6))
plt.plot(z_list, q_result_radius, label="Ray tracing by q-parameter (theory)")
plt.plot(z_list_MC, mean_r_1e2, label="Ray tracing by Monte Carlo Method (This work)", linestyle="--")
plt.xlabel("z (mm)")
plt.ylabel("ビーム半径 (mm)")
plt.legend()
plt.grid()
display(gcf())
近傍においては基本ガウシアンビームに対して広がり角が大きくなっていることがしっかり確認され、遠方においてはコリメート後は基本ガウシアンビームに対してほとんど広がり角が変わらないことが確認された。
コメント