토카막 평형의 자기 좌표계#

이 튜토리얼은 자기핵융합에서 쓰이는 네 가지 자속 좌표계를 비교합니다: PEST, Boozer, Hamada, Equal-arc.

자속 좌표란 무엇인가?#

토카막에서 자기력선은 자속면이라고 부르는 중첩된 토로이달 곡면 위에 놓입니다. 자속 좌표계 \((\psi, \theta, \varphi)\)는 다음 조건을 갖습니다.

  • \(\psi\)는 자속면을 라벨링한다. 예를 들어 \(\psi_{\rm norm}=0\)은 축, \(1\)은 LCFS이다.

  • \(\varphi\)는 표준 토로이달 각이다.

  • \(\theta\)는 유용한 성질을 갖도록 정의한 폴로이달 각이다.

네 좌표계는 \(\theta\)의 선택만 다릅니다.

좌표계

폴로이달 각 \(\theta\) 정의

핵심 성질

PEST

\(\partial(\mathbf{B}\cdot\nabla\theta^*) / \partial\theta^* = 0\): straight field lines

가장 단순한 직선 자기력선 좌표

Boozer

PEST + Jacobian이 자속 함수: \(\sqrt{g}=\sqrt{g}(\psi)\)

\(\mathbf{J}\cdot\nabla\varphi={\rm const}(\psi)\); 스텔러레이터와 신고전 계산에서 널리 쓰임

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\) drift가 순수하게 반경 방향이 되어 신고전 수송 이론과 운동 방정식이 단순해집니다.

  • Hamada: 전류 유선도 직선이 되므로 일부 MHD 안정성 코드에 필요합니다.

  • Equal-arc: 수치 격자에 실용적이며 플라즈마 경계를 잘 해상합니다.

설정: Solov’ev 해석 평형#

Cerfon & Freidberg (2010)의 Solov’ev 평형을 EAST 규모(R0 ≈ 1.86 m)로 조정해 사용합니다.

[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 straight-field-line 좌표#

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에 더해 Jacobian \(\sqrt{g}\)가 자속면 양, 즉 \(\psi\)에만 의존하고 \(\theta\)에는 의존하지 않을 것을 요구합니다. 이 조건은 \(\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. 비교 패널 – 네 좌표계#

여기서는 같은 자속면 집합에 대해 네 가지 격자를 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 분해가 깨끗한 모드 구조를 가집니다.

반대로 기하학적 각 또는 equal-arc 각을 쓰면, 단일 \((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 좌표계 사이의 관계#

네 좌표계는 모두 \(\theta_{\rm new}=\theta^*+f(\psi,\theta^*)\) 형태의 각도 변환으로 연결됩니다.

  • Equal-arc → 등호길이 매개변수화, 순수하게 기하학적인 변환

  • PEST → 자기력선 감김을 균일화하며, \(B_R, B_Z\) 적분이 필요

  • Boozer → PEST에 Jacobian을 평탄하게 만드는 주기적 보정을 더함

  • Hamada → 등면적 변환이며, 둘러싼 면적을 사용하고 \(\sqrt{g}\)를 전체 곡면 면적에 비례하게 만드는 것과 동등

토카막에서는 Boozer와 Hamada가 자주 매우 비슷합니다. 두 경우 모두 자속면 평균된 양, 즉 Jacobian과 면적이 같은 압력 평형에 의해 지배되기 때문입니다.

[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.