地球の上空約60～1000 kmには、太陽からの紫外線やX線によって大気が電離し、自由電子とイオンが存在する領域があります。この領域を 電離層（ionosphere） と呼びます。
電離層は電波の伝搬に大きな影響を与えます。短波（HF: 3～30 MHz）帯の電波は電離層で屈折・反射され、地球の裏側まで到達できます。これが短波通信や国際放送の原理です。一方、衛星通信で使われるマイクロ波帯（GHz帯）の電波は電離層を透過しますが、ファラデー回転や群遅延などの影響を受けます。
本記事の内容
- 電離層の構成（D/E/F層）と各層の特性
- プラズマ周波数の導出
- 屈折率の導出（電子の運動方程式→分極→誘電率）
- 臨界周波数と最大使用周波数（MUF）
- スネルの法則による電波屈折経路
- 衛星通信への電離層影響（ファラデー回転、群遅延）
- Pythonでのレイトレーシングシミュレーション
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
電離層の構成
電離層の層構造
電離層は電子密度の高さプロファイルに基づいて、いくつかの層に分類されます。
| 層 | 高度 [km] | 電子密度 [m$^{-3}$] | 特徴 |
|---|---|---|---|
| D層 | 60～90 | $10^8$ ～ $10^9$ | 昼間のみ存在。HF帯の吸収が大きい |
| E層 | 90～150 | $\sim 10^{11}$ | 比較的安定。スポラディックE（Es）による異常伝搬 |
| F1層 | 150～250 | $\sim 2 \times 10^{11}$ | 昼間のみ。夜間はF2層に統合 |
| F2層 | 250～500 | $\sim 10^{12}$ | 最も電子密度が高い。HF通信の主要反射層 |
電子密度プロファイルの日変化
電離層の電子密度は太陽活動によって大きく変動します。
- 昼間: 太陽光による電離が活発。D層・F1層が形成され、全体の電子密度が増加
- 夜間: 電離源がなくなり、D層とF1層は消滅。F2層は残存（再結合が遅い）
- 太陽活動極大期: 電子密度が増加し、臨界周波数が上昇
- 季節変動: 夏期と冬期で電子密度分布が変化
チャップマン関数
単純化された電離層のモデルとして、チャップマン関数 で電子密度プロファイルを表すことがあります。
$$ N_e(h) = N_0 \exp\left[\frac{1}{2}\left(1 – z – \sec\chi \cdot e^{-z}\right)\right] $$
ここで $z = (h – h_0)/H$ は正規化高度（$h_0$ は最大電子密度の高度、$H$ はスケールハイト）、$\chi$ は太陽天頂角です。
プラズマ周波数の導出
電子の運動方程式
電離層中の自由電子が電磁波の電場 $\bm{E} = E_0 e^{j\omega t}\hat{\bm{x}}$ によって駆動される場合を考えます。磁場の影響と衝突を無視すると、電子の運動方程式は、
$$ m_e \frac{d^2 \bm{x}}{dt^2} = -e\bm{E} $$
ここで $m_e = 9.109 \times 10^{-31}$ kg は電子質量、$e = 1.602 \times 10^{-19}$ C は電気素量です。
変位の解
電場が $E_0 e^{j\omega t}$ で振動するとき、電子の変位 $\bm{x}$ も同じ周波数で振動すると仮定します。
$$ \bm{x}(t) = \bm{x}_0 e^{j\omega t} $$
運動方程式に代入すると、
$$ m_e (j\omega)^2 \bm{x}_0 e^{j\omega t} = -e E_0 e^{j\omega t} \hat{\bm{x}} $$
$$ -m_e \omega^2 \bm{x}_0 = -e E_0 \hat{\bm{x}} $$
$$ \begin{equation} \bm{x}_0 = \frac{eE_0}{m_e \omega^2}\hat{\bm{x}} \end{equation} $$
分極ベクトル
単位体積あたり $N_e$ 個の電子が変位 $\bm{x}_0$ すると、分極ベクトル（電気双極子モーメント密度）は、
$$ \bm{P} = -N_e e \bm{x}_0 = -\frac{N_e e^2}{m_e \omega^2} E_0 e^{j\omega t}\hat{\bm{x}} $$
比誘電率の導出
電束密度 $\bm{D}$ と電場 $\bm{E}$ の関係は、
$$ \bm{D} = \varepsilon_0 \bm{E} + \bm{P} = \varepsilon_0 \bm{E} – \frac{N_e e^2}{m_e \omega^2}\bm{E} $$
$$ \bm{D} = \varepsilon_0\left(1 – \frac{N_e e^2}{\varepsilon_0 m_e \omega^2}\right)\bm{E} $$
したがって、比誘電率は、
$$ \varepsilon_r = 1 – \frac{N_e e^2}{\varepsilon_0 m_e \omega^2} $$
プラズマ周波数の定義
プラズマ角周波数 $\omega_p$ を以下で定義します。
$$ \omega_p^2 = \frac{N_e e^2}{\varepsilon_0 m_e} $$
プラズマ周波数 $f_p$ は、
$$ \begin{equation} \boxed{f_p = \frac{1}{2\pi}\sqrt{\frac{N_e e^2}{\varepsilon_0 m_e}} \approx 9\sqrt{N_e} \text{ [Hz]}} \end{equation} $$
ここで $N_e$ の単位は [m$^{-3}$] です。近似値の係数を確認します。
$$ \frac{1}{2\pi}\sqrt{\frac{e^2}{\varepsilon_0 m_e}} = \frac{1}{2\pi}\sqrt{\frac{(1.602 \times 10^{-19})^2}{8.854 \times 10^{-12} \times 9.109 \times 10^{-31}}} $$
$$ = \frac{1}{2\pi}\sqrt{\frac{2.566 \times 10^{-38}}{8.072 \times 10^{-42}}} = \frac{1}{2\pi}\sqrt{3178.6} = \frac{56.38}{2\pi} \approx 8.98 \approx 9 $$
数値例
F2層の電子密度 $N_e = 10^{12}$ m$^{-3}$ のとき、
$$ f_p = 9\sqrt{10^{12}} = 9 \times 10^6 = 9 \text{ MHz} $$
屈折率の導出
屈折率とプラズマ周波数
電離層の屈折率 $n$ は $n^2 = \varepsilon_r$（非磁性媒質なので $\mu_r = 1$）から、
$$ \begin{equation} n^2 = 1 – \frac{f_p^2}{f^2} = 1 – \frac{\omega_p^2}{\omega^2} = 1 – \frac{N_e e^2}{\varepsilon_0 m_e \omega^2} \end{equation} $$
これを アプルトン・ハートリーの式（磁場と衝突を無視した簡略版）と呼びます。
屈折率の周波数依存性
$f > f_p$ のとき $n^2 > 0$ であり、$n < 1$ です。電磁波は伝搬しますが、位相速度は $v_p = c/n > c$ となります。
$f < f_p$ のとき $n^2 < 0$ であり、$n$ は純虚数です。電磁波は指数的に減衰し、伝搬できません。これは カットオフ と呼ばれます。
$f = f_p$ のとき $n = 0$ であり、波は完全に反射されます。
位相速度と群速度
電離層中の位相速度と群速度は以下のようになります。
位相速度:
$$ v_p = \frac{c}{n} = \frac{c}{\sqrt{1 – f_p^2/f^2}} > c $$
群速度:
$$ v_g = cn = c\sqrt{1 – \frac{f_p^2}{f^2}} < c $$
この関係を確認します。分散関係 $\omega^2 = \omega_p^2 + c^2 k^2$ から、
$$ v_g = \frac{d\omega}{dk} = \frac{c^2 k}{\omega} = \frac{c^2 k}{\omega} = c \cdot \frac{ck}{\omega} = cn $$
また、$v_p \cdot v_g = c^2$ が成り立ちます。
群速度は情報・エネルギーの伝達速度であり、常に $c$ 以下です。位相速度が $c$ を超えても相対性理論と矛盾しません。
臨界周波数と最大使用周波数
臨界周波数
垂直入射（真上に向けて電波を発射した場合）で電離層から反射される最高周波数を 臨界周波数 $f_c$ と呼びます。
垂直入射では、電子密度の最大値 $N_{\max}$ で $n = 0$ になる条件から、
$$ \begin{equation} f_c = f_p(N_{\max}) = 9\sqrt{N_{\max}} \text{ [Hz]} \end{equation} $$
$f_c$ より高い周波数の垂直入射波は電離層を突き抜けます。
最大使用周波数（MUF）
斜め入射では、スネルの法則により垂直入射よりも高い周波数でも反射が可能です。
スネルの法則による導出
成層大気（水平に均一な層が重なった構造）を仰角 $\theta_0$（水平面からの角度）で入射する電波を考えます。地上での入射角（法線からの角度）を $\phi_0 = 90° – \theta_0$ とすると、スネルの法則は、
$$ n_0 \sin\phi_0 = n(h)\sin\phi(h) $$
地上では $n_0 = 1$ なので、
$$ \sin\phi_0 = n(h)\sin\phi(h) $$
反射が起こる条件は $\phi(h) = 90°$（電波が水平になる）、すなわち、
$$ \sin\phi_0 = n(h_{\text{ref}}) $$
$n^2 = 1 – f_p^2/f^2$ を代入し、$\sin^2\phi_0 = 1 – f_p^2(h_{\text{ref}})/f^2$ で反射が起こるので、
$$ \cos^2\phi_0 = \frac{f_p^2(h_{\text{ref}})}{f^2} $$
最大電子密度の高さで反射が起こる場合（ぎりぎり反射）、
$$ f = \frac{f_c}{\cos\phi_0} = \frac{f_c}{\sin\theta_0} $$
$$ \begin{equation} \boxed{f_{\text{MUF}} = \frac{f_c}{\cos\phi_0} = \frac{f_c}{\sin\theta_0}} \end{equation} $$
ここで $\theta_0$ は仰角（水平面からの角度）です。
セカント法則
上の関係は セカント法則 とも呼ばれます。
$$ f_{\text{MUF}} = f_c \sec\phi_0 $$
ただし、これは平面地球モデルでの近似です。地球の曲率を考慮すると、修正が必要です。
数値例
$f_c = 10$ MHz、仰角 $\theta_0 = 10°$（入射角 $\phi_0 = 80°$）のとき、
$$ f_{\text{MUF}} = \frac{10}{\cos 80°} = \frac{10}{0.1736} = 57.6 \text{ MHz} $$
低い仰角ほどMUFが高くなり、より高い周波数でも電離層反射が利用できます。
跳躍距離
電離層で反射された電波が地上に戻る距離（跳躍距離）は、仰角が低いほど長くなります。1回のホップで数百～数千 kmの距離をカバーでき、複数回のホップで地球の裏側まで到達可能です。
衛星通信への電離層影響
ファラデー回転
地磁場が存在する電離層中では、直線偏波の電磁波の偏波面が回転します。これを ファラデー回転 と呼びます。
ファラデー回転角 $\Omega$ は、
$$ \begin{equation} \Omega = \frac{e^3}{8\pi^2 \varepsilon_0 m_e^2 c} \cdot \frac{1}{f^2} \int_{\text{path}} N_e(s) B_{\parallel}(s) \, ds \end{equation} $$
ここで $B_{\parallel}$ は伝搬方向の地磁場成分です。
ファラデー回転の特徴：
- 周波数の2乗に反比例: $\Omega \propto 1/f^2$
- TEC（全電子数）に比例
- 1 GHz以下では数十度の回転がありうる
- 対策: 円偏波 の使用（円偏波はファラデー回転の影響を受けない）
群遅延（電離層遅延）
電離層中では群速度が $v_g = c\sqrt{1 – f_p^2/f^2} < c$ なので、電波は真空中より遅れて到達します。
群遅延の余分な時間は、
$$ \Delta t = \int_{\text{path}}\left(\frac{1}{v_g} – \frac{1}{c}\right)ds $$
$f \gg f_p$ の近似（$\sqrt{1 – x} \approx 1 – x/2$ for $x \ll 1$）を使うと、
$$ \frac{1}{v_g} = \frac{1}{c\sqrt{1 – f_p^2/f^2}} \approx \frac{1}{c}\left(1 + \frac{f_p^2}{2f^2}\right) $$
$$ \frac{1}{v_g} – \frac{1}{c} \approx \frac{f_p^2}{2cf^2} = \frac{N_e e^2}{8\pi^2 \varepsilon_0 m_e c f^2} $$
したがって、
$$ \begin{equation} \Delta t \approx \frac{e^2}{8\pi^2 \varepsilon_0 m_e c f^2} \int_{\text{path}} N_e \, ds = \frac{40.3}{c f^2} \cdot \text{TEC} \end{equation} $$
ここで TEC（Total Electron Content） は経路に沿った全電子数 [el/m$^2$] です。
$$ \text{TEC} = \int_{\text{path}} N_e \, ds $$
TEC の単位は TECU（$1 \text{ TECU} = 10^{16}$ el/m$^2$）がよく使われます。
数値例：GPS信号の電離層遅延
GPS L1周波数 $f = 1575.42$ MHz、典型的な垂直TEC = 30 TECU のとき、
$$ \Delta t = \frac{40.3}{3 \times 10^8 \times (1575.42 \times 10^6)^2} \times 30 \times 10^{16} $$
$$ = \frac{40.3 \times 30 \times 10^{16}}{3 \times 10^8 \times 2.482 \times 10^{18}} = \frac{1.209 \times 10^{19}}{7.446 \times 10^{26}} \approx 16.2 \text{ ns} $$
距離換算で $\Delta d = c \Delta t \approx 4.9$ m の測位誤差に相当します。これがGPSの電離層補正が重要な理由です。
シンチレーション
電離層の不規則構造（電子密度のゆらぎ）により、電波の振幅と位相が変動する現象を 電離層シンチレーション と呼びます。赤道域と極域で顕著であり、衛星通信やGNSSの信号品質を劣化させます。
Pythonでの実装
電離層の電子密度プロファイル
import numpy as np
import matplotlib.pyplot as plt
# --- 電離層の電子密度プロファイル ---
def chapman_layer(h, h_max, N_max, H):
"""
チャップマン関数による電離層モデル
Parameters
----------
h : array
高度 [km]
h_max : float
電子密度最大の高度 [km]
N_max : float
最大電子密度 [m^-3]
H : float
スケールハイト [km]
Returns
-------
Ne : array
電子密度 [m^-3]
"""
z = (h - h_max) / H
return N_max * np.exp(0.5 * (1 - z - np.exp(-z)))
# 各層のパラメータ（昼間モデル）
layers_day = {
'D': {'h_max': 80, 'N_max': 1e9, 'H': 10},
'E': {'h_max': 110, 'N_max': 1e11, 'H': 10},
'F1': {'h_max': 200, 'N_max': 2e11, 'H': 25},
'F2': {'h_max': 300, 'N_max': 1e12, 'H': 50},
}
# 夜間モデル（D層、F1層消滅）
layers_night = {
'E': {'h_max': 110, 'N_max': 5e9, 'H': 10},
'F2': {'h_max': 350, 'N_max': 5e11, 'H': 70},
}
h = np.linspace(50, 800, 1000)
fig, axes = plt.subplots(1, 3, figsize=(18, 7))
# (1) 昼間の各層と合成プロファイル
ax = axes[0]
Ne_total_day = np.zeros_like(h)
for name, params in layers_day.items():
Ne = chapman_layer(h, **params)
Ne_total_day += Ne
ax.plot(Ne, h, '--', linewidth=1.5, label=f'{name} layer', alpha=0.7)
ax.plot(Ne_total_day, h, 'k-', linewidth=2.5, label='Total (daytime)')
ax.set_xlabel('Electron Density [m$^{-3}$]', fontsize=12)
ax.set_ylabel('Altitude [km]', fontsize=12)
ax.set_title('Daytime Ionospheric Profile', fontsize=13)
ax.set_xscale('log')
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3, which='both')
ax.set_xlim([1e7, 2e12])
ax.set_ylim([50, 800])
# (2) 昼夜比較
ax = axes[1]
Ne_total_night = np.zeros_like(h)
for name, params in layers_night.items():
Ne = chapman_layer(h, **params)
Ne_total_night += Ne
ax.plot(Ne_total_day, h, 'r-', linewidth=2.5, label='Daytime')
ax.plot(Ne_total_night, h, 'b-', linewidth=2.5, label='Nighttime')
ax.fill_betweenx(h, Ne_total_night, Ne_total_day, alpha=0.1, color='orange')
ax.set_xlabel('Electron Density [m$^{-3}$]', fontsize=12)
ax.set_ylabel('Altitude [km]', fontsize=12)
ax.set_title('Day vs Night Comparison', fontsize=13)
ax.set_xscale('log')
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3, which='both')
ax.set_xlim([1e7, 2e12])
ax.set_ylim([50, 800])
# (3) プラズマ周波数のプロファイル
ax = axes[2]
fp_day = 9 * np.sqrt(Ne_total_day) / 1e6 # [MHz]
fp_night = 9 * np.sqrt(Ne_total_night) / 1e6 # [MHz]
ax.plot(fp_day, h, 'r-', linewidth=2.5, label='Daytime')
ax.plot(fp_night, h, 'b-', linewidth=2.5, label='Nighttime')
ax.axvline(x=9 * np.sqrt(1e12) / 1e6, color='gray', linestyle=':',
linewidth=1.5, label=f'$f_c$ (day) = {9*np.sqrt(1e12)/1e6:.1f} MHz')
ax.set_xlabel('Plasma Frequency [MHz]', fontsize=12)
ax.set_ylabel('Altitude [km]', fontsize=12)
ax.set_title('Plasma Frequency Profile', fontsize=13)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
ax.set_xlim([0, 12])
ax.set_ylim([50, 800])
plt.tight_layout()
plt.savefig('ionosphere_profile.png', dpi=150, bbox_inches='tight')
plt.show()
print("=== 各層の臨界周波数 ===")
for name, params in layers_day.items():
fc = 9 * np.sqrt(params['N_max']) / 1e6
print(f" {name}層: N_max = {params['N_max']:.0e} m^-3, "
f"f_c = {fc:.2f} MHz")
屈折率と位相速度・群速度
import numpy as np
import matplotlib.pyplot as plt
# --- 屈折率の周波数依存性 ---
# プラズマ周波数
fp = 9.0 # MHz（F2層の臨界周波数）
f = np.linspace(0.1, 50, 1000) # MHz
# 屈折率 n² = 1 - (fp/f)²
n_sq = 1 - (fp / f)**2
# n (実数部のみ、n² > 0の領域)
n = np.where(n_sq > 0, np.sqrt(n_sq), np.nan)
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# (1) n² vs 周波数
ax = axes[0, 0]
ax.plot(f, n_sq, 'b-', linewidth=2.5)
ax.axhline(y=0, color='red', linestyle='--', linewidth=1)
ax.axvline(x=fp, color='green', linestyle=':', linewidth=1.5,
label=f'$f_p$ = {fp} MHz')
ax.fill_between(f, n_sq, 0, where=n_sq < 0, color='red', alpha=0.15,
label='Evanescent (cutoff)')
ax.set_xlabel('Frequency [MHz]', fontsize=12)
ax.set_ylabel('$n^2$', fontsize=12)
ax.set_title('Refractive Index Squared', fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_ylim([-2, 1.2])
# (2) 屈折率 n vs 周波数
ax = axes[0, 1]
ax.plot(f, n, 'b-', linewidth=2.5, label='$n = \\sqrt{1 - (f_p/f)^2}$')
ax.axhline(y=1, color='gray', linestyle='--', linewidth=1)
ax.axvline(x=fp, color='green', linestyle=':', linewidth=1.5,
label=f'$f_p$ = {fp} MHz')
ax.set_xlabel('Frequency [MHz]', fontsize=12)
ax.set_ylabel('Refractive index $n$', fontsize=12)
ax.set_title('Refractive Index', fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_ylim([0, 1.1])
# (3) 位相速度と群速度
ax = axes[1, 0]
c_norm = 1 # 光速で正規化
f_prop = f[f > fp] # 伝搬可能な周波数のみ
n_prop = np.sqrt(1 - (fp / f_prop)**2)
vp = c_norm / n_prop
vg = c_norm * n_prop
ax.plot(f_prop, vp, 'r-', linewidth=2, label='$v_p/c = 1/n$')
ax.plot(f_prop, vg, 'b-', linewidth=2, label='$v_g/c = n$')
ax.axhline(y=1, color='gray', linestyle='--', linewidth=1, label='$c$')
ax.axvline(x=fp, color='green', linestyle=':', linewidth=1.5)
ax.set_xlabel('Frequency [MHz]', fontsize=12)
ax.set_ylabel('Velocity / $c$', fontsize=12)
ax.set_title('Phase and Group Velocity', fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
ax.set_ylim([0, 5])
ax.set_xlim([0, 50])
# (4) MUF vs 仰角
ax = axes[1, 1]
fc_values = [5, 7, 9, 12, 15] # 臨界周波数 [MHz]
elevation = np.linspace(1, 90, 200) # 仰角 [度]
for fc in fc_values:
muf = fc / np.sin(np.radians(elevation))
muf_clipped = np.clip(muf, 0, 100)
ax.plot(elevation, muf_clipped, linewidth=2, label=f'$f_c$ = {fc} MHz')
ax.set_xlabel('Elevation Angle [deg]', fontsize=12)
ax.set_ylabel('MUF [MHz]', fontsize=12)
ax.set_title('Maximum Usable Frequency vs Elevation', fontsize=13)
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3)
ax.set_ylim([0, 100])
plt.tight_layout()
plt.savefig('ionosphere_refraction.png', dpi=150, bbox_inches='tight')
plt.show()
電離層レイトレーシング
import numpy as np
import matplotlib.pyplot as plt
# --- 電離層レイトレーシング ---
def electron_density_profile(h, N_max=1e12, h_max=300, H=50):
"""チャップマン層の電子密度"""
z = (h - h_max) / H
Ne = N_max * np.exp(0.5 * (1 - z - np.exp(-z)))
return np.where(h > 60, Ne, 0) # 60 km以下は電離層なし
def plasma_freq(Ne):
"""プラズマ周波数 [Hz]"""
e = 1.602e-19
eps0 = 8.854e-12
me = 9.109e-31
return np.sqrt(Ne * e**2 / (eps0 * me)) / (2 * np.pi)
def refractive_index(f, Ne):
"""屈折率"""
fp = plasma_freq(Ne)
n_sq = 1 - (fp / f)**2
return np.where(n_sq > 0, np.sqrt(n_sq), 0)
def ray_trace(f, elevation_deg, N_max=1e12, h_max=300, H=50,
R_earth=6371e3, dh=500):
"""
電離層レイトレーシング（球面地球モデル）
Parameters
----------
f : float
周波数 [Hz]
elevation_deg : float
初期仰角 [度]
N_max : float
最大電子密度 [m^-3]
h_max : float
最大電子密度の高度 [km]
H : float
スケールハイト [km]
R_earth : float
地球半径 [m]
dh : float
高度ステップ [m]
Returns
-------
x_path, y_path : arrays
レイパスの座標 [m]
reflected : bool
反射が起こったか
"""
elevation_rad = np.radians(elevation_deg)
# 初期条件
r = R_earth # 地表面の距離
theta = 0 # 地心角
# レイの方向（極座標で）
# 地表面での入射角（法線＝動径方向からの角度）
phi = np.pi / 2 - elevation_rad # 法線からの角度
# スネルの法則の保存量（球面）: n * r * sin(phi) = const
# ボウヘメットの法則
n0 = 1.0 # 地表面の屈折率
snell_const = n0 * r * np.sin(phi)
x_path = [r * np.sin(theta)]
y_path = [r * np.cos(theta)]
going_up = True
max_h = 0
reflected = False
for _ in range(100000):
h_km = (r - R_earth) / 1e3 # 高度 [km]
if h_km > 800:
break # 電離層を通過
if h_km < 0 and not going_up:
break # 地表に到達
# 電子密度と屈折率
Ne = electron_density_profile(h_km, N_max, h_max, H)
n = refractive_index(f, Ne)
if n < 1e-6:
# カットオフ → 反射
going_up = False
reflected = True
n = 1e-6
# ボウヘメットの法則: sin(phi) = snell_const / (n * r)
sin_phi = snell_const / (n * r)
if sin_phi > 1:
# 全反射
going_up = False
reflected = True
sin_phi = 1.0
cos_phi = np.sqrt(max(1 - sin_phi**2, 0))
if not going_up:
cos_phi = -cos_phi # 下向き
# 微小ステップ ds に沿った移動
ds = dh if going_up else -dh
# 座標の更新（極座標）
dr = ds * cos_phi
r_dtheta = ds * sin_phi
r_new = r + dr
theta_new = theta + r_dtheta / r
if r_new < R_earth:
# 地表面に到達
r_new = R_earth
x_path.append(r_new * np.sin(theta_new))
y_path.append(r_new * np.cos(theta_new))
break
r = r_new
theta = theta_new
max_h = max(max_h, (r - R_earth) / 1e3)
x_path.append(r * np.sin(theta))
y_path.append(r * np.cos(theta))
# 反射後に地表面近くまで降りてきたら上昇に切り替わる判定は不要
# （1ホップのみ）
return np.array(x_path), np.array(y_path), reflected
# --- レイトレーシング実行 ---
fig, axes = plt.subplots(1, 2, figsize=(16, 8))
R_earth = 6371e3
N_max = 1e12
h_max = 300 # km
H = 50
# (1) 周波数を変えたレイパス（仰角固定）
ax = axes[0]
# 地球の描画
theta_earth = np.linspace(-0.1, 0.15, 300)
x_earth = R_earth * np.sin(theta_earth)
y_earth = R_earth * np.cos(theta_earth)
ax.fill_between(x_earth / 1e3, 0, y_earth / 1e3, color='lightblue', alpha=0.3)
ax.plot(x_earth / 1e3, y_earth / 1e3, 'b-', linewidth=1)
# 電離層の描画（高度60-800 km）
for h_layer in [60, 300, 800]:
r_layer = R_earth + h_layer * 1e3
ax.plot(r_layer * np.sin(theta_earth) / 1e3,
r_layer * np.cos(theta_earth) / 1e3,
'g--', linewidth=0.5, alpha=0.5)
elevation = 20 # 仰角 20度
freqs_mhz = [3, 5, 7, 9, 12, 15, 20, 30]
colors = plt.cm.plasma(np.linspace(0.1, 0.9, len(freqs_mhz)))
for f_mhz, color in zip(freqs_mhz, colors):
f_hz = f_mhz * 1e6
x_path, y_path, reflected = ray_trace(
f_hz, elevation, N_max=N_max, h_max=h_max, H=H, dh=200
)
label = f'{f_mhz} MHz'
if reflected:
label += ' (reflected)'
ls = '-' if reflected else '--'
ax.plot(x_path / 1e3, y_path / 1e3, color=color, linewidth=2,
linestyle=ls, label=label)
ax.set_xlabel('Horizontal Distance [km]', fontsize=12)
ax.set_ylabel('Vertical [km]', fontsize=12)
ax.set_title(f'Ray Tracing: Varying Frequency (elevation = {elevation}°)',
fontsize=13)
ax.legend(fontsize=8, loc='upper right')
ax.grid(True, alpha=0.3)
ax.set_aspect('equal')
# 表示範囲の調整
ax.set_xlim([-100, 2500])
ax.set_ylim([R_earth / 1e3 - 50, R_earth / 1e3 + 900])
# (2) 仰角を変えたレイパス（周波数固定）
ax = axes[1]
ax.fill_between(x_earth / 1e3, 0, y_earth / 1e3, color='lightblue', alpha=0.3)
ax.plot(x_earth / 1e3, y_earth / 1e3, 'b-', linewidth=1)
for h_layer in [60, 300, 800]:
r_layer = R_earth + h_layer * 1e3
ax.plot(r_layer * np.sin(theta_earth) / 1e3,
r_layer * np.cos(theta_earth) / 1e3,
'g--', linewidth=0.5, alpha=0.5)
f_fixed = 7e6 # 7 MHz
elevations = [5, 10, 15, 20, 30, 45, 60, 80]
colors = plt.cm.viridis(np.linspace(0.1, 0.9, len(elevations)))
for elev, color in zip(elevations, colors):
x_path, y_path, reflected = ray_trace(
f_fixed, elev, N_max=N_max, h_max=h_max, H=H, dh=200
)
label = f'{elev}°'
if reflected:
label += ' (ref.)'
ls = '-' if reflected else '--'
ax.plot(x_path / 1e3, y_path / 1e3, color=color, linewidth=2,
linestyle=ls, label=label)
ax.set_xlabel('Horizontal Distance [km]', fontsize=12)
ax.set_ylabel('Vertical [km]', fontsize=12)
ax.set_title(f'Ray Tracing: Varying Elevation (f = {f_fixed/1e6:.0f} MHz)',
fontsize=13)
ax.legend(fontsize=8, loc='upper right')
ax.grid(True, alpha=0.3)
ax.set_aspect('equal')
ax.set_xlim([-100, 2500])
ax.set_ylim([R_earth / 1e3 - 50, R_earth / 1e3 + 900])
plt.tight_layout()
plt.savefig('ionosphere_raytracing.png', dpi=150, bbox_inches='tight')
plt.show()
# 臨界周波数と MUF の確認
fc = 9 * np.sqrt(N_max) / 1e6
print(f"\n=== 電離層パラメータ ===")
print(f" 最大電子密度: N_max = {N_max:.0e} m^-3")
print(f" 臨界周波数: f_c = {fc:.1f} MHz")
print(f"\n=== MUF（セカント法則） ===")
for elev in [5, 10, 20, 30, 45, 90]:
muf = fc / np.sin(np.radians(elev))
print(f" 仰角 {elev:>2d}° → MUF = {muf:>7.1f} MHz")
ファラデー回転と群遅延の計算
import numpy as np
import matplotlib.pyplot as plt
# --- ファラデー回転と群遅延 ---
# 定数
e = 1.602e-19 # 電気素量 [C]
me = 9.109e-31 # 電子質量 [kg]
eps0 = 8.854e-12 # 真空誘電率 [F/m]
c = 3e8 # 光速 [m/s]
# TEC (Total Electron Content)
TEC_TECU = np.linspace(1, 100, 200) # [TECU]
TEC = TEC_TECU * 1e16 # [el/m²]
# 群遅延の計算: Δt = 40.3 * TEC / (c * f²)
def group_delay_ns(TEC, f_hz):
"""群遅延 [ns]"""
return 40.3 * TEC / (c * f_hz**2) * 1e9
def range_error_m(TEC, f_hz):
"""測距誤差 [m]"""
return 40.3 * TEC / f_hz**2
# ファラデー回転角の計算（簡易モデル）
# Ω ≈ 1.355e-5 * B_parallel * TEC / f² [rad]
# B_parallel ≈ 4e-5 T (地磁場の典型値)
B_parallel = 4e-5 # [T]
def faraday_rotation_deg(TEC, f_hz, B=4e-5):
"""ファラデー回転角 [度]"""
omega = (e**3 * B * TEC) / (8 * np.pi**2 * eps0 * me**2 * c * f_hz**2)
return np.degrees(omega)
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# (1) 群遅延 vs TEC（各周波数）
ax = axes[0, 0]
frequencies = {
'150 MHz (VHF)': 150e6,
'400 MHz (UHF)': 400e6,
'1.2 GHz (GPS L2)': 1227.6e6,
'1.6 GHz (GPS L1)': 1575.42e6,
'2.3 GHz (S-band)': 2.3e9,
'8 GHz (X-band)': 8e9,
}
for name, f_val in frequencies.items():
dt = group_delay_ns(TEC, f_val)
ax.plot(TEC_TECU, dt, linewidth=2, label=name)
ax.set_xlabel('Vertical TEC [TECU]', fontsize=12)
ax.set_ylabel('Group Delay [ns]', fontsize=12)
ax.set_title('Ionospheric Group Delay', fontsize=13)
ax.legend(fontsize=8)
ax.grid(True, alpha=0.3)
ax.set_yscale('log')
ax.set_ylim([0.01, 10000])
# (2) 測距誤差 vs 周波数（TEC固定）
ax = axes[0, 1]
f_range = np.logspace(np.log10(100e6), np.log10(30e9), 500)
TEC_values = [10, 30, 50, 100] # TECU
for tec_val in TEC_values:
tec = tec_val * 1e16
dr = range_error_m(tec, f_range)
ax.plot(f_range / 1e9, dr, linewidth=2, label=f'TEC = {tec_val} TECU')
ax.axhline(y=1, color='gray', linestyle=':', linewidth=1, label='1 m error')
ax.axhline(y=0.1, color='gray', linestyle='--', linewidth=1, label='10 cm error')
ax.set_xlabel('Frequency [GHz]', fontsize=12)
ax.set_ylabel('Range Error [m]', fontsize=12)
ax.set_title('Ionospheric Range Error', fontsize=13)
ax.set_xscale('log')
ax.set_yscale('log')
ax.legend(fontsize=9)
ax.grid(True, alpha=0.3, which='both')
ax.set_ylim([1e-4, 1e3])
# (3) ファラデー回転 vs 周波数
ax = axes[1, 0]
tec_fixed = 30e16 # 30 TECU
fr_deg = faraday_rotation_deg(tec_fixed, f_range)
ax.plot(f_range / 1e9, fr_deg, 'b-', linewidth=2.5)
ax.axhline(y=90, color='red', linestyle=':', linewidth=1.5,
label='90° (critical for linear pol.)')
ax.axhline(y=10, color='orange', linestyle='--', linewidth=1,
label='10°')
ax.set_xlabel('Frequency [GHz]', fontsize=12)
ax.set_ylabel('Faraday Rotation [deg]', fontsize=12)
ax.set_title(f'Faraday Rotation (TEC = 30 TECU, B = {B_parallel*1e5:.0f}$\\mu$T)',
fontsize=13)
ax.set_xscale('log')
ax.set_yscale('log')
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3, which='both')
ax.set_ylim([1e-3, 1e4])
# (4) GPSの2周波補正
ax = axes[1, 1]
f_L1 = 1575.42e6
f_L2 = 1227.60e6
# 各周波数での遅延
delay_L1 = range_error_m(TEC, f_L1)
delay_L2 = range_error_m(TEC, f_L2)
# 2周波補正: Δρ = (f1² * ρ1 - f2² * ρ2) / (f1² - f2²) で電離層遅延を除去
# 残留誤差は高次項のみ（ここでは理論上ゼロ）
# 代わりに各周波数の遅延を比較
ax.plot(TEC_TECU, delay_L1, 'b-', linewidth=2, label=f'L1 ({f_L1/1e6:.1f} MHz)')
ax.plot(TEC_TECU, delay_L2, 'r-', linewidth=2, label=f'L2 ({f_L2/1e6:.1f} MHz)')
ax.plot(TEC_TECU, delay_L2 - delay_L1, 'g--', linewidth=2,
label='$\\Delta$(L2 - L1)')
# L1の代表的な値に注釈
tec_30 = 30e16
dl1_30 = range_error_m(tec_30, f_L1)
ax.annotate(f'TEC=30: {dl1_30:.1f} m', xy=(30, dl1_30),
xytext=(50, dl1_30 + 3), fontsize=10,
arrowprops=dict(arrowstyle='->', color='blue'))
ax.set_xlabel('Vertical TEC [TECU]', fontsize=12)
ax.set_ylabel('Range Error [m]', fontsize=12)
ax.set_title('GPS Ionospheric Delay', fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('ionosphere_effects.png', dpi=150, bbox_inches='tight')
plt.show()
# 数値サマリー
print("=== GPS電離層遅延（TEC = 30 TECU） ===")
tec = 30e16
for name, f_val in [('L1', 1575.42e6), ('L2', 1227.60e6), ('L5', 1176.45e6)]:
delay = range_error_m(tec, f_val)
dt = group_delay_ns(tec, f_val)
fr = faraday_rotation_deg(tec, f_val)
print(f" {name}: f = {f_val/1e6:.2f} MHz, "
f"遅延 = {dt:.2f} ns, "
f"測距誤差 = {delay:.2f} m, "
f"Faraday = {fr:.1f}°")
まとめ
本記事では、電離層伝搬と電波の屈折について基礎から丁寧に解説しました。
- 電離層の構成: D/E/F1/F2層の4層構造で、F2層が最も電子密度が高く短波通信の主要反射層です
- プラズマ周波数: $f_p = 9\sqrt{N_e}$ Hz（$N_e$ [m$^{-3}$]）で、電離層の電波特性を決定する基本量です
- 屈折率: $n^2 = 1 – (f_p/f)^2$ を導出しました。$f < f_p$ ではカットオフ（電波が伝搬不能）です
- 位相速度と群速度: $v_p = c/n > c$, $v_g = cn < c$ で、$v_p v_g = c^2$ の関係があります
- 臨界周波数とMUF: 垂直入射の臨界周波数 $f_c$ とセカント法則 $f_{\text{MUF}} = f_c/\cos\phi_0$ を導出しました
- ファラデー回転: 偏波面の回転は $\propto 1/f^2$ で、低周波ほど影響が大きいため衛星通信では円偏波が使われます
- 群遅延: $\Delta t \propto \text{TEC}/f^2$ で、GPS L1帯で数メートルの測位誤差に相当します
次のステップとして、以下の記事も参考にしてください。