import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

# --- 設定パラメータ ---
N_SOURCES = 16
WAVELENGTH = 1.0
K = 2 * np.pi / WAVELENGTH
SPACING = WAVELENGTH / 2
STEER_ANGLE = 30  # 30度方向にビームを向ける
SPEED_OF_LIGHT = 1.0 # シミュレーション上の波の速さ(c)

# 計算領域
X_RANGE = np.linspace(-10, 10, 300)
Y_RANGE = np.linspace(0, 15, 300)
X, Y = np.meshgrid(X_RANGE, Y_RANGE)

# 波源の配置
source_xs = np.linspace(-(N_SOURCES-1)*SPACING/2, (N_SOURCES-1)*SPACING/2, N_SOURCES)

def get_initial_phase(i):
    """ビーム走査のための各波源の初期位相オフセット"""
    theta_rad = np.radians(STEER_ANGLE)
    return -K * source_xs[i] * np.sin(theta_rad)

# プロットの初期化
fig, ax = plt.subplots(figsize=(8, 6))
# 最初のフレームを描画
initial_field = np.zeros_like(X)
cont = ax.contourf(X, Y, initial_field, levels=50, cmap='RdBu_r', vmin=-1.5, vmax=1.5)
ax.scatter(source_xs, np.zeros_like(source_xs), color='black', s=5, zorder=10)
ax.set_title(f"Phased Array Beamforming (Angle: {STEER_ANGLE}°)")
ax.set_aspect('equal')

def update(frame):
    """各フレームでの波の状態を計算"""
    ax.clear() # 前のフレームをクリア
    t = frame * 0.1 # 時間の進み
    omega = K * SPEED_OF_LIGHT
    
    total_field = np.zeros_like(X, dtype=complex)
    for i, sx in enumerate(source_xs):
        r = np.sqrt((X - sx)**2 + Y**2)
        phi_0 = get_initial_phase(i)
        # 複素形式: exp(j * (ωt - kr + φ0))
        total_field += np.exp(1j * (omega * t - K * r + phi_0)) / np.sqrt(r + 0.5)
    
    # 再描画
    field_real = np.real(total_field)
    im = ax.contourf(X, Y, field_real, levels=50, cmap='RdBu_r', vmin=-1.5, vmax=1.5)
    ax.scatter(source_xs, np.zeros_like(source_xs), color='black', s=5)
    ax.set_title(f"Phased Array Beamforming (Angle: {STEER_ANGLE}°)")
    ax.set_aspect('equal')
    return im,

# アニメーションの作成 (framesでループ回数を指定)
ani = FuncAnimation(fig, update, frames=50, interval=50, blit=False)

plt.show()