地表波伝搬 — 中波・長波の地表面回折メカニズム

import numpy as np
import matplotlib.pyplot as plt

def fresnel_reflection(psi_deg, eps_r, sigma, freq_hz):
    """
    フレネル反射係数を計算する
    psi_deg: グレージング角 [度]
    eps_r: 比誘電率
    sigma: 導電率 [S/m]
    freq_hz: 周波数 [Hz]
    返り値: (Gamma_V, Gamma_H) 複素反射係数
    """
    psi = np.radians(psi_deg)
    omega = 2 * np.pi * freq_hz
    eps0 = 8.854e-12

    # 複素比誘電率
    eps_c = eps_r - 1j * sigma / (omega * eps0)

    sin_psi = np.sin(psi)
    cos_psi = np.cos(psi)
    sqrt_term = np.sqrt(eps_c - cos_psi**2)

    # 垂直偏波（TM）
    Gamma_V = (eps_c * sin_psi - sqrt_term) / (eps_c * sin_psi + sqrt_term)

    # 水平偏波（TE）
    Gamma_H = (sin_psi - sqrt_term) / (sin_psi + sqrt_term)

    return Gamma_V, Gamma_H

# グレージング角
psi = np.linspace(0.1, 90, 500)

# 異なる地表面条件
surfaces = [
    ("Sea water", 80, 4.0, '#1565C0'),
    ("Wet soil", 25, 0.02, '#2E7D32'),
    ("Average ground", 15, 0.005, '#F57F17'),
    ("Dry soil", 7, 0.001, '#BF360C'),
    ("Desert sand", 4, 0.0003, '#6D4C41'),
]

freq = 1e6  # 1 MHz（中波帯）

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 垂直偏波 |Gamma_V|
for name, eps_r, sigma, color in surfaces:
    Gamma_V, Gamma_H = fresnel_reflection(psi, eps_r, sigma, freq)
    axes[0, 0].plot(psi, np.abs(Gamma_V), linewidth=2, color=color, label=name)
axes[0, 0].set_xlabel('Grazing Angle [deg]', fontsize=12)
axes[0, 0].set_ylabel('|$\\Gamma_V$|', fontsize=12)
axes[0, 0].set_title('Vertical Polarization Reflection (1 MHz)', fontsize=13)
axes[0, 0].legend(fontsize=9)
axes[0, 0].grid(True, alpha=0.3)
axes[0, 0].set_xlim([0, 90])
axes[0, 0].set_ylim([0, 1])

# 水平偏波 |Gamma_H|
for name, eps_r, sigma, color in surfaces:
    Gamma_V, Gamma_H = fresnel_reflection(psi, eps_r, sigma, freq)
    axes[0, 1].plot(psi, np.abs(Gamma_H), linewidth=2, color=color, label=name)
axes[0, 1].set_xlabel('Grazing Angle [deg]', fontsize=12)
axes[0, 1].set_ylabel('|$\\Gamma_H$|', fontsize=12)
axes[0, 1].set_title('Horizontal Polarization Reflection (1 MHz)', fontsize=13)
axes[0, 1].legend(fontsize=9)
axes[0, 1].grid(True, alpha=0.3)
axes[0, 1].set_xlim([0, 90])
axes[0, 1].set_ylim([0, 1])

# 垂直偏波の位相
for name, eps_r, sigma, color in surfaces:
    Gamma_V, _ = fresnel_reflection(psi, eps_r, sigma, freq)
    axes[1, 0].plot(psi, np.degrees(np.angle(Gamma_V)),
                     linewidth=2, color=color, label=name)
axes[1, 0].set_xlabel('Grazing Angle [deg]', fontsize=12)
axes[1, 0].set_ylabel('Phase of $\\Gamma_V$ [deg]', fontsize=12)
axes[1, 0].set_title('Phase of Vertical Reflection (1 MHz)', fontsize=13)
axes[1, 0].legend(fontsize=9)
axes[1, 0].grid(True, alpha=0.3)
axes[1, 0].set_xlim([0, 90])

# 周波数による変化（平均的な大地）
freqs_sweep = [0.1e6, 0.3e6, 1e6, 3e6, 10e6, 30e6]  # Hz
eps_r_avg, sigma_avg = 15, 0.005

for f in freqs_sweep:
    Gamma_V, _ = fresnel_reflection(psi, eps_r_avg, sigma_avg, f)
    label = f'{f/1e6:.1f} MHz' if f < 1e6 else f'{f/1e6:.0f} MHz'
    axes[1, 1].plot(psi, np.abs(Gamma_V), linewidth=2, label=label)
axes[1, 1].set_xlabel('Grazing Angle [deg]', fontsize=12)
axes[1, 1].set_ylabel('|$\\Gamma_V$|', fontsize=12)
axes[1, 1].set_title('Frequency Dependence of $\\Gamma_V$ (avg. ground)', fontsize=13)
axes[1, 1].legend(fontsize=9)
axes[1, 1].grid(True, alpha=0.3)
axes[1, 1].set_xlim([0, 90])
axes[1, 1].set_ylim([0, 1])

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

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import erfc

def ground_wave_field(d_km, freq_mhz, eps_r, sigma, Pt_kw=1.0):
    """
    地表波の電界強度を計算する（平面大地近似）
    d_km: 距離 [km]
    freq_mhz: 周波数 [MHz]
    eps_r: 大地の比誘電率
    sigma: 大地の導電率 [S/m]
    Pt_kw: 送信電力 [kW]
    返り値: 電界強度 [dB(uV/m)]
    """
    d = d_km * 1000  # m
    f = freq_mhz * 1e6  # Hz
    lam = 3e8 / f  # 波長 [m]
    k = 2 * np.pi / lam

    # 複素比誘電率
    eps0 = 8.854e-12
    omega = 2 * np.pi * f
    eps_c = eps_r - 1j * sigma / (omega * eps0)

    # 正規化表面インピーダンス
    Delta = 1.0 / np.sqrt(eps_c)

    # 数値距離
    w = -1j * np.pi * d / lam * (Delta / (Delta + 1))**2

    # 減衰関数の近似（ノートンの公式）
    sqrt_w = np.sqrt(w)
    # 大きい|w|では近似を使用
    F = np.where(np.abs(w) < 100,
                  1 - 1j * np.sqrt(np.pi) * sqrt_w * np.exp(-w) * erfc(1j * sqrt_w),
                  1.0 / (2 * w))

    # 自由空間の電界強度（垂直ダイポール、1 kW）
    # E = sqrt(30 * Pt * Gt) / d [V/m]
    # 垂直短ダイポールの利得 Gt ≈ 1.5 (1.76 dBi)
    Gt = 1.5
    E_free = np.sqrt(30 * Pt_kw * 1000 * Gt) / d  # V/m

    # 地表波の電界強度（垂直偏波なので係数2を考慮）
    E = E_free * np.abs(F) * 2  # 直接波+像の寄与

    # dB(uV/m)に変換
    E_dbuvm = 20 * np.log10(E * 1e6 + 1e-30)  # +1e-30でlog(0)回避

    return E_dbuvm, np.abs(F)

# 距離配列
d = np.logspace(-1, 3, 500)  # 0.1〜1000 km

# 周波数を変えた場合（平均的な大地）
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

freqs = [0.1, 0.3, 1.0, 3.0, 10.0]  # MHz
eps_r, sigma = 15, 0.005  # 平均的な大地

for f in freqs:
    E_dbuvm, F_abs = ground_wave_field(d, f, eps_r, sigma)
    axes[0].plot(d, E_dbuvm, linewidth=2, label=f'{f} MHz')

# 自由空間参考線
E_fs = 20 * np.log10(np.sqrt(30 * 1000 * 1.5) / (d * 1000) * 1e6)
axes[0].plot(d, E_fs, 'k--', linewidth=1, alpha=0.5, label='Free space')

axes[0].set_xlabel('Distance [km]', fontsize=12)
axes[0].set_ylabel('Field Strength [dB($\\mu$V/m)]', fontsize=12)
axes[0].set_title('Ground Wave: Frequency Dependence\n(Average Ground, 1 kW)', fontsize=13)
axes[0].legend(fontsize=10)
axes[0].grid(True, alpha=0.3)
axes[0].set_xscale('log')
axes[0].set_xlim([0.1, 1000])
axes[0].set_ylim([0, 140])

# 地表面の種類を変えた場合（1 MHz固定）
surfaces = [
    ("Sea water", 80, 4.0, '#1565C0'),
    ("Wet soil", 25, 0.02, '#2E7D32'),
    ("Average ground", 15, 0.005, '#F57F17'),
    ("Dry soil", 7, 0.001, '#BF360C'),
    ("Desert sand", 4, 0.0003, '#6D4C41'),
]

for name, eps_r, sigma, color in surfaces:
    E_dbuvm, _ = ground_wave_field(d, 1.0, eps_r, sigma)
    axes[1].plot(d, E_dbuvm, linewidth=2, color=color, label=name)

axes[1].set_xlabel('Distance [km]', fontsize=12)
axes[1].set_ylabel('Field Strength [dB($\\mu$V/m)]', fontsize=12)
axes[1].set_title('Ground Wave: Surface Dependence\n(1 MHz, 1 kW)', fontsize=13)
axes[1].legend(fontsize=10)
axes[1].grid(True, alpha=0.3)
axes[1].set_xscale('log')
axes[1].set_xlim([0.1, 1000])
axes[1].set_ylim([0, 140])

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

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import erfc

def attenuation_function(w_abs):
    """
    減衰関数 |F(w)| の計算（数値距離の絶対値から）
    簡易版: w を実数として近似
    """
    w = w_abs * np.exp(1j * np.pi / 4)  # 位相 pi/4 の近似
    sqrt_w = np.sqrt(w)
    F = np.where(np.abs(w) < 200,
                  1 - 1j * np.sqrt(np.pi) * sqrt_w * np.exp(-w) * erfc(1j * sqrt_w),
                  1.0 / (2 * w))
    return np.abs(F)

# 数値距離
w_abs = np.logspace(-3, 3, 1000)

F_abs = attenuation_function(w_abs)

fig, ax = plt.subplots(figsize=(10, 6))
ax.loglog(w_abs, F_abs, 'b-', linewidth=2.5, label='|F(w)|')

# 近似曲線
# 近距離: F ≈ 1
ax.axhline(y=1, color='gray', linestyle=':', alpha=0.5, label='F = 1 (free space)')
# 遠距離: F ≈ 1/(2|w|)
w_far = np.logspace(0.5, 3, 100)
ax.loglog(w_far, 1 / (2 * w_far), 'r--', linewidth=2,
           label='F $\\approx$ 1/(2|w|) (far field)')

ax.set_xlabel('Numerical Distance |w|', fontsize=13)
ax.set_ylabel('|F(w)|', fontsize=13)
ax.set_title('Ground Wave Attenuation Function', fontsize=14)
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3, which='both')
ax.set_xlim([1e-3, 1e3])
ax.set_ylim([1e-4, 2])

# 領域のアノテーション
ax.text(0.01, 1.2, 'Near field\n(free space)', fontsize=11,
         ha='center', color='green')
ax.text(3, 0.05, 'Transition', fontsize=11, ha='center', color='orange')
ax.text(100, 0.002, 'Far field\n($\\propto 1/d^2$)', fontsize=11,
         ha='center', color='red')

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

import numpy as np
import matplotlib.pyplot as plt

def two_ray_model(d_km, ht, hr, freq_mhz, eps_r, sigma):
    """
    2波モデルによる受信電力の相対値
    d_km: 距離 [km]
    ht, hr: 送受信アンテナ高さ [m]
    freq_mhz: 周波数 [MHz]
    eps_r, sigma: 大地パラメータ
    """
    d = d_km * 1000  # m
    f = freq_mhz * 1e6
    lam = 3e8 / f
    k = 2 * np.pi / lam

    # 直接経路長と反射経路長
    d1 = np.sqrt(d**2 + (ht - hr)**2)
    d2 = np.sqrt(d**2 + (ht + hr)**2)

    # グレージング角
    psi = np.arctan((ht + hr) / d)

    # フレネル反射係数（垂直偏波）
    eps0 = 8.854e-12
    omega = 2 * np.pi * f
    eps_c = eps_r - 1j * sigma / (omega * eps0)
    sin_psi = np.sin(psi)
    cos_psi = np.cos(psi)
    sqrt_term = np.sqrt(eps_c - cos_psi**2)
    Gamma_V = (eps_c * sin_psi - sqrt_term) / (eps_c * sin_psi + sqrt_term)

    # 合成電界（自由空間基準で正規化）
    E = np.abs(d[0] / d1 * np.exp(-1j * k * d1) +
               Gamma_V * d[0] / d2 * np.exp(-1j * k * d2))

    # dB変換（自由空間の d[0] 地点を基準）
    P_rel = 20 * np.log10(E + 1e-30) - 20 * np.log10(2)  # 補正

    return P_rel

d = np.logspace(-1, 2, 2000)  # 0.1〜100 km
freq = 100  # 100 MHz（VHF帯）
eps_r, sigma = 15, 0.005

fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# アンテナ高さを変えた場合
ht_values = [10, 30, 50, 100]
hr = 2  # 受信 2 m

for ht in ht_values:
    P = two_ray_model(d, ht, hr, freq, eps_r, sigma)
    axes[0].plot(d, P, linewidth=1.5, label=f'$h_t$ = {ht} m')

# 自由空間参考
P_fs = -20 * np.log10(d / d[0])
axes[0].plot(d, P_fs, 'k--', linewidth=1, alpha=0.5, label='Free space ($1/d^2$)')

# d^4 参考
P_d4 = -40 * np.log10(d / 1) + 10  # 適当なオフセット
axes[0].plot(d, P_d4, 'r:', linewidth=1, alpha=0.5, label='Two-ray ($1/d^4$)')

axes[0].set_xlabel('Distance [km]', fontsize=12)
axes[0].set_ylabel('Relative Power [dB]', fontsize=12)
axes[0].set_title(f'Two-Ray Model ({freq} MHz, $h_r$ = {hr} m)', fontsize=13)
axes[0].legend(fontsize=9)
axes[0].grid(True, alpha=0.3)
axes[0].set_xscale('log')
axes[0].set_xlim([0.1, 100])
axes[0].set_ylim([-80, 10])

# 高度パターン（距離固定、受信アンテナ高さを変化）
d_fixed = [1, 5, 10, 30]  # km
hr_array = np.linspace(0.5, 50, 500)

for d_val in d_fixed:
    d_arr = np.array([d_val])
    P_list = []
    for h in hr_array:
        P = two_ray_model(d_arr, 50, h, freq, eps_r, sigma)
        P_list.append(P[0])
    axes[1].plot(np.array(P_list), hr_array, linewidth=1.5,
                  label=f'd = {d_val} km')

axes[1].set_xlabel('Relative Power [dB]', fontsize=12)
axes[1].set_ylabel('Receiver Height [m]', fontsize=12)
axes[1].set_title(f'Height-Gain Pattern ({freq} MHz, $h_t$ = 50 m)', fontsize=13)
axes[1].legend(fontsize=10)
axes[1].grid(True, alpha=0.3)

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