反射器  放射器  導波器1  導波器2  導波器3  ...  → 最大放射方向
  |       |       |        |        |
  ←  d_r →← d_1 →← d_2  →← d_3  →

import numpy as np
import matplotlib.pyplot as plt

def array_factor(N, d_lambda, phi, alpha=None):
    """
    N素子等間隔等振幅アレイのアレイファクターを計算する。

    Parameters
    ----------
    N : int
        素子数
    d_lambda : float
        素子間隔（波長比 d/λ）
    phi : array
        観測角 [rad]（0がエンドファイア方向）
    alpha : float or None
        隣接素子間の位相差 [rad]。Noneならエンドファイア条件を適用

    Returns
    -------
    AF_norm : array
        正規化電力パターン（最大値1）
    """
    k = 2 * np.pi  # k*λ = 2π なので k*d = 2π*(d/λ)
    kd = k * d_lambda

    if alpha is None:
        # エンドファイア条件: α = -kd (Hansen-Woodyard改良なし)
        alpha = -kd

    psi = kd * np.cos(phi) + alpha

    # sin(N*psi/2) / sin(psi/2) の計算（ゼロ除算回避）
    numerator = np.sin(N * psi / 2)
    denominator = np.sin(psi / 2)

    # |psi/2| が非常に小さい場合の処理
    AF = np.where(
        np.abs(denominator) < 1e-10,
        N,  # ロピタルの定理から極限値はN
        numerator / denominator
    )

    AF_power = np.abs(AF)**2
    AF_norm = AF_power / np.max(AF_power)
    return AF_norm

# 角度の定義
phi = np.linspace(0, 2 * np.pi, 3601)

# 素子間隔
d_lambda = 0.25  # λ/4

# 素子数を変えてパターンを比較
N_values = [2, 3, 5, 7, 10, 15]

fig, axes = plt.subplots(2, 3, figsize=(18, 12),
                         subplot_kw=dict(projection='polar'))

for idx, N in enumerate(N_values):
    ax = axes[idx // 3][idx % 3]
    AF = array_factor(N, d_lambda, phi)

    # dBスケール
    AF_dB = np.where(AF > 1e-10, 10 * np.log10(AF), -40)
    AF_dB_plot = np.clip(AF_dB, -40, 0) + 40  # 0〜40の範囲に変換

    ax.plot(phi, AF_dB_plot, 'b-', linewidth=1.5)
    ax.set_title(f'N = {N} elements\n$d = {d_lambda}\\lambda$', pad=15)
    ax.set_theta_zero_location('E')  # 0°を右に

plt.suptitle('Endfire Array Factor (dB, floor=-40dB)', fontsize=14, y=1.02)
plt.tight_layout()
plt.savefig('yagi_array_patterns.png', dpi=150, bbox_inches='tight')
plt.show()

import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate

def compute_endfire_directivity(N, d_lambda):
    """
    N素子エンドファイアアレイの指向性を数値計算する。
    素子パターンは等方性（アレイファクターのみ考慮）。

    Parameters
    ----------
    N : int
        素子数
    d_lambda : float
        素子間隔（波長比）

    Returns
    -------
    D : float
        指向性（リニア値）
    """
    kd = 2 * np.pi * d_lambda
    alpha = -kd  # エンドファイア条件

    def AF_power(theta):
        """球面座標theta方向のアレイファクター電力パターン"""
        psi = kd * np.cos(theta) + alpha
        denom = np.sin(psi / 2)
        if np.abs(denom) < 1e-12:
            return float(N**2)
        return (np.sin(N * psi / 2) / denom)**2

    # 放射パターンは軸対称と仮定（φ依存なし）
    # 最大放射強度 (theta=0)
    U_max = AF_power(0.0)

    # 全放射電力の計算: ∫ AF^2(theta) sin(theta) dtheta × 2π
    def integrand(theta):
        return AF_power(theta) * np.sin(theta)

    P_total, _ = integrate.quad(integrand, 0, np.pi, limit=200)
    P_total *= 2 * np.pi

    # 指向性
    D = 4 * np.pi * U_max / P_total
    return D

# 素子数を変えて指向性を計算
d_lambda = 0.25
N_list = np.arange(2, 21)
D_list = []

for N in N_list:
    D = compute_endfire_directivity(N, d_lambda)
    D_list.append(D)

D_array = np.array(D_list)
G_dBi = 10 * np.log10(D_array)

# プロット
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# リニアスケール
ax = axes[0]
ax.plot(N_list, D_array, 'bo-', markersize=5)
ax.set_xlabel('Number of elements $N$', fontsize=12)
ax.set_ylabel('Directivity $D$', fontsize=12)
ax.set_title('Directivity vs Number of Elements (Linear)', fontsize=12)
ax.grid(True, alpha=0.3)

# dBスケール
ax = axes[1]
ax.plot(N_list, G_dBi, 'ro-', markersize=5)
ax.set_xlabel('Number of elements $N$', fontsize=12)
ax.set_ylabel('Directivity [dBi]', fontsize=12)
ax.set_title('Directivity vs Number of Elements (dB)', fontsize=12)
ax.grid(True, alpha=0.3)

# 表形式で結果を表示
print(f"{'N':>4} {'D':>10} {'D [dBi]':>10}")
print("-" * 30)
for N, D, GdB in zip(N_list, D_array, G_dBi):
    print(f"{N:4d} {D:10.2f} {GdB:10.2f}")

plt.tight_layout()
plt.savefig('yagi_directivity_vs_elements.png', dpi=150, bbox_inches='tight')
plt.show()

import numpy as np
import matplotlib.pyplot as plt

def yagi_pattern(phi, reflector_pos, driven_pos, director_pos_list,
                 reflector_phase, director_phases, reflector_amp=0.9,
                 director_amps=None):
    """
    簡易八木・宇田アンテナの放射パターンを計算する。

    Parameters
    ----------
    phi : array
        観測角 [rad]
    reflector_pos : float
        反射器のx座標（波長単位）
    driven_pos : float
        放射器のx座標（波長単位）
    director_pos_list : list of float
        各導波器のx座標（波長単位）
    reflector_phase : float
        反射器の電流位相 [rad]
    director_phases : list of float
        各導波器の電流位相 [rad]
    reflector_amp : float
        反射器の電流振幅比
    director_amps : list of float or None
        各導波器の電流振幅比

    Returns
    -------
    F_norm : array
        正規化電力パターン
    """
    k = 2 * np.pi  # k = 2π/λ, 位置は波長単位

    if director_amps is None:
        # 導波器は遠いほど振幅が減衰すると仮定
        director_amps = [0.8 * (0.9)**i for i in range(len(director_pos_list))]

    # 複素電界の重ね合わせ
    E = np.zeros(len(phi), dtype=complex)

    # 放射器（位相基準、振幅1）
    E += 1.0 * np.exp(1j * k * driven_pos * np.cos(phi))

    # 反射器
    E += reflector_amp * np.exp(1j * reflector_phase) * \
         np.exp(1j * k * reflector_pos * np.cos(phi))

    # 導波器
    for pos, phase, amp in zip(director_pos_list, director_phases,
                                director_amps):
        E += amp * np.exp(1j * phase) * np.exp(1j * k * pos * np.cos(phi))

    F_power = np.abs(E)**2
    F_norm = F_power / np.max(F_power)
    return F_norm

# 角度
phi = np.linspace(0, 2 * np.pi, 3601)

# --- 設計1: 3素子（反射器+放射器+導波器1本）---
driven_pos = 0.0
reflector_pos = -0.20  # 反射器は後方0.20λ
director_pos_1 = [0.25]  # 導波器は前方0.25λ
reflector_phase_1 = 0.5  # 遅れ位相
director_phases_1 = [-0.5]  # 進み位相

F_3elem = yagi_pattern(phi, reflector_pos, driven_pos, director_pos_1,
                       reflector_phase_1, director_phases_1)

# --- 設計2: 5素子 ---
director_pos_5 = [0.25, 0.55, 0.85]
director_phases_5 = [-0.5, -1.0, -1.5]
F_5elem = yagi_pattern(phi, reflector_pos, driven_pos, director_pos_5,
                       reflector_phase_1, director_phases_5)

# --- 設計3: 10素子 ---
director_pos_10 = [0.25 + 0.30 * i for i in range(8)]
director_phases_10 = [-0.5 - 0.5 * i for i in range(8)]
F_10elem = yagi_pattern(phi, reflector_pos, driven_pos, director_pos_10,
                        reflector_phase_1, director_phases_10)

# プロット
fig, axes = plt.subplots(1, 3, figsize=(18, 6),
                         subplot_kw=dict(projection='polar'))

configs = [
    (F_3elem, '3 elements\n(Refl + Driven + 1 Dir)'),
    (F_5elem, '5 elements\n(Refl + Driven + 3 Dir)'),
    (F_10elem, '10 elements\n(Refl + Driven + 8 Dir)')
]

for ax, (F, title) in zip(axes, configs):
    # dBスケール
    F_dB = np.where(F > 1e-10, 10 * np.log10(F), -40)
    F_plot = np.clip(F_dB, -40, 0) + 40

    ax.plot(phi, F_plot, 'b-', linewidth=1.5)
    ax.set_title(title, pad=15, fontsize=11)
    ax.set_theta_zero_location('E')

    # 前方/後方比(F/B ratio)の計算
    idx_front = np.argmin(np.abs(phi - 0))
    idx_back = np.argmin(np.abs(phi - np.pi))
    fb_ratio = 10 * np.log10(F[idx_front] / max(F[idx_back], 1e-10))
    ax.annotate(f'F/B = {fb_ratio:.1f} dB', xy=(0.05, 0.95),
                xycoords='axes fraction', fontsize=9,
                bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))

plt.suptitle('Yagi-Uda Antenna Radiation Patterns (dB)', fontsize=14)
plt.tight_layout()
plt.savefig('yagi_uda_patterns.png', dpi=150, bbox_inches='tight')
plt.show()

import numpy as np
import matplotlib.pyplot as plt

def compute_fb_ratio(reflector_phase, director_phase, d_r=0.20, d_d=0.25):
    """
    3素子八木・宇田アンテナの前後比(F/B ratio)を計算する。

    Parameters
    ----------
    reflector_phase : float
        反射器の電流位相 [rad]
    director_phase : float
        導波器の電流位相 [rad]
    d_r : float
        反射器〜放射器の間隔（波長単位）
    d_d : float
        放射器〜導波器の間隔（波長単位）

    Returns
    -------
    fb_dB : float
        前後比 [dB]
    """
    k = 2 * np.pi
    # 前方 (phi=0) の電界
    E_front = (1.0 +
               0.9 * np.exp(1j * reflector_phase) * np.exp(-1j * k * d_r) +
               0.8 * np.exp(1j * director_phase) * np.exp(1j * k * d_d))

    # 後方 (phi=π) の電界
    E_back = (1.0 +
              0.9 * np.exp(1j * reflector_phase) * np.exp(1j * k * d_r) +
              0.8 * np.exp(1j * director_phase) * np.exp(-1j * k * d_d))

    P_front = np.abs(E_front)**2
    P_back = np.abs(E_back)**2

    if P_back < 1e-15:
        return 40.0
    return 10 * np.log10(P_front / P_back)

# 反射器位相と導波器位相を走査
refl_phases = np.linspace(-1.0, 2.0, 200)
dir_phases = np.linspace(-2.0, 1.0, 200)
REFL, DIR = np.meshgrid(refl_phases, dir_phases)

FB = np.zeros_like(REFL)
for i in range(len(dir_phases)):
    for j in range(len(refl_phases)):
        FB[i, j] = compute_fb_ratio(REFL[i, j], DIR[i, j])

# プロット
fig, ax = plt.subplots(figsize=(10, 7))
levels = np.arange(-10, 25, 2)
cs = ax.contourf(np.degrees(REFL) * 180 / np.pi / 180,
                 np.degrees(DIR) * 180 / np.pi / 180,
                 FB, levels=levels, cmap='RdYlGn')
# 角度をラジアン単位で表示
cs = ax.contourf(REFL, DIR, FB, levels=levels, cmap='RdYlGn')
plt.colorbar(cs, ax=ax, label='F/B Ratio [dB]')
ax.set_xlabel('Reflector current phase [rad]', fontsize=12)
ax.set_ylabel('Director current phase [rad]', fontsize=12)
ax.set_title('Front-to-Back Ratio vs Element Phases\n'
             '(3-element Yagi-Uda, $d_r=0.20\\lambda$, $d_d=0.25\\lambda$)',
             fontsize=12)
ax.grid(True, alpha=0.3)

# 最適点をマーク
opt_idx = np.unravel_index(np.argmax(FB), FB.shape)
ax.plot(REFL[opt_idx], DIR[opt_idx], 'k*', markersize=15,
        label=f'Optimum: F/B = {FB[opt_idx]:.1f} dB')
ax.legend(fontsize=10)

plt.tight_layout()
plt.savefig('yagi_fb_ratio_sensitivity.png', dpi=150, bbox_inches='tight')
plt.show()

print(f"最適反射器位相: {REFL[opt_idx]:.2f} rad ({np.degrees(REFL[opt_idx]):.1f}°)")
print(f"最適導波器位相: {DIR[opt_idx]:.2f} rad ({np.degrees(DIR[opt_idx]):.1f}°)")
print(f"最大F/B比: {FB[opt_idx]:.1f} dB")