光線行列によるJuliaとPythonを用いた光線追跡の計算

import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib

f = 0.9 # lens focal length in mm
rin = 0 # input ray radius in mm : 8um
theta = 20 # angle in degrees
theta = np.deg2rad(theta) # convert to radians
lam = 940e-6 # wavelength in mm

input_matrix = np.array([[rin],[theta]]) # input ray matrix
lens_matrix = np.array([[1, 0], [-1/f, 1]]) # lens matrix

z_list = np.arange(0,10,0.01) # distance from input ray in mm
r_result = []
theta_result = []

for z in z_list:
    if z <= 0.9:  # before lens
        M = np.array([[1, z], [0, 1]])

    else:   # after lens
        before_lens = np.array([[1, 0.9], [0, 1]])
        z_matrix = np.array([[1, z-0.9], [0, 1]])
        M = np.linalg.multi_dot([z_matrix, lens_matrix, before_lens])

    output_matrix = np.dot(M, input_matrix)

    r_out= output_matrix[0][0]
    theta_out = output_matrix[1][0]

    r_result.append(r_out)
    theta_result.append(np.rad2deg(theta_out)) # 角度[deg]で保存

plt.figure(figsize=(16, 5))
plt.subplot(1, 2, 1)
plt.plot(z_list, r_result, label='Ray Radius [mm]')
plt.xlabel("z (mm)")
plt.ylabel("ビーム半径 (mm)")
plt.legend()
plt.grid()

plt.subplot(1, 2, 2)
plt.plot(z_list, theta_result, label='Ray Angle [deg]')
plt.xlabel("z (mm)")
plt.ylabel("広がり角 (deg)")
plt.legend()
plt.grid()
plt.show()

# パラメータ設定
f = 0.9  # レンズの焦点距離 (mm)
rin = 8e-3  # 初期ビームウエスト半径 (mm) : w0
lam = 940e-6  # 波長 (mm)

# 初期qパラメータの計算
q_in = 1 / (1j * lam / (np.pi * rin**2)) 

# ABCD行列の定義
input_matrix = np.array([[q_in], [1]], dtype=complex)  # input ray matrix (q, 1)
lens_matrix = np.array([[1, 0], [-1 / f, 1]], dtype=complex)  # lens matrix

z_list = np.arange(0, 10, 0.01)  
q_result_radius = [] 

# レンズの位置
z_lens = 0.9 # mm

for z_val in z_list:
    if z_val <= z_lens:  # レンズの手前
        # 自由空間伝播
        M_prop = np.array([[1, z_val], [0, 1]], dtype=complex)
        M_total = M_prop
    else:  # レンズの後
        # レンズまでの伝播
        M_before_lens = np.array([[1, z_lens], [0, 1]], dtype=complex)
        # レンズ通過後の残りの距離
        z_after_lens = z_val - z_lens
        M_after_lens_prop = np.array([[1, z_after_lens], [0, 1]], dtype=complex)
        
        # システム全体のABCD行列
        M_total = np.linalg.multi_dot([M_after_lens_prop, lens_matrix, M_before_lens])


    A, B, C, D = M_total[0,0], M_total[0,1], M_total[1,0], M_total[1,1]
    q_out = (A * q_in + B) / (C * q_in + D)
    q_out = 1 / q_out
    #ここまでq_outの計算終わり

    # qパラメータを使ったビーム半径の計算
    q_out_im = np.imag(q_out)
    beam_radius = np.sqrt(lam / (np.pi * np.abs(q_out_im)))
    q_result_radius.append(beam_radius)

# --- 漸近線の計算 ---

# 後ろ側焦点距離 -> 近似線は後ろ側焦点から広がっていく
after_focus_postion = z_lens + f
w0_after_focus = f * lam / (np.pi * rin) 
theta_ff_after_focus = lam / (np.pi * w0_after_focus) 
# 漸近線の計算
z_asymptote = np.linspace(after_focus_postion, z_list[-1], 200)
asymptote_radius = theta_ff_after_focus * np.abs(z_asymptote - after_focus_postion)

# プロット
plt.figure(figsize=(10, 6))
plt.plot(z_list, q_result_radius, label="計算されたビーム半径 (mm)")
plt.plot(z_asymptote, asymptote_radius, '--', label=f"遠方場での漸近線 (θ ≈ {theta_ff_after_focus:.4f} rad)")
# 新しいビームウエストの位置とサイズをプロット上に示す
plt.plot(after_focus_postion, w0_after_focus, 'ro', label=f"後ろ側焦点距離 (z={after_focus_postion:.2f}mm, w0={w0_after_focus*1000:.1f}um)")


plt.xlabel("z (mm)")
plt.ylabel("ビーム半径 (mm)") # Y軸ラベルを修正
plt.legend()
plt.grid(True)
plt.ylim(0, max(q_result_radius)*1.1) # y軸の範囲を調整

plt.axvline(z_lens, color='gray', linestyle=':', label=f"レンズ位置 (z={z_lens}mm)")
plt.legend() # axvline のラベルを表示するために再度呼び出す
plt.show()