Übersetzung möglicherweise veraltet Die englische Quellseite ist neuer; diese deutsche Seite ist möglicherweise noch nicht vollständig synchronisiert.
Englische Quelle öffnen

Magnetische Koordinatensysteme in Tokamak-Gleichgewichten#

Dieses Tutorial vergleicht vier Flusskoordinatensysteme, die in der magnetischen Fusion verwendet werden: PEST, Boozer, Hamada und Equal-Arc.

Was sind Flusskoordinaten?#

In einem Tokamak liegen magnetische Feldlinien auf verschachtelten toroidalen Flächen, den Flussflächen. Ein Flusskoordinatensystem \((\psi, \theta, \varphi)\) erfüllt:

  • \(\psi\) beschriftet Flussflächen, z. B. \(\psi_{\rm norm} = 0\) auf der Achse und \(1\) an der LCFS

  • \(\varphi\) ist der übliche toroidale Winkel

  • \(\theta\) ist ein poloidaler Winkel mit nützlichen Eigenschaften

Die vier Systeme unterscheiden sich nur in der Wahl von \(\theta\):

System

Definition des poloidalen Winkels \(\theta\)

Haupteigenschaft

PEST

\(\partial(\mathbf{B}\cdot\nabla\theta^*) / \partial\theta^* = 0\): gerade Feldlinien

Einfachste Gerade-Feldlinien-Koordinaten

Boozer

PEST + Jacobi-Determinante ist Flussfunktion: \(\sqrt{g} = \sqrt{g}(\psi)\)

\(\mathbf{J}\cdot\nabla\varphi = {\rm const}(\psi)\); genutzt in Stellarator-/Tokamak-Codeausgaben

Hamada

Gleiche Fläche ab der magnetischen Achse

\(\mathbf{J}\cdot\nabla\theta = {\rm const}(\psi)\); vereinfacht MHD-Stabilitätsmatrizen

Equal-Arc

Bogenlänge entlang der Flussfläche ist gleichmäßig in \(\theta_{\rm ea}\)

Numerisch robust; einfache Konstruktion

Warum sind diese Wahlmöglichkeiten wichtig?#

  • Gerade Feldlinien (PEST, Boozer, Hamada): Ein Modus mit toroidaler Zahl \(n\) und poloidaler Zahl \(m\) erfüllt bei Resonanz \(q = m/n\). Das vereinfacht die Fourier-Analyse von MHD-Moden.

  • Boozer: Die \(1/B^2\)-Drift ist rein radial, was die neoklassische Transporttheorie und die kinetische Gleichung vereinfacht.

  • Hamada: Stromlinien sind ebenfalls gerade; das wird in einigen MHD-Stabilitätscodes benötigt.

  • Equal-Arc: Praktisch für numerische Gitter; löst die Plasmagrenze gut auf.

Einrichtung: analytisches Solov’ev-Gleichgewicht#

Wir verwenden das Solov’ev-Gleichgewicht nach Cerfon & Freidberg (2010), skaliert auf eine Referenz-Tokamakgröße (\(R_0 \approx 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 → referenz-tokamakgroß, R0≈1.86 m)
eq = solovev_iter_like(scale=0.3)
print(f"Gleichgewicht: 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"Magnetische Achse: R={Rmaxis:.3f} m, Z={Zmaxis:.3f} m")

Gleichgewicht: R0=1.86 m, a=0.60 m, B0=5.3 T
κ=1.70, δ=0.33, q0=1.50
Magnetische Achse: 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 Flächen, ntheta=181 poloidale Punkte
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-Gitter aufgebaut: {ns} Flächen, {ntheta} poloidale Punkte")
print(f"S-Bereich: [{S[1]:.3f}, {S[-1]:.3f}]")
print(f"q-Bereich: [{q_iS[1]:.2f}, {q_iS[-1]:.2f}]")

PEST-Gitter aufgebaut: 40 Flächen, 181 poloidale Punkte
S-Bereich: [0.028, 0.972]
q-Bereich: [1.49, 2.39]

Abschnitt 3: Equal-Arc-Koordinaten#

Der Equal-Arc-Winkel \(\theta_{\rm ea}\) parametrisiert jede Flussfläche so, dass die Bogenlänge entlang der Fläche proportional zu \(\theta_{\rm ea}\) ist. Dies ist die einfachste nichttriviale Koordinatentransformation.

[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-Gitter aufgebaut: 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"Bogenlängensegmente auf Fläche {i_mid}: Mittel={ds.mean()*100:.2f} cm, Std={ds.std()*100:.3f} cm")
print(f"  Gleichmäßigkeit (Std/Mittel): {ds.std()/ds.mean():.4f}")

Equal-Arc-Gitter aufgebaut: shape (40, 181)
Bogenlängensegmente auf Fläche 20: Mittel=1.27 cm, Std=0.000 cm
  Gleichmäßigkeit (Std/Mittel): 0.0001

[4]:

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

# Left: R-Z view of flux Flächen 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='Achse')
ax.set_xlabel('R (m)')
ax.set_ylabel('Z (m)')
ax.set_title('Equal-Arc-Gitter (R-Z-Ansicht)')
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 (gleichmäßig)')
ax.set_xlabel(r'$\theta_{\rm ea}$ (rad)')
ax.set_ylabel('Bogenlängeninkrement / Mittelwert')
ax.set_title('Bogenlängen-Gleichmäßigkeit in Equal-Arc-Koordinaten')
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: Bogenlängensegmente sind bis auf ~1 % gleichmäßig (Interpolationsdiskretisierung)")

Equal-Arc: Bogenlängensegmente sind bis auf ~1 % gleichmäßig (Interpolationsdiskretisierung)

Abschnitt 4: PEST-Koordinaten mit geraden Feldlinien#

In PEST-Koordinaten erfüllt der Sicherheitsfaktor \(q(\psi)\):

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

Das bedeutet, dass Feldlinien in der \((\theta^*, \varphi)\)-Ebene als Geraden mit Steigung \(1/q\) erscheinen.

[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='Analytisches q(S)')
ax.set_xlabel('S = √ψ_norm')
ax.set_ylabel('Sicherheitsfaktor q')
ax.set_title('Sicherheitsfaktorprofil: PEST gegen analytisch')
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-Gitter (R-Z-Ansicht)')
ax.set_aspect('equal')

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

[6]:

# Demonstrate straight Feldlinien 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('Feldlinien in PEST $(\\theta^*, \\varphi)$ space — Geraden mit Steigung $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()

Abschnitt 5: Boozer-Koordinaten#

Boozer-Koordinaten erweitern PEST durch die zusätzliche Forderung, dass die Jacobi-Determinante \(\sqrt{g}\) eine Flussflächen-Größe ist, also nur von \(\psi\) und nicht von \(\theta\) abhängt. Dadurch gilt \(\mathbf{J}\cdot\nabla\varphi = {\rm const}(\psi)\).

Konstruktion: \(\theta_B = \theta^* + \lambda(\psi, \theta^*)\) mit

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

wobei \(\langle\cdot\rangle\) den Flussflächen-Mittelwert \(\frac{1}{2\pi}\int_0^{2\pi}\cdot\,d\theta^*\) bezeichnet.

[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-Gitter aufgebaut: shape {R_B.shape}")
print(f"Max. |λ|-Korrektur: {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 Flächen
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-Winkelkorrektur $\\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-Gitter (R-Z-Ansicht)')
ax.set_aspect('equal')

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

Boozer-Gitter aufgebaut: shape (40, 181)
Max. |λ|-Korrektur: 0.5429 rad

Abschnitt 6: Hamada-Koordinaten#

Hamada-Koordinaten verlangen, dass sowohl Feldlinien als auch Stromdichtelinien gerade sind:

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

Für achsensymmetrische Gleichgewichte reduziert sich das darauf, den poloidalen Winkel flächenäquidistant zu machen: Die von der magnetischen Achse bis zum Winkel \(\theta_H\) überstrichene Fläche ist proportional zu \(\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-Gitter aufgebaut: 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 Flächen
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('Kumulative Fläche / Gesamtfläche')
ax.set_title('Hamada: kumulative Fläche ist 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-Gitter (R-Z-Ansicht)')
ax.set_aspect('equal')

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

Hamada-Gitter aufgebaut: shape (40, 181)

Abschnitt 7: Vergleichspanel - alle vier Koordinatensysteme#

Hier überlagern wir alle vier Gitter in der R-Z-Ebene für denselben Satz von Flussflächen.

[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 Flächen (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} Koordinaten', fontsize=13, fontweight='bold')
    ax.set_aspect('equal')

    # Add annotation
    ann = {'PEST': 'Gerade B-Feldlinien',
           'Equal-arc': 'Gleichmäßige Bogenlänge',
           '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('Vergleich magnetischer Koordinatensysteme\n(farbige Linien = Flussflächen, dünne Linien = poloidales Gitter)',
             fontsize=14, y=1.01)
plt.savefig('coords_comparison.png', dpi=100, bbox_inches='tight')
plt.show()

Abschnitt 8: Physikalische Interpretation#

8.1 Fourier-Modenkopplung#

In Gerade-Feldlinien-Koordinaten (PEST, Boozer, Hamada) resoniert eine Störung mit einer einzelnen toroidalen Modenzahl \(n\) auf der Fläche, auf der \(q = m/n\) gilt. Die Fourier-Zerlegung in \(\theta\) besitzt dort eine klare Modenstruktur.

Im Gegensatz dazu erscheint bei geometrischem Winkel oder Equal-Arc ein einzelner \((m,n)\)-Modus als mehrere Fourier-Komponenten und koppelt verschiedene \(m\)-Werte - das sogenannte Fourier-Kopplungsproblem.

[10]:

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

m_show = 12   # show Moduss 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 = Modus m
    log_data = np.log10(data + 1e-6)
    im = ax.imshow(
        log_data.T,          # shape (m_show, ns): Modus 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'Normalisiertes Flusslabel $S$', fontsize=10)
    ax.set_ylabel('Poloidaler Modus m', fontsize=10)
    ax.set_yticks(np.arange(0, m_show, 2))
    ax.set_title(f'R(θ) Fourier-Spektrum: {name}', fontsize=11)

# Annotation: PEST sollte sein 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-Inhalt: PEST gegen Equal-Arc-Koordinaten', fontsize=12)
plt.tight_layout()
plt.show()

8.2 Übersichtstabelle#

[11]:

summary = """
╔══════════════╦═══════════════════════╦══════════════════════════════╦══════════════════════════╗
║ Koordinate   ║ Definition von θ       ║ Hauptvorteil               ║ Verwendung                  ║
╠══════════════╬═══════════════════════╬══════════════════════════════╬══════════════════════════╣
║ PEST         ║ Gerade B-Feldlinien    ║ Minimale Fourier-Kopplung     ║ generische MHD-Codes     ║
║              ║ B·∇θ*/B·∇φ = q(ψ)   ║ q = m/n bei Resonanz         ║                          ║
╠══════════════╬═══════════════════════╬══════════════════════════════╬══════════════════════════╣
║ Boozer       ║ PEST + √g = √g(ψ)    ║ 1/B²-Drift ist rein radial   ║ Codeausgaben             ║
║              ║ (gleichmäßiger Jac.)   ║ klarere neoklassische Theorie ║                          ║
╠══════════════╬═══════════════════════╬══════════════════════════════╬══════════════════════════╣
║ Hamada       ║ gleiche Fläche ab Achse║ J·∇θ = const, MHD-Stabilität ║ Stabilitätscodes         ║
║              ║ ∝ überstrichene Fläche ║ Matrizen vereinfacht          ║                          ║
╠══════════════╬═══════════════════════╬══════════════════════════════╬══════════════════════════╣
║ Equal-arc    ║ gleichmäßige Bogenlänge║ Einfache Konstruktion         ║ numerische Gitter, FEM   ║
║              ║ ds/dθ_ea = const     ║ löst Rand gut auf             ║                          ║
╚══════════════╩═══════════════════════╩══════════════════════════════╩══════════════════════════╝
"""
print(summary)


╔══════════════╦═══════════════════════╦══════════════════════════════╦══════════════════════════╗
║ Koordinate   ║ Definition von θ       ║ Hauptvorteil               ║ Verwendung                  ║
╠══════════════╬═══════════════════════╬══════════════════════════════╬══════════════════════════╣
║ PEST         ║ Gerade B-Feldlinien    ║ Minimale Fourier-Kopplung     ║ generische MHD-Codes     ║
║              ║ B·∇θ*/B·∇φ = q(ψ)   ║ q = m/n bei Resonanz         ║                          ║
╠══════════════╬═══════════════════════╬══════════════════════════════╬══════════════════════════╣
║ Boozer       ║ PEST + √g = √g(ψ)    ║ 1/B²-Drift ist rein radial   ║ Codeausgaben             ║
║              ║ (gleichmäßiger Jac.)   ║ klarere neoklassische Theorie ║                          ║
╠══════════════╬═══════════════════════╬══════════════════════════════╬══════════════════════════╣
║ Hamada       ║ gleiche Fläche ab Achse║ J·∇θ = const, MHD-Stabilität ║ Stabilitätscodes         ║
║              ║ ∝ überstrichene Fläche ║ Matrizen vereinfacht          ║                          ║
╠══════════════╬═══════════════════════╬══════════════════════════════╬══════════════════════════╣
║ Equal-arc    ║ gleichmäßige Bogenlänge║ Einfache Konstruktion         ║ numerische Gitter, FEM   ║
║              ║ ds/dθ_ea = const     ║ löst Rand gut auf             ║                          ║
╚══════════════╩═══════════════════════╩══════════════════════════════╩══════════════════════════╝

8.3 Beziehung zwischen den Systemen#

Alle vier Systeme sind durch Winkeltransformationen der Form \(\theta_{\rm new} = \theta^* + f(\psi, \theta^*)\) verbunden:

  • Equal-Arc -> Parametrisierung mit gleicher Bogenlänge (rein geometrisch)

  • PEST -> gleiche Feldlinien-Windung (erfordert Integration von \(B_R, B_Z\))

  • Boozer -> PEST + periodische Korrektur zur Glättung der Jacobi-Determinante

  • Hamada -> Flächentransformation (bezieht die eingeschlossene Fläche ein und entspricht einer \(\sqrt{g}\)-Skalierung proportional zur Gesamtfläche)

Für einen Tokamak sind Boozer und Hamada oft sehr ähnlich, weil beide flussflächen-gemittelten Größen, Jacobi-Determinante und Fläche, durch dieselbe Druckbilanz gesteuert werden.

[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"Vergleiche Darstellungen des poloidalen Winkels auf Fläche S={S[i_surf]:.3f}")
print(f"  Sicherheitsfaktor 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='Identität (geom = coord)')
ax.set_xlabel(r'Koordinatenwinkel $\theta / \pi$')
ax.set_ylabel(r'Geometrischer Winkel $\theta_{\rm geom} / \pi$')
ax.set_title('Geometrischer Winkel gegen Koordinatenwinkel für jedes 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("Alle vier Systeme unterscheiden sich nur darin, wie θ um die Flussfläche verteilt ist.")

Vergleiche Darstellungen des poloidalen Winkels auf Fläche S=0.474
  Sicherheitsfaktor q = 1.638
Alle vier Systeme unterscheiden sich nur darin, wie θ um die Flussfläche verteilt ist.