"""Plot and explicitly save neoclassical transport comparisons."""
from __future__ import annotations
from pathlib import Path
import numpy as np
from .neoclassical import E_CHARGE, NeoResult, PLOT_GROUPS, XGC_UNCERTAINTY
[docs]
def plot_neoclassical_transport(result: NeoResult, quantity: str, *, xlim=None, ylim=None, show_uncertainty: bool = True, uncertainty_sigma: float = 1.0):
"""Plot XGC and analytical curves with optional temporal-spread bands.
The band is the selected multiple of the standard deviation across time
samples, as in the IDL plots. It is not a standard error or a confidence
interval for the time-averaged value. The bootstrap band covers variation
in XGC current only, not uncertainty in the modeled Ohmic subtraction.
Parameters
----------
result : NeoResult
Computed transport comparison.
quantity : str
One key from the plotted transport quantities.
xlim, ylim : pair, optional
Axis bounds; ``None`` leaves Matplotlib autoscaling in effect.
show_uncertainty : bool
Whether to shade temporal spread of the XGC curve.
uncertainty_sigma : float
Number of temporal standard deviations in the shaded band.
Returns
-------
matplotlib.figure.Figure
Figure ready for display or SVG export.
"""
from matplotlib import pyplot as plt
if quantity not in PLOT_GROUPS:
raise KeyError(f"Unknown neoclassical plot quantity: {quantity}")
if not np.isfinite(uncertainty_sigma) or uncertainty_sigma < 0:
raise ValueError("Uncertainty width must be a finite nonnegative sigma multiple")
title, unit, variables = PLOT_GROUPS[quantity]
fig, ax = plt.subplots(figsize=(8.5, 6.375))
psi = result.profiles["psi_pol_norm"]
labels = {"jpar_plasma": "XGC", "jpar_efit": "EFIT", "jboot_sauter": "Sauter", "sauter_total": "Sauter + Ohmic", "xgc_bootstrap_estimate": "XGC bootstrap estimate", "hh_bootstrap": "Hinton–Hazeltine bootstrap", "jpar_plasma_hh": "Hinton–Hazeltine", "ptl_flux_HH": "Hinton–Hazeltine", "heatflux_i_CH": "Chang–Hinton", "heatflux_e_HH": "Hinton–Hazeltine", "energy_flux_hh": "Hinton–Hazeltine", "D_HH": "Hinton–Hazeltine", "chi_i_ch": "Chang–Hinton", "chi_e_hh": "Hinton–Hazeltine"}
derived = {
"sauter_total": result.profiles["jboot_sauter"] + result.profiles["jpar_ohmic"],
"xgc_bootstrap_estimate": result.profiles["jpar_plasma"] - result.profiles["jpar_ohmic"],
"hh_bootstrap": result.profiles["jpar_plasma_hh"] - result.profiles["jpar_ohmic_hh"],
"energy_flux_hh": result.profiles["heatflux_e_HH"] + 2.5*E_CHARGE*result.profiles["T_e"]*result.profiles["ptl_flux_HH"],
}
for name in variables:
if name in result.profiles or name in derived:
values = result.profiles[name] if name in result.profiles else derived[name]
line, = ax.plot(psi, values, label=labels.get(name, "XGC"), linewidth=2)
sigma_name = XGC_UNCERTAINTY.get(name)
if show_uncertainty and uncertainty_sigma > 0 and result.sample_count > 1 and sigma_name in result.profiles:
width = uncertainty_sigma*np.asarray(result.profiles[sigma_name])
ax.fill_between(psi, values-width, values+width, color=line.get_color(), alpha=0.18, linewidth=0)
average_label = " | line FSA" if result.fsa == "line" else ""
ax.set(xlabel=r"$\psi_N$", ylabel=unit, title=f"{title} | gsteps {result.gstep_range[0]}–{result.gstep_range[1]}{average_label}")
if xlim is not None:
ax.set_xlim(*xlim)
if ylim is not None:
ax.set_ylim(*ylim)
ax.grid(alpha=.25)
ax.legend()
fig.tight_layout()
return fig
[docs]
def save_neoclassical_svg_plots(result: NeoResult, directory: str | Path, *, quantities=None, show_uncertainty: bool = True, uncertainty_sigma: float = 1.0) -> list[Path]:
"""Write selected comparison plots as standalone SVG files.
Parameters
----------
result
Previously computed analysis result.
directory
Directory where SVG files are written.
quantities
Optional subset of plot quantity names; all are used by default.
show_uncertainty
Shade the temporal spread when plotting.
uncertainty_sigma
Width of the shaded band in standard deviations.
Returns
-------
list[Path]
paths of written SVG plots.
"""
from matplotlib import pyplot as plt
directory = Path(directory)
directory.mkdir(parents=True, exist_ok=True)
output = []
for quantity in quantities or PLOT_GROUPS:
fig = plot_neoclassical_transport(result, quantity, show_uncertainty=show_uncertainty, uncertainty_sigma=uncertainty_sigma)
path = directory / f"neoclassical_{quantity}.svg"
fig.savefig(path, format="svg")
plt.close(fig)
output.append(path)
return output