import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# --- 自由空間損失の計算関数 ---
def fspl_db(d_km, f_mhz):
    """
    自由空間損失 [dB]

    Parameters
    ----------
    d_km : float or array
        距離 [km]
    f_mhz : float or array
        周波数 [MHz]

    Returns
    -------
    loss_db : float or array
        自由空間損失 [dB]
    """
    return 20 * np.log10(d_km) + 20 * np.log10(f_mhz) + 32.44

# --- 3Dプロット ---
d_km = np.logspace(-2, 5, 200)     # 10 m ～ 100,000 km
f_mhz = np.logspace(1, 5, 200)     # 10 MHz ～ 100 GHz

D, F = np.meshgrid(d_km, f_mhz)
L = fspl_db(D, F)

fig = plt.figure(figsize=(14, 6))

# 3Dサーフェス
ax1 = fig.add_subplot(121, projection='3d')
surf = ax1.plot_surface(np.log10(D), np.log10(F), L,
                        cmap='viridis', alpha=0.85, edgecolor='none')

ax1.set_xlabel('log$_{10}$(d [km])', fontsize=11)
ax1.set_ylabel('log$_{10}$(f [MHz])', fontsize=11)
ax1.set_zlabel('FSPL [dB]', fontsize=11)
ax1.set_title('Free Space Path Loss', fontsize=13)
fig.colorbar(surf, ax=ax1, shrink=0.5, aspect=10, label='FSPL [dB]')
ax1.view_init(elev=25, azim=-130)

# 2D等高線
ax2 = fig.add_subplot(122)
levels = np.arange(20, 300, 20)
cs = ax2.contourf(np.log10(D), np.log10(F), L,
                  levels=levels, cmap='viridis')
ax2.contour(np.log10(D), np.log10(F), L,
            levels=levels, colors='white', linewidths=0.5, alpha=0.5)
fig.colorbar(cs, ax=ax2, label='FSPL [dB]')

# 代表的なシナリオをプロット
scenarios = {
    'Wi-Fi': (-2, np.log10(5000)),
    'Cellular': (0, np.log10(2000)),
    'GEO Sat': (np.log10(36000), np.log10(12000)),
}
for name, (logd, logf) in scenarios.items():
    loss = fspl_db(10**logd, 10**logf)
    ax2.plot(logd, logf, 'r*', markersize=15)
    ax2.annotate(f'{name}\n{loss:.0f} dB', xy=(logd, logf),
                 xytext=(logd + 0.3, logf + 0.2),
                 fontsize=9, color='red',
                 arrowprops=dict(arrowstyle='->', color='red'))

ax2.set_xlabel('log$_{10}$(d [km])', fontsize=11)
ax2.set_ylabel('log$_{10}$(f [MHz])', fontsize=11)
ax2.set_title('FSPL Contour Map', fontsize=13)

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

import numpy as np
import matplotlib.pyplot as plt

def fspl_db(d_km, f_mhz):
    """自由空間損失 [dB]"""
    return 20 * np.log10(d_km) + 20 * np.log10(f_mhz) + 32.44

# --- 周波数ごとの距離-損失カーブ ---
d_km = np.logspace(-3, 6, 500)  # 1 m ～ 1,000,000 km

frequencies = {
    '900 MHz (LTE)': 900,
    '2.4 GHz (Wi-Fi)': 2400,
    '5 GHz (Wi-Fi)': 5000,
    '12 GHz (Ku-band)': 12000,
    '28 GHz (5G mmWave)': 28000,
    '60 GHz (WiGig)': 60000,
}

fig, ax = plt.subplots(figsize=(12, 7))

colors = plt.cm.tab10(np.linspace(0, 1, len(frequencies)))
for (name, f_mhz), color in zip(frequencies.items(), colors):
    L = fspl_db(d_km, f_mhz)
    ax.plot(d_km, L, color=color, linewidth=2, label=name)

# 距離マーカー
distances = {
    '10 m': 0.01,
    '100 m': 0.1,
    '1 km': 1,
    '100 km': 100,
    'GEO': 36000,
    'Moon': 384400,
}
for name, d in distances.items():
    ax.axvline(x=d, color='gray', linestyle=':', linewidth=0.8, alpha=0.5)
    ax.text(d, 30, name, rotation=90, fontsize=8, color='gray',
            ha='right', va='bottom')

ax.set_xlabel('Distance [km]', fontsize=12)
ax.set_ylabel('FSPL [dB]', fontsize=12)
ax.set_title('Free Space Path Loss vs Distance', fontsize=14)
ax.set_xscale('log')
ax.legend(fontsize=9, loc='lower right')
ax.grid(True, alpha=0.3, which='both')
ax.set_ylim([20, 300])
ax.set_xlim([1e-3, 1e6])

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

import numpy as np

# --- リンクバジェット計算 ---

def link_budget(P_t_dBW, G_t_dBi, G_r_dBi, d_km, f_mhz,
                L_other_dB=0, T_sys_K=290, bandwidth_hz=1e6):
    """
    リンクバジェット計算

    Parameters
    ----------
    P_t_dBW : float
        送信電力 [dBW]
    G_t_dBi : float
        送信アンテナ利得 [dBi]
    G_r_dBi : float
        受信アンテナ利得 [dBi]
    d_km : float
        距離 [km]
    f_mhz : float
        周波数 [MHz]
    L_other_dB : float
        その他の損失 [dB]
    T_sys_K : float
        システム雑音温度 [K]
    bandwidth_hz : float
        帯域幅 [Hz]

    Returns
    -------
    results : dict
        計算結果の辞書
    """
    # 自由空間損失
    L_fs_dB = 20 * np.log10(d_km) + 20 * np.log10(f_mhz) + 32.44

    # EIRP
    EIRP_dBW = P_t_dBW + G_t_dBi

    # 受信電力
    P_r_dBW = EIRP_dBW + G_r_dBi - L_fs_dB - L_other_dB

    # 雑音電力
    k_B = 1.38e-23  # ボルツマン定数 [J/K]
    k_B_dBW = 10 * np.log10(k_B)  # -228.6 dBW/K/Hz
    N_dBW = k_B_dBW + 10 * np.log10(T_sys_K) + 10 * np.log10(bandwidth_hz)

    # C/N
    CN_dB = P_r_dBW - N_dBW

    # C/N0
    CN0_dBHz = P_r_dBW - k_B_dBW - 10 * np.log10(T_sys_K)

    return {
        'EIRP [dBW]': EIRP_dBW,
        'FSPL [dB]': L_fs_dB,
        'Other Loss [dB]': L_other_dB,
        'Rx Power [dBW]': P_r_dBW,
        'Rx Power [dBm]': P_r_dBW + 30,
        'Noise Power [dBW]': N_dBW,
        'C/N [dB]': CN_dB,
        'C/N0 [dBHz]': CN0_dBHz,
    }


# === シナリオ1: 静止衛星ダウンリンク ===
print("=" * 55)
print("シナリオ1: 静止衛星ダウンリンク (Ku-band)")
print("=" * 55)
result1 = link_budget(
    P_t_dBW=10 * np.log10(100),   # 100 W = 20 dBW
    G_t_dBi=35,                    # 衛星アンテナ利得
    G_r_dBi=40,                    # 地上パラボラアンテナ（直径1.2m）
    d_km=36000,                    # 静止軌道
    f_mhz=12000,                   # 12 GHz
    L_other_dB=2,                  # 大気損失・指向誤差等
    T_sys_K=150,                   # 受信システム雑音温度
    bandwidth_hz=36e6              # 36 MHz帯域
)
for key, val in result1.items():
    print(f"  {key:>20s} = {val:>10.2f}")

# === シナリオ2: Wi-Fi (5 GHz) ===
print("\n" + "=" * 55)
print("シナリオ2: Wi-Fi (5 GHz, 屋内10m)")
print("=" * 55)
result2 = link_budget(
    P_t_dBW=10 * np.log10(0.1),   # 100 mW = -10 dBW
    G_t_dBi=3,                     # Wi-Fiアンテナ
    G_r_dBi=3,                     # 端末アンテナ
    d_km=0.01,                     # 10 m
    f_mhz=5000,                    # 5 GHz
    L_other_dB=5,                  # 壁・家具等の損失
    T_sys_K=600,                   # 室内雑音環境
    bandwidth_hz=80e6              # 80 MHz帯域
)
for key, val in result2.items():
    print(f"  {key:>20s} = {val:>10.2f}")

# === シナリオ3: 携帯基地局 ===
print("\n" + "=" * 55)
print("シナリオ3: 携帯基地局 (2 GHz, 1 km)")
print("=" * 55)
result3 = link_budget(
    P_t_dBW=10 * np.log10(20),    # 20 W = 13 dBW
    G_t_dBi=18,                    # セクタアンテナ
    G_r_dBi=0,                     # スマートフォン
    d_km=1,                        # 1 km
    f_mhz=2000,                    # 2 GHz
    L_other_dB=10,                 # 建物遮蔽・フェージング余裕
    T_sys_K=290,                   # 室温
    bandwidth_hz=10e6              # 10 MHz帯域
)
for key, val in result3.items():
    print(f"  {key:>20s} = {val:>10.2f}")

import numpy as np
import matplotlib.pyplot as plt

# --- リンクバジェットのウォーターフォールチャート ---

# 静止衛星ダウンリンクのリンクバジェット
items = [
    ('Tx Power\n(100W)', 20.0, 'gain'),
    ('Tx Antenna\nGain', 35.0, 'gain'),
    ('FSPL\n(36000km, 12GHz)', -205.15, 'loss'),
    ('Atm. Loss', -2.0, 'loss'),
    ('Rx Antenna\nGain', 40.0, 'gain'),
]

fig, ax = plt.subplots(figsize=(12, 7))

cumulative = 0
x_positions = []
heights = []
bottoms = []
colors_bar = []

for i, (label, value, type_) in enumerate(items):
    if type_ == 'gain':
        bottom = cumulative
        color = '#2196F3'
    else:
        bottom = cumulative + value
        color = '#F44336'

    x_positions.append(i)
    heights.append(abs(value))
    bottoms.append(min(cumulative, cumulative + value))
    colors_bar.append(color)
    cumulative += value

# 受信電力（最終結果）
items.append(('Rx Power', cumulative, 'result'))
x_positions.append(len(items) - 1)
heights.append(abs(cumulative))
bottoms.append(min(0, cumulative))
colors_bar.append('#4CAF50')

labels = [item[0] for item in items]

ax.bar(x_positions, heights, bottom=bottoms, color=colors_bar,
       width=0.6, edgecolor='black', linewidth=0.8)

# 値のラベル
running = 0
for i, (label, value, type_) in enumerate(items[:-1]):
    running += value
    # バーの上に値を表示
    y_pos = bottoms[i] + heights[i] + 2
    ax.text(i, y_pos, f'{value:+.1f} dB', ha='center', fontsize=9,
            fontweight='bold')
    # 累積値を線で接続
    if i < len(items) - 2:
        ax.plot([i + 0.3, i + 0.7], [running, running], 'k--',
                linewidth=0.8, alpha=0.5)

# 最終結果
ax.text(len(items)-1, bottoms[-1] + heights[-1] + 2,
        f'{cumulative:.1f} dBW\n({cumulative+30:.1f} dBm)',
        ha='center', fontsize=10, fontweight='bold', color='green')

ax.set_xticks(x_positions)
ax.set_xticklabels(labels, fontsize=9)
ax.set_ylabel('Power Level [dBW]', fontsize=12)
ax.set_title('Link Budget Waterfall: GEO Satellite Downlink (Ku-band)', fontsize=14)
ax.grid(True, alpha=0.3, axis='y')
ax.axhline(y=0, color='gray', linewidth=0.8)

# 凡例
from matplotlib.patches import Patch
legend_elements = [
    Patch(facecolor='#2196F3', label='Gain'),
    Patch(facecolor='#F44336', label='Loss'),
    Patch(facecolor='#4CAF50', label='Result'),
]
ax.legend(handles=legend_elements, fontsize=11, loc='upper right')

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

import numpy as np
import matplotlib.pyplot as plt

# --- 電力密度の距離減衰 ---
d = np.logspace(0, 6, 500)  # 1 m ～ 1000 km

# 各シナリオのEIRP [W]
scenarios = {
    'Wi-Fi AP (100mW, 3dBi)': 0.1 * 10**(3/10),
    'Cell Tower (20W, 18dBi)': 20 * 10**(18/10),
    'Radar (1MW, 40dBi)': 1e6 * 10**(40/10),
    'GEO Sat (100W, 35dBi)': 100 * 10**(35/10),
}

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

for name, eirp in scenarios.items():
    S = eirp / (4 * np.pi * d**2)  # 電力密度 [W/m²]
    ax1.plot(d / 1e3, S, linewidth=2, label=name)

ax1.set_xlabel('Distance [km]', fontsize=12)
ax1.set_ylabel('Power Density [W/m$^2$]', fontsize=12)
ax1.set_title('Power Density vs Distance', fontsize=13)
ax1.set_xscale('log')
ax1.set_yscale('log')
ax1.legend(fontsize=9)
ax1.grid(True, alpha=0.3, which='both')
ax1.set_ylim([1e-20, 1e5])

# 逆2乗則の確認
d_norm = np.logspace(0, 4, 100)
S_norm = 1 / d_norm**2

ax2.plot(d_norm, S_norm, 'b-', linewidth=2, label='$S \\propto 1/d^2$')
ax2.plot(d_norm, 1 / d_norm**3, 'r--', linewidth=1.5, label='$\\propto 1/d^3$ (reference)')
ax2.plot(d_norm, 1 / d_norm, 'g--', linewidth=1.5, label='$\\propto 1/d$ (reference)')

# 距離2倍→電力1/4のマーカー
d1, d2 = 10, 20
ax2.plot([d1, d1], [1/d1**2, 1e-10], 'k:', linewidth=1)
ax2.plot([d2, d2], [1/d2**2, 1e-10], 'k:', linewidth=1)
ax2.annotate(f'd={d1}: S={1/d1**2:.4f}', xy=(d1, 1/d1**2),
             xytext=(d1*2, 1/d1**2*3), fontsize=9,
             arrowprops=dict(arrowstyle='->', color='black'))
ax2.annotate(f'd={d2}: S={1/d2**2:.5f}\n(= 1/4 of above)', xy=(d2, 1/d2**2),
             xytext=(d2*2, 1/d2**2*5), fontsize=9,
             arrowprops=dict(arrowstyle='->', color='black'))

ax2.set_xlabel('Distance (normalized)', fontsize=12)
ax2.set_ylabel('Power density (normalized)', fontsize=12)
ax2.set_title('Inverse Square Law', fontsize=13)
ax2.set_xscale('log')
ax2.set_yscale('log')
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3, which='both')

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