"""
波の干渉シミュレーションモジュール。
このモジュールは、2つの光源から発生する波の干渉を計算し、3Dアニメーションとして描画します。
関連リンク:
:doc:`interfere_usage`
"""
import argparse
import os
import numpy as np
import multiprocessing
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.animation import FuncAnimation
# 波の干渉計算
[ドキュメント]
def interference_wave(x, y, t, k, omega, source1, source2):
"""
2つの光源から発生する波の干渉を計算する。
指定された空間座標 (x, y) および時間 t における、2つの光源からの波の重ね合わせを計算します。
それぞれの波の複素振幅を求め、それらの和を干渉結果として返します。
:param x: (numpy.ndarray) 空間のX座標を表すメッシュグリッド。
:param y: (numpy.ndarray) 空間のY座標を表すメッシュグリッド。
:param t: (float) 現在の時間。
:param k: (float) 波数。
:param omega: (float) 角周波数。
:param source1: (numpy.ndarray) 光源1の位置座標 [x, y]。
:param source2: (numpy.ndarray) 光源2の位置座標 [x, y]。
:returns: (numpy.ndarray) 2つの波の干渉結果である複素振幅の配列。
"""
# 光源1からの波
r1 = np.sqrt((x - source1[0])**2 + (y - source1[1])**2)
wave1 = np.exp(1j * (k * r1 - omega * t))
# 光源2からの波
r2 = np.sqrt((x - source2[0])**2 + (y - source2[1])**2)
wave2 = np.exp(1j * (k * r2 - omega * t))
# 干渉結果
total_wave = wave1 + wave2
return total_wave
[ドキュメント]
def parse_args():
parser = argparse.ArgumentParser(
description="2つの光源から発生する波の干渉を3Dアニメーションとして描画します。"
)
# パラメータ設定
parser.add_argument("--k", type=float, default=2 * np.pi,
help="波数 (単位: m^-1)")
parser.add_argument("--omega", type=float, default=2 * np.pi,
help="角周波数 (単位: rad/s)")
parser.add_argument("--c", type=float, default=1.0,
help="波動の伝播速度 (単位: m/s)")
parser.add_argument("--t_max", type=float, default=10.0,
help="最大時間")
parser.add_argument("--t_steps", type=int, default=200,
help="時間ステップ数")
parser.add_argument("--xmin", type=float, default=-5.0,
help="空間の範囲（x軸）の最小値")
parser.add_argument("--xmax", type=float, default=5.0,
help="空間の範囲（x軸）の最大値")
parser.add_argument("--nx", type=int, default=100,
help="空間の範囲（x軸）の分割数")
parser.add_argument("--ymin", type=float, default=0.0,
help="空間の範囲（y軸）の最小値")
parser.add_argument("--ymax", type=float, default=5.0,
help="空間の範囲（y軸）の最大値")
parser.add_argument("--ny", type=int, default=100,
help="空間の範囲（y軸）の分割数")
parser.add_argument("--source1_x", type=float, default=-1.0,
help="光源1のx座標")
parser.add_argument("--source1_y", type=float, default=0.0,
help="光源1のy座標")
parser.add_argument("--source2_x", type=float, default=1.0,
help="光源2のx座標")
parser.add_argument("--source2_y", type=float, default=0.0,
help="光源2のy座標")
parser.add_argument("--interval", type=int, default=100,
help="アニメーションのフレーム間隔 (ms)")
parser.add_argument("--zmin", type=float, default=-5.0,
help="z軸の最小値")
parser.add_argument("--zmax", type=float, default=5.0,
help="z軸の最大値")
parser.add_argument("--contour_y", type=float, default=5.0,
help="強度等高線を描画するy座標")
return parser.parse_args()
[ドキュメント]
def main():
args = parse_args()
# パラメータ設定
k = args.k # 波数 (単位: m^-1)
omega = args.omega # 角周波数 (単位: rad/s)
c = args.c # 波動の伝播速度 (単位: m/s)
t_max = args.t_max # 最大時間
t_steps = args.t_steps # 時間ステップ数
x_range = np.linspace(args.xmin, args.xmax, args.nx) # 空間の範囲（x軸）
y_range = np.linspace(args.ymin, args.ymax, args.ny) # 空間の範囲（y軸）y >= 0に制限
# 空間メッシュグリッドの生成
x, y = np.meshgrid(x_range, y_range)
# 光源の位置 (光源1, 光源2)
source1 = np.array([args.source1_x, args.source1_y]) # 光源1の位置 (-1, 0)
source2 = np.array([args.source2_x, args.source2_y]) # 光源2の位置 (1, 0)
# プロット設定
fig = plt.figure(figsize=(10, 7))
ax = fig.add_subplot(111, projection='3d')
# 光源位置に球を描画
ax.scatter(source1[0], source1[1], 0, color="red", s=100, label="Source 1") # 光源1
ax.scatter(source2[0], source2[1], 0, color="blue", s=100, label="Source 2") # 光源2
# 時間ステップにおける波の干渉をプロット
t = 0 # 初期時間
X, Y = np.meshgrid(x_range, y_range)
# アニメーションの作成
def animate(t):
"""
アニメーションの各フレームを更新・描画する。
与えられた時間 t に基づいて干渉波とその強度を計算し、3Dプロット上の
波の振幅サーフェスおよび特定のY座標（y=5）での強度等高線を更新します。
:param t: (float) 現在の時間（アニメーションのフレームに応じた時間）。
:returns: (None) 戻り値なし。
"""
# 干渉波を計算
total_wave = interference_wave(X, Y, t, k, omega, source1, source2)
# 干渉波の強度（|ψ|^2）を計算
intensity = np.abs(total_wave)**2
# y = 5の位置での強度を取得
# y = 5 に対応するインデックスを取得
y_index = np.abs(y_range - args.contour_y).argmin()
# y = 5 の位置での強度を抽出し、2D配列として生成
Z_contour = np.zeros_like(intensity)
Z_contour[y_index, :] = intensity[y_index, :] # y = 5の位置だけ強度を設定
# y = 5 での等高線を描画
ax.clear()
ax.scatter(source1[0], source1[1], 0, color="red", s=100, label="Source 1") # 光源1
ax.scatter(source2[0], source2[1], 0, color="blue", s=100, label="Source 2") # 光源2
ax.plot_surface(X, Y, np.real(total_wave), cmap='Blues', edgecolor='none')
ax.contour(X, Y, Z_contour, 20, cmap='hot') # 等高線の追加
ax.set_xlim([args.xmin, args.xmax])
ax.set_ylim([args.ymin, args.ymax])
ax.set_zlim([args.zmin, args.zmax])
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Wave Amplitude')
ax.set_title(f'Wave Interference at time t={t:.2f}s')
ax.legend()
# プロット
ani = FuncAnimation(
fig,
animate,
frames=np.linspace(0, t_max, t_steps),
interval=args.interval,
)
plt.show()
if __name__ == "__main__":
main()