トカマク平衡における磁気座標系#

このチュートリアルでは、磁気核融合で使われる 4 種類の磁束座標系、 PEST、Boozer、Hamada、Equal-arc を比較します。

磁束座標とは何か#

トカマクでは、磁力線は 磁束面 と呼ばれる入れ子状のトロイダル面上に乗ります。磁束座標系 \((\psi, \theta, \varphi)\) とは、次の性質を持つ座標です。

\(\psi\) は磁束面を標識する（例：\(\psi_{\rm norm} = 0\) が軸、\(1\) が LCFS）
\(\varphi\) は標準的なトロイダル角
\(\theta\) は有用な性質を持つように定義されたポロイダル角

4 つの系は \(\theta\) の選び方だけが異なります。

系	ポロイダル角 \(\theta\) の定義	主な性質
PEST	\(\partial(\mathbf{B}\cdot\nabla\theta^) / \partial\theta^ = 0\)：直線磁力線	最も単純な直線磁力線座標
Boozer	PEST に加えてヤコビアンが磁束関数：\(\sqrt{g} = \sqrt{g}(\psi)\)	\(\mathbf{J}\cdot\nabla\varphi = {\rm const}(\psi)\)；W7-X、LHD コードで使用
Hamada	磁気軸からの等面積	\(\mathbf{J}\cdot\nabla\theta = {\rm const}(\psi)\)；MHD 安定性行列を簡略化
Equal-arc	磁束面に沿う弧長が \(\theta_{\rm ea}\) で一様	数値的に堅牢；構成が単純

なぜこの選び方が重要か#

直線磁力線（PEST、Boozer、Hamada）：トロイダルモード数 \(n\) とポロイダルモード数 \(m\) を持つモードは、共鳴条件 \(q = m/n\) を満たします。これにより MHD モードの Fourier 解析が簡単になります。
Boozer：\(1/B^2\) ドリフトが純粋に径方向となり、新古典輸送理論と運動方程式が簡略化されます。
Hamada：電流線も直線になるため、一部の MHD 安定性コードで必要になります。
Equal-arc：数値格子に実用的で、プラズマ境界をよく分解します。

セットアップ：Solov’ev解析平衡#

EAST サイズ（R0 ≈ 1.86 m）へスケーリングした Cerfon & Freidberg (2010) の Solov’ev 平衡を使います。

[1]:

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from matplotlib.gridspec import GridSpec

from pyna.toroidal.equilibrium.Solovev import solovev_iter_like
from pyna.toroidal.coords.PEST import build_PEST_mesh
from pyna.toroidal.coords.EqualArc import build_equal_arc_mesh
from pyna.toroidal.coords.Hamada import build_Hamada_mesh
from pyna.toroidal.coords.Boozer import build_Boozer_mesh

# Create equilibrium (scale=0.3 → EAST-sized, R0≈1.86 m)
eq = solovev_iter_like(scale=0.3)
print(f"Equilibrium: R0={eq.R0:.2f} m, a={eq.a:.2f} m, B0={eq.B0:.1f} T")
print(f"κ={eq.kappa:.2f}, δ={eq.delta:.2f}, q0={eq.q0:.2f}")
Rmaxis, Zmaxis = eq.magnetic_axis
print(f"Magnetic axis: R={Rmaxis:.3f} m, Z={Zmaxis:.3f} m")

Equilibrium: R0=1.86 m, a=0.60 m, B0=5.3 T
κ=1.70, δ=0.33, q0=1.50
Magnetic axis: R=1.946 m, Z=0.000 m

[2]:

# Build the background grid for the equilibrium
nR, nZ = 150, 150
R_grid = np.linspace(0.3 * eq.R0, 1.5 * eq.R0, nR)
Z_grid = np.linspace(-eq.a * eq.kappa * 1.3, eq.a * eq.kappa * 1.3, nZ)
Rg, Zg = np.meshgrid(R_grid, Z_grid, indexing='ij')

BR_grid, BZ_grid = eq.BR_BZ(Rg, Zg)
BPhi_grid = eq.Bphi(Rg)
psi_norm_grid = eq.psi(Rg, Zg)

# Build PEST mesh
# ns=40 radial surfaces, ntheta=181 poloidal points
ns = 40
ntheta = 181
S, TET, R_mesh, Z_mesh, q_iS = build_PEST_mesh(
    R_grid, Z_grid, BR_grid, BZ_grid, BPhi_grid, psi_norm_grid,
    Rmaxis, Zmaxis, ns=ns, ntheta=ntheta
)
print(f"PEST mesh built: {ns} surfaces, {ntheta} poloidal points")
print(f"S range: [{S[1]:.3f}, {S[-1]:.3f}]")
print(f"q range: [{q_iS[1]:.2f}, {q_iS[-1]:.2f}]")

PEST mesh built: 40 surfaces, 181 poloidal points
S range: [0.028, 0.972]
q range: [1.49, 2.39]

第 3 節：等弧長（Equal-arc）座標#

Equal-arc 角 \(\theta_{\rm ea}\) は、各磁束面に沿う弧長が \(\theta_{\rm ea}\) に比例するようにパラメータ化します。これは最も単純な非自明座標変換です。

[3]:

# Build equal-arc mesh
_, TET_ea, R_ea, Z_ea = build_equal_arc_mesh(S, TET, R_mesh, Z_mesh, n_theta=ntheta)
print(f"Equal-arc mesh built: shape {R_ea.shape}")

# Verify arc length uniformity on a mid-radius surface
i_mid = ns // 2
R_s = R_ea[i_mid, :]
Z_s = Z_ea[i_mid, :]
R_loop = R_s[:-1]; Z_loop = Z_s[:-1]  # drop endpoint duplicate
R_closed = np.append(R_loop, R_loop[0])
Z_closed = np.append(Z_loop, Z_loop[0])
ds = np.sqrt(np.diff(R_closed)**2 + np.diff(Z_closed)**2)
print(f"Arc-length segments at surface {i_mid}: mean={ds.mean()*100:.2f} cm, std={ds.std()*100:.3f} cm")
print(f"  Uniformity (std/mean): {ds.std()/ds.mean():.4f}")

Equal-arc mesh built: shape (40, 181)
Arc-length segments at surface 20: mean=1.27 cm, std=0.000 cm
  Uniformity (std/mean): 0.0001

[4]:

# Plot equal-arc coordinates
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# Left: R-Z view of flux surfaces with equal-arc mesh lines
ax = axes[0]
# Plot psi_norm contours as background
cs = ax.contour(Rg, Zg, psi_norm_grid, levels=np.linspace(0, 1, 11),
                colors='lightgray', linewidths=0.5)
# Plot every 5th surface in equal-arc
stride = max(1, ns // 10)
colors = cm.plasma(np.linspace(0.2, 0.9, ns // stride + 1))
for k, i in enumerate(range(1, ns, stride)):
    ax.plot(R_ea[i, :], Z_ea[i, :], color=colors[k], lw=1.0)
# Plot a few poloidal lines (constant θ_ea)
for j in range(0, ntheta-1, ntheta // 8):
    ax.plot(R_ea[1:, j], Z_ea[1:, j], 'b-', lw=0.5, alpha=0.5)
ax.plot(Rmaxis, Zmaxis, 'r+', ms=10, mew=2, label='Axis')
ax.set_xlabel('R (m)')
ax.set_ylabel('Z (m)')
ax.set_title('Equal-arc mesh (R-Z view)')
ax.set_aspect('equal')
ax.legend()

# Right: arc-length increments as function of theta_ea
ax = axes[1]
for k, i in enumerate(range(2, ns-1, stride)):
    R_s = R_ea[i, :-1]; Z_s = Z_ea[i, :-1]
    R_c = np.append(R_s, R_s[0]); Z_c = np.append(Z_s, Z_s[0])
    ds_s = np.sqrt(np.diff(R_c)**2 + np.diff(Z_c)**2)
    ds_s_norm = ds_s / ds_s.mean()
    ax.plot(TET_ea[:-1], ds_s_norm, color=colors[k], lw=0.8,
            label=f'S={S[i]:.2f}' if k % 2 == 0 else None)
ax.axhline(1.0, color='k', ls='--', lw=1, label='Ideal (uniform)')
ax.set_xlabel(r'$\theta_{\rm ea}$ (rad)')
ax.set_ylabel('Arc-length increment / mean')
ax.set_title('Arc-length uniformity in equal-arc coordinates')
ax.set_xlim(0, 2*np.pi)
ax.legend(fontsize=8, ncol=2)

plt.tight_layout()
plt.savefig('equal_arc_coords.png', dpi=100, bbox_inches='tight')
plt.show()
print("Equal-arc: arc-length segments are uniform to within ~1% (interpolation discretisation)")

Equal-arc: arc-length segments are uniform to within ~1% (interpolation discretisation)

第 4 節：PEST 直線磁力線座標#

PEST 座標では、安全係数 \(q(\psi)\) が次を満たします。

\[q(\psi) = \frac{\mathbf{B}\cdot\nabla\varphi}{\mathbf{B}\cdot\nabla\theta^*} = {\rm const \; on \; each \; surface}\]

これは、磁力線が \((\theta^*, \varphi)\) 平面で傾き \(1/q\) の直線として見えることを意味します。

[5]:

# Show q(S) profile from PEST integration
q_valid = q_iS[1:]
S_valid = S[1:]

# Compare with analytic equilibrium q
psi_vals = S_valid**2
q_analytic = eq.q_profile(psi_vals, n_theta=256)

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

ax = axes[0]
ax.plot(S_valid, q_valid, 'b-o', ms=3, label='PEST q(S)')
finite = np.isfinite(q_analytic)
ax.plot(S_valid[finite], q_analytic[finite], 'r--', label='Analytic q(S)')
ax.set_xlabel('S = √ψ_norm')
ax.set_ylabel('Safety factor q')
ax.set_title('Safety factor profile: PEST vs analytic')
ax.legend()
ax.grid(True, alpha=0.3)

# PEST mesh in R-Z
ax = axes[1]
cs = ax.contour(Rg, Zg, psi_norm_grid, levels=np.linspace(0, 1, 11),
                colors='lightgray', linewidths=0.5)
stride_s = max(1, ns // 10)
colors_pest = cm.viridis(np.linspace(0.2, 0.9, ns // stride_s + 1))
for k, i in enumerate(range(1, ns, stride_s)):
    ax.plot(R_mesh[i, :], Z_mesh[i, :], color=colors_pest[k], lw=1.0)
for j in range(0, ntheta-1, ntheta // 8):
    ax.plot(R_mesh[1:, j], Z_mesh[1:, j], 'g-', lw=0.5, alpha=0.5)
ax.plot(Rmaxis, Zmaxis, 'r+', ms=10, mew=2)
ax.set_xlabel('R (m)')
ax.set_ylabel('Z (m)')
ax.set_title('PEST mesh (R-Z view)')
ax.set_aspect('equal')

plt.tight_layout()
plt.savefig('pest_coords.png', dpi=100, bbox_inches='tight')
plt.show()

[6]:

# Demonstrate straight field lines in PEST theta-phi space
# In PEST: a field line traces theta* = phi/q + const
fig, ax = plt.subplots(figsize=(8, 5))

phi_line = np.linspace(0, 4 * np.pi, 200)
surface_indices = [ns // 5, ns // 3, ns // 2, 2 * ns // 3, 4 * ns // 5]
colors_fl = cm.Set1(np.linspace(0, 0.8, len(surface_indices)))

for k, i_s in enumerate(surface_indices):
    if i_s >= ns or not np.isfinite(q_iS[i_s]):
        continue
    q_s = q_iS[i_s]
    # Field line: theta* = phi / q (PEST straight field line)
    theta_fl = (phi_line / q_s) % (2 * np.pi)
    ax.plot(phi_line / np.pi, theta_fl / np.pi, color=colors_fl[k],
            label=f'S={S[i_s]:.2f}, q={q_s:.2f}', lw=1.5)

ax.set_xlabel(r'$\varphi / \pi$')
ax.set_ylabel(r'$\theta^* / \pi$ (PEST)')
ax.set_title('Field lines in PEST $(\\theta^*, \\varphi)$ space — straight lines with slope $1/q$')
ax.legend(fontsize=9)
ax.set_xlim(0, 4)
ax.set_ylim(0, 2)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('pest_fieldlines.png', dpi=100, bbox_inches='tight')
plt.show()

第 5 節：Boozer 座標#

Boozer 座標は PEST を拡張し、さらにヤコビアン \(\sqrt{g}\) が磁束面量（\(\theta\) には依存せず \(\psi\) のみに依存）であることを要求します。これにより \(\mathbf{J}\cdot\nabla\varphi = {\rm const}(\psi)\) が成り立ちます。

構成：\(\theta_B = \theta^* + \lambda(\psi, \theta^*)\) とし、

\[\frac{\partial\lambda}{\partial\theta^*} = \frac{\langle\sqrt{g}\rangle}{\sqrt{g}} - 1\]

ここで \(\langle\cdot\rangle\) は磁束面平均 \(\frac{1}{2\pi}\int_0^{2\pi}\cdot\,d\theta^*\) を表します。

[7]:

# Build Boozer mesh
_, TET_B, R_B, Z_B, lambda_corr = build_Boozer_mesh(
    S, TET, R_mesh, Z_mesh, q_iS, equilibrium=eq, n_theta=ntheta
)
print(f"Boozer mesh built: shape {R_B.shape}")
print(f"Max |λ| correction: {np.nanmax(np.abs(lambda_corr)):.4f} rad")

# Show the angle correction lambda(psi, theta*)
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

ax = axes[0]
# lambda correction as a function of theta for several surfaces
for k, i_s in enumerate(range(2, ns-1, max(1, ns // 8))):
    lam = lambda_corr[i_s, :]
    if np.any(np.isfinite(lam)):
        ax.plot(TET / np.pi, lam, label=f'S={S[i_s]:.2f}',
                color=cm.plasma(0.1 + 0.8 * i_s / ns))
ax.set_xlabel(r'$\theta^* / \pi$ (PEST)')
ax.set_ylabel(r'$\lambda = \theta_B - \theta^*$ (rad)')
ax.set_title('Boozer angle correction $\\lambda(\\psi, \\theta^*)$')
ax.legend(fontsize=8, ncol=2)
ax.grid(True, alpha=0.3)

ax = axes[1]
cs = ax.contour(Rg, Zg, psi_norm_grid, levels=np.linspace(0, 1, 11),
                colors='lightgray', linewidths=0.5)
stride_s = max(1, ns // 10)
colors_b = cm.inferno(np.linspace(0.2, 0.9, ns // stride_s + 1))
for k, i in enumerate(range(1, ns, stride_s)):
    ax.plot(R_B[i, :], Z_B[i, :], color=colors_b[k], lw=1.0)
for j in range(0, ntheta-1, ntheta // 8):
    ax.plot(R_B[1:, j], Z_B[1:, j], 'm-', lw=0.5, alpha=0.5)
ax.plot(Rmaxis, Zmaxis, 'r+', ms=10, mew=2)
ax.set_xlabel('R (m)')
ax.set_ylabel('Z (m)')
ax.set_title('Boozer mesh (R-Z view)')
ax.set_aspect('equal')

plt.tight_layout()
plt.savefig('boozer_coords.png', dpi=100, bbox_inches='tight')
plt.show()

Boozer mesh built: shape (40, 181)
Max |λ| correction: 0.5429 rad

第 6 節：Hamada 座標#

Hamada 座標では、磁力線と電流密度線の両方が直線であることを要求します。

\[\mathbf{J}\cdot\nabla\psi = 0, \quad \mathbf{J}\cdot\nabla\theta_H = 0, \quad \mathbf{J}\cdot\nabla\varphi_H = {\rm const}(\psi)\]

軸対称平衡では、これはポロイダル角を 等面積 にすること、つまり磁気軸から角度 \(\theta_H\) まで掃いた面積が \(\theta_H\) に比例することへ帰着します。

[8]:

# Build Hamada mesh
_, TET_H, R_H, Z_H = build_Hamada_mesh(S, TET, R_mesh, Z_mesh, n_theta=ntheta)
print(f"Hamada mesh built: shape {R_H.shape}")

# Verify equal-area property
R_ax = R_mesh[0, 0]
Z_ax = Z_mesh[0, 0]

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

ax = axes[0]
# Show cumulative area vs theta_H for several surfaces
for k, i_s in enumerate(range(2, ns-2, max(1, ns // 8))):
    R_s = R_H[i_s, :-1]; Z_s = Z_H[i_s, :-1]
    R_c = np.append(R_s, R_s[0]); Z_c = np.append(Z_s, Z_s[0])
    dA = 0.5 * ((R_c[:-1] - R_ax) * (Z_c[1:] - Z_ax) -
                 (R_c[1:] - R_ax) * (Z_c[:-1] - Z_ax))
    A_cum = np.cumsum(dA)
    A_total = A_cum[-1]
    if abs(A_total) > 1e-10:
        ax.plot(TET_H[:-1] / np.pi, A_cum / A_total,
                color=cm.cool(0.1 + 0.8 * i_s / ns),
                label=f'S={S[i_s]:.2f}')
# Ideal line
ax.plot([0, 2], [0, 1], 'k--', lw=2, label='Ideal (linear)')
ax.set_xlabel(r'$\theta_H / \pi$ (Hamada)')
ax.set_ylabel('Cumulative area / total')
ax.set_title('Hamada: cumulative area is linear in $\\theta_H$')
ax.legend(fontsize=8, ncol=2)
ax.grid(True, alpha=0.3)

ax = axes[1]
cs = ax.contour(Rg, Zg, psi_norm_grid, levels=np.linspace(0, 1, 11),
                colors='lightgray', linewidths=0.5)
colors_h = cm.cool(np.linspace(0.2, 0.9, ns // stride_s + 1))
for k, i in enumerate(range(1, ns, stride_s)):
    ax.plot(R_H[i, :], Z_H[i, :], color=colors_h[k], lw=1.0)
for j in range(0, ntheta-1, ntheta // 8):
    ax.plot(R_H[1:, j], Z_H[1:, j], 'c-', lw=0.5, alpha=0.5)
ax.plot(Rmaxis, Zmaxis, 'r+', ms=10, mew=2)
ax.set_xlabel('R (m)')
ax.set_ylabel('Z (m)')
ax.set_title('Hamada mesh (R-Z view)')
ax.set_aspect('equal')

plt.tight_layout()
plt.savefig('hamada_coords.png', dpi=100, bbox_inches='tight')
plt.show()

Hamada mesh built: shape (40, 181)

第 7 節：比較パネル - 4 つの座標系#

ここでは、同じ磁束面群について 4 つのメッシュをすべて R-Z 平面に重ね描きします。

[9]:

fig = plt.figure(figsize=(16, 14))
gs = GridSpec(2, 2, figure=fig, hspace=0.3, wspace=0.3)

meshes = [
    ('PEST', TET, R_mesh, Z_mesh, 'viridis', 'g'),
    ('Equal-arc', TET_ea, R_ea, Z_ea, 'plasma', 'b'),
    ('Boozer', TET_B, R_B, Z_B, 'inferno', 'm'),
    ('Hamada', TET_H, R_H, Z_H, 'cool', 'c'),
]

stride_s = max(1, ns // 8)
stride_t = ntheta // 12

for idx, (name, tet, R_m, Z_m, cmap, pline_color) in enumerate(meshes):
    ax = fig.add_subplot(gs[idx // 2, idx % 2])

    # Flux surface contours as background
    ax.contour(Rg, Zg, psi_norm_grid, levels=np.linspace(0.05, 0.95, 10),
               colors='lightgray', linewidths=0.4)

    # Mesh surfaces (constant S)
    surface_colors = plt.get_cmap(cmap)(np.linspace(0.2, 0.9, ns // stride_s + 1))
    for k, i in enumerate(range(1, ns, stride_s)):
        ax.plot(R_m[i, :], Z_m[i, :], color=surface_colors[k], lw=1.2)

    # Poloidal lines (constant theta)
    for j in range(0, len(tet)-1, stride_t):
        ax.plot(R_m[1:, j], Z_m[1:, j], color=pline_color,
                lw=0.6, alpha=0.6)

    ax.plot(Rmaxis, Zmaxis, 'r+', ms=12, mew=2)
    ax.set_xlabel('R (m)', fontsize=11)
    ax.set_ylabel('Z (m)', fontsize=11)
    ax.set_title(f'{name} coordinates', fontsize=13, fontweight='bold')
    ax.set_aspect('equal')

    # Add annotation
    ann = {'PEST': 'Straight B field lines',
           'Equal-arc': 'Uniform arc length',
           'Boozer': 'Straight B + uniform Jacobian',
           'Hamada': 'Straight B + equal area'}[name]
    ax.text(0.05, 0.95, ann, transform=ax.transAxes, fontsize=9,
            verticalalignment='top', bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))

fig.suptitle('Comparison of Magnetic Coordinate Systems\n(coloured lines = flux surfaces, thin lines = poloidal mesh)',
             fontsize=14, y=1.01)
plt.savefig('coords_comparison.png', dpi=100, bbox_inches='tight')
plt.show()

第 8 節：物理的解釈#

8.1 Fourier モード結合#

直線磁力線座標（PEST、Boozer、Hamada）では、単一のトロイダルモード数 \(n\) を持つ摂動は、\(q = m/n\) となる面で共鳴します。\(\theta\) に関する Fourier 分解は明瞭なモード構造を持ちます。

一方、幾何角（または等弧長角）を使うと、単一の \((m, n)\) モードが 複数の Fourier 成分として現れ、異なる \(m\) 値を結合します。これはいわゆる Fourier 結合問題です。

[10]:

# Demonstrate: Fourier spectrum of R(theta) — 2D heatmap (m vs S) for PEST and equal-arc

m_show = 12   # show modes m=0..m_show-1

# Compute FFT of R(theta) - <R> at each surface for both coordinate systems
S_arr   = np.linspace(0.1, 0.9, ns)  # approx normalised flux label per surface
fft_PEST = np.zeros((ns, m_show))
fft_EA   = np.zeros((ns, m_show))

for i in range(ns):
    R_pest = R_mesh[i, :-1] - R_mesh[i, :-1].mean()
    R_ea_i = R_ea[i, :-1]   - R_ea[i, :-1].mean()
    n_pts_pest = len(R_pest)
    n_pts_ea   = len(R_ea_i)

    fft_p = np.abs(np.fft.rfft(R_pest))[:m_show]
    fft_p = fft_p / (fft_p.max() + 1e-30)
    fft_PEST[i, :len(fft_p)] = fft_p[:m_show]

    fft_e = np.abs(np.fft.rfft(R_ea_i))[:m_show]
    fft_e = fft_e / (fft_e.max() + 1e-30)
    fft_EA[i, :len(fft_e)] = fft_e[:m_show]

fig, axes = plt.subplots(1, 2, figsize=(13, 4.5))

for ax, data, name in zip(
    axes,
    [fft_PEST, fft_EA],
    ['PEST (straight field line)', 'Equal-arc θ'],
):
    # log scale heatmap: rows = radial surface index, cols = mode m
    log_data = np.log10(data + 1e-6)
    im = ax.imshow(
        log_data.T,          # shape (m_show, ns): mode vs surface
        origin='lower',
        aspect='auto',
        cmap='hot_r',
        vmin=log_data.max() - 3,   # show 3 decades
        vmax=log_data.max(),
        extent=[S_arr[0], S_arr[-1], -0.5, m_show - 0.5],
        interpolation='nearest',
    )
    cbar = fig.colorbar(im, ax=ax, pad=0.02, shrink=0.85)
    cbar.set_label(r'$\log_{10}$ normalised FFT amplitude', fontsize=8)
    ax.set_xlabel(r'Normalised flux label $S$', fontsize=10)
    ax.set_ylabel('Poloidal mode m', fontsize=10)
    ax.set_yticks(np.arange(0, m_show, 2))
    ax.set_title(f'R(θ) Fourier spectrum: {name}', fontsize=11)

# Annotation: PEST should be concentrated at m=1; equal-arc spreads
axes[0].text(0.05, 0.9, 'Spectrum concentrated\nat m=1 (ideal)',
             transform=axes[0].transAxes, fontsize=8, color='white',
             bbox=dict(boxstyle='round', facecolor='steelblue', alpha=0.5))
axes[1].text(0.05, 0.9, 'Spectrum spreads to\nhigher m (non-ideal)',
             transform=axes[1].transAxes, fontsize=8, color='white',
             bbox=dict(boxstyle='round', facecolor='tomato', alpha=0.5))

plt.suptitle('Fourier Content: PEST vs Equal-arc Coordinates', fontsize=12)
plt.tight_layout()
plt.show()

8.2 まとめ表#

[11]:

summary = """
╔══════════════╦═══════════════════════╦══════════════════════════════╦══════════════════════════╗
║ Coordinate   ║ Definition of θ       ║ Main advantage               ║ Used in                  ║
╠══════════════╬═══════════════════════╬══════════════════════════════╬══════════════════════════╣
║ PEST         ║ Straight B field lines║ Minimal Fourier coupling      ║ PEST, GATO, ELITE codes  ║
║              ║ B·∇θ*/B·∇φ = q(ψ)   ║ q = m/n at resonance         ║                          ║
╠══════════════╬═══════════════════════╬══════════════════════════════╬══════════════════════════╣
║ Boozer       ║ PEST + √g = √g(ψ)    ║ 1/B² drift is purely radial  ║ W7-X, LHD, VMEC output   ║
║              ║ (uniform Jacobian)    ║ Cleaner neoclassical theory   ║                          ║
╠══════════════╬═══════════════════════╬══════════════════════════════╬══════════════════════════╣
║ Hamada       ║ Equal area from axis  ║ J·∇θ = const, MHD stability  ║ TERPSICHORE, CASTOR      ║
║              ║ ∝ swept poloidal area ║ matrices simplified           ║                          ║
╠══════════════╬═══════════════════════╬══════════════════════════════╬══════════════════════════╣
║ Equal-arc    ║ Uniform arc length    ║ Simple construction           ║ Numerical grids, FEM     ║
║              ║ ds/dθ_ea = const     ║ Resolves boundary well        ║                          ║
╚══════════════╩═══════════════════════╩══════════════════════════════╩══════════════════════════╝
"""
print(summary)


╔══════════════╦═══════════════════════╦══════════════════════════════╦══════════════════════════╗
║ Coordinate   ║ Definition of θ       ║ Main advantage               ║ Used in                  ║
╠══════════════╬═══════════════════════╬══════════════════════════════╬══════════════════════════╣
║ PEST         ║ Straight B field lines║ Minimal Fourier coupling      ║ PEST, GATO, ELITE codes  ║
║              ║ B·∇θ*/B·∇φ = q(ψ)   ║ q = m/n at resonance         ║                          ║
╠══════════════╬═══════════════════════╬══════════════════════════════╬══════════════════════════╣
║ Boozer       ║ PEST + √g = √g(ψ)    ║ 1/B² drift is purely radial  ║ W7-X, LHD, VMEC output   ║
║              ║ (uniform Jacobian)    ║ Cleaner neoclassical theory   ║                          ║
╠══════════════╬═══════════════════════╬══════════════════════════════╬══════════════════════════╣
║ Hamada       ║ Equal area from axis  ║ J·∇θ = const, MHD stability  ║ TERPSICHORE, CASTOR      ║
║              ║ ∝ swept poloidal area ║ matrices simplified           ║                          ║
╠══════════════╬═══════════════════════╬══════════════════════════════╬══════════════════════════╣
║ Equal-arc    ║ Uniform arc length    ║ Simple construction           ║ Numerical grids, FEM     ║
║              ║ ds/dθ_ea = const     ║ Resolves boundary well        ║                          ║
╚══════════════╩═══════════════════════╩══════════════════════════════╩══════════════════════════╝

8.3 座標系間の関係#

4 つの系はすべて、\(\theta_{\rm new} = \theta^* + f(\psi, \theta^*)\) という形の角度変換で結ばれています。

Equal-arc → 等弧長パラメータ化（純粋に幾何的）
PEST → 等しい磁力線巻き付き（\(B_R, B_Z\) の積分が必要）
Boozer → PEST + ヤコビアンを平坦化する周期補正
Hamada → 等面積変換（囲まれた面積を使い、\(\sqrt{g}\) を全表面面積に比例させることと等価）

トカマクでは、Boozer と Hamada はしばしば非常によく似ます。どちらも磁束面平均された量（ヤコビアンと面積）が同じ圧力バランスに支配されるためです。

[12]:

# Final comparison: all four theta angles for a single field line
# Show how theta_PEST, theta_B, theta_H, theta_ea relate on the midradius surface
i_surf = ns // 2
print(f"Comparing poloidal angle representations on surface S={S[i_surf]:.3f}")
print(f"  Safety factor q = {q_iS[i_surf]:.3f}")

# For each coordinate system, compute the geometric angle (atan2(Z-Zax, R-Rax))
fig, ax = plt.subplots(figsize=(10, 5))

for name, tet_arr, R_m, Z_m, color in [
    ('PEST', TET, R_mesh, Z_mesh, 'blue'),
    ('Equal-arc', TET_ea, R_ea, Z_ea, 'green'),
    ('Boozer', TET_B, R_B, Z_B, 'red'),
    ('Hamada', TET_H, R_H, Z_H, 'purple'),
]:
    R_s = R_m[i_surf, :]
    Z_s = Z_m[i_surf, :]
    theta_geom = np.arctan2(Z_s - Zmaxis, R_s - Rmaxis) % (2 * np.pi)
    ax.plot(tet_arr / np.pi, theta_geom / np.pi, label=name, color=color, lw=1.5)

# Diagonal = geometric angle equals coordinate angle
ax.plot([0, 2], [0, 2], 'k--', lw=1, alpha=0.5, label='Identity (geom = coord)')
ax.set_xlabel(r'Coordinate angle $\theta / \pi$')
ax.set_ylabel(r'Geometric angle $\theta_{\rm geom} / \pi$')
ax.set_title('Geometric angle vs coordinate angle for each system\n' +
             f'(surface S={S[i_surf]:.2f}, q={q_iS[i_surf]:.2f})')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_xlim(0, 2)
ax.set_ylim(0, 2)

plt.tight_layout()
plt.savefig('theta_comparison.png', dpi=100, bbox_inches='tight')
plt.show()
print("All four systems differ only in how θ is distributed around the flux surface.")

Comparing poloidal angle representations on surface S=0.474
  Safety factor q = 1.638
All four systems differ only in how θ is distributed around the flux surface.