diff --git a/.gitattributes b/.gitattributes new file mode 100644 index 0000000..3695e37 --- /dev/null +++ b/.gitattributes @@ -0,0 +1,32 @@ +articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/d2cross_Ee0100eV-10MeV-80log_Eph100eV-10MeV-81log_ntheta61_EH.npz filter=lfs diff=lfs merge=lfs -text +articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/d2cross_Ee01eV-100MeV-80log_Eph1eV-100MeV-81log_ntheta61_EH.npz filter=lfs diff=lfs merge=lfs -text +articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/d2cross_phi_Ee01eV-100MeV-80log_Eph1eV-100MeV-81log_nthetaph61_nthetae060_EH.npz filter=lfs diff=lfs merge=lfs -text +articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/responsivities.npz filter=lfs diff=lfs merge=lfs -text +articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/d2cross_phi_Ee01eV-100MeV-240log_Eph1eV-100MeV-241log_nthetaph61_nthetae060_EH.npz filter=lfs diff=lfs merge=lfs -text +articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/260616_SCRAM86_Ar_data.npz filter=lfs diff=lfs merge=lfs -text +articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/260616_SCRAM86_Fe_data.npz filter=lfs diff=lfs merge=lfs -text +articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/260616_SCRAM86_Kr_data.npz filter=lfs diff=lfs merge=lfs -text +articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/260616_SCRAM86_O_data.npz filter=lfs diff=lfs merge=lfs -text +articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/260616_SCRAM86_W_data.npz filter=lfs diff=lfs merge=lfs -text +articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/260625_FLYCHK_H_data.npz filter=lfs diff=lfs merge=lfs -text +articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/260625_FLYCHK_He_data.npz filter=lfs diff=lfs merge=lfs -text +articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/Inversion_HXRVA_ptcam_dvezinet_20260709-204104.npz filter=lfs diff=lfs merge=lfs -text +articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/TFG_Config_ExpSPARC_SPARC-V2_sh00000_Vers1.8.18.npz filter=lfs diff=lfs merge=lfs -text +articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/Inversion_BOLO_DMS_ptcam_dvezinet_20260709-205535.npz filter=lfs diff=lfs merge=lfs -text +articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/HXR_MPP200.png filter=lfs diff=lfs merge=lfs -text +articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/HXR_MPP200_zoom.png filter=lfs diff=lfs merge=lfs -text +articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig03_validate_ff.py filter=lfs diff=lfs merge=lfs -text +articles/RSI_2026_RunawayBremsstrahlungDetection/figures/fig02_dist_type.png filter=lfs diff=lfs merge=lfs -text +articles/RSI_2026_RunawayBremsstrahlungDetection/figures/fig02_RE_ff.png filter=lfs diff=lfs merge=lfs -text +articles/RSI_2026_RunawayBremsstrahlungDetection/figures/fig03_validate_ff.png filter=lfs diff=lfs merge=lfs -text +articles/RSI_2026_RunawayBremsstrahlungDetection/figures/fig08_perfs_single.png filter=lfs diff=lfs merge=lfs -text +articles/RSI_2026_RunawayBremsstrahlungDetection/figures/fig07_responsivities.png filter=lfs diff=lfs merge=lfs -text +articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/260715_SCRAM86_H_data.npz filter=lfs diff=lfs merge=lfs -text +articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/260715_SCRAM86_He_data.npz filter=lfs diff=lfs merge=lfs -text +articles/RSI_2026_RunawayBremsstrahlungDetection/figures/HXR_MPP200_zoom.png filter=lfs diff=lfs merge=lfs -text +articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/OMP000U.png filter=lfs diff=lfs merge=lfs -text +articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/OMP000U_Zoom.png filter=lfs diff=lfs merge=lfs -text +articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/OMP000U_Zoom_inside.png filter=lfs diff=lfs merge=lfs -text +articles/RSI_2026_RunawayBremsstrahlungDetection/figures/fig14_bolo_cad.png filter=lfs diff=lfs merge=lfs -text +articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/Inversion_HXRVA_dvezinet_20260713-154229.npz filter=lfs diff=lfs merge=lfs -text +articles/RSI_2026_RunawayBremsstrahlungDetection/figures/fig11_hxrva_cad.png filter=lfs diff=lfs merge=lfs -text diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/.gitattributes b/articles/RSI_2026_RunawayBremsstrahlungDetection/.gitattributes new file mode 100644 index 0000000..75ea5b1 --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/.gitattributes @@ -0,0 +1,5 @@ +figures/fig01_cross.png filter=lfs diff=lfs merge=lfs -text +figures/fig02_dist.png filter=lfs diff=lfs merge=lfs -text +figures/fig04_bremsstrahlung.png filter=lfs diff=lfs merge=lfs -text +figures/fig05_tokamak.png filter=lfs diff=lfs merge=lfs -text +figures/fig06_responsivities.png filter=lfs diff=lfs merge=lfs -text diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/__init__.py b/articles/RSI_2026_RunawayBremsstrahlungDetection/__init__.py new file mode 100644 index 0000000..c4f97d0 --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/__init__.py @@ -0,0 +1 @@ +from .code import * diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/code/__init__.py b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/__init__.py new file mode 100644 index 0000000..0528c14 --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/__init__.py @@ -0,0 +1,21 @@ +from ._main import main +from ._load_spect import main as load_spect +from ._load_spect_anis import main as load_spect_anis +from ._perfs import main as perfs +from ._singlept_sensor import main as singlept_sensor + +from ._fig01_cross import main as fig01_cross +from ._fig02_dist import main as fig02_dist +from ._fig02_dist_type import main as fig02_dist_type +from ._fig02_RE_ff import main as fig02_RE_ff +from ._fig03_validate_ff import main as fig03_validate_ff +from ._fig04_bremsstrahlung import main as fig04_bremsstrahlung +from ._fig05_emiss import main as fig05_emiss +from ._fig06_tokamak import main as fig06_tokamak +from ._fig07_responsivities import main as fig07_responsivities +from ._fig08_perfs_single import main as fig08_perfs_single +from ._fig09_perfs_scans import main as fig09_perfs_scans +from ._fig11_hxrva_cad import main as fig11_hxrva_cad +from ._fig12_hxrva_images import main as fig12_hxrva_images +from ._fig13_hxrva_metrics import main as fig13_hxrva_metrics +from ._fig14_bolo_cad import main as fig14_bolo_cad diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig01_cross.py b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig01_cross.py new file mode 100644 index 0000000..fb132ed --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig01_cross.py @@ -0,0 +1,315 @@ + + +import os + + +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.gridspec as gridspec +import tofu as tf + + +from ._savefig import main as savefig + + +tfphysemis = tf.physics_tools.electrons.emission + + +# ##################################################### +# ##################################################### +# DEFAULTS +# ##################################################### + + +_PATH_HERE = os.path.dirname(__file__) +_PATH_PAPER = os.path.dirname(_PATH_HERE) + + +_DPFE_DCROSS = { + 'EH0': os.path.join( + _PATH_PAPER, + 'inputs', + 'd2cross_Ee01eV-100MeV-80log_Eph1eV-100MeV-81log_ntheta61_EH.npz', + ), + 'EH1': os.path.join( + _PATH_PAPER, + 'inputs', + 'd2cross_Ee0100eV-10MeV-80log_Eph100eV-10MeV-81log_ntheta61_EH.npz', + ), + 'BHE': os.path.join( + _PATH_PAPER, + 'inputs', + 'd2cross_Ee01eV-100MeV-400log_Eph1eV-100MeV-401log_ntheta181_BHE.npz', + ), +} + + +# ##################################################### +# ##################################################### +# Fig 1 - cross-section +# ##################################################### + + +def main( + figsize=(12, 6), + pfe_cross='EH0', + version='EH', + Eph_eV=np.r_[1e3, 10e3, 500e3], + Ee0_eV=np.r_[20e3, 1e6], + fontsize=14, + # save + path_save=None, + pfe_save=None, + # unused + **kwdargs, +): + """ Plot EH cross-section over a wide range of Eph and Ee0 + + Add 4 hand-picked cases vs angle of emission theta_ph + + """ + + # -------------- + # load + # -------------- + + pfe = _DPFE_DCROSS[pfe_cross] + + dout = { + kk: vv.tolist() + for kk, vv in dict(np.load(pfe, allow_pickle=True)).items() + } + units = dout['cross'][version]['units'] + Z = dout.get('Z', {'data': 1})['data'] + + # -------------- + # prepare axes + # -------------- + + dmargin_theta = { + 'left': 0.07, 'right': 0.46, + 'bottom': 0.09, 'top': 0.92, + 'wspace': 0.08, 'hspace': 0.05, + } + dmargin_map = { + 'left': 0.54, 'right': 0.99, + 'bottom': 0.09, 'top': 0.92, + 'wspace': 0.10, 'hspace': 0.10, + } + + fig = plt.figure(figsize=figsize) + + gs_theta = gridspec.GridSpec(ncols=2, nrows=2, **dmargin_theta) + gs_map = gridspec.GridSpec(ncols=1, nrows=1, **dmargin_map) + dax = {} + + # -------------- + # prepare axes + # -------------- + + # -------------- + # ax - isolines + + ax = fig.add_subplot( + gs_map[0, 0], + xscale='log', + yscale='log', + aspect='equal', + ) + ax.set_xlabel( + r"$E_{e,0}$ (keV)", + size=fontsize, + fontweight='bold', + ) + ax.set_ylabel( + r"$E_{ph}$ (keV)", + size=fontsize, + fontweight='bold', + ) + ax.set_title( + r"$d^2\sigma(E_{e0}, E_{ph}, \theta_{ph}, Z)$" + + f"\n Z = {Z} - version = {version}", + size=fontsize, + fontweight='bold', + ) + + # store + dax['map'] = {'handle': ax, 'type': 'isolines'} + + # -------------- + # ax - theta_norm + + # theta_norm0 + ax = fig.add_subplot( + gs_theta[0, 0], + xscale='linear', + ) + ax.set_ylabel( + "normalized cross-section (adim.)", + size=fontsize, + fontweight='bold', + ) + ax.set_title( + r"$E_{e0}$" + f" = {Ee0_eV[0]*1e-3:2.0f} keV", + size=fontsize, + fontweight='bold', + ) + ax.tick_params(labelbottom=False) + + # store + dax['theta_norm0'] = {'handle': ax, 'type': 'isolines'} + + # theta_norm1 + ax = fig.add_subplot( + gs_theta[0, 1], + sharex=dax['theta_norm0']['handle'], + sharey=dax['theta_norm0']['handle'], + ) + ax.set_title( + r"$E_{e0}$" + f" = {Ee0_eV[1]*1e-6:2.0f} MeV", + size=fontsize, + fontweight='bold', + ) + ax.tick_params(labelbottom=False) + ax.tick_params(labelleft=False) + ax.set_xticks(np.r_[0, 45, 90, 135, 180]) + + # store + dax['theta_norm1'] = {'handle': ax, 'type': 'isolines'} + + # -------------- + # ax - theta_abs + + # theta_abs0 + ax = fig.add_subplot( + gs_theta[1, 0], + sharex=dax['theta_norm0']['handle'], + ) + ax.set_xlabel( + r"$\theta_{ph}$ (deg)", + size=fontsize, + fontweight='bold', + ) + ax.set_ylabel( + r"$d\sigma$" + f"({units})", + size=fontsize, + fontweight='bold', + ) + ax.set_yticks(np.logspace(-12, -5, 8)) + + # store + dax['theta_abs0'] = {'handle': ax, 'type': 'isolines'} + + # theta_abs1 + ax = fig.add_subplot( + gs_theta[1, 1], + sharex=dax['theta_norm0']['handle'], + sharey=dax['theta_abs0']['handle'], + ) + ax.set_xlabel( + r"$\theta_{ph}$ (deg)", + size=fontsize, + fontweight='bold', + ) + ax.tick_params(labelleft=False) + + # store + dax['theta_abs1'] = {'handle': ax, 'type': 'isolines'} + + # ticklabels size + for kax, vax in dax.items(): + dax[kax]['handle'].tick_params( + axis='both', + which='major', + labelsize=fontsize - 1, + ) + + # ------------------ + # call built-in + # ------------------ + + # cases 100 keV + lc = ['r', 'g', 'b'] + for i0, e0 in enumerate(Ee0_eV): + iphok = Eph_eV < e0 + dcases = { + i1: { + 'E_e0_eV': e0, + 'E_ph_eV': eph, + 'color': lc[i1], + 'label': f"Eph = {eph*1e-3:3.0f} keV", + } + for i1, eph in enumerate(Eph_eV[iphok]) + } + _dax, _d2cross = tfphysemis.plot_xray_thin_d2cross_ei_anisotropy( + d2cross=pfe, + dcases=dcases, + dax={ + 'map': dax['map']['handle'], + 'theta_norm': dax[f'theta_norm{i0}']['handle'], + 'theta_abs': dax[f'theta_abs{i0}']['handle'], + }, + ) + + # remove contour plots + if i0 == 0: + for cc in dax['map']['handle'].get_children(): + if cc.__class__.__name__ == 'QuadContourSet': + cc.remove() + + # remove legend + dax['theta_norm0']['handle'].get_legend().remove() + + # --------------------- + # Adjust map x/y scales + # --------------------- + + dax['map']['handle'].set_xlim(0.1, 10e3) + dax['map']['handle'].set_ylim(0.1, 10e3) + dax['theta_abs1']['handle'].set_ylabel('') + dax['theta_norm1']['handle'].legend(loc='lower right') + + # -------------- + # add a, b, c, d, e + # -------------- + + dabc = { + 'theta_norm0': '(a)', + 'theta_norm1': '(c)', + 'theta_abs0': '(b)', + 'theta_abs1': '(d)', + } + for kax, abc in dabc.items(): + dax[kax]['handle'].grid(visible=True, which='major', axis='both') + dax[kax]['handle'].text( + 0.95, 0.95, + abc, + fontsize=fontsize, + fontweight='bold', + horizontalalignment='right', + verticalalignment='top', + transform=dax[kax]['handle'].transAxes, + ) + + dax['map']['handle'].text( + 0., 1.02, + "(e)", + fontsize=fontsize, + fontweight='bold', + horizontalalignment='left', + verticalalignment='bottom', + transform=dax['map']['handle'].transAxes, + ) + + # -------------- + # save + # -------------- + + savefig( + fig=fig, + pfe_save=pfe_save, + path_save=path_save, + file=__file__, + ) + + return dax diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig02_RE_ff.py b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig02_RE_ff.py new file mode 100644 index 0000000..7b49e58 --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig02_RE_ff.py @@ -0,0 +1,405 @@ +import os + + +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.gridspec as gridspec +import datastock as ds + + +from . import _load_spect_anis +from ._fig02_dist_type import _DDIST, _DDIST_PLOT +from ._savefig import main as savefig + + +# ##################################################### +# ##################################################### +# DEFAULTS +# ##################################################### + + +# PATHS +_PATH_HERE = os.path.dirname(__file__) +_PATH_PAPER = os.path.dirname(_PATH_HERE) + + +# ##################################################### +# ##################################################### +# Main +# ##################################################### + + +def main( + dmix='H', + # plotting temp + Te_eV=1e3, + # d2cross + d2cross_phi=None, + # Eph + Eph=np.r_[0.05, 5, 50]*1e3, + # plot + figsize=(5, 7), + fontsize=12, + # save + path_save=None, + pfe_save=None, + # unused + **kwdargs, +): + """ Validate Bremsstrahlung vs SCRAM and FLYCHK + + """ + + # -------------- + # inputs + # -------------- + + # -------------- + # load elements + # -------------- + + demiss = {} + for kdist, vdist in _DDIST['RE'].items(): + + kwd = { + kk: vv for kk, vv in _DDIST['maxwell'].items() + if 'Te' not in kk + } + kwd.update(**vdist) + # kwd.update(**_DDIST['coords']) + + demiss[kdist], ddist, dmix, data_source = _load_spect_anis.main( + dmix=dmix, + # d2cross + d2cross_phi=d2cross_phi, + # dist + **kwd, + ) + + # -------------- + # extract + # -------------- + + # Te, ne, units + kdomref = list(_DDIST['RE'].keys())[0] + Teu = np.unique(ddist['plasma']['Te_eV']['data']) + units = demiss[kdomref]['emiss']['maxwell']['ff']['units'] + + indTe = np.argmin(np.abs(Teu - Te_eV)) + Te_eV = ddist['plasma']['Te_eV']['data'][(0, indTe, 0)] + ne = np.unique(ddist['plasma']['ne_m3']['data'])[0] + jp = np.unique(ddist['plasma']['jp_Am2']['data'])[0] + jp_frac = np.unique(ddist['plasma']['jp_fraction_re']['data'])[0] + + nE = Eph.size + + # vmax, vmin + vmax = np.nanmax(demiss[kdomref]['emiss']['maxwell']['ff']['data']) + vmax_log10 = np.ceil(np.log10(vmax)) + vmax = 10**vmax_log10 + vmin = 10**(vmax_log10 - 21) + + # title + tit = ( + r"$n_e$" + f" = {ne:1.0e}" + r"$/m^3$, " + + r"$j_P$" + f" = {jp*1e-6:1.0f}" + r"$MA/m^2$" + "\n" + + r"$T_e$" + f" = {Te_eV*1e-3:1.0f} keV, " + + r"$F_{RE}$" + f" = {jp_frac:2.1f}" + ) + + # -------------- + # prepare axes + # -------------- + + dmargin_spect = { + 'left': 0.13, 'right': 0.97, + 'bottom': 0.06, 'top': 0.94, + 'wspace': 0.25, 'hspace': 0.20, + } + dmargin_theta = { + 'left': 0.13, 'right': 0.97, + 'bottom': 0.07, 'top': 0.55, + 'wspace': 0.15, 'hspace': 0.10, + } + + fig = plt.figure(figsize=figsize) + + gs_spect = gridspec.GridSpec(ncols=1, nrows=3, **dmargin_spect) + gs_theta = gridspec.GridSpec(ncols=nE, nrows=2, **dmargin_theta) + dax = {} + + # ---------------- + # ax - spectra + # ---------------- + + ax = fig.add_subplot( + gs_spect[0, :], + aspect='auto', + ) + ax.set_ylabel( + r"$\epsilon_{ff}$" + f' ({units})', + fontsize=fontsize, + fontweight='bold', + ) + ax.set_xlabel( + r"$E_{ph}$" + ' (keV)', + fontsize=fontsize, + fontweight='bold', + ) + ax.set_title( + tit, + fontsize=fontsize, + fontweight='bold', + ) + ax.text( + 0.95, + 0.95, + '(a)', + horizontalalignment='right', + verticalalignment='top', + transform=ax.transAxes, + ) + + dax['spectra'] = ax + + # ---------------- + # ax - theta_rel + # ---------------- + + ax0_rel = None + for ie, ee in enumerate(Eph): + + if ee < 1e3: + estr = f"{ee:3.0f} eV" + elif ee < 1e6: + estr = f"{ee*1e-3:3.0f} keV" + else: + estr = f"{ee*1e-6:3.0f} MeV" + + ax = fig.add_subplot( + gs_theta[0, ie], + aspect='auto', + sharex=ax0_rel, + sharey=ax0_rel, + ) + ax.set_title( + r"$E_{ph}$" + f" = {estr}", + fontsize=fontsize, + fontweight='bold', + ) + if ie == 0: + ax0_rel = ax + ax.set_xlim(0, 180) + ax.set_ylim(0, 1) + ax.set_xticks([0, 45, 90, 135, 180]) + ax.set_ylabel( + r"$\epsilon_{ff} / max(\epsilon_{ff})$", + fontsize=fontsize, + fontweight='bold', + ) + else: + ax.tick_params(labelleft=False) + + ax.text( + 0.95, + 0.95, + ['(b)', '(c)', '(d)'][ie], + horizontalalignment='right', + verticalalignment='top', + transform=ax.transAxes, + ) + ax.tick_params(labelbottom=False) + dax[f'theta_rel_{ie}'] = ax + + # ---------------- + # ax - theta_abs + # ---------------- + + ax0_abs = None + for ie, ee in enumerate(Eph): + + ax = fig.add_subplot( + gs_theta[1, ie], + aspect='auto', + sharex=ax0_rel, + sharey=ax0_abs, + ) + ax.set_xlabel( + r"$\theta_{ph,B}$" + " (deg)", + fontsize=fontsize, + fontweight='bold', + ) + if ie == 0: + ax0_abs = ax + ax.set_ylabel( + r"$\epsilon_{ff}$" + f' ({units})', + fontsize=fontsize, + fontweight='bold', + ) + ax.set_ylim(1e6, 1e13) + else: + ax.tick_params(labelleft=False) + + ax.text( + 0.95, + 0.95, + ['(e)', '(f)', '(g)'][ie], + horizontalalignment='right', + verticalalignment='top', + transform=ax.transAxes, + ) + dax[f'theta_abs_{ie}'] = ax + + # ---------------- + # check dax format + + dax = ds._generic_check._check_dax(dax) + + # ticklabels size + for kax, vax in dax.items(): + dax[kax]['handle'].tick_params( + axis='both', + which='major', + labelsize=fontsize - 1, + ) + + # -------------- + # plot - spectra + # -------------- + + kax = 'spectra' + if dax.get(kax) is not None: + ax = dax[kax]['handle'] + + # slice + sli = (0, indTe, 0, slice(None), slice(None)) + + # ---------- + # Maxwell + + data = demiss[kdomref]['emiss']['maxwell']['ff']['data'][sli] + ax.loglog( + demiss[kdomref]['E_ph']['data']*1e-3, + np.mean(data, axis=-1), + c='k', + ls='-', + label='Maxwellian', + ) + + # ---------- + # RE + + for kdist in demiss.keys(): + emiss_E = demiss[kdist]['emiss']['RE']['ff']['data'][sli] + + ax.fill_between( + demiss[kdomref]['E_ph']['data']*1e-3, + np.nanmin(emiss_E, axis=-1), + np.nanmax(emiss_E, axis=-1), + hatch=_DDIST_PLOT[kdist]['hatch'], + facecolor='None', + alpha=0.5, + edgecolor=_DDIST_PLOT[kdist]['color'], + ls='-', + label=kdist, + ) + + # -------- + # vlines + + for ie, ee in enumerate(Eph): + ax.axvline( + ee*1e-3, + c='k', + ls='--', + lw=1, + # label=f"E_ph = {ee*1e-3:3.1f} keV", + ) + + ax.set_xscale('log') + ax.set_yscale('log') + ax.set_ylim(vmin, vmax) + ax.set_xlim(1e-3, 1e5) + ax.set_ylim(1e-3, 1e18) + ax.set_yticks(np.logspace(-3, 18, 8)) + ax.grid(True) + ax.legend(loc='lower left', fontsize=fontsize - 2) + + # -------------- + # plot vs theta + # -------------- + + for ie, ee in enumerate(Eph): + + # indTe + indTe = np.argmin(np.abs(Teu - Te_eV)) + Te_eV = ddist['plasma']['Te_eV']['data'][(0, indTe, 0)] + + # indE + indE = np.argmin(np.abs(demiss[kdomref]['E_ph']['data'] - ee)) + ee = demiss[kdomref]['E_ph']['data'][indE] + sli = (0, indTe, 0, indE, slice(None)) + + # ----------- + # theta_rel + + kax = f"theta_rel_{ie}" + if dax.get(kax) is not None: + ax = dax[kax]['handle'] + + # plot + for kdist in demiss.keys(): + + # data + emiss_RE = demiss[kdist]['emiss']['RE']['ff']['data'][sli] + + # plot + ax.plot( + demiss[kdomref]['theta_ph_vsB']['data']*180/np.pi, + emiss_RE / np.nanmax(emiss_RE), + c=_DDIST_PLOT[kdist]['color'], + ls='-', + lw=1, + label=kdist, + ) + ax.grid(True) + + # ----------- + # theta_abs + + kax = f"theta_abs_{ie}" + if dax.get(kax) is not None: + ax = dax[kax]['handle'] + + # plot + for kdist in demiss.keys(): + + # data + emiss_RE = demiss[kdist]['emiss']['RE']['ff']['data'][sli] + + # plot + ax.semilogy( + demiss[kdomref]['theta_ph_vsB']['data']*180/np.pi, + emiss_RE, + c=_DDIST_PLOT[kdist]['color'], + ls='-', + lw=1, + label=kdist, + ) + + if ie == 0: + ax.set_yticks(np.logspace(6, 13, 6)) + ax.grid(True) + + # -------------- + # save + # -------------- + + savefig( + fig=fig, + pfe_save=pfe_save, + path_save=path_save, + file=__file__, + ) + + return dax, demiss diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig02_dist.py b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig02_dist.py new file mode 100644 index 0000000..6fdcb01 --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig02_dist.py @@ -0,0 +1,410 @@ +import os + +import numpy as np +import scipy.integrate as scpinteg +import astropy.units as asunits +import matplotlib.pyplot as plt +import matplotlib.gridspec as gridspec +import datastock as ds +import tofu as tf + + +from ._savefig import main as savefig + + +tfphysdist = tf.physics_tools.electrons.distribution + + +# ##################################################### +# ##################################################### +# DEFAULTS +# ##################################################### + + +_PATH_HERE = os.path.dirname(__file__) + + +_DDIST = { + # maxwell + 'Te_eV': np.r_[0.1e3, 2e3, 0.1e3, 2e3], + 'ne_m3': 1e19, + 'jp_Am2': 1e6, + # RE + 'jp_fraction_re': np.r_[0.1, 0.1, 0.9, 0.9], + 'dominant': 'bump', + 'pnormW': np.r_[0.1, 4, 0.1, 4], + 'Ekin_max_eV': np.r_[100e3, 10e6, 100e3, 10e6], + 'Ekin_min_eV': 100, + 'step': 0.1, + 'theta_width': 20*np.pi/180, + # coords + 'E_eV': np.logspace(0, 8, 80), + 'theta': np.linspace(0, 180, 181) * np.pi / 180, +} + + +# ##################################################### +# ##################################################### +# main +# ##################################################### + + +def main( + # coords + E_eV=None, + theta=None, + # Maxwell + ne_m3=None, + jp_Am2=None, + # RE + dominant=None, + jp_fraction_re=None, + Ekin_max_eV=None, + Ekin_min_eV=None, + step=None, + pnormW=None, + # plot + figsize=(5, 7), + fontsize=12, + # save + path_save=None, + pfe_save=None, + # unused + **kwdargs, +): + """ Plot a Maxwellian + RE distribution + + Includes: + - a 2d (E, pitch) contour plot + - a 1d (E,) plot + + """ + + # ------------ + # inputs + # ------------ + + din = locals() + din = { + kk: vv if din.get(kk) is None else din[kk] + for kk, vv in _DDIST.items() + } + + # ------------ + # compute + # ------------ + + # ddist = {'dist': dict, 'plasma': dist, 'coords': dist} + ddist = tfphysdist.get_distribution(**din) + + # units + units2d = asunits.Unit(ddist['dist']['RE']['dist']['units']) + units1d = units2d * asunits.Unit(ddist['coords']['x1']['units']) + + # ------------ + # Derive 1d data + # ------------ + + dataRE = scpinteg.trapezoid( + ddist['dist']['RE']['dist']['data'], + x=ddist['coords']['x1']['data'], + axis=-1, + ) + dataMax = scpinteg.trapezoid( + ddist['dist']['maxwell']['dist']['data'], + x=ddist['coords']['x1']['data'], + axis=-1, + ) + + # ------------ + # Derive levels, vmin, vmax + # ------------ + + Ekin_max = ddist['plasma']['Ekin_max_eV']['data'] + vminRE_2d = np.inf + vminRE_1d = np.inf + for ind in np.ndindex(dataRE.shape[:-1]): + indE = np.argmin(np.abs(ddist['coords']['x0']['data'] - Ekin_max[ind])) + sli = ind + (indE, slice(None)) + vmaxRE_2d = np.nanmax(ddist['dist']['RE']['dist']['data'][sli]) + vminRE_2d = min(vminRE_2d, vmaxRE_2d) + sli = ind + (indE,) + vmaxRE_1d = dataRE[sli] + vminRE_1d = min(vminRE_1d, vmaxRE_1d) + vmaxRE_2d = np.nanmax(ddist['dist']['RE']['dist']['data']) + vmaxRE_1d = np.nanmax(dataRE) + vmaxMax_2d = np.nanmax(ddist['dist']['maxwell']['dist']['data']) + vmaxMax_1d = np.nanmax(dataMax) + + # 1d + vmaxlog10_1d = np.log10(max(vmaxRE_1d, vmaxMax_1d)) + dlog10_1d = vmaxlog10_1d - np.log10(vminRE_1d) + vmaxlog10_1d = np.ceil(vmaxlog10_1d) + vminlog10_1d = np.floor(vmaxlog10_1d - 1 - 1.2*dlog10_1d) + vmax_1d = 10**vmaxlog10_1d + vmin_1d = 10**vminlog10_1d + + # 2d + vmaxlog10_2d = np.log10(max(vmaxRE_2d, vmaxMax_2d)) + dlog10_2d = vmaxlog10_2d - np.log10(vminRE_2d) + vmaxlog10_2d = np.ceil(vmaxlog10_2d) + vminlog10_2d = np.floor(vmaxlog10_2d - 1 - 1.2*dlog10_2d) + levels_2d = np.logspace(vminlog10_2d, vmaxlog10_2d - 1, 6) + + # -------------- + # labels + # -------------- + + dlabel = {} + for ind in np.ndindex(ddist['dist']['RE']['dist']['data'].shape[:-2]): + Te = ddist['plasma']['Te_eV']['data'][ind] * 1e-3 + jpf = ddist['plasma']['jp_fraction_re']['data'][ind] + Ek = ddist['plasma']['Ekin_max_eV']['data'][ind] + Ek = f"{Ek*1e-3:3.0f} keV" if np.log10(Ek) <= 6 else f"{Ek*1e-6:2.0f} Mev" + + dlabel[ind] = f"{jpf:2.1f} , {Te:2.1f} keV, {Ek}" + + # title + ne = np.unique(ddist['plasma']['ne_m3']['data']) + assert ne.size == 1 + jp = np.unique(ddist['plasma']['jp_Am2']['data']) + assert jp.size == 1 + tit = ( + r"$n_e$" + + f" = {ne[0]:2.1e}, " + + r"$j_{P,tot}$" + + f" = {jp[0]*1e-6:2.1f} MA/m2" + ) + + # -------------- + # print + # -------------- + + _print(ddist) + + # -------------- + # prepare axes + # -------------- + + dmargin = { + 'left': 0.12, 'right': 0.98, + 'bottom': 0.06, 'top': 0.97, + 'wspace': 0.25, 'hspace': 0.10, + } + + fig = plt.figure(figsize=figsize) + + gs = gridspec.GridSpec(ncols=1, nrows=2, **dmargin) + dax = {} + + # -------------- + # axes - 2d + # -------------- + + ax = fig.add_subplot(gs[0, 0], aspect='auto', xscale='log') + ax.set_ylabel( + r'$\theta_{e_0,B}$ (deg)', + fontsize=fontsize, + fontweight='bold', + ) + ax.set_title(tit, fontsize=fontsize, fontweight='bold') + ax.text( + 0.01, + 0.99, + '(a)', + horizontalalignment='left', + verticalalignment='top', + fontsize=fontsize, + fontweight='bold', + transform=ax.transAxes, + ) + + dax['2d'] = ax + + # -------------- + # axes - 1d + # -------------- + + ax = fig.add_subplot(gs[1, 0], aspect='auto', sharex=ax, yscale='log') + ax.set_xlabel('E (keV)', fontsize=fontsize, fontweight='bold') + ax.set_ylabel( + f" ({units1d})", + fontsize=fontsize, + fontweight='bold', + ) + ax.text( + 0.01, + 0.99, + '(b)', + horizontalalignment='left', + verticalalignment='top', + fontsize=fontsize, + fontweight='bold', + transform=ax.transAxes, + ) + + dax['1d'] = ax + + dax = ds._generic_check._check_dax(dax) + + # ------------ + # plot 1d + # ------------ + + kax = '1d' + dcolor = {} + if dax.get(kax) is not None: + ax = dax[kax]['handle'] + + # for legend + ax.plot( + [], [], + c='w', + ls='-', + lw=1, + label=r"$j_{frac}$, $T_e$, $E_{e_0,max}$", + ) + + # loop plot + for ind in np.ndindex(dataRE.shape[:-1]): + sli = ind + (slice(None),) + + # Max + l0, = ax.plot( + ddist['coords']['x0']['data']*1e-3, + dataMax[sli], + ls='-', + lw=1, + ) + dcolor[ind] = l0.get_color() + + # RE + ax.plot( + ddist['coords']['x0']['data']*1e-3, + dataRE[sli], + ls='--', + lw=1, + color=dcolor[ind], + ) + + # Total + ax.plot( + ddist['coords']['x0']['data']*1e-3, + dataMax[sli] + dataRE[sli], + ls='-', + lw=2, + color=dcolor[ind], + label=dlabel[ind], + ) + + # Add critical energy + Ec = tf.physics_tools.electrons.convert_momentum_velocity_energy( + momentum_normalized=ddist['dist']['RE']['p_crit']['data'], + )['energy_kinetic_eV']['data'] + for ec in np.unique(Ec): + ax.axvline(ec*1e-3, c='k', lw=1, ls='--') + + # decorate + ax.legend(loc="upper right") + ax.grid(True) + ax.set_ylim(vmin_1d, vmax_1d) + + # ------------ + # plot 2d + # ------------ + + kax = '2d' + if dax.get(kax) is not None: + ax = dax[kax]['handle'] + + # loop plot + for ind in np.ndindex(dataRE.shape[:-1]): + sli = ind + (slice(None), slice(None)) + + # data + data = ( + ddist['dist']['maxwell']['dist']['data'][sli] + + ddist['dist']['RE']['dist']['data'][sli] + ) + + # contour + ax.contour( + ddist['coords']['x0']['data']*1e-3, + ddist['coords']['x1']['data']*180/np.pi, + data.T, + levels_2d, + colors=dcolor[ind], + ) + + # decorate + ax.grid(True) + + # -------------- + # save + # -------------- + + savefig( + fig=fig, + pfe_save=pfe_save, + path_save=path_save, + file=__file__, + ) + + return dax, ddist + + +# ##################################################### +# ##################################################### +# _print +# ##################################################### + + +def _print(ddist, sep=' '): + + # ----------- + # header + + head = [ + 'ind', + 'Te (keV)', + 'ne (1e20/m3)', 'Max / RE', + 'jp (MA/m2)', 'Max / RE', + ] + lmax = np.max([len(ss) for ss in head]) + + # ----------- + # header + + lc = [] + for ind in np.ndindex(ddist['dist']['RE']['dist']['data'].shape[:-2]): + Te = ddist['plasma']['Te_eV']['data'][ind]*1e-3 + ne = ddist['plasma']['ne_m3']['data'][ind]*1e-20 + jp = ddist['plasma']['jp_Am2']['data'][ind]*1e-6 + ne_max = ddist['dist']['maxwell']['integ_ne']['data'][ind]*1e-20 + ne_RE = ddist['dist']['RE']['integ_ne']['data'][ind]*1e-20 + jp_max = ddist['dist']['maxwell']['integ_jp']['data'][ind]*1e-6 + jp_RE = ddist['dist']['RE']['integ_jp']['data'][ind]*1e-6 + + cc = [ + str(ind), + f'{Te:2.1f}', + f'{ne:2.2f}', f"{ne_max:2.2f} / {ne_RE:2.2f}", + f'{jp:2.2f}', f"{jp_max:2.2f} / {jp_RE:2.2f}", + ] + lc.append(cc) + lmax = max(lmax, np.max([len(ss) for ss in cc])) + + # ---------------- + # concatenate + + line = sep.join(['-'*lmax for ss in head]) + head = sep.join([ss.ljust(lmax) for ss in head]) + lc = [ + sep.join([ss.ljust(lmax) for ss in cc]) + for cc in lc + ] + + msg = '\n'.join([head, line] + lc) + print(msg) + + return diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig02_dist_type.py b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig02_dist_type.py new file mode 100644 index 0000000..a44cb55 --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig02_dist_type.py @@ -0,0 +1,398 @@ +import os + +import numpy as np +import scipy.integrate as scpinteg +import astropy.units as asunits +import matplotlib.pyplot as plt +import matplotlib.gridspec as gridspec +import datastock as ds +import tofu as tf + + +from ._savefig import main as savefig + + +tfphysdist = tf.physics_tools.electrons.distribution + + +# ##################################################### +# ##################################################### +# DEFAULTS +# ##################################################### + + +_PATH_HERE = os.path.dirname(__file__) + + +_DDIST = { + "maxwell": { + 'Te_eV': 1e3, + 'ne_m3': 1e19, + 'jp_Am2': 1e6, + }, + 'RE': { + 'dreicer': { + 'dominant': 'dreicer', + 'jp_fraction_re': 0.9, + 'Efield_par_Vm': 0.5, + 'Ekin_min_eV': 100, + # used for legacy + 'pnormW': 0, + 'Ekin_max_eV': 0, + }, + 'avalanche 100 keV': { + 'dominant': 'avalanche', + 'jp_fraction_re': 0.9, + 'Efield_par_Vm': 0.5, + 'Ekin_max_eV': 100e3, + 'Ekin_min_eV': 100, + 'sigmap': 0.1, + # used for legacy + 'pnormW': 0, + }, + 'avalanche 10 MeV': { + 'dominant': 'avalanche', + 'jp_fraction_re': 0.9, + 'Efield_par_Vm': 0.5, + 'Ekin_max_eV': 10e6, + 'Ekin_min_eV': 100, + 'sigmap': 10, + # used for legacy + 'pnormW': 0, + }, + }, + # coords + 'coords': { + 'E_eV': np.logspace(0, 8, 80), + 'theta': np.linspace(0, 180, 181) * np.pi / 180, + } +} + + +_DDIST_PLOT = { + 'maxwell': {'color': 'k', 'hatch': 'None'}, + 'dreicer': {'color': 'tab:blue', 'hatch': '||'}, + 'avalanche 100 keV': {'color': 'tab:orange', 'hatch': '--'}, + 'avalanche 10 MeV': {'color': 'tab:green', 'hatch': 'o'}, +} + + +# ##################################################### +# ##################################################### +# main +# ##################################################### + + +def main( + # maxwell + jp_Am2=None, + # plot + figsize=(5, 7), + fontsize=14, + # save + path_save=None, + pfe_save=None, + # unused + **kwdargs, +): + """ Plot a Maxwellian + RE distribution + + Includes: + - a 2d (E, pitch) contour plot + - a 1d (E,) plot + + """ + + # ------------ + # inputs + # ------------ + + # ------------ + # compute + # ------------ + + ddist = {} + for kdist, vdist in _DDIST['RE'].items(): + + kwd = dict(_DDIST['maxwell']) + kwd.update(**vdist) + kwd.update(**_DDIST['coords']) + + if jp_Am2 is not None: + kwd['jp_Am2'] = jp_Am2 + + ddist[kdist] = tfphysdist.get_distribution(**kwd) + + # sanity check - jp_RE + lip = [v0['dist']['RE']['integ_jp']['data'] for v0 in ddist.values()] + assert np.allclose(lip, lip[0]) + + # sanity check - Maxwell all the same + lMax = np.array([ + v0['dist']['maxwell']['dist']['data'] for v0 in ddist.values() + ]) + err_rel = (lMax - lMax[0:1, ...]) + iok = lMax[0] > 0 + sli = (slice(None), iok) + err_rel[sli] = 100 * err_rel[sli] / lMax[(slice(0, 1), iok)] + assert np.all(err_rel < 0.5), np.max(err_rel) + + # units + kdomref = list(_DDIST['RE'].keys())[0] + units2d = asunits.Unit(ddist[kdomref]['dist']['RE']['dist']['units']) + units1d = units2d * asunits.Unit(ddist[kdomref]['coords']['x1']['units']) + + # ------------ + # Derive 1d data + # ------------ + + d1d = {} + for dom in ddist.keys(): + d1d[dom] = scpinteg.trapezoid( + ddist[dom]['dist']['RE']['dist']['data'], + x=ddist[dom]['coords']['x1']['data'], + axis=-1, + ) + + dataMax = scpinteg.trapezoid( + ddist[kdomref]['dist']['maxwell']['dist']['data'], + x=ddist[kdomref]['coords']['x1']['data'], + axis=-1, + ) + + # ------------ + # Derive levels, vmin, vmax + # ------------ + + # 1d + vmaxlog10_1d = np.ceil(np.log10(np.nanmax(dataMax))) + vmax_1d = 10**(vmaxlog10_1d) + vmin_1d = 10**(vmaxlog10_1d - 8) + + # 2d + dvminmaxlog10_2d = {} + for dom in ddist.keys(): + + Ekin_max = ddist[dom]['plasma']['Ekin_max_eV']['data'] + if np.allclose(Ekin_max, 0.): + Ekin_max = np.full(Ekin_max.shape, 1e7) + vminRE_2d = np.inf + for ind in np.ndindex(d1d[dom].shape[:-1]): + indE = np.argmin(np.abs( + ddist[dom]['coords']['x0']['data'] - Ekin_max[ind] + )) + sli = ind + (indE, slice(None)) + vmaxRE_2d = np.nanmax( + ddist[dom]['dist']['RE']['dist']['data'][sli] + ) + vminRE_2d = min(vminRE_2d, vmaxRE_2d) + sli = ind + (indE,) + vmaxRE_2d = np.nanmax(ddist[dom]['dist']['RE']['dist']['data']) + vmaxMax_2d = np.nanmax( + ddist[dom]['dist']['maxwell']['dist']['data'] + ) + + # 2d + vmaxlog10_2d = np.log10(max(vmaxRE_2d, vmaxMax_2d)) + dlog10_2d = vmaxlog10_2d - np.log10(vminRE_2d) + dvminmaxlog10_2d[dom] = { + 'max': np.ceil(vmaxlog10_2d), + 'min': np.floor(vmaxlog10_2d - 1 - 1.2*dlog10_2d), + } + + levels_2d = np.logspace( + np.min([vv['min'] for vv in dvminmaxlog10_2d.values()]), + np.max([vv['max'] for vv in dvminmaxlog10_2d.values()]) - 1, + 6, + ) + + # -------------- + # labels + # -------------- + + dlabel = {} + for dom in ddist.keys(): + dlabel[dom] = f"{dom}" + + # title + Te = np.unique(ddist[kdomref]['plasma']['Te_eV']['data']) + assert Te.size == 1 + ne = np.unique(ddist[kdomref]['plasma']['ne_m3']['data']) + assert ne.size == 1 + jp = np.unique(ddist[kdomref]['plasma']['jp_Am2']['data']) + assert jp.size == 1 + jp_fraction_re = np.unique(ddist[kdomref]['plasma']['jp_fraction_re']['data']) + assert jp_fraction_re.size == 1 + Ekin_max_eV = np.unique(ddist[kdomref]['plasma']['Ekin_max_eV']['data']) + assert Ekin_max_eV.size == 1 + sigmap = np.unique(ddist[kdomref]['plasma']['sigmap']['data']) + assert sigmap.size == 1 + vloop = np.unique(ddist[kdomref]['plasma']['Efield_par_Vm']['data']) + assert vloop.size == 1 + tit = ( + r"$n_e$" + f" = {ne[0]:2.1e} " + r"$/m^3$" "+, " + + r"$T_e$" + f" = {Te[0]*1e-3:2.1f} keV, " + + r"$V_{loop}$" + f" = {vloop[0]:2.1f} V/m\n" + + r"$j_{P,tot}$" + f" = {jp[0]*1e-6:2.1f} MA/m2, " + + r"$F_{RE}$" + f" = {jp_fraction_re[0]:2.1f}" + ) + + # -------------- + # prepare axes + # -------------- + + dmargin = { + 'left': 0.13, 'right': 0.97, + 'bottom': 0.07, 'top': 0.93, + 'wspace': 0.25, 'hspace': 0.12, + } + + fig = plt.figure(figsize=figsize) + + gs = gridspec.GridSpec(ncols=1, nrows=2, **dmargin) + dax = {} + + # -------------- + # axes - 2d + # -------------- + + ax = fig.add_subplot(gs[0, 0], aspect='auto', xscale='log') + ax.set_ylabel( + r'$\theta_{e_0,B}$ (deg)', + fontsize=fontsize, + fontweight='bold', + ) + ax.set_title(tit, fontsize=fontsize, fontweight='bold') + ax.text( + 0.01, + 0.99, + '(a)', + horizontalalignment='left', + verticalalignment='top', + fontsize=fontsize, + fontweight='bold', + transform=ax.transAxes, + ) + ax.set_yticks([0, 30, 60, 90, 120, 150, 180]) + ax.set_xticks(np.logspace(-3, 5, 9)) + + dax['2d'] = ax + + # -------------- + # axes - 1d + # -------------- + + ax = fig.add_subplot(gs[1, 0], aspect='auto', sharex=ax, yscale='log') + ax.set_xlabel('E (keV)', fontsize=fontsize, fontweight='bold') + ax.set_ylabel( + f" ({units1d})", + fontsize=fontsize, + fontweight='bold', + ) + ax.text( + 0.01, + 0.99, + '(b)', + horizontalalignment='left', + verticalalignment='top', + fontsize=fontsize, + fontweight='bold', + transform=ax.transAxes, + ) + + ax.set_yticks(np.power(10, vmaxlog10_1d + np.arange(-8, 1))) + dax['1d'] = ax + + dax = ds._generic_check._check_dax(dax) + + # ticklabels size + for kax, vax in dax.items(): + dax[kax]['handle'].tick_params( + axis='both', + which='major', + labelsize=fontsize - 1, + ) + + # ------------ + # plot 1d + # ------------ + + kax = '1d' + if dax.get(kax) is not None: + ax = dax[kax]['handle'] + + # Maxwell + ax.plot( + ddist[kdomref]['coords']['x0']['data']*1e-3, + dataMax[0, ...], + ls='-', + lw=1, + c='k', + label='Maxwellian', + ) + + # loop plot + for dom in ddist.keys(): + + # RE + l0, = ax.plot( + ddist[dom]['coords']['x0']['data']*1e-3, + d1d[dom][0, ...], + ls='-', + lw=1, + color=_DDIST_PLOT[dom]['color'], + label=dlabel[dom], + ) + + # decorate + ax.legend(loc="upper right") + ax.grid(True) + ax.set_ylim(vmin_1d, vmax_1d) + + # ------------ + # plot 2d + # ------------ + + kax = '2d' + if dax.get(kax) is not None: + ax = dax[kax]['handle'] + + # Maxwell + ax.contour( + ddist[kdomref]['coords']['x0']['data']*1e-3, + ddist[kdomref]['coords']['x1']['data']*180/np.pi, + ddist[kdomref]['dist']['maxwell']['dist']['data'][0, ...].T, + levels_2d, + colors='k', + ) + + # loop plot + for dom in ddist.keys(): + + # data + data = ddist[dom]['dist']['RE']['dist']['data'][0, ...] + + # contour + ax.contour( + ddist[dom]['coords']['x0']['data']*1e-3, + ddist[dom]['coords']['x1']['data']*180/np.pi, + data.T, + levels_2d, + colors=_DDIST_PLOT[dom]['color'], + ) + + # decorate + ax.grid(True) + + # -------------- + # save + # -------------- + + savefig( + fig=fig, + pfe_save=pfe_save, + path_save=path_save, + file=__file__, + ) + + return dax, ddist diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig03_validate_ff.py b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig03_validate_ff.py new file mode 100644 index 0000000..5ffc00c --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig03_validate_ff.py @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:8a45ef1dab424101ec3ee218115c2ba3fca98e9d97b7505f6d7da3119ca42778 +size 7677 diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig04_bremsstrahlung.py b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig04_bremsstrahlung.py new file mode 100644 index 0000000..202ddb4 --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig04_bremsstrahlung.py @@ -0,0 +1,459 @@ + + +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.gridspec as gridspec +import matplotlib.transforms as transforms +import datastock as ds +import tofu as tf + + +from . import _load_spect_anis +from ._savefig import main as savefig + + +tfphysemis = tf.physics_tools.electrons.emission + + +# ##################################################### +# ##################################################### +# DEFAULTS +# ##################################################### + + +_CASES = { + 'case': { + '0': { + 'Te': 0.1e3, + 'jp_frac': 0.9, + 'Ekin_max_eV': 100e3, + 'color': 'r', + 'hatch': '//', + }, + '1': { + 'Te': 2.0e3, + 'jp_frac': 0.1, + 'Ekin_max_eV': 10e6, + 'color': 'b', + 'hatch': "\\", + }, + }, + 'theta_ph_vsB': { + 'val': np.r_[0, 0.5, 1]*np.pi, + 'ls': ['-', '--', ':'], + }, + 'E_ph_eV': { + 'val': np.r_[0.1, 2, 20]*1e3, + 'ls': ['-', '-', '-'], + }, +} + + +# ###################################### +# ###################################### +# main +# ###################################### + + +def main( + d2cross_phi=None, + # dist + ne_m3=None, + pnormW=None, + Ekin_max_eV=None, + # Te_eV=1e3 * np.linspace(0.1, 2.5, 25), + Te_eV=None, + # jp_fraction_re=np.linspace(0.025, 0.975, 39), + jp_fraction_re=None, + # emiss + E_ph_eV=None, + # cases + cases=None, + # plot + figsize=(15, 8), + fontsize=14, + # save + path_save=None, + pfe_save=None, + # unused + **kwdargs, +): + """ Plot a selection of Maxwellian and RE free-free spectra + + Includes selected angular amissivity + Includes a 2d contour plot of photon energy threshold for RE dominance + + """ + + # -------------- + # cases + # -------------- + + if cases is None: + cases = _CASES + + # -------------- + # demiss + # -------------- + + demiss, ddist, dmix, data_source = _load_spect_anis.main( + dmix='H', + # d2cross + d2cross_phi=d2cross_phi, + # dist + ne_m3=ne_m3, + pnormW=pnormW, + Ekin_max_eV=Ekin_max_eV, + jp_fraction_re=jp_fraction_re, + ) + + units = demiss['emiss']['RE']['ff']['units'] + ne = np.unique(ddist['plasma']['ne_m3']['data'])[0] + jp = np.unique(ddist['plasma']['jp_Am2']['data'])[0] + nEkin = demiss['emiss']['maxwell']['ff']['data'].shape[0] + + tit = ( + r"$n_e$" + f" = {ne:1.0e} /m3\n" + + r"$j_{P,RE}$" + f" = {jp*1e-6:2.1f} MA/m2" + ) + + # -------------- + # Elim + # -------------- + + shape = ddist['plasma']['Te_eV']['data'].shape + Elim = np.full(shape, np.nan) + for ii, ind in enumerate(np.ndindex(shape)): + + sli_emiss = ind + (slice(None), 0) + emiss_RE = demiss['emiss']['RE']['ff']['data'][sli_emiss] + emiss_max = demiss['emiss']['maxwell']['ff']['data'][sli_emiss] + + ilim = (emiss_RE > emiss_max) + if np.any(ilim): + Elim[ind] = np.min(demiss['E_ph']['data'][ilim]) + + # -------------- + # prepare axes + # -------------- + + dmargin = { + 'left': 0.06, 'right': 0.98, + 'bottom': 0.06, 'top': 0.93, + 'wspace': 0.25, 'hspace': 0.30, + } + dmargin_theta = { + 'left': 0.06, 'right': 0.98, + 'bottom': 0.06, 'top': 0.60, + 'wspace': 0.25, 'hspace': 0.10, + } + dmargin_map = dict(dmargin) + dmargin_map['hspace'] = 0.10 + + fig = plt.figure(figsize=figsize) + + nE = len(cases['E_ph_eV']['val']) + gs = gridspec.GridSpec(ncols=nE + 2, nrows=3, **dmargin) + gs_theta = gridspec.GridSpec(ncols=nE + 2, nrows=2, **dmargin_theta) + gs_map = gridspec.GridSpec(ncols=nE + 2, nrows=nEkin, **dmargin_map) + dax = {} + + # ---------------- + # ax - spectra + # ---------------- + + ax = fig.add_subplot(gs[0, :nE], aspect='auto') + ax.set_xlabel( + r"$E_{ph}$" + ' (keV)', + fontsize=fontsize, + fontweight='bold', + ) + ax.set_ylabel( + r"$\epsilon$" + f' ({units})', + fontsize=fontsize, + fontweight='bold', + ) + dax['spectra'] = ax + + # ---------------- + # ax - theta - lin + # ---------------- + + ax0 = None + for ii in range(nE): + ax = fig.add_subplot( + gs_theta[0, ii], + aspect='auto', + sharex=ax0, + sharey=ax0, + ) + if ii == 0: + ax.set_ylabel( + 'emiss (norm.)', + fontsize=fontsize, + fontweight='bold', + ) + ax.set_xlim(0, 180) + ax.set_ylim(0, 1) + ax0 = ax + + ax.set_xticks(np.r_[0, 45, 90, 135, 180]) + dax[f'theta_{ii}_lin'] = ax + + # ---------------- + # ax - theta - log + # ---------------- + + ax0 = None + for ii in range(nE): + ax = fig.add_subplot( + gs_theta[1, ii], + aspect='auto', + sharex=ax0, + sharey=ax0, + ) + ax.set_xlabel( + r'$\theta_{ph,B}$' + ' (deg)', + fontsize=fontsize, + fontweight='bold', + ) + if ii == 0: + ax.set_ylabel( + f'{units}', + fontsize=fontsize, + fontweight='bold', + ) + ax.set_xlim(0, 180) + # ax.set_ylim(0, 1) + ax0 = ax + + dax[f'theta_{ii}_log'] = ax + + # ---------------- + # ax - Elim + # ---------------- + + ax0 = None + for ii in range(nEkin): + ax = fig.add_subplot( + gs_map[ii, nE:], + aspect='auto', + sharex=ax0, + sharey=ax0, + ) + if ii == 0: + ax0 = ax + ax.set_title( + tit, + fontsize=fontsize, + fontweight='bold', + ) + elif ii == nEkin - 1: + ax.set_xlabel('Te (keV)', fontsize=fontsize, fontweight='bold') + ax.set_ylabel('jp_frac', fontsize=fontsize, fontweight='bold') + + dax[f'Elim_{ii}'] = ax + + dax = ds._generic_check._check_dax(dax) + + # -------------- + # plot - cases + # -------------- + + Teu = np.unique(ddist['plasma']['Te_eV']['data']) + jp_fracu = np.unique(ddist['plasma']['jp_fraction_re']['data']) + Ekinu = np.unique(ddist['plasma']['Ekin_max_eV']['data']) + + # loop on cases + for i0, (k0, v0) in enumerate(cases['case'].items()): + + # slice + Te = Teu[np.argmin(np.abs(Teu - v0['Te']))] + jpf = jp_fracu[np.argmin(np.abs(jp_fracu - v0['jp_frac']))] + Ekin = Ekinu[np.argmin(np.abs(Ekinu - v0['Ekin_max_eV']))] + ic = ( + (ddist['plasma']['jp_fraction_re']['data'] == jpf) + & (ddist['plasma']['Te_eV']['data'] == Te) + & (ddist['plasma']['Ekin_max_eV']['data'] == Ekin) + ) + assert ic.sum() == 1 + ic = tuple([cc[0] for cc in ic.nonzero()]) + + # -------- + # spectra + + kax = 'spectra' + if dax.get(kax) is not None: + ax = dax[kax]['handle'] + + sli = ic + (slice(None), slice(None)) + for kdist in demiss['emiss'].keys(): + emiss_E = demiss['emiss'][kdist]['ff']['data'][sli] + + # plot + ax.fill_between( + demiss['E_ph']['data']*1e-3, + np.nanmin(emiss_E, axis=-1), + np.nanmax(emiss_E, axis=-1), + hatch=v0['hatch'], + facecolor='None', + alpha=0.5, + edgecolor=v0['color'], + ls='--' if kdist == 'RE' else '-', + label=f'{kdist}_{ic}', + ) + + ax.set_xscale('log') + ax.set_yscale('log') + ax.grid(True) + + # vlines + for i1, cc in enumerate(cases['E_ph_eV']['val']): + ax.axvline( + cc*1e-3, + c='k', + ls='--', + lw=1, + label=f"E_ph = {cc*1e-3:3.1f} keV", + ) + + # Elim + iE = np.argmin(np.abs(demiss['E_ph']['data'] - Elim[ic])) + ax.plot( + np.r_[Elim[ic], Elim[ic]] * 1e-3, + [emiss_E[iE, 0], 1e16], + color=v0['color'], + ls='--', + lw=1, + label=f"E_lim = {Elim[ic]*1e-3:3.1f} keV", + ) + + # text + trans = transforms.blended_transform_factory( + ax.transData, + ax.transAxes, + ) + ax.text( + Elim[ic] * 1e-3, + 1, + r"$E_{ph,lim}$" + f"\n = {Elim[ic] * 1e-3:2.1f} keV", + horizontalalignment='center', + verticalalignment='bottom', + fontsize=fontsize, + fontweight='bold', + color=v0['color'], + transform=trans, + ) + + ax.set_ylim(1e0, 1e16) + ax.set_xlim(1e-2, 2e4) + + # -------- + # theta + + for i1, cc in enumerate(cases['E_ph_eV']['val']): + + kax = f'theta_{i1}_lin' + if dax.get(kax) is not None: + ax = dax[kax]['handle'] + + iE = np.argmin(np.abs(demiss['E_ph']['data'] - cc)) + sli = ic + (iE, slice(None)) + + for kdist in demiss['emiss'].keys(): + emiss_theta = demiss['emiss'][kdist]['ff']['data'][sli] + + # plot + ax.plot( + demiss['theta_ph_vsB']['data']*180/np.pi, + emiss_theta / emiss_theta.max(), + # ls=cases['E_ph_eV']['ls'][i1], + ls='--' if kdist == 'RE' else '-', + lw=1, + marker='None', + color=v0['color'], + label=f'{kdist}_{ic}_{cc*1e-3:3.1f}keV', + ) + + # deco + ax.set_title( + r"$E_{ph,B}$" + f" = {cc*1e-3:3.1f} keV", + fontsize=fontsize, + fontweight='bold', + ) + + # --------- + # theta - log + + kax = f'theta_{i1}_log' + if dax.get(kax) is not None: + ax = dax[kax]['handle'] + + iE = np.argmin(np.abs(demiss['E_ph']['data'] - cc)) + sli = ic + (iE, slice(None)) + + for kdist in demiss['emiss'].keys(): + emiss_theta = demiss['emiss'][kdist]['ff']['data'][sli] + + # plot + ax.semilogy( + demiss['theta_ph_vsB']['data']*180/np.pi, + emiss_theta, + # ls=cases['E_ph_eV']['ls'][i1], + ls='--' if kdist == 'RE' else '-', + lw=1, + marker='None', + color=v0['color'], + label=f'{kdist}_{ic}_{cc*1e-3:3.1f}keV', + ) + + ax.set_ylim(2e4, 6e14) + + # -------------- + # plot - Elim + # -------------- + + for ii in range(nEkin): + kax = f'Elim_{ii}' + if dax.get(kax) is not None: + ax = dax[kax]['handle'] + + sli = (ii, slice(None), slice(None)) + cs = ax.contour( + ddist['plasma']['Te_eV']['data'][sli] * 1e-3, + ddist['plasma']['jp_fraction_re']['data'][sli], + Elim[sli] * 1e-3, + cmap=plt.cm.viridis, + levels=np.r_[1, 2, 5, 7.5, 10, 15], + vmin=0.1, + vmax=20, + ) + + # cases + for i0, (k0, v0) in enumerate(cases['case'].items()): + ax.plot( + v0['Te']*1e-3, + v0['jp_frac'], + marker='*', + markersize=8, + markerfacecolor=v0['color'], + color=v0['color'], + ) + + ax.clabel(cs, cs.levels, fontsize=12) + + ax.set_xlim(0, ddist['plasma']['Te_eV']['data'].max()*1e-3) + ax.set_ylim(0, 1) + + # -------------- + # save + # -------------- + + savefig( + fig=fig, + pfe_save=pfe_save, + path_save=path_save, + file=__file__, + ) + + return dax, demiss, ddist, d2cross_phi diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig05_emiss.py b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig05_emiss.py new file mode 100644 index 0000000..2f5e456 --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig05_emiss.py @@ -0,0 +1,416 @@ +import os + + +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.gridspec as gridspec +import datastock as ds + + +from . import _load_spect_anis +from ._fig02_dist_type import _DDIST +from ._fig04_bremsstrahlung import _CASES +from ._savefig import main as savefig + + +# ##################################################### +# ##################################################### +# DEFAULTS +# ##################################################### + + +# PATHS +_PATH_HERE = os.path.dirname(__file__) +_PATH_PAPER = os.path.dirname(_PATH_HERE) + + +_RE = 'avalanche 10 MeV' +_JP_FRAC = np.linspace(0.025, 0.975, 39) + +_DDMIX = { + 0: 'H', + 1: {'H': 0.90, 'O': 0.08, 'Fe': 0.02}, +} + + +# ##################################################### +# ##################################################### +# Main +# ##################################################### + + +def main( + dmix=None, + # cases + cases=None, + # d2cross + d2cross_phi=None, + # dist + ne_m3=None, + jp_Am2=None, + # RE + re=None, + dominant=None, + jp_fraction_re=None, + Efield_par_Vm=None, + Ekin_min_eV=None, + Ekin_max_eV=None, + sigmap=None, + pnormW=0, + # jp_fraction_re=np.linspace(0.025, 0.975, 39), + # plot + figsize=(15, 8), + fontsize=14, + # save + path_save=None, + pfe_save=None, + # unused + **kwdargs, +): + """ Validate Bremsstrahlung vs SCRAM and FLYCHK + + """ + + # -------------- + # inputs + # -------------- + + if dmix is None: + dmix = _DDMIX + nmix = len(dmix) + + if jp_fraction_re is None: + jp_fraction_re = _JP_FRAC + + # Maxwell + kwd_max = {'ne_m3': ne_m3, 'jp_Am2': jp_Am2} + + # RE + if re is None: + re = _RE + assert re in _DDIST['RE'].keys() + lRE = [ + 'dominant', 'jp_fraction_re', 'Efield_par_Vm', + 'Efield_par_Vm', 'Ekin_max_eV', 'Ekin_min_eV', + 'sigmap', 'pnormW' + ] + kwd_RE = { + kk: _DDIST['RE'][re].get(kk) if vv is None else vv + for kk, vv in locals().items() + if kk in lRE + } + + if cases is None: + cases = _CASES + + # -------------- + # load elements + # -------------- + + demiss = {} + kwd = dict(kwd_max) + kwd.update(**kwd_RE) + for ii in dmix.keys(): + demiss[ii], ddist, dmix[ii], data_source = _load_spect_anis.main( + dmix=dmix[ii], + # d2cross + d2cross_phi=d2cross_phi, + # dist + **kwd, + ) + + # extract + # E_ph = demiss[0]['common']['E_photon']['data'] + Te = ddist['plasma']['Te_eV']['data'] + ne = np.unique(ddist['plasma']['ne_m3']['data'])[0] + jp = np.unique(ddist['plasma']['jp_Am2']['data'])[0] + units = demiss[0]['emiss']['maxwell']['ff']['units'] + + # -------------- + # Elim + # -------------- + + shape = ddist['plasma']['Te_eV']['data'].shape + Elim = {kk: np.full(shape, np.nan) for kk in dmix.keys()} + for kk in dmix.keys(): + for ind in np.ndindex(shape): + + sli_anis = ind + (slice(None), 0) + sli_iso = ind[:-1] + (0, slice(None), 0) + + emiss_max = ( + demiss[kk]['emiss']['maxwell']['ff']['data'][sli_anis] + + demiss[kk]['emiss']['maxwell']['fb']['data'][sli_iso] + + demiss[kk]['emiss']['maxwell']['bb']['data'][sli_iso] + ) + + emiss_RE = demiss[kk]['emiss']['RE']['ff']['data'][sli_anis] + + ilim = (emiss_RE > emiss_max) + if np.any(ilim): + Elim[kk][ind] = np.min(demiss[kk]['E_ph']['data'][ilim]) + + # -------------- + # prepare axes + # -------------- + + dmargin = { + 'left': 0.05, 'right': 0.99, + 'bottom': 0.06, 'top': 0.93, + 'wspace': 0.20, 'hspace': 0.20, + } + + fig = plt.figure(figsize=figsize) + + gs = gridspec.GridSpec(ncols=3, nrows=nmix, **dmargin) + dax = {} + + # ---------------- + # ax - spectra + # ---------------- + + tit0 = ( + r"$n_e$" + f" = {ne:1.1e} " + r"$/m^3$, " + + r"$j_P$" + f" = {jp:1.1e} " + r"$A/m^2$" + "\n" + + f"{re} - " + ) + + ax0_spect = None + ax0_map = None + for im, mix in enumerate(sorted(dmix.keys())): + + lk = list(dmix[mix].keys()) + lc = [np.unique(dmix[mix][kk])[0] for kk in lk] + inds = np.argsort(lc)[::-1] + lstr = [f"{lk[ss]} {lc[ss]*100:3.1f} \\%" for ss in inds] + tit = ", ".join(lstr) + if im == 0: + tit = tit0 + tit + + ax = fig.add_subplot( + gs[im, :2], + sharex=ax0_spect, + sharey=ax0_spect, + aspect='auto', + ) + ax.set_title( + tit, + fontsize=fontsize, + fontweight='bold', + ) + ax.set_ylabel( + r"$\epsilon$" + f' ({units})', + fontsize=fontsize, + fontweight='bold', + ) + if im == 0: + ax0_spect = ax + elif im == nmix - 1: + ax.set_xlabel( + r"$E_{ph}$" + ' (keV)', + fontsize=fontsize, + fontweight='bold', + ) + ax.text( + 0., + 1.05, + ['(a)', '(b)'][im], + horizontalalignment='left', + verticalalignment='bottom', + transform=ax.transAxes, + ) + + dax[f'spect_{im}'] = ax + + # ---------------- + # ax - Elim + # ---------------- + + ax = fig.add_subplot( + gs[im, 2:], + aspect='auto', + sharex=ax0_map, + sharey=ax0_map, + ) + if im == 0: + ax0_map = ax + elif im == nmix - 1: + ax.set_xlabel( + r"$T_e$" + ' (keV)', + fontsize=fontsize, + fontweight='bold', + ) + ax.set_ylabel(r"$F_{RE}$", fontsize=fontsize, fontweight='bold') + + ax.text( + 0., + 1.05, + ['(c)', '(d)'][im], + horizontalalignment='left', + verticalalignment='bottom', + transform=ax.transAxes, + ) + dax[f'Elim_{im}'] = ax + + # ---------------- + # check dax format + + dax = ds._generic_check._check_dax(dax) + + # -------------- + # plot - cases + # -------------- + + Teu = np.unique(ddist['plasma']['Te_eV']['data']) + jp_fracu = np.unique(ddist['plasma']['jp_fraction_re']['data']) + Ekinu = np.unique(ddist['plasma']['Ekin_max_eV']['data']) + + # loop on cases + for i0, (k0, v0) in enumerate(cases['case'].items()): + + # slice + Te = Teu[np.argmin(np.abs(Teu - v0['Te']))] + jpf = jp_fracu[np.argmin(np.abs(jp_fracu - v0['jp_frac']))] + Ekin = Ekinu[np.argmin(np.abs(Ekinu - v0['Ekin_max_eV']))] + ic = ( + (ddist['plasma']['jp_fraction_re']['data'] == jpf) + & (ddist['plasma']['Te_eV']['data'] == Te) + & (ddist['plasma']['Ekin_max_eV']['data'] == Ekin) + ) + assert ic.sum() == 1 + ic = tuple([cc[0] for cc in ic.nonzero()]) + print(i0, k0, ic) + + # -------- + # spectra + + # loop on dmix + for ii in dmix.keys(): + + kax = f'spect_{ii}' + if dax.get(kax) is not None: + ax = dax[kax]['handle'] + + sli_anis = ic + (slice(None), slice(None)) + sli_iso = ic[:-1] + (0, slice(None), slice(None)) + + # ----------- + # total -maxwell + + emiss_max = ( + demiss[ii]['emiss']['maxwell']['ff']['data'][sli_anis] + + demiss[ii]['emiss']['maxwell']['fb']['data'][sli_iso] + + demiss[ii]['emiss']['maxwell']['bb']['data'][sli_iso] + ) + + # plot + ax.fill_between( + demiss[ii]['E_ph']['data']*1e-3, + np.nanmin(emiss_max, axis=-1), + np.nanmax(emiss_max, axis=-1), + hatch=v0['hatch'], + facecolor='None', + alpha=0.5, + edgecolor=v0['color'], + ls='-', + ) + + # ----------- + # ff -maxwell + + ff_max = demiss[ii]['emiss']['maxwell']['ff']['data'][sli_anis] + + # plot + ax.plot( + demiss[ii]['E_ph']['data']*1e-3, + np.nanmean(ff_max, axis=-1), + color=v0['color'], + ls='--', + ) + + # ----------- + # ff - RE + + emiss_RE = demiss[ii]['emiss']['RE']['ff']['data'][sli_anis] + + # plot + ax.fill_between( + demiss[ii]['E_ph']['data']*1e-3, + np.nanmin(emiss_RE, axis=-1), + np.nanmax(emiss_RE, axis=-1), + hatch=v0['hatch'], + facecolor='None', + alpha=0.5, + edgecolor=v0['color'], + ls='--', + label=r'$\epsilon_{ff}^{RE}$' if i0 == 0 else None, + ) + + # decorate + if i0 == 0: + vmax_log10 = np.ceil(np.log10(np.nanmax(emiss_max))) + + ax.set_xscale('log') + ax.set_yscale('log') + ax.set_ylim(10**(vmax_log10 - 18), 10**vmax_log10) + ax.set_xlim(left=demiss[ii]['E_ph']['data'][0]*1e-3) + ax.grid(True) + + # legend + lab = ( + r"$\epsilon_{ff}^{Max} + \epsilon_{fb}^{Max} + $" + r"$\epsilon_{bb}^{Max}$" + ) + ax.plot([], [], c='k', ls='-', label=lab) + lab = r"$\epsilon_{ff}^{Max}$" + ax.plot([], [], c='k', ls='--', label=lab) + ax.legend(loc='upper right') + + # -------------- + # plot - Elim + # -------------- + + for kk in dmix.keys(): + kax = f'Elim_{kk}' + if dax.get(kax) is not None: + ax = dax[kax]['handle'] + + sli = (0, slice(None), slice(None)) + cs = ax.contour( + ddist['plasma']['Te_eV']['data'][sli] * 1e-3, + ddist['plasma']['jp_fraction_re']['data'][sli], + Elim[kk][sli] * 1e-3, + cmap=plt.cm.viridis, + levels=np.r_[1, 2, 5, 7.5, 10, 15, 20], + vmin=0.1, + vmax=20, + ) + + # cases + for i0, (k0, v0) in enumerate(cases['case'].items()): + ax.plot( + v0['Te']*1e-3, + v0['jp_frac'], + marker='*', + markersize=8, + markerfacecolor=v0['color'], + color=v0['color'], + ) + + ax.clabel(cs, cs.levels, fontsize=12) + + ax.set_xlim(0, ddist['plasma']['Te_eV']['data'].max()*1e-3) + ax.set_ylim(0, 1) + + # -------------- + # save + # -------------- + + savefig( + fig=fig, + pfe_save=pfe_save, + path_save=path_save, + file=__file__, + ) + + return dax, demiss, Elim, dmix diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig06_tokamak.py b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig06_tokamak.py new file mode 100644 index 0000000..9665de1 --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig06_tokamak.py @@ -0,0 +1,1223 @@ +import os +import copy + + +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.gridspec as gridspec +import matplotlib.path as mpath +import matplotlib.patches as mpatches +import matplotlib.transforms as transforms +import datastock as ds +import tofu as tf + + +from ._savefig import main as savefig + + +# ##################################################### +# ##################################################### +# DEFAULTS +# ##################################################### + + +_PATH_HERE = os.path.dirname(__file__) +_PATH_SAVE = os.path.join(os.path.dirname(_PATH_HERE), 'figures') + + +_DR = { + 'R0': 1.8, + 'rplasma': 0.60, + 'RVes': [1.2, 2.66], + 'Rcryo': 4.6, + '' + 'PP_R': np.r_[2.50, 4.7], # 4.2 + 'PP_width': 0.47, + 'PP_phi': np.r_[-180, -20, 0, 20] * np.pi/180, +} + + +_DSENSORS = { + 'in': { + 'pp': 0, + 'R': 2.55, + 'cw': False, + 'rplasma_ratio': 0.7, + 'color': 'b', + 'marker': '.', + 'ms': 2, + 'alpha': 0.6, + 'shielding': True, + 'dogleg': { + 'frac_h': 0.5, + 'frac_v': 0.75, + }, + 'text': { + 'str': "(0)", + }, + }, + 'ex1': { + 'pp': 1, + 'R': 6, + 'cw': False, + 'color': 'g', + 'width': 0.10, + 'dist': 4, + 'marker': '.', + 'ms': 2, + 'alpha': 0.6, + 'shielding': True, + 'opening': True, + 'wall': True, + 'beamdump': True, + 'neutrons_length': 1, + 'neutrons_width': 0.2, + 'reflector': { + 'angle': 15*np.pi/180, + 'R': _DR['PP_R'][1] + 0.1, + 'width': 0.25, + 'color': 'k', + 'lw': 2, + }, + 'text': { + 'str': "(2)", + }, + }, + 'ex2': { + 'pp': 2, + 'R': 6, + 'cw': False, + 'color': 'g', + 'width': 0.10, + 'dist': 4, + 'marker': '.', + 'ms': 2, + 'alpha': 0.6, + 'shielding': True, + 'opening': True, + 'wall': True, + 'beamdump': True, + 'neutrons_length': 1, + 'neutrons_width': 0.2, + 'text': { + 'str': "(1)", + }, + }, + 'ex3': { + 'pp': 3, + 'R': 6, + 'cw': False, + 'color': (0.1, 0.9, 0.1), + 'width': 0.10, + 'dist': 4, + 'marker': '.', + 'ms': 2, + 'alpha': 0.6, + 'shielding': True, + 'opening': True, + 'wall': True, + 'beamdump': True, + 'neutrons_length': 1, + 'neutrons_width': 0.2, + 'reflector': { + 'angle': -15*np.pi/180, + 'R': _DR['PP_R'][0], + 'width': 0.25, + 'color': 'k', + 'lw': 2, + }, + 'text': { + 'str': "(3)", + }, + }, +} + + +# ##################################################### +# ##################################################### +# main +# ##################################################### + + +def main( + # tokamak + R0=None, + rplasma=None, + RVes=None, + Rcryo=None, + # port plug + PP_R=None, + PP_width=None, + PP_phi=None, + # sensors + res=None, + dsensors=None, + # plot + figsize=(5, 7), + fontsize=12, + # save + path_save=None, + pfe_save=None, + # unused + **kwdargs, +): + """ Plot a top view of a SPARC-like tokamak with diagnostics + + Includes: + - a in-vessel view + - 3 ex-cryostat view with or without reflectors + + Also plots the associated observation angles vs the local magnetic field + + """ + + # -------------- + # Load SPARC + # -------------- + + config, dinput = _fig02_check(**locals()) + + phi = np.pi * np.linspace(-1, 1, 181) + cos = np.cos(phi) + sin = np.sin(phi) + + # -------------- + # prepare axes + # -------------- + + dmargin = { + 'left': 0.11, 'right': 0.94, + 'bottom': 0.06, 'top': 0.99, + 'wspace': 0.25, 'hspace': 0.20, + } + + fig = plt.figure(figsize=figsize) + + gs = gridspec.GridSpec(ncols=1, nrows=2, **dmargin) + dax = {} + + # -------------- + # axes - hor + # -------------- + + ax = fig.add_subplot(gs[0, 0], aspect='equal') + ax.set_xlabel('X (m)', fontsize=fontsize, fontweight='bold') + ax.set_ylabel('Y (m)', fontsize=fontsize, fontweight='bold') + ax.text( + 0.01, + 0.99, + '(a)', + horizontalalignment='left', + verticalalignment='top', + fontsize=fontsize, + fontweight='bold', + transform=ax.transAxes, + ) + + dax['hor'] = ax + + # -------------- + # axes - ang vs rplasma + # -------------- + + ax = fig.add_subplot(gs[1, 0], aspect='auto') + ax.set_xlabel('r / a', fontsize=fontsize, fontweight='bold') + ax.set_ylabel( + r'$\theta_{ph,B}$ (deg)', + fontsize=fontsize, + fontweight='bold', + ) + ax.text( + 0.01, + 0.99, + '(b)', + horizontalalignment='left', + verticalalignment='top', + fontsize=fontsize, + fontweight='bold', + transform=ax.transAxes, + ) + + dax['theta_vs_B'] = ax + + dax = ds._generic_check._check_dax(dax) + + # -------------- + # plot hor + # -------------- + + kax = 'hor' + if dax.get(kax) is not None: + ax = dax[kax]['handle'] + + # -------------- + # plot R + + lk = [kk for kk in dinput.keys() if kk[0] == 'R'] + for k0 in lk: + + if 'plasma' in k0: + inner = dinput[k0]['data'][0] * np.array([cos, sin]).T + outer = dinput[k0]['data'][1] * np.array([cos, sin]).T + vertices = np.concatenate((inner, outer[::-1]), axis=0) + + codes = np.ones( + len(inner), + dtype=mpath.Path.code_type, + ) * mpath.Path.LINETO + codes[0] = mpath.Path.MOVETO + all_codes = np.concatenate((codes, codes)) + + path = mpath.Path(vertices, all_codes) + patch = mpatches.PathPatch( + path, + facecolor='r', + alpha=0.1, + edgecolor='r', + ) + ax.add_patch(patch) + + else: + for v1 in dinput[k0]['data']: + ax.plot( + v1*cos, + v1*sin, + **dinput[k0]['prop'], + ) + + # -------------- + # add arrows + + R = dinput['R0']['data'][0] + 0.5 * dinput['rplasma']['data'][0] + phi = np.r_[130, 170] * np.pi / 180 + # dist = R * np.hypot( + # np.cos(phi[0]) - np.cos(phi[1]), + # np.sin(phi[0]) - np.sin(phi[1]), + # ) + # rad = (R * (1 - np.cos(np.abs(np.diff(phi)/2))) / dist)[0] + rad = -0.3 + ax.annotate( + "", + xy=(R*np.cos(phi[0]), R*np.sin(phi[0])), + xycoords='data', + xytext=(R*np.cos(phi[1]), R*np.sin(phi[1])), + textcoords='data', + color='r', + fontweight='bold', + fontsize=fontsize, + horizontalalignment='center', + verticalalignment='center', + arrowprops=dict( + arrowstyle="->", + lw=1.5, + color='r', + shrinkA=5, shrinkB=5, + patchA=None, patchB=None, + connectionstyle=f'arc3,rad={rad}', + ), + ) + ax.text( + R*np.cos(np.mean(phi)), + R*np.sin(np.mean(phi)), + "RE", + color='r', + horizontalalignment='left', + verticalalignment='top', + fontweight='bold', + fontsize=fontsize, + ) + + # -------------- + # plot port plug + + for ii, phi in enumerate(dinput['PP_phi']['data']): + + # edges + width = dinput['PP_width']['data'][0] + ppR0 = dinput['PP_R']['data'][0] + ppR1 = dinput['PP_R']['data'][1] + length = ppR1 - ppR0 + cent = 0.5 * (ppR0+ppR1) * np.r_[np.cos(phi), np.sin(phi)] + xy = ( + cent[0] - 0.5 * length, + cent[1] - 0.5 * width, + ) + + # patch white + patch = mpatches.Rectangle( + xy, + length, + width, + angle=phi*180/np.pi, + rotation_point='center', + facecolor='w', + alpha=1., + edgecolor='None', + zorder=10, + ) + ax.add_patch(patch) + + # central line + ppR0 = dinput['PP_R']['data'][0] + ppR1 = dinput['PP_R']['data'][1] + + centx = np.r_[ppR0, ppR1] * np.cos(phi) + centy = np.r_[ppR0, ppR1] * np.sin(phi) + + ax.plot( + centx, + centy, + c='k', + lw=1, + ls='--', + alpha=0.3, + zorder=15, + ) + + # patch dark + patch = mpatches.Rectangle( + (xy[0], cent[1] - 0.4*width), + length, + width*0.80, + angle=phi*180/np.pi, + rotation_point='center', + facecolor='k', + alpha=0.5, + edgecolor='None', + zorder=20, + ) + ax.add_patch(patch) + + # edges + ephi = np.r_[-np.sin(phi), np.cos(phi)] + edgex = ( + centx[None, :] + + 0.5 * width * ephi[0] * np.r_[1, np.nan, -1][:, None] + ).ravel() + edgey = ( + centy[None, :] + + 0.5 * width * ephi[1] * np.r_[1, np.nan, -1][:, None] + ).ravel() + + ax.plot( + edgex, + edgey, + c='k', + lw=1, + ls='-', + zorder=20, + ) + + # -------------- + # add sensors + + dsensors = _sensors( + res=res, + dsensors=dsensors, + dinput=dinput, + ) + + for k0, v0 in dsensors.items(): + + # FOV + patch = mpatches.PathPatch( + v0['path'], + facecolor=v0['color'], + alpha=v0['alpha'], + zorder=50, + edgecolor=v0['color'], + ) + ax.add_patch(patch) + + # sensor + ax.plot( + [v0['cent'][0]], + [v0['cent'][1]], + c=v0['color'], + ls='None', + lw=2, + marker=v0['marker'], + zorder=30, + label=v0.get('label', k0), + ) + + # FOV sampling + ax.plot( + v0['ptsx'], + v0['ptsy'], + c=v0['color'], + marker=v0.get('marker', '.'), + ls='None', + zorder=60, + ms=v0.get('ms', 4), + ) + + # patch + if v0.get('patch') is not None: + ax.add_patch( + v0['patch'], + ) + + # reflectors + if v0.get('reflector') is not None: + ax.plot( + v0['reflector']['x'], + v0['reflector']['y'], + c=v0['reflector']['color'], + lw=v0['reflector']['lw'], + label='reflector', + zorder=100, + ) + + # dogleg + if v0.get('dogleg') is not None: + for ii, (cc, lw) in enumerate([('w', 3), (v0['color'], 1)]): + ax.plot( + v0['dogleg']['x'], + v0['dogleg']['y'], + ls='-', + lw=lw, + color=cc, + zorder=100 + 10*ii, + ) + + # neutrons + if v0.get('neutrons') is not None: + ax.fill( + v0['neutrons']['x'], + v0['neutrons']['y'], + v0['neutrons']['color'], + alpha=0.5, + ) + + # text + if v0.get('text') is not None: + ax.text( + v0['text']['x'], + v0['text']['y'], + v0['text']['str'], + horizontalalignment=v0['text']['horizontalalignment'], + verticalalignment=v0['text']['verticalalignment'], + fontsize=fontsize, + fontweight='bold', + color=v0['color'], + transform=ax.transData, + ) + + # overall text in-vessel + lkin = [k0 for k0 in dsensors.keys() if k0.startswith('in')] + xx = np.mean([dsensors[k0]['text']['x'] for k0 in lkin]) + yy = np.max([dsensors[k0]['text']['y'] for k0 in lkin]) + ax.text( + xx, + yy + 0.5, + 'In-vessel\nsensors', + horizontalalignment='center', + verticalalignment='bottom', + fontsize=fontsize, + fontweight='bold', + color=dsensors[lkin[0]]['color'], + transform=ax.transData, + ) + + # overall text ex-cryostat + lkex = [k0 for k0 in dsensors.keys() if k0.startswith('ex')] + xx = np.mean([dsensors[k0]['text']['x'] for k0 in lkex]) + yy = np.max([dsensors[k0]['text']['y'] for k0 in lkex]) + ax.text( + xx, + yy + 0.5, + 'Ex-cryostat\nsensors', + horizontalalignment='center', + verticalalignment='bottom', + fontsize=fontsize, + fontweight='bold', + color=dsensors['ex1']['color'], + transform=ax.transData, + ) + + # overall text neutrons + ind = np.argmin([ + np.min(dsensors[k0]['neutrons']['y']) for k0 in lkex + ]) + xx = np.max(dsensors[lkex[ind]]['neutrons']['x']) + yy = np.min(dsensors[lkex[ind]]['neutrons']['y']) + ax.text( + xx, + yy - 0.5, + 'neutrons', + horizontalalignment='center', + verticalalignment='bottom', + fontsize=fontsize, + fontweight='bold', + color=dsensors[lkex[ind]]['neutrons']['color'], + transform=ax.transData, + ) + + # ---------------- + # plot theta_vs_B + # ---------------- + + kax = 'theta_vs_B' + if dax.get(kax) is not None: + ax = dax[kax]['handle'] + + for k0, v0 in dsensors.items(): + + ax.plot( + v0['rplasma_norm'], + v0['theta_vs_B'] * 180 / np.pi, + marker=v0.get('marker', '.'), + ms=v0.get('ms', 6), + color=v0['color'], + ls=v0.get('ls', 'None'), + label=v0.get('label', k0), + ) + + # text + if v0.get('text') is not None: + trans = transforms.blended_transform_factory( + ax.transAxes, ax.transData, + ) + ind = v0['rplasma_norm'] > 0.9 + yy = v0['theta_vs_B'][ind] + if k0 == 'ex3': + yy = np.mean(yy[yy < np.pi/4]) + else: + yy = np.mean(yy) + if k0 == 'ex1': + txt = '(1,2)' + else: + txt = v0['text']['str'] + + if k0 != 'ex2': + ax.text( + 1., + yy * 180 / np.pi, + txt, + horizontalalignment='left', + verticalalignment='center', + fontsize=fontsize, + fontweight='bold', + color=v0['color'], + transform=trans, + ) + + ax.axhline(90, c='k', ls='--', lw=1) + ax.set_xlim(-1, 1) + ax.set_ylim(0, 180) + ax.grid(True) + + # comments + ax.text( + -0.97, + 55, + 'forward', + horizontalalignment='left', + verticalalignment='top', + rotation=90, + fontsize=fontsize, + fontweight='bold', + transform=ax.transData, + ) + ax.text( + -0.97, + 110, + 'backward', + horizontalalignment='left', + verticalalignment='bottom', + rotation=90, + fontsize=fontsize, + fontweight='bold', + transform=ax.transData, + ) + + # -------------- + # save + # -------------- + + savefig( + fig=fig, + pfe_save=pfe_save, + path_save=path_save, + file=__file__, + ) + + return dax + + +def _fig02_check( + config=None, + # tokamak + R0=None, + rplasma=None, + RVes=None, + Rcryo=None, + PP_R=None, + PP_width=None, + PP_phi=None, + # unused + **kwdargs, +): + + # ------------------ + # config + # ------------------ + + if config is None: + config = 'SPARC' + + if isinstance(config, str): + config = tf.load_config(config) + + # ------------------ + # Geometry - R + # ------------------ + + lk = ['R0', 'rplasma', 'RVes', 'Rcryo', 'PP_R', 'PP_width', 'PP_phi'] + dinput = { + kk: { + 'data': None, + 'prop': {'color': 'k', 'ls': '-', 'lw': 1, 'label': None} + } + for kk in lk + } + + # R0 + dinput['R0']['data'] = np.r_[float(ds._generic_check._check_var( + R0, 'R0', + types=(float, int), + sign='>0', + default=_DR['R0'], + ))] + dinput['R0']['prop']['ls'] = '--' + + # rplasma + dinput['rplasma']['data'] = np.r_[float(ds._generic_check._check_var( + rplasma, 'rplasma', + types=(float, int), + sign='>0', + default=_DR['rplasma'], + ))] + assert dinput['rplasma']['data'] < dinput['R0']['data'] + + dinput['Rplasma'] = { + 'data': dinput['R0']['data'] + dinput['rplasma']['data']*np.r_[-1, 1], + 'prop': { + 'color': 'r', + 'lw': 1, + 'ls': '-', + 'label': 'plasma', + } + } + + # RVes + if RVes is None: + RVes = _DR['RVes'] + dinput['RVes']['data'] = ds._generic_check._check_flat1darray( + RVes, 'RVes', + dtype=float, + sign='>0', + unique=True, + size=2, + ) + Rlim = dinput['R0']['data'] - dinput['rplasma']['data'] + assert dinput['RVes']['data'][0] < Rlim + Rlim = dinput['R0']['data'] + dinput['rplasma']['data'] + assert dinput['RVes']['data'][1] > Rlim + dinput['RVes']['prop']['lw'] = 2 + + # rplasma + dinput['Rcryo']['data'] = np.r_[float(ds._generic_check._check_var( + Rcryo, 'Rcryo', + types=(float, int), + sign='>0', + default=_DR['Rcryo'], + ))] + assert dinput['Rcryo']['data'] > dinput['RVes']['data'][1] + dinput['Rcryo']['prop']['lw'] = 2 + + # ------------------ + # Geometry - Port plug + # ------------------ + + # PP_R + if PP_R is None: + PP_R = _DR['PP_R'] + dinput['PP_R']['data'] = ds._generic_check._check_flat1darray( + PP_R, 'PP_R', + dtype=float, + sign='>0', + unique=True, + size=2, + ) + Rin = dinput['R0']['data'] + dinput['rplasma']['data'] + assert dinput['PP_R']['data'][0] > Rin + assert dinput['PP_R']['data'][1] > dinput['Rcryo']['data'] + + # PP_width + dinput['PP_width']['data'] = np.r_[float(ds._generic_check._check_var( + PP_width, 'PP_width', + types=(float, int), + sign='>0', + default=_DR['PP_width'], + ))] + + # PP_phi + if PP_phi is None: + PP_phi = _DR['PP_phi'] + PP_phi = ds._generic_check._check_flat1darray( + PP_phi, 'PP_phi', + dtype=float, + unique=True, + ) + dinput['PP_phi']['data'] = np.arctan2(np.sin(PP_phi), np.cos(PP_phi)) + + return config, dinput + + +def _sensors( + res=None, + dsensors=None, + # dinput + dinput=None, + # unused + **kwdargs, +): + + # -------------- + # inputs + # -------------- + + res = ds._generic_check._check_var( + res, 'res', + types=float, + sign='>0', + default=0.01, + ) + + # -------------- + # initialize + # -------------- + + if dsensors is None: + dsensors = copy.deepcopy(_DSENSORS) + + # -------------- + # check + # -------------- + + for k0, v0 in dsensors.items(): + for k1, v1 in v0.items(): + if v1 is None: + dsensors[k0][k1] = _DSENSORS[k0].get(k1) + for k1, v1 in _DSENSORS[k0].items(): + if dsensors[k0].get(k1) is None: + dsensors[k0][k1] = v1 + + # -------------- + # Derive + # -------------- + + for k0, v0 in dsensors.items(): + + phi = dinput['PP_phi']['data'][v0['pp']] + eR = np.r_[np.cos(phi), np.sin(phi)] + ephi = np.r_[-np.sin(phi), np.cos(phi)] + sign = v0["cw"] * 2 - 1 + width = dinput['PP_width']['data'] + length = dinput['PP_R']['data'][1] - dinput['PP_R']['data'][0] + ppc = np.mean(dinput['PP_R']['data']) * eR + + # ---------- + # cent + + if k0 == 'in': + cent = v0['R'] * eR + sign * 0.5 * width * ephi + else: + dphi = np.arctan2(width - v0['width'], length) + eRs = eR * np.cos(dphi) + sign * ephi * np.sin(dphi) + cent = ppc + v0["dist"] * eRs + + # store + dsensors[k0]["dphi"] = dphi + dsensors[k0]["ppc"] = ppc + dsensors[k0]["eRs"] = eRs + + # ---------- + # FOV + + if k0 == 'in': + + R0 = dinput['R0']['data'][0] + rplasma = dinput['rplasma']['data'][0] + + # out + R = R0 + rplasma * v0["rplasma_ratio"] + vect_out = _tangent(cent, R, sign) + + # in + R = R0 - rplasma * v0["rplasma_ratio"] + vect_in = _tangent(cent, R, sign) + + else: + ephis = np.r_[-eRs[1], eRs[0]] + vect_out = ( + (length + v0["dist"]) * (-eRs) + 0.5 * v0['width'] * ephis + ) + vect_in = ( + (length + v0["dist"]) * (-eRs) - 0.5 * v0['width'] * ephis + ) + vect_out = vect_out / np.linalg.norm(vect_out) + vect_in = vect_in / np.linalg.norm(vect_in) + + # ------------ + # reflector in + + if dsensors[k0].get('reflector') is not None: + + ref_angle = dsensors[k0]['reflector']['angle'] + ref_R = dsensors[k0]['reflector']['R'] + ref_width = dsensors[k0]['reflector']['width'] + + ref_kk, ref_isout = _intersect(cent, -eRs, ref_R) + ref_kk = ref_kk[~ref_isout] + + ref_cent = cent + ref_kk * (-eRs) + ref_epar = eRs * np.cos(ref_angle) + ephis * np.sin(ref_angle) + ref_nin = np.r_[-ref_epar[1], ref_epar[0]] + + refx = ref_cent[0] + 0.5 * ref_width * np.r_[-1, 1] * ref_epar[0] + refy = ref_cent[1] + 0.5 * ref_width * np.r_[-1, 1] * ref_epar[1] + + dsensors[k0]["reflector"]['x'] = refx + dsensors[k0]["reflector"]['y'] = refy + + cent_out, vect_out2 = _intersect_line( + cent, vect_out, ref_cent, ref_nin, + ) + cent_in, vect_in2 = _intersect_line( + cent, vect_in, ref_cent, ref_nin, + ) + + # if ref_R > 3 => update cent + if ref_R > 3: + + cent_new = ( + ref_cent + + np.sum((cent - ref_cent) * ref_epar) * ref_epar + - np.sum((cent - ref_cent) * ref_nin) * ref_nin + ) + vect_out2 = vect_out + vect_in2 = vect_in + else: + cent_new = cent + + if False: + msg = ( + f"\nReflector {k0}:\n" + f"\t- ref_R: {ref_R} vs {np.linalg.norm(ref_cent)}\n" + f"\t- ref_angle: {ref_angle*180/np.pi:3.1f} deg\n" + f"\t- ref_width: {ref_width}\n" + f"\t- ref_cent: {ref_cent}\n" + f"\t- ref_nin: {ref_nin}\n" + f"\t- ref_epar: {ref_epar}\n" + f"\t- ref_isout: {ref_isout}\n" + f"\t- ref_kk: {ref_kk}\n" + f"\t- cent: {cent}\n" + f"\t- cent_new: {cent_new}\n" + f"\t- cent_out: {cent_out}\n" + f"\t- cent_in: {cent_in}\n" + f"\t- eRs: {eRs}\n" + f"\t- vect_out: {vect_out}\n" + f"\t- vect_out2: {vect_out2}\n" + f"\t- vect_in: {vect_in}\n" + f"\t- vect_in2: {vect_in2}\n" + ) + print(msg) + else: + cent_out = cent + cent_in = cent + cent_new = None + vect_out2 = vect_out + vect_in2 = vect_in + + dsensors[k0]['cent'] = cent if cent_new is None else cent_new + + # FOV + xx, yy = _FOV( + cent_out, vect_out2, + cent_in, vect_in2, + R0, + rplasma, + cent_new, + ) + + path = mpath.Path(np.array([xx, yy]).T) + + # Sample FOV + DX = np.max(xx) - np.min(xx) + DY = np.max(yy) - np.min(yy) + nptsx = int(DX / res) + nptsy = int(DY / res) + ptsx = np.linspace(np.min(xx), np.max(xx), nptsx) + ptsy = np.linspace(np.min(yy), np.max(yy), nptsy) + ptsx = np.repeat(ptsx[:, None], nptsy, axis=1).ravel() + ptsy = np.repeat(ptsy[None, :], nptsx, axis=0).ravel() + iok = ( + path.contains_points(np.array([ptsx, ptsy]).T) + & (np.hypot(ptsx, ptsy) >= R0 - rplasma) + & (np.hypot(ptsx, ptsy) <= R0 + rplasma) + ) + ptsx = ptsx[iok] + ptsy = ptsy[iok] + + # ------------ + # Angle + + pts_phi = np.arctan2(ptsy, ptsx) + pts_ephi0 = -np.sin(pts_phi) + pts_ephi1 = np.cos(pts_phi) + vect0 = ptsx - cent[0] + vect1 = ptsy - cent[1] + vectn = np.sqrt(vect0**2 + vect1**2) + vect0 = vect0 / vectn + vect1 = vect1 / vectn + theta_vs_B = np.arccos(vect0 * pts_ephi0 + vect1 * pts_ephi1) + + # ------------ + # store + + dsensors[k0]["ptsx"] = ptsx + dsensors[k0]["ptsy"] = ptsy + dsensors[k0]["rplasma_norm"] = ( + (np.hypot(ptsx, ptsy) - R0) / dinput['rplasma']['data'] + ) + dsensors[k0]["theta_vs_B"] = theta_vs_B + dsensors[k0]["path"] = path + + # patch + if k0.startswith('ex'): + xpath = ( + ppc[0] + + 0.5 * v0['width'] * np.r_[-1, -1, 1, 1] * ephis[0] + + 0.55 * length * np.r_[-1, 1, 1, -1] * eRs[0] + ) + ypath = ( + ppc[1] + + 0.5 * v0['width'] * np.r_[-1, -1, 1, 1] * ephis[1] + + 0.55 * length * np.r_[-1, 1, 1, -1] * eRs[1] + ) + + path_patch = mpath.Path(np.array([xpath, ypath]).T) + dsensors[k0]['patch'] = mpatches.PathPatch( + path_patch, + facecolor='w', + alpha=1., + zorder=30, + edgecolor='w', + ) + + # ---------- + # dogleg + + if v0.get('dogleg') is not None: + # frac_h = v0['dogleg']['frac_h'] + frac_v = v0['dogleg']['frac_v'] + xx = np.r_[ + cent[0], + ppc[0] - ephi[0] * frac_v * width/2 - eR[0] * length/2, + ppc[0] - ephi[0] * frac_v * width/2, + ppc[0] + ephi[0] * frac_v * width/2, + ppc[0] + ephi[0] * frac_v * width/2 + eR[0] * 2*length/3, + ] + yy = np.r_[ + cent[1], + ppc[1] - ephi[1] * frac_v * width/2 - eR[1] * length/2, + ppc[1] - ephi[1] * frac_v * width/2, + ppc[1] + ephi[1] * frac_v * width/2, + ppc[1] + ephi[1] * frac_v * width/2 + eR[1] * 2*length/3, + ] + + dsensors[k0]['dogleg']['x'] = xx + dsensors[k0]['dogleg']['y'] = yy + + # ---------- + # neutrons + + if dsensors[k0].get('neutrons_length') is not None: + + poutpp = ppc + 0.5 * np.diff(dinput['PP_R']['data']) * eRs + neutrons_length = dsensors[k0]['neutrons_length'] * eRs + rad = dsensors[k0]['neutrons_width'] + theta = np.linspace(-1, 1, 101) * np.pi / 2 + cos = np.cos(theta) + sin = np.sin(theta) + xx = poutpp[0] + np.r_[ + 0, + neutrons_length[0] + rad * (cos * eRs[0] + sin * ephis[0]), + 0, + ] + yy = poutpp[1] + + np.r_[ + 0, + neutrons_length[1] + rad * (cos * eRs[1] + sin * ephis[1]), + 0, + ] + + dsensors[k0]['neutrons'] = { + 'x': xx, + 'y': yy, + 'color': 'r', + } + + # ---------- + # text + + if v0.get('text') is not None: + if k0.startswith('in'): + dsensors[k0]['text']['x'] = ppc[0] + length * 0.7 * eR[0] + dsensors[k0]['text']['y'] = ppc[1] + length * 0.7 * eR[1] + dsensors[k0]['text']['horizontalalignment'] = 'right' + dsensors[k0]['text']['verticalalignment'] = 'center' + else: + if cent_new is None: + cent_new = cent + dsensors[k0]['text']['x'] = cent_new[0] + length/10 + dsensors[k0]['text']['y'] = cent_new[1] + dsensors[k0]['text']['horizontalalignment'] = 'left' + dsensors[k0]['text']['verticalalignment'] = 'center' + + return dsensors + + +def _tangent( + cent=None, + R=None, + sign=None, +): + + # --------- + # + + phi = np.arctan2(cent[1], cent[0]) + eR = np.r_[np.cos(phi), np.sin(phi)] + ephi = np.r_[-np.sin(phi), np.cos(phi)] + + ang = np.arcsin(R / np.hypot(*cent)) + vect = (-eR) * np.cos(ang) - sign * ephi * np.sin(ang) + vect = vect / np.linalg.norm(vect) + + return vect + + +def _FOV( + cent_out, vect_out, + cent_in, vect_in, + R0, + rplasma, + cent=None, +): + + # ------------------ + # intersect vect_out + # ------------------ + + kk_out_out, isout_out_out = _intersect(cent_out, vect_out, R0 + rplasma) + kk_out_in, isout_out_in = _intersect(cent_out, vect_out, R0 - rplasma) + + kk_out = np.r_[kk_out_out, kk_out_in] + iok_out = np.r_[isout_out_out, ~isout_out_in] + iok = np.isfinite(kk_out) & iok_out + assert iok.sum() >= 1 + kout = np.min(kk_out[iok]) + pt_out = cent_out + kout * vect_out + + # ------------------ + # intersect vect_in + # ------------------ + + kk_in_out, isout_in_out = _intersect(cent_in, vect_in, R0 + rplasma) + kk_in_in, isout_in_in = _intersect(cent_in, vect_in, R0 - rplasma) + + kk_in = np.r_[kk_in_out, kk_in_in] + iok_in = np.r_[isout_in_out, ~isout_in_in] + iok = np.isfinite(kk_in) & iok_in + assert iok.sum() >= 1 + kin = np.min(kk_in[iok]) + pt_in = cent_in + kin * vect_in + + assert np.allclose(np.linalg.norm(pt_out), np.linalg.norm(pt_in)) + Rpts = np.linalg.norm(pt_out) + + # ------------------ + # polyx, polyy + # ------------------ + + # ang + ang_out = np.arctan2(pt_out[1], pt_out[0]) + ang_in = np.arctan2(pt_in[1], pt_in[0]) + ang_min = min(ang_out, ang_in) + ang_max = max(ang_out, ang_in) + if np.abs(ang_min - ang_max) > np.pi: + ang_min, ang_max = ang_max, ang_min + 2*np.pi + ang = np.linspace(ang_min, ang_max, 31)[::-1] + + polyx = np.r_[cent_out[0], Rpts * np.cos(ang), cent_in[0]] + polyy = np.r_[cent_out[1], Rpts * np.sin(ang), cent_in[1]] + + # cent + if cent is not None: + polyx = np.r_[polyx, cent[0], cent_out[0]] + polyy = np.r_[polyy, cent[1], cent_out[1]] + + return polyx, polyy + + +def _intersect(cent, vect, R): + + # ---------- + # kk + + # AM = ku + # R^2 = (OA + AM)^2 = RA^2 + k^2 + 2 ku OA + a = 1 + b = 2 * np.sum(vect * cent) + c = np.sum(cent**2) - R**2 + delta = b**2 - 4 * a * c + + kk = np.full((2,), np.nan) + if delta == 0: + kk[0] = -b / (2*a) + elif delta > 0: + kk = (-b + np.r_[1, -1] * np.sqrt(delta)) / (2 * a) + + # ---------- + # isout + + xx = cent[0] + kk * vect[0] + yy = cent[1] + kk * vect[1] + phi = np.arctan2(yy, xx) + isout = (vect[0] * np.cos(phi) + vect[1] * np.sin(phi)) > 0. + + assert isout.sum() <= 1 + + return kk, isout + + +def _intersect_line(cent, vect, ref, nin): + + # ---------- + # kk + + # cent_M = kk * vect + # (ref_cent + cent_M).nin = 0 + # (ref_cent + kk * vect).nin = 0 + kk = - np.sum((cent - ref) * nin) / np.sum(vect * nin) + + # ---------- + # reflected vector + + vect2 = vect - 2 * np.sum(vect * nin) * nin + vect2 = vect2 / np.linalg.norm(vect2) + + return cent + kk * vect, vect2 diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig07_responsivities.py b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig07_responsivities.py new file mode 100644 index 0000000..058d65b --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig07_responsivities.py @@ -0,0 +1,158 @@ +import os + + +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.gridspec as gridspec +import datastock as ds + + +from ._savefig import main as savefig + + +# ##################################################### +# ##################################################### +# DEFAULTS +# ##################################################### + + +_PATH_HERE = os.path.dirname(__file__) +_PATH_PAPER = os.path.dirname(_PATH_HERE) +_PATH_SAVE = os.path.join(_PATH_PAPER, 'figures') +_PFE_RESPONSIVITIES = os.path.join( + os.path.join(_PATH_PAPER, 'inputs'), + 'responsivities.npz', +) + + +_LRESP = [ + 'bolo', + 'cvd_bare', 'cvd_filter', + 'spectro', + 'mesxr_11_keV', + 'mehxr_60_keV', +] + + +# ##################################################### +# ##################################################### +# main +# ##################################################### + + +def main( + pfe=None, + lw=2, + figsize=(7, 4), + fontsize=14, + # save + path_save=None, + pfe_save=None, + # unused + **kwdargs, +): + """ Plot the spectral responsivities of a series of sensors + + """ + + # -------------- + # load + # -------------- + + if pfe is None: + pfe = _PFE_RESPONSIVITIES + + dresp = { + k0: v0.tolist() + for k0, v0 in np.load(pfe, allow_pickle=True).items() + } + + # -------------- + # prepare axes + # -------------- + + dmargin = { + 'left': 0.11, 'right': 0.97, + 'bottom': 0.12, 'top': 0.98, + 'wspace': 0.25, 'hspace': 0.20, + } + + fig = plt.figure(figsize=figsize) + + gs = gridspec.GridSpec(ncols=1, nrows=1, **dmargin) + dax = {} + + # -------------- + # axes - resp + # -------------- + + ax = fig.add_subplot(gs[0, 0], aspect='auto') + ax.set_xlabel('E (keV)', fontsize=fontsize, fontweight='bold') + ax.set_ylabel('responsivity', fontsize=fontsize, fontweight='bold') + + dax['resp'] = ax + + dax = ds._generic_check._check_dax(dax) + + # -------------- + # plot - resp + # -------------- + + kax = 'resp' + if dax.get(kax) is not None: + ax = dax[kax]['handle'] + + # --------------- + # loop on sensors + + dme = {'mesxr': False, 'mehxr': False, 'cvd': False} + for k0 in _LRESP: + v0 = dresp[k0] + + # resp, color, lab + lk = [kk for kk in dme.keys() if kk in k0] + if len(lk) == 1: + kk = lk[0] + if dme[kk] is False: + dme[kk] = v0.get('color', 'k') + lab = f"{kk} - {v0['responsivity']['units']}" + else: + v0['color'] = dme[kk] + v0['ls'] = '--' + lab = None + else: + lab = f"{k0} - {v0['responsivity']['units']}" + + # lw + if lw is None: + lwi = v0.get('lw', 1) + else: + lwi = lw + + # plot + ax.loglog( + v0['E_eV']['data']*1e-3, + v0['responsivity']['data'], + ls=v0.get('ls', '-'), + lw=lwi, + c=v0.get('color', 'k'), + marker=v0.get('marker', 'None'), + label=lab, + ) + + ax.set_ylim(1e-5, 2) + ax.grid(True) + ax.legend() + + # -------------- + # save + # -------------- + + savefig( + fig=fig, + pfe_save=pfe_save, + path_save=path_save, + file=__file__, + ) + + return dax, dresp diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig08_perfs_single.py b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig08_perfs_single.py new file mode 100644 index 0000000..47b7af9 --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig08_perfs_single.py @@ -0,0 +1,524 @@ +import string + + +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.gridspec as gridspec +import matplotlib.transforms as transforms +import datastock as ds + + +from . import _perfs +from ._savefig import main as savefig + + +# ####################################### +# ####################################### +# DEFAULTS +# ####################################### + + +_DCASES = { + 0: { + 'Te_eV': 1e3, + 'jp_fraction_re': 0.3, + 're': 'avalanche 100 keV', + }, + 1: { + 'Te_eV': 0.2e3, + 'jp_fraction_re': 0.5, + 're': 'avalanche 10 MeV', + }, + 2: { + 'Te_eV': 2e3, + 'jp_fraction_re': 0.8, + 're': 'dreicer', + }, +} + + +_DMARKER = { + 'bolo': { + 'marker': 'None', + 'ls': '-', + }, + 'cvd_bare': { + 'marker': 'None', + 'ls': '--', + }, + 'cvd_filter': { + 'marker': 'None', + 'ls': '-.', + }, + 'spectro': { + 'marker': 'None', + 'ls': '-', + }, + 'mesxr_11_keV': { + 'marker': 'None', + 'ls': '--', + }, + 'mehxr_60_keV': { + 'marker': 'None', + 'ls': '-.', + }, +} + + +_LCOMP = ['maxwell bb', 'maxwell fb', 'maxwell ff', 'RE ff'] +_DCOLOR = { + 'maxwell bb': { + 'label': r'$\epsilon_{bb}^{Max}$', + 'color': 'b', + }, + 'maxwell fb': { + 'label': r'$\epsilon_{fb}^{Max}$', + 'color': 'g', + }, + 'maxwell ff': { + 'label': r'$\epsilon_{ff}^{Max}$', + 'color': 'm', + }, + 'RE ff': { + 'label': r'$\epsilon_{ff}^{RE}$', + 'color': 'r', + }, +} + + +# ####################################### +# ####################################### +# Main +# ####################################### + + +def main( + dmix=None, + # cases + dcases=None, + # d2cross + d2cross_phi=None, + # dist + ne_m3=None, + jp_Am2=None, + # jp_fraction_re=np.linspace(0.025, 0.975, 39), + # plot + figsize=(14, 10), + # figsize=(8, 10), + fontsize=14, + # save + path_save=None, + pfe_save=None, + # unused + **kwdargs, +): + + # ----------- + # inputs + # ----------- + + # dcases + if dcases is None: + dcases = _DCASES + + # -------------- + # compute + # -------------- + + dout = {} + for kcase, vcase in dcases.items(): + ( + demiss_integ, dsignal, ddist, dmix, + total_headon, diff_RE, diff_max, + dang, theta, + lresp, ldist, + ) = _perfs.main( + dmix=dmix, + d2cross_pphi=d2cross_phi, + Te_eV=vcase['Te_eV'], + jp_fraction_re=vcase['jp_fraction_re'], + re=vcase['re'], + ) + + dout[kcase] = { + 'Te_eV': vcase['Te_eV'], + 're': vcase['re'], + 'ddist': ddist, + 'demiss_integ': demiss_integ, + 'dsignal': dsignal, + 'total_headon': total_headon, + 'diff_RE': diff_RE, + 'diff_max': diff_max, + } + + # -------------- + # prepare plot params + # -------------- + + ne = np.unique(ddist['plasma']['ne_m3']['data'])[0] + jp = np.unique(ddist['plasma']['jp_Am2']['data'])[0] + + width = 0.2 + + lk = list(dmix.keys()) + lc = [np.unique(dmix[kk])[0] for kk in lk] + inds = np.argsort(lc)[::-1] + lstr = [f"{lk[ss]} {lc[ss]*100:3.1f} \\%" for ss in inds] + tit_mix = ", ".join(lstr) + + # -------------- + # prepare axes + # -------------- + + dmargin_theta = { + 'left': 0.05, 'right': 0.99, + 'bottom': 0.06, 'top': 0.93, + 'wspace': 0.18, 'hspace': 0.20, + } + dmargin_perf = { + 'left': 0.05, 'right': 0.99, + 'bottom': 0.08, 'top': 0.60, + 'wspace': 0.18, 'hspace': 0.08, + } + + fig = plt.figure(figsize=figsize) + + gs_theta = gridspec.GridSpec(ncols=len(dout), nrows=3, **dmargin_theta) + gs_perf = gridspec.GridSpec(ncols=len(dout), nrows=2, **dmargin_perf) + dax = {} + + # ---------------------- + # loop on cases for axes + # ---------------------- + + ax0_theta = None + ax0_abs = None + ax0_xi = None + for kcase, vcase in dout.items(): + + Teu = np.unique(vcase['ddist']['plasma']['Te_eV']['data']) + indTe = np.argmin(np.abs(Teu - vcase['Te_eV'])) + Te_eV = vcase['ddist']['plasma']['Te_eV']['data'][(0, indTe, 0)] + jp_frac = np.unique( + vcase['ddist']['plasma']['jp_fraction_re']['data'] + )[0] + + # tit + tit = ( + r"$n_e$" + f" = {ne:1.0e}" + r"$/m^3$, " + + r"$j_P$" + f" = {jp*1e-6:1.0f}" + r"$MA/m^2$" + "\n" + + r"$T_e$" + f" = {Te_eV*1e-3:2.1f} keV, " + + r"$F_{RE}$" + f" = {jp_frac:2.1f}\n" + + f"{tit_mix}\n" + + vcase['re'] + ) + + # -------------- + # axes - vs theta + + ax = fig.add_subplot( + gs_theta[0, kcase], + aspect='auto', + sharex=ax0_theta, + sharey=ax0_theta, + ) + ax.set_title( + tit, + fontsize=fontsize, + fontweight='bold', + ) + ax.set_xlabel( + r'$\theta_{ph,B}$ (deg)', + fontsize=fontsize, + fontweight='bold', + ) + if kcase == 0: + ax.set_ylabel( + r"$\epsilon_{ff}^{RE,\eta} / " + r"\max\left(\epsilon_{ff}^{RE,\eta}\right)$", + fontsize=fontsize, + fontweight='bold', + ) + ax0_theta = ax + + ax.text( + 0.01, + 0.99, + f'({string.ascii_lowercase[kcase]})', + horizontalalignment='left', + verticalalignment='top', + fontsize=fontsize, + fontweight='bold', + transform=ax.transAxes, + ) + + dax[f'theta_{kcase}'] = ax + + # -------------- + # axes - Absolute + # -------------- + + ax = fig.add_subplot( + gs_perf[0, kcase], + aspect='auto', + sharex=ax0_abs, + sharey=ax0_abs, + ) + if kcase == 0: + ax.set_ylabel( + r"$\epsilon^{\eta} / \max\left(\epsilon^{\eta}\right)$", + fontsize=fontsize, + fontweight='bold', + ) + ax.set_xlim(0, len(lresp) + 2) + ax0_abs = ax + + ax.text( + 0.01, + 0.99, + f'({string.ascii_lowercase[kcase + len(dout)]})', + horizontalalignment='left', + verticalalignment='top', + fontsize=fontsize, + fontweight='bold', + transform=ax.transAxes, + ) + ax.tick_params(labelbottom=False) + + dax[f'abs_{kcase}'] = ax + + # -------------- + # axes - rel + # -------------- + + ax = fig.add_subplot( + gs_perf[1, kcase], + sharex=ax0_abs, + sharey=ax0_xi, + aspect='auto', + ) + if kcase == 0: + ax.set_ylabel( + r"$\xi = \Delta_{i,ff}^{RE} / M_{i}^{head-on}$", + fontsize=fontsize, + fontweight='bold', + ) + ax0_xi = ax + + ax.text( + 0.01, + 0.99, + f'({string.ascii_lowercase[kcase + 2*len(dout)]})', + horizontalalignment='left', + verticalalignment='top', + fontsize=fontsize, + fontweight='bold', + transform=ax.transAxes, + ) + + dax[f'xi_{kcase}'] = ax + + dax = ds._generic_check._check_dax(dax) + + # ticklabels size + for kax, vax in dax.items(): + dax[kax]['handle'].tick_params( + axis='both', + which='major', + labelsize=fontsize - 1, + ) + + # -------------- + # plot vs theta + # -------------- + + for kcase, vcase in dout.items(): + + kax = f'theta_{kcase}' + if dax.get(kax) is not None: + ax = dax[kax]['handle'] + + # ------------ + # loop on resp + + for kresp in _perfs._LRESP: + + data = vcase['demiss_integ'][kresp]['RE']['ff']['data'] + l0, = ax.plot( + theta*180/np.pi, + data / data.max(), + ls=_DMARKER[kresp]['ls'], + marker=_DMARKER[kresp]['marker'], + color=_perfs._DANGLES[dang[kresp]]['color'], + lw=1, + label=kresp, + ) + + # ------------ + # loop on angles + + for kang, vang in _perfs._DANGLES.items(): + for kk in ['head-on', 'back']: + ax.axvspan( + vang[kk][0]*180/np.pi, + vang[kk][1]*180/np.pi, + facecolor=vang['color'], + alpha=vang['alpha'], + ) + + # ------------ + # decorate + + if kcase == 0: + ax.set_xlim(0, 180) + ax.set_ylim(0, 1) + ax.set_xticks([0, 45, 90, 135, 180]) + ax.grid(True) + elif '10 MeV' in vcase['re']: + ax.legend(loc='upper right', fontsize=fontsize - 2) + + # -------------- + # plot abs + # -------------- + + kax = f'abs_{kcase}' + if dax.get(kax) is not None: + ax = dax[kax]['handle'] + + # ------------ + # loop on resp + + weight_counts = {'head-on': {}, 'back': {}} + + # loop + shape = vcase['dsignal']['bolo']['RE']['ff']['back']['data'].shape + zeros = np.zeros(shape, dtype=float) + for kdir in weight_counts.keys(): + + # detail + for kdist in ldist: + lk = sorted(vcase['dsignal'][lresp[0]][kdist].keys()) + for kemiss in lk: + + data = np.array([ + zeros if (kresp == 'spectro' and kemiss == 'bb') + else vcase['dsignal'][kresp][kdist][kemiss][kdir]['data'].squeeze() + for kresp in lresp + ]) + weight_counts[kdir][f"{kdist} {kemiss}"] = data + + # normalize + for kk, vv in weight_counts[kdir].items(): + weight_counts[kdir][kk] = vv / vcase['total_headon'] + + # -------- + # plot + + for idir, kdir in enumerate(weight_counts.keys()): + bottom = 0 + for kk in _LCOMP: + ax.bar( + np.arange(len(lresp)) + 1 - 0.15 + 0.3*idir, + weight_counts[kdir][kk], + width, + label=_DCOLOR[kk]['label'] if idir == 0 else None, + bottom=bottom, + color=_DCOLOR[kk]['color'], + ) + bottom += weight_counts[kdir][kk] + + # decorate + ax.set_ylim(0, 1) + ax.grid(True) + ax.legend(loc='upper right', fontsize=fontsize - 2) + + # -------------- + # plot rel + # -------------- + + kax = f'xi_{kcase}' + if dax.get(kax) is not None: + ax = dax[kax]['handle'] + + # ---------- + # plot + + # RE + ax.semilogy( + np.arange(len(lresp)) + 1, + vcase['diff_RE'] / vcase['total_headon'], + marker='o', + ms=8, + ls='None', + c='k', + label='RE', + ) + + # maxwell + ax.semilogy( + np.arange(len(lresp)) + 1, + vcase['diff_max'] / vcase['total_headon'], + marker='x', + ms=8, + ls='None', + c='k', + label='maxwellian', + ) + + # vmin + iok = vcase['diff_max'] > 0 + vmin_log10 = np.floor(np.log10(min( + np.min(vcase['diff_RE'] / vcase['total_headon']), + np.min(vcase['diff_max'][iok] / vcase['total_headon'][iok]), + ))) + + # digitizer + trans = transforms.blended_transform_factory( + ax.transAxes, ax.transData, + ) + ddigit = {16: 'r', 12: 'b'} + for nbits, cc in ddigit.items(): + yy = 1/2**nbits + if yy > 10**np.floor(vmin_log10): + ax.axhline(yy, c=cc, ls='--') + ax.text( + 1, + yy, + f'{nbits} bits', + ha='right', + va='bottom', + color=cc, + fontsize=fontsize - 2, + transform=trans, + ) + + # decorate + if kcase == 0: + ax.set_ylim(10**np.floor(vmin_log10), 1) + ax.set_xticks(np.arange(len(lresp)) + 1) + ax.set_xticklabels( + lresp, + rotation=30, + ha="right", + ) + ax.legend(loc='center right') + ax.grid(True) + + # -------------- + # save + # -------------- + + savefig( + fig=fig, + pfe_save=pfe_save, + path_save=path_save, + file=__file__, + ) + + return dax, demiss_integ, dsignal + + +# ####################################### +# ####################################### +# Subroutine +# ####################################### diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig09_perfs_scans.py b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig09_perfs_scans.py new file mode 100644 index 0000000..3d7f0d1 --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig09_perfs_scans.py @@ -0,0 +1,390 @@ +import string + + +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.lines as mlines +import matplotlib.gridspec as gridspec +import datastock as ds + + +from ._load_spect_anis import _JP_FRAC +from . import _perfs +from ._fig08_perfs_single import _DCASES +from ._savefig import main as savefig + + +# ####################################### +# ####################################### +# DEFAULTS +# ####################################### + + +_RE = ['dreicer', 'avalanche 100 keV', 'avalanche 10 MeV'] + + +_DLEVELS = { + 'bolo': { + 'xi': np.r_[1e-6, 1e-5, 5e-5, 1e-4, 1e-3], + 'kappa': np.r_[10, 20, 30, 50, 70, 80], + 'total': np.r_[0.1, 0.3, 0.5, 0.7, 0.9], + }, + 'cvd_bare': { + 'xi': np.r_[1e-5, 5e-5, 1e-4, 5e-4, 1e-3], + 'kappa': np.r_[10, 20, 30, 50, 70, 80, 90], + 'total': np.r_[0.1, 0.3, 0.5, 0.7, 0.9], + }, + 'cvd_filter': { + 'xi': np.r_[1e-4, 1e-3, 1e-2, 1e-1, 0.5], + 'kappa': np.r_[10, 20, 30, 50, 70, 80, 90], + 'total': np.r_[0.1, 0.3, 0.5, 0.7, 0.9], + }, + 'spectro': { + 'xi': 10, + 'kappa': 10, + 'total': np.r_[0.1, 0.3, 0.5, 0.7, 0.9], + }, + 'mesxr_11_keV': { + 'xi': np.r_[0.01, 0.1, 0.2, 0.3], + 'kappa': 10, + 'total': np.r_[0.1, 0.3, 0.5, 0.7, 0.9], + }, + 'mehxr_60_keV': { + 'xi': 10, + 'kappa': 10, + 'total': np.r_[0.1, 0.3, 0.5, 0.7, 0.9], + }, +} + + +# ####################################### +# ####################################### +# Main +# ####################################### + + +def main( + dmix=None, + # d2cross + d2cross_phi=None, + # dist + ne_m3=None, + jp_Am2=None, + # RE + re=None, + dominant=None, + jp_fraction_re=None, + Efield_par_Vm=None, + Ekin_min_eV=None, + Ekin_max_eV=None, + sigmap=None, + pnormW=0, + # plot + figsize=(10, 14), + fontsize=14, + # save + path_save=None, + pfe_save=None, + # unused + **kwdargs, +): + + # ----------- + # inputs + # ----------- + + if jp_fraction_re is None: + jp_fraction_re = _JP_FRAC + + if re is None: + re = _RE + + # -------------- + # compute + # -------------- + + dout = {} + for ii, rei in enumerate(re): + ( + demiss_integ, dsignal, ddist, dmix, + total_headon, diff_RE, diff_max, + dang, theta, + lresp, ldist, + ) = _perfs.main( + dmix=dmix, + ne_m3=ne_m3, + jp_Am2=jp_Am2, + d2cross_phi=d2cross_phi, + re=rei, + jp_fraction_re=jp_fraction_re, + ) + + dout[rei] = { + 'ii': ii, + 'demiss_integ': demiss_integ, + 'dsignal': dsignal, + 'ddist': ddist, + 'total_headon': total_headon, + 'diff_RE': diff_RE, + 'diff_max': diff_max, + } + + # -------------- + # extract + # -------------- + + nresp = len(lresp) + nre = len(re) + + # -------------- + # prepare + # -------------- + + for rei, vout in dout.items(): + # xi range + dout[rei]['xi'] = vout['diff_RE'] / vout['total_headon'] + + # RE vs Maxwell + difftot = vout['diff_max'] + vout['diff_RE'] + dout[rei]['kappa'] = vout['diff_RE'] / difftot + + # -------------- + # prepare Te, F cases + # -------------- + + lkcase = sorted(_DCASES.keys()) + Te_case = np.array([_DCASES[kk]['Te_eV'] for kk in lkcase]) + Fre_case = np.array([_DCASES[kk]['jp_fraction_re'] for kk in lkcase]) + + ne = np.unique(ddist['plasma']['ne_m3']['data'])[0] + jp = np.unique(ddist['plasma']['jp_Am2']['data'])[0] + + lk = list(dmix.keys()) + lc = [np.unique(dmix[kk])[0] for kk in lk] + inds = np.argsort(lc)[::-1] + lstr = [f"{lk[ss]} {lc[ss]*100:3.1f} \\%" for ss in inds] + tit_mix = ", ".join(lstr) + + tit = ( + r"$n_e$" + f" = {ne:1.0e}" + r"$/m^3$, " + + r"$j_P$" + f" = {jp*1e-6:1.0f}" + r"$MA/m^2$" + "\n" + + tit_mix + ) + + # -------------- + # prepare axes + # -------------- + + dmargin = { + 'left': 0.08, 'right': 0.90, + 'bottom': 0.04, 'top': 0.98, + 'wspace': 0.10, 'hspace': 0.10, + } + + fig = plt.figure(figsize=figsize) + fig.suptitle(tit, fontsize=fontsize, fontweight='bold') + + gs = gridspec.GridSpec(ncols=nre, nrows=nresp, **dmargin) + dax = {} + + # -------------- + # axes + # -------------- + + ax0 = None + for ie, rei in enumerate(re): + for iresp, kresp in enumerate(lresp): + + ax = fig.add_subplot( + gs[iresp, ie], + aspect='auto', + sharex=ax0, + sharey=ax0, + ) + if ie == iresp == 0: + ax0 = ax + + # title + if iresp == 0: + ax.set_title( + rei, + fontsize=fontsize, + fontweight='bold', + ) + + # xlabel + if iresp == len(lresp) - 1: + ax.set_xlabel( + r'$T_e$ (keV)', + fontsize=fontsize, + fontweight='bold', + ) + else: + ax.tick_params(labelbottom=False) + + # ylabel + if ie == 0: + ax.set_ylabel( + f"{kresp}\n" + r"$F_{RE}$", + fontsize=fontsize, + fontweight='bold', + ) + + # char + ax.text( + 0.01, + 0.99, + f"({string.ascii_lowercase[iresp + ie * nre]})", + horizontalalignment='left', + verticalalignment='top', + fontsize=fontsize, + fontweight='bold', + transform=ax.transAxes, + ) + + dax[f"{rei}_{kresp}"] = ax + + # check + dax = ds._generic_check._check_dax(dax) + + # -------------- + # plot vs theta + # -------------- + + dcolor = { + 'xi': { + 'color': 'k', + 'label': r'$\xi_{RE} = \frac{\Delta_{ff}^{RE}}{M_i}$', + }, + 'kappa': { + 'color': 'b', + 'label': r"$\kappa = \frac{\Delta_{ff}^{RE}}{\Delta_{ff}^{Max}}$", + }, + 'tot': { + 'color': 'r', + 'label': r"$M_i$", + }, + } + for ie, rei in enumerate(re): + for iresp, kresp in enumerate(lresp): + + kax = f"{rei}_{kresp}" + if dax.get(kax) is not None: + ax = dax[kax]['handle'] + + # ------------- + # xi range + + # set levels xi + if _DLEVELS.get(kresp, {}).get('xi') is not None: + levels = _DLEVELS[kresp]['xi'] + + # plot xi range + sli = (0, slice(None), slice(None)) + cs = ax.contour( + ddist['plasma']['Te_eV']['data'][sli] * 1e-3, + ddist['plasma']['jp_fraction_re']['data'][sli], + dout[rei]['xi'][iresp], + colors=dcolor['xi'], + linestyles='-', + levels=levels, + label='xi', + ) + ax.clabel( + cs, + cs.levels, + fmt=lambda vv: f"{vv:2.1e}", + fontsize=12, + ) + + # ------------- + # kappa range + + # set levels kappa + if _DLEVELS.get(kresp, {}).get('kappa') is not None: + levels = _DLEVELS[kresp]['kappa'] + + # plot kappa + cs = ax.contour( + ddist['plasma']['Te_eV']['data'][sli] * 1e-3, + ddist['plasma']['jp_fraction_re']['data'][sli], + dout[rei]['kappa'][iresp], + colors=dcolor['kappa'], + linestyles='-', + levels=levels, + label='kappa', + ) + ax.clabel(cs, cs.levels, fontsize=12) + + # ------------- + # total range + + # set levels total + if _DLEVELS.get(kresp, {}).get('total') is not None: + levels = _DLEVELS[kresp]['total'] + + # plot total head-on + total = dout[rei]['total_headon'][iresp] + cs = ax.contour( + ddist['plasma']['Te_eV']['data'][sli] * 1e-3, + ddist['plasma']['jp_fraction_re']['data'][sli], + total / np.nanmax(total), + colors=dcolor['total'], + linestyles='-', + levels=levels, + label='total_headon_norm', + ) + ax.clabel(cs, cs.levels, fontsize=12) + + # add case + ax.plot( + Te_case*1e-3, + Fre_case, + ls='None', + marker='x', + ms=6, + color='k', + ) + + # ------------ + # decorate + + if ie == iresp == 0: + ax.set_xlim(0, 2.5) + ax.set_ylim(0, 1) + ax.grid(True) + + if ie == len(re) - 1: + lh = [ + mlines.Line2D( + [], [], + ls='-', + c=cc, + label=dcolor[kk]['label'], + ) + for kk, cc in dcolor.items() + ] + ax.legend( + handles=lh, + loc='upper right', + bbox_to_anchor=(1.2, 1.), + ) + + # -------------- + # save + # -------------- + + savefig( + fig=fig, + pfe_save=pfe_save, + path_save=path_save, + file=__file__, + ) + + return dax, dout + + +# ####################################### +# ####################################### +# Subroutine +# ####################################### diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig11_hxrva_cad.py b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig11_hxrva_cad.py new file mode 100644 index 0000000..c29b7b2 --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig11_hxrva_cad.py @@ -0,0 +1,373 @@ +import os +from PIL import Image + + +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.gridspec as gridspec +import datastock as ds +import tofu as tf + + +from ._load_spect import _PATH_INPUTS +from ._savefig import main as savefig + + +# ####################################### +# ####################################### +# DEFAULTS +# ####################################### + + +_PFE_CAD = os.path.join( + _PATH_INPUTS, + 'HXR_MPP200_zoom.png', +) +_PFE_CONFIG = os.path.join( + _PATH_INPUTS, + 'TFG_Config_ExpSPARC_SPARC-V2_sh00000_Vers1.8.18.npz', +) +_PFE_COLL = os.path.join( + _PATH_INPUTS, + 'Inversion_HXRVA_dvezinet_20260713-154229.npz', +) + + +# ####################################### +# ####################################### +# Main +# ####################################### + + +def main( + coll=None, + key_cam=None, + config=None, + # ptcam + angle0=None, + angle1=None, + # plot + figsize=(5, 4), + # figsize=(8, 10), + fontsize=14, + # save + path_save=None, + pfe_save=None, + # unused + **kwdargs, +): + + # ----------- + # inputs + # ----------- + + # load CAD image + image = Image.open(_PFE_CAD) + + # config + if config is None: + config = _PFE_CONFIG + if isinstance(config, str): + config = tf.load(config) + + # load coll + if coll is None: + coll = _PFE_COLL + if isinstance(coll, str): + coll = tf.data.load(coll) + + dkrays, angle0, angle1 = _add_ptcam( + coll=coll, + key_cam=key_cam, + angle0=angle0, + angle1=angle1, + config=config, + ) + + # -------------- + # compute + # -------------- + + # data + data, _, ref = coll.get_rays_quantity( + key='LEFT_pt', + quantity='alpha', + ) + + # extent + dang0 = angle0[1] - angle0[0] + dang1 = angle1[1] - angle1[0] + extent = ( + (angle0[0] - 0.5*dang0) * 180/np.pi, + (angle0[-1] + 0.5*dang0) * 180/np.pi, + (angle1[0] - 0.5*dang1) * 180/np.pi, + (angle1[-1] + 0.5*dang1) * 180/np.pi, + ) + + # -------------- + # prepare figure + # -------------- + + dmargin = { + 'left': 0.03, 'right': 0.97, + 'bottom': 0.06, 'top': 0.97, + 'wspace': 0.30, 'hspace': 0.30, + } + + fig = plt.figure(figsize=figsize) + gs = gridspec.GridSpec(ncols=2, nrows=15, **dmargin) + dax = {} + + # ---------------------- + # axes + # ---------------------- + + # -------------- + # axes - CAD + + ax = fig.add_subplot(gs[:, 0], aspect='auto') + + ax.text( + 0.0, + 1.0, + '(a)', + horizontalalignment='left', + verticalalignment='bottom', + fontsize=fontsize, + fontweight='bold', + transform=ax.transAxes, + ) + + dax['CAD'] = ax + + # -------------- + # axes - View + + ax = fig.add_subplot(gs[:-3, 1], aspect='equal') + + ax.set_title( + "point-sphere camera view", + fontsize=fontsize, + fontweight='bold', + ) + ax.set_xlabel( + r"$\theta_0$", + fontsize=fontsize, + fontweight='bold', + labelpad=0, + ) + ax.set_ylabel( + r"$\theta_1$", + fontsize=fontsize, + fontweight='bold', + labelpad=-12, + ) + + ax.text( + -0.05, + 1.01, + '(b)', + horizontalalignment='right', + verticalalignment='bottom', + fontsize=fontsize, + fontweight='bold', + transform=ax.transAxes, + ) + + dax['view'] = ax + + # -------------- + # axes - colorbar + + ax = fig.add_subplot(gs[-1, 1], aspect='auto') + ax.set_title( + 'Incidence angle on PFC (deg)', + fontsize=fontsize, + fontweight='bold', + ) + + dax['cbar'] = ax + + dax = ds._generic_check._check_dax(dax) + + # ticklabels size + for kax, vax in dax.items(): + dax[kax]['handle'].tick_params( + axis='both', + which='major', + labelsize=fontsize - 1, + ) + + # -------------- + # plot CAD + # -------------- + + kax = "CAD" + if dax.get(kax) is not None: + ax = dax[kax]['handle'] + + ax.imshow( + image, + ) + + ax.set_axis_off() + + # -------------- + # plot view + # -------------- + + kax = "view" + if dax.get(kax) is not None: + ax = dax[kax]['handle'] + + im = ax.imshow( + data.T*180/np.pi, + extent=extent, + origin='lower', + aspect='equal', + cmap=plt.cm.gray_r, + interpolation='nearest', + vmin=0, + vmax=90, + ) + + plt.colorbar( + im, + cax=dax['cbar']['handle'], + orientation='horizontal', + ) + dax['cbar']['handle'].set_xticks([0, 30, 60, 90]) + + # ---------- + # save + # ---------- + + savefig( + fig=fig, + pfe_save=pfe_save, + path_save=path_save, + file=__file__, + ) + + return dax, coll + + +# ############################################ +# ############################################ +# add ptcam +# ############################################ + + +def _add_ptcam( + coll=None, + key_cam=None, + angle0=None, + angle1=None, + config=None, +): + + # ------- + # inputs + # ------- + + # angle0 + if angle0 is None: + angle0 = (30*np.pi/180) * np.linspace(-1, 1, 180) + + angle0 = ds._generic_check._check_flat1darray( + angle0, "angle0", + dtype=float, + unique=True, + ) + + # angle1 + if angle1 is None: + angle1 = (40*np.pi/180) * np.linspace(-1, 1, 240) + + angle1 = ds._generic_check._check_flat1darray( + angle1, "angle1", + dtype=float, + unique=True, + ) + # ------- + # angles ref + # ------- + + ref_rays = ('nangle0', 'nangle1') + nrays = (angle0.size, angle1.size) + coll.add_ref(ref_rays[0], size=nrays[0]) + coll.add_ref(ref_rays[1], size=nrays[1]) + + # --------- + # angles + # --------- + + coll.add_data( + 'angle0', + data=angle0*180/np.pi, + units='deg', + ref=ref_rays[0], + ) + + coll.add_data( + 'angle1', + data=angle1*180/np.pi, + units='deg', + ref=ref_rays[1], + ) + + # --------------- + # add single points + # --------------- + + if key_cam is None: + key_cam = sorted(coll.dobj['camera'].keys()) + + dkrays = {} + for kcam in key_cam: + + # ------------- + # cent, vect + + cent = np.mean(coll.get_camera_cents_xyz(kcam), axis=1) + phi0 = np.arctan2(cent[1], cent[0]) + ephi0 = np.r_[-np.sin(phi0), np.cos(phi0), 0] + + vect = {} + dvect = coll.get_camera_unit_vectors(kcam) + ls = ['x', 'y', 'z'] + for kv in ['nin', 'e0', 'e1']: + vect[kv] = np.array([dvect[f'{kv}_{ss}'] for ss in ls]) + if vect[kv].ndim > 1: + laxis = range(1, vect[kv].ndim) + vect[kv] = np.mean(vect[kv], axis=tuple(laxis)) + + # adjust e0 + e1 = np.r_[0, 0, 1.] + e0 = np.cross(vect['nin'], e1) + e0 = e0 / np.linalg.norm(e0) + if np.sum(e0 * ephi0) < 0.: + e0 = -e0 + vect['e0'] = e0 + e1 = np.cross(vect['nin'], e0) + e1 = e1 / np.linalg.norm(e1) + if e1[2] < 0: + e1 = -e1 + vect['e1'] = e1 + + # -------------------- + # using tofu built-in + + kray = f"{kcam}_pt" + coll.add_single_point_camera2d( + key=kray, + cent=cent, + angle0='angle0', + angle1='angle1', + config=config, + **vect, + ) + + dkrays[kcam] = kray + + return dkrays, angle0, angle1 diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig12_hxrva_images.py b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig12_hxrva_images.py new file mode 100644 index 0000000..f9795b3 --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig12_hxrva_images.py @@ -0,0 +1,410 @@ +import string + + +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.gridspec as gridspec +import datastock as ds + + +from . import _singlept_sensor +from ._savefig import main as savefig + + +# ####################################### +# ####################################### +# DEFAULTS +# ####################################### + + +_RE = ['avalanche 10 MeV', 'dreicer'] +_JP_FRAC = 0.9 +_TE = 0.5e3 +_T = 5 # 5 + + +_CASES = { + 0: { + 'helicity': False, + 'pitch': True, + 'ne': False, + }, + 1: { + 'helicity': True, + 'pitch': True, + 'ne': False, + }, + 2: { + 'helicity': True, + 'pitch': True, + 'ne': True, + }, +} + + +_FONTSIZE = 14 + + +# ####################################### +# ####################################### +# Main +# ####################################### + + +def main( + # coll + coll=None, + key_cam=None, + config=None, + res=None, + # responsivity + key_resp=None, + # equilibrium + t=None, + # dmix + dmix=None, + # d2cross + d2cross_phi=None, + # dist + ne_m3=None, + jp_Am2=None, + Te_eV=None, + jp_fraction_re=None, + # RE + re=None, + # assumptions + cases=None, + # ptcam + angle0=None, + angle1=None, + # plot params + dvmax=None, + # plot + dmargin=None, + figsize=None, + fontsize=None, + # saving + pfe_save=None, + path_save=None, + # unused + **kwdargs, +): + + # ----------- + # inputs + # ----------- + + if t is None: + t = _T + + if Te_eV is None: + Te_eV = _TE + + if jp_fraction_re is None: + jp_fraction_re = _JP_FRAC + + if re is None: + re = _RE + if isinstance(re, str): + re = [re] + + if cases is None: + cases = _CASES + + if key_resp is None: + key_resp = 'cvd_filter' + + # -------------- + # compute + # -------------- + + ( + coll, config, + dangles, dsig_los, dmetrics, + krays_max, krays_min, + angle0, angle1, + ne, jp, jp_frac, Te_eV, + ) = _singlept_sensor.main( + # coll + coll=coll, + key_cam=key_cam, + config=config, + res=res, + # responsivity + key_resp=key_resp, + # equilibrium + t=t, + # dmix + dmix=dmix, + # d2cross + d2cross_phi=d2cross_phi, + # dist + ne_m3=ne_m3, + jp_Am2=jp_Am2, + Te_eV=Te_eV, + jp_fraction_re=jp_fraction_re, + # RE + re=re, + # assumptions + cases=cases, + # ptcam + angle0=angle0, + angle1=angle1, + ) + + # ------------- + # inputs + # ------------- + + if fontsize is None: + fontsize = _FONTSIZE + + if dvmax is None: + dvmax = {} + + for rei in re: + if dvmax.get(rei) is None: + dvmax[rei] = np.nanmax([ + dsig_los[kcase][krays_max][rei]['RE']['ff']['data'] + for kcase in cases.keys() + ]) + + # ------------- + # prepare + # ------------- + + dang0 = angle0[1] - angle0[0] + dang1 = angle1[1] - angle1[0] + extent = ( + (angle0[0] - 0.5*dang0) * 180/np.pi, + (angle0[-1] + 0.5*dang0) * 180/np.pi, + (angle1[0] - 0.5*dang1) * 180/np.pi, + (angle1[-1] + 0.5*dang1) * 180/np.pi, + ) + + kcase = sorted(cases.keys())[0] + units = dsig_los[kcase][krays_max][re[0]]['RE']['ff']['units'] + + # ------------- + # prepare fig + # ------------- + + if dmargin is None: + dmargin = { + 'left': 0.05, 'right': 0.90, + 'bottom': 0.06, 'top': 0.93, + 'wspace': 0.18, 'hspace': 0.20, + } + dmargin_cbar = { + 'left': 0.93, 'right': 0.97, + 'bottom': 0.06, 'top': 0.93, + 'wspace': 0.18, 'hspace': 0.20, + } + + fig = plt.figure(figsize=figsize) + + nc = np.max([kk for kk in cases.keys()]) + 2 + gs = gridspec.GridSpec(ncols=nc, nrows=len(re), **dmargin) + gs_cbar = gridspec.GridSpec(ncols=1, nrows=len(re), **dmargin_cbar) + dax = {} + + # ---------------------- + # fig title + # ---------------------- + + # tit + tit = ( + r"$n_e$" + f" = {ne:1.0e}" + r"$/m^3$, " + + r"$j_P$" + f" = {jp*1e-6:1.0f}" + r"$MA/m^2$" + "\n" + + r"$T_e$" + f" = {Te_eV*1e-3:2.1f} keV, " + + r"$F_{RE}$" + f" = {jp_frac:2.1f}\n" + + key_resp + ) + fig.suptitle( + tit, + fontsize=fontsize, + fontweight='bold', + ) + + # ---------------------- + # loop on cases for axes + # ---------------------- + + ax0 = None + for ire, rei in enumerate(re): + for icase, (kcase, vcase) in enumerate(cases.items()): + + # tit + hel = 'helicity' if vcase['helicity'] else 'no helicity' + pitch = 'pitch angle distrib.' if vcase['pitch'] else 'no pitch' + nep = 'ne profile' if vcase['ne'] else 'flat profile' + tit = ( + f"{hel}\n" + f"{pitch}\n" + f"{nep}\n" + ) + + # -------------- + # axes - image + + ax = fig.add_subplot( + gs[ire, kcase], + aspect='equal', + sharex=ax0, + sharey=ax0, + ) + if ire == 0: + ax.set_title( + tit, + fontsize=fontsize, + fontweight='bold', + ) + ax.set_xlabel( + r'$\theta_0$ (deg)', + fontsize=fontsize, + fontweight='bold', + ) + if icase == 0: + ax.set_ylabel( + f"{rei}\n" + r"$\theta_1$ (deg)", + fontsize=fontsize, + fontweight='bold', + ) + ax0 = ax + + ax.text( + 0.01, + 0.99, + f'({string.ascii_lowercase[icase + nc*ire]})', + horizontalalignment='left', + verticalalignment='top', + fontsize=fontsize, + fontweight='bold', + transform=ax.transAxes, + ) + + dax[f'{rei}_{kcase}'] = ax + + # --------- + # diff_RE + + ax = fig.add_subplot( + gs[ire, -1], + aspect='equal', + sharex=ax0, + sharey=ax0, + ) + if ire == 0: + ax.set_title( + tit, + fontsize=fontsize, + fontweight='bold', + ) + ax.set_xlabel( + r'$\theta_0$ (deg)', + fontsize=fontsize, + fontweight='bold', + ) + + ax.text( + 0.01, + 0.99, + f'({string.ascii_lowercase[nc - 1 + nc*ire]})', + horizontalalignment='left', + verticalalignment='top', + fontsize=fontsize, + fontweight='bold', + transform=ax.transAxes, + ) + + dax[f'{rei}_diff'] = ax + + # --------- + # colorbar + + ax = fig.add_subplot( + gs_cbar[ire, 0], + aspect='auto', + ) + ax.set_title( + str(units), + fontsize=fontsize, + fontweight='bold', + ) + + dax[f'{rei}_cbar'] = ax + + dax = ds._generic_check._check_dax(dax) + + # ticklabels size + for kax, vax in dax.items(): + dax[kax]['handle'].tick_params( + axis='both', + which='major', + labelsize=fontsize - 1, + ) + + # -------------- + # plot vs theta + # -------------- + + for ire, rei in enumerate(re): + for icase, (kcase, vcase) in enumerate(cases.items()): + + kax = f"{rei}_{kcase}" + if dax.get(kax) is not None: + ax = dax[kax]['handle'] + + data = dsig_los[kcase][krays_max][rei]['RE']['ff']['data'] + + im = ax.imshow( + data.T, + extent=extent, + origin='lower', + cmap=plt.cm.viridis, + interpolation='nearest', + vmin=0, + vmax=dvmax[rei], + ) + + if icase == len(cases) - 1: + ax = dax[f'{rei}_cbar']['handle'] + plt.colorbar(im, cax=ax) + + # --------- + # diff + + kax = f"{rei}_diff" + if dax.get(kax) is not None: + ax = dax[kax]['handle'] + + kcase = [ + k0 for k0, v0 in cases.items() + if all([v1 for v1 in v0.values()]) + ][0] + data = dsig_los[kcase][krays_max][rei]['RE']['ff']['data'] + data_back = dsig_los[kcase][krays_min][rei]['RE']['ff']['data'] + diff = data - data_back[::-1, :] + + im = ax.imshow( + diff.T, + extent=extent, + origin='lower', + cmap=plt.cm.viridis, + interpolation='nearest', + vmin=0, + vmax=dvmax[rei], + ) + + # ---------- + # save + # ---------- + + savefig( + fig=fig, + pfe_save=pfe_save, + path_save=path_save, + file=__file__, + ) + + return dax diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig13_hxrva_metrics.py b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig13_hxrva_metrics.py new file mode 100644 index 0000000..4247cec --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig13_hxrva_metrics.py @@ -0,0 +1,351 @@ +import string + + +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.gridspec as gridspec +import datastock as ds + + +from . import _singlept_sensor +from ._savefig import main as savefig + + +# ####################################### +# ####################################### +# DEFAULTS +# ####################################### + + +_RE = ['avalanche 10 MeV', 'dreicer'] +_JP_FRAC = 0.9 +_TE = 0.5e3 +_T = 5 # 5 + + +_FONTSIZE = 14 + + +# ####################################### +# ####################################### +# Main +# ####################################### + + +def main( + # coll + coll=None, + key_cam=None, + config=None, + res=None, + # responsivity + key_resp=None, + # equilibrium + t=None, + # dmix + dmix=None, + # d2cross + d2cross_phi=None, + # dist + ne_m3=None, + jp_Am2=None, + Te_eV=None, + jp_fraction_re=None, + # RE + re=None, + # ptcam + angle0=None, + angle1=None, + # plot params + dvmax=None, + dvmin=None, + # plot + dmargin=None, + figsize=None, + fontsize=None, + # saving + pfe_save=None, + path_save=None, + # unused + **kwdargs, +): + + # ----------- + # inputs + # ----------- + + if t is None: + t = _T + + if Te_eV is None: + Te_eV = _TE + + if jp_fraction_re is None: + jp_fraction_re = _JP_FRAC + + if re is None: + re = _RE + if isinstance(re, str): + re = [re] + + if key_resp is None: + key_resp = 'cvd_filter' + + # -------------- + # compute + # -------------- + + cases = {0: {'helicity': True, 'pitch': True, 'ne': True}} + ( + coll, config, + dangles, dsig_los, dmetrics, + krays_max, krays_min, + angle0, angle1, + ne, jp, jp_frac, Te_eV, + ) = _singlept_sensor.main( + # coll + coll=coll, + key_cam=key_cam, + config=config, + res=res, + # responsivity + key_resp=key_resp, + # equilibrium + t=t, + # dmix + dmix=dmix, + # d2cross + d2cross_phi=d2cross_phi, + # dist + ne_m3=ne_m3, + jp_Am2=jp_Am2, + Te_eV=Te_eV, + jp_fraction_re=jp_fraction_re, + # RE + re=re, + # assumptions + cases=cases, + # ptcam + angle0=angle0, + angle1=angle1, + ) + + # ------------- + # inputs + # ------------- + + if fontsize is None: + fontsize = _FONTSIZE + + # ------------- + # prepare + # ------------- + + dang0 = angle0[1] - angle0[0] + dang1 = angle1[1] - angle1[0] + extent = ( + (angle0[0] - 0.5*dang0) * 180/np.pi, + (angle0[-1] + 0.5*dang0) * 180/np.pi, + (angle1[0] - 0.5*dang1) * 180/np.pi, + (angle1[-1] + 0.5*dang1) * 180/np.pi, + ) + + kcase = sorted(cases.keys())[0] + units = dsig_los[kcase][krays_max][re[0]]['RE']['ff']['units'] + + # ---------------------- + # loop on cases for axes + # ---------------------- + + dmetrics_plot = { + 0: { + 'key': 'xi', + 'tit': 'xi', + }, + 1: { + 'key': 'kappa', + 'tit': 'kappa', + }, + 2: { + 'key': 'meas_RE_diff', + 'tit': r'$\delta \epsilon_{ff}^{RE}$', + }, + 3: { + 'key': 'meas_headon', + 'tit': r"$\epsilon^{head-on}$", + }, + } + + # ----------- + # vmin, vmax + + if dvmax is None: + dvmax = {} + if dvmin is None: + dvmin = {} + + lm = sorted(dmetrics_plot.keys()) + for km, vm in dmetrics_plot.items(): + if dvmax.get(km) is None: + dvmax[km] = np.nanmax([ + dmetrics[rei][vm['key']]['data'] for rei in re + ]) + if dvmin.get(km) is None: + dvmin[km] = dvmax[km] / 1000 + + # ------------- + # prepare fig + # ------------- + + if dmargin is None: + dmargin = { + 'left': 0.05, 'right': 0.90, + 'bottom': 0.10, 'top': 0.93, + 'wspace': 0.18, 'hspace': 0.20, + } + dmargin_cbar = { + 'left': 0.05, 'right': 0.90, + 'bottom': 0.05, 'top': 0.10, + 'wspace': 0.18, 'hspace': 0.20, + } + + fig = plt.figure(figsize=figsize) + + gs = gridspec.GridSpec(ncols=4, nrows=len(re), **dmargin) + gs_cbar = gridspec.GridSpec(ncols=4, nrows=1, **dmargin_cbar) + dax = {} + + # ---------------------- + # fig title + # ---------------------- + + # tit + tit = ( + r"$n_e$" + f" = {ne:1.0e}" + r"$/m^3$, " + + r"$j_P$" + f" = {jp*1e-6:1.0f}" + r"$MA/m^2$" + "\n" + + r"$T_e$" + f" = {Te_eV*1e-3:2.1f} keV, " + + r"$F_{RE}$" + f" = {jp_frac:2.1f}\n" + + key_resp + ) + fig.suptitle( + tit, + fontsize=fontsize, + fontweight='bold', + ) + + # ---------------------- + # loop on cases for axes + # ---------------------- + + ax0 = None + for ire, rei in enumerate(re): + for im, km in enumerate(lm): + + # -------------- + # axes - image + + ax = fig.add_subplot( + gs[ire, im], + aspect='equal', + sharex=ax0, + sharey=ax0, + ) + if ire == 0: + ax.set_title( + dmetrics_plot[km]['tit'], + fontsize=fontsize, + fontweight='bold', + ) + if ire == len(re) - 1: + ax.set_xlabel( + r'$\theta_0$ (deg)', + fontsize=fontsize, + fontweight='bold', + ) + if im == 0: + ax.set_ylabel( + f"{rei}\n" + r"$\theta_1$ (deg)", + fontsize=fontsize, + fontweight='bold', + ) + ax0 = ax + + ax.text( + 0.01, + 0.99, + f'({string.ascii_lowercase[im + 3*ire]})', + horizontalalignment='left', + verticalalignment='top', + fontsize=fontsize, + fontweight='bold', + transform=ax.transAxes, + ) + + dax[f'{rei}_{km}'] = ax + + # --------- + # colorbar + + if ire == 0: + ax = fig.add_subplot( + gs_cbar[0, im], + aspect='auto', + ) + ax.set_title( + str(units), + fontsize=fontsize, + fontweight='bold', + ) + + dax[f'{km}_cbar'] = ax + + dax = ds._generic_check._check_dax(dax) + + # ticklabels size + for kax, vax in dax.items(): + dax[kax]['handle'].tick_params( + axis='both', + which='major', + labelsize=fontsize - 1, + ) + + # -------------- + # plot vs theta + # -------------- + + for ire, rei in enumerate(re): + for im, km in enumerate(lm): + + kax = f"{rei}_{km}" + if dax.get(kax) is not None: + ax = dax[kax]['handle'] + + data = dmetrics[rei][dmetrics_plot[km]['key']]['data'] + + im = ax.contour( + angle0, + angle1, + data.T, + cmap=plt.cm.viridis, + vmin=dvmin[km], + vmax=dvmax[km], + levels=10, + ) + + # colorbar + if ire == 0: + ax = dax[f'{km}_cbar']['handle'] + plt.colorbar(im, cax=ax) + + # ---------- + # save + # ---------- + + savefig( + fig=fig, + pfe_save=pfe_save, + path_save=path_save, + file=__file__, + ) + + return dax, dmetrics diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig14_bolo_cad.py b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig14_bolo_cad.py new file mode 100644 index 0000000..987f619 --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig14_bolo_cad.py @@ -0,0 +1,527 @@ +import os +from PIL import Image + + +import numpy as np +import matplotlib.pyplot as plt +import matplotlib.gridspec as gridspec +import datastock as ds +import tofu as tf + + +from ._load_spect import _PATH_INPUTS +from ._savefig import main as savefig + + +# ####################################### +# ####################################### +# DEFAULTS +# ####################################### + + +_PFE_CAD0 = os.path.join( + _PATH_INPUTS, + 'OMP000U.png', +) +_PFE_CAD1 = os.path.join( + _PATH_INPUTS, + 'OMP000U_Zoom.png', +) +_PFE_CAD2 = os.path.join( + _PATH_INPUTS, + 'OMP000U_Zoom_inside.png', +) +_PFE_CONFIG = os.path.join( + _PATH_INPUTS, + 'TFG_Config_ExpSPARC_SPARC-V2_sh00000_Vers1.8.18.npz', +) +_PFE_COLL = os.path.join( + _PATH_INPUTS, + 'Inversion_BOLO_DMS_ptcam_dvezinet_20260709-205535.npz', +) + + +_DSAMPLING = { + 'dedge': {'res': 'min'}, + 'dsurface': {'nb': 7}, +} + + +# ####################################### +# ####################################### +# Main +# ####################################### + + +def main( + coll=None, + key_cam=None, + config=None, + # ptcam + angle0=None, + angle1=None, + # fov + dsampling=None, + # plot + figsize=(5, 7), + # figsize=(8, 10), + fontsize=14, + # save + path_save=None, + pfe_save=None, + # unused + **kwdargs, +): + + # ----------- + # inputs + # ----------- + + if key_cam is None: + key_cam = 'IDIV_000U' + + if dsampling is None: + dsampling = _DSAMPLING + + # load CAD image + image0 = Image.open(_PFE_CAD0) + image1 = Image.open(_PFE_CAD1) + image2 = Image.open(_PFE_CAD2) + + # config + if config is None: + config = _PFE_CONFIG + if isinstance(config, str): + config = tf.load(config) + + # load coll + if coll is None: + coll = _PFE_COLL + if isinstance(coll, str): + coll = tf.data.load(coll) + + dkrays, angle0, angle1 = _add_ptcam( + coll=coll, + key_cam=key_cam, + angle0=angle0, + angle1=angle1, + config=config, + ) + + wdiag = coll._which_diagnostic + kdiag = [kd for kd in coll.dobj[wdiag].keys() if 'bolo' in kd.lower()][0] + + # ------------- + # add FOV + # ------------- + + # make rays + kcam = list(dkrays.keys())[0] + coll.add_rays_from_diagnostic( + key=kdiag, + dsampling_pixel=dsampling, + dsampling_optics=dsampling, + optics=-1, + config=config, + store=True, + strict=False, + key_rays=None, + overwrite=None, + ) + key_rays = f"{kcam.replace('_los', '').replace('BOLO_', '')}_rays" + + # FOV + key = list(dkrays.values())[0] + wrays = coll._which_rays + shape = coll.dobj[wrays][key]['shape'][1:] + ind_fov = np.zeros(shape, dtype=bool) + dout_rays = coll.get_rays_angles_from_single_point_camera2d( + key_single_pt_cam=key, + key_rays=key_rays, + return_indices=True, + convex_axis=(-1, -2), + ) + for k1, v1 in dout_rays['hull'].items(): + ind_fov = ind_fov | v1 + + # -------------- + # compute + # -------------- + + # data + assert len(dkrays) == 1 + kray = list(dkrays.values())[0] + data, _, ref = coll.get_rays_quantity( + key=kray, + quantity='alpha', + ) + + # colormap + colormap = np.zeros(data.shape[1:] + (4,), dtype=float) + colormap[..., -1] = data[0, ...] / (np.pi/2) + colormap[(ind_fov, 0)] = 1. + + # extent + dang0 = angle0[1] - angle0[0] + dang1 = angle1[1] - angle1[0] + extent = ( + (angle0[0] - 0.5*dang0) * 180/np.pi, + (angle0[-1] + 0.5*dang0) * 180/np.pi, + (angle1[0] - 0.5*dang1) * 180/np.pi, + (angle1[-1] + 0.5*dang1) * 180/np.pi, + ) + + # -------------- + # prepare figure + # -------------- + + dmargin = { + 'left': 0.04, 'right': 0.90, + 'bottom': 0.07, 'top': 0.99, + 'wspace': 0.18, 'hspace': 0.20, + } + + n0 = 5 + n1 = 2*n0 - 2 + nc = 8 + fig = plt.figure(figsize=figsize) + gs = gridspec.GridSpec(ncols=2*nc, nrows=(2*n0 + 1 + n1), **dmargin) + dax = {} + + # ---------------------- + # axes + # ---------------------- + + # -------------- + # axes - CAD0 + + ax = fig.add_subplot(gs[:2*n0, :nc], aspect='auto') + + ax.text( + 0., + 1., + '(a)', + horizontalalignment='left', + verticalalignment='bottom', + fontsize=fontsize, + fontweight='bold', + transform=ax.transAxes, + ) + + dax['CAD0'] = ax + + # -------------- + # axes - CAD1 + + ax = fig.add_subplot(gs[:n0, nc:], aspect='auto') + + ax.text( + 0., + 1., + '(b)', + horizontalalignment='left', + verticalalignment='bottom', + fontsize=fontsize, + fontweight='bold', + transform=ax.transAxes, + ) + + dax['CAD1'] = ax + + # -------------- + # axes - CAD2 + + ax = fig.add_subplot(gs[n0:2*n0, nc:], aspect='auto') + + ax.text( + 0., + 1., + '(c)', + horizontalalignment='left', + verticalalignment='bottom', + fontsize=fontsize, + fontweight='bold', + transform=ax.transAxes, + ) + + dax['CAD2'] = ax + + # -------------- + # axes - View + + ax = fig.add_subplot(gs[2*n0+1:, :-2], aspect='equal') + + ax.set_title( + "view from a point-sphere camera model", + fontsize=fontsize, + fontweight='bold', + ) + ax.set_xlabel( + r"$\theta_0$", + fontsize=fontsize, + fontweight='bold', + ) + ax.set_ylabel( + r"$\theta_1$", + fontsize=fontsize, + fontweight='bold', + ) + + ax.text( + -0.15, + 1.0, + '(d)', + horizontalalignment='right', + verticalalignment='bottom', + fontsize=fontsize, + fontweight='bold', + transform=ax.transAxes, + ) + + dax['view'] = ax + + # -------------- + # axes - colorbar + + ax = fig.add_subplot(gs[2*n0+1:, -2], aspect='auto') + dax['cbar'] = ax + + dax = ds._generic_check._check_dax(dax) + + # ticklabels size + for kax, vax in dax.items(): + dax[kax]['handle'].tick_params( + axis='both', + which='major', + labelsize=fontsize - 1, + ) + + # -------------- + # plot CAD + # -------------- + + kax = "CAD0" + if dax.get(kax) is not None: + ax = dax[kax]['handle'] + + ax.imshow(image0) + ax.set_axis_off() + + kax = "CAD1" + if dax.get(kax) is not None: + ax = dax[kax]['handle'] + + ax.imshow(image1) + ax.set_axis_off() + + kax = "CAD2" + if dax.get(kax) is not None: + ax = dax[kax]['handle'] + + ax.imshow(image2) + ax.set_axis_off() + + # -------------- + # plot view + # -------------- + + kax = "view" + if dax.get(kax) is not None: + ax = dax[kax]['handle'] + + im = ax.imshow( + np.swapaxes(colormap, 0, 1), + extent=extent, + origin='lower', + aspect='equal', + cmap=plt.cm.gray_r, + # interpolation='nearest', + # vmin=0, + # vmax=1, + ) + + # cax + cax = dax['cbar']['handle'] + plt.colorbar( + im, + cax=cax, + orientation='vertical', + ) + + + deg_ticks = np.r_[0, 30, 60, 90] + cax.set_yticks(deg_ticks/90.) + cax.set_yticklabels(deg_ticks) + cax.set_ylabel( + 'Incidence angle on PFC (deg)', + fontsize=fontsize, + fontweight='bold', + ) + cax.yaxis.set_label_position("right") + cax.yaxis.tick_right() + + # ---------- + # save + # ---------- + + savefig( + fig=fig, + pfe_save=pfe_save, + path_save=path_save, + file=__file__, + ) + + return dax, coll + + +# ############################################ +# ############################################ +# add ptcam +# ############################################ + + +def _add_ptcam( + coll=None, + key_cam=None, + angle0=None, + angle1=None, + config=None, +): + + # ------- + # inputs + # ------- + + # angle0 + if angle0 is None: + angle0 = (np.pi/180) * np.linspace(-50, 30, 240) + + angle0 = ds._generic_check._check_flat1darray( + angle0, "angle0", + dtype=float, + unique=True, + ) + + # angle1 + if angle1 is None: + angle1 = (40*np.pi/180) * np.linspace(-1, 1, 240) + + angle1 = ds._generic_check._check_flat1darray( + angle1, "angle1", + dtype=float, + unique=True, + ) + + # isBOLO + wdiag = coll._which_diagnostic + isbolo = any(['bolo' in kd.lower() for kd in coll.dobj[wdiag].keys()]) + + # ------- + # angles ref + # ------- + + ref_rays = ('nangle0', 'nangle1') + nrays = (angle0.size, angle1.size) + coll.add_ref(ref_rays[0], size=nrays[0]) + coll.add_ref(ref_rays[1], size=nrays[1]) + + # --------- + # angles + # --------- + + coll.add_data( + 'angle0', + data=angle0*180/np.pi, + units='deg', + ref=ref_rays[0], + ) + + coll.add_data( + 'angle1', + data=angle1*180/np.pi, + units='deg', + ref=ref_rays[1], + ) + + # --------------- + # add single points + # --------------- + + wcam = coll._which_cam + if key_cam is None: + key_cam = 'IDIV_000U' + if isinstance(key_cam, str): + key_cam = [key_cam] + + + dkrays = {} + for kcam in key_cam: + + # -------------------- + # HXRVA => cent, vect + + if isbolo is True: + + wrays = coll._which_rays + if kcam not in coll.dobj[wrays].keys(): + if kcam in coll.dobj[wcam].keys(): + kdiag = [ + kd for kd, vd in coll.dobj[wdiag].items() + if kcam in vd[wcam] + ][0] + kcam = coll.dobj[wdiag][kdiag]['doptics'][kcam]['los'] + + kray = f"{kcam}_pt" + coll.add_single_point_camera2d( + key=kray, + key_rays=kcam, + angle0='angle0', + angle1='angle1', + e1=np.r_[0, 0, 1], + config=config, + ) + + else: + + cent = np.mean(coll.get_camera_cents_xyz(kcam), axis=1) + phi0 = np.arctan2(cent[1], cent[0]) + ephi0 = np.r_[-np.sin(phi0), np.cos(phi0), 0] + + vect = {} + dvect = coll.get_camera_unit_vectors(kcam) + ls = ['x', 'y', 'z'] + for kv in ['nin', 'e0', 'e1']: + vect[kv] = np.array([dvect[f'{kv}_{ss}'] for ss in ls]) + if vect[kv].ndim > 1: + laxis = range(1, vect[kv].ndim) + vect[kv] = np.mean(vect[kv], axis=tuple(laxis)) + + # adjust e0 + e1 = np.r_[0, 0, 1.] + e0 = np.cross(vect['nin'], e1) + e0 = e0 / np.linalg.norm(e0) + if np.sum(e0 * ephi0) < 0.: + e0 = -e0 + vect['e0'] = e0 + e1 = np.cross(vect['nin'], e0) + e1 = e1 / np.linalg.norm(e1) + if e1[2] < 0: + e1 = -e1 + vect['e1'] = e1 + + # -------------------- + # using tofu built-in + + kray = f"{kcam}_pt" + coll.add_single_point_camera2d( + key=kray, + cent=cent, + angle0='angle0', + angle1='angle1', + config=config, + **vect, + ) + + dkrays[kcam] = kray + + return dkrays, angle0, angle1 diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_load_spect.py b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_load_spect.py new file mode 100644 index 0000000..589d5e1 --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_load_spect.py @@ -0,0 +1,397 @@ +import os + + +import numpy as np +import scipy.constants as scpct +import datastock as ds + + +# ##################################################### +# ##################################################### +# DEFAULTS +# ##################################################### + + +# PATHS +_PATH_HERE = os.path.dirname(__file__) +_PATH_PAPER = os.path.dirname(_PATH_HERE) +_PATH_INPUTS = os.path.join(_PATH_PAPER, 'inputs') + + +# SPECTRAL MODELLING FILES +_LTYPES = ['SCRAM86', 'CHIANTI', 'FLYCHK'] + + +# PLASMA +_NE = 1e19 # /m3 + + +# ##################################################### +# ##################################################### +# Main +# ##################################################### + + +def main( + dmix=None, + ne_m3=None, + data_source=None, +): + """ Extract ff, fb, bb emissivities from SCRAM / FLYCHK / CHIANTY data + + Uses user-provided ne_m3 and concentrations (via dmix) + Uses the Te values foudn in the source files + + Returns dplasma: + { + 'concentration': {'H': float / array, 'Ni': float / array, ...}, + 'emiss': { + 'H': { + 'ff': {'data': array, 'units': str}}, + 'fb': {'data': array, 'units': str}}, + 'bb': {'data': array, 'units': str}}, + }, + ... + 'Ni': { + 'ff': {'data': array, 'units': str}}, + 'fb': {'data': array, 'units': str}}, + 'bb': {'data': array, 'units': str}}, + }, + }, + 'emiss_tot': { + 'ff': {'data': array, 'units': str}}, + 'fb': {'data': array, 'units': str}}, + 'bb': {'data': array, 'units': str}}, + }, + 'common': { + 'Te': {'data': array, 'units': str}, + 'ne': {'data': array, 'units': str}, + 'E_photon': {'data': array, 'units': str}, + }, + } + + """ + + # -------------- + # inputs + # -------------- + + # data_source + lok = ['SCRAM86', 'CHIANTI', 'FLYCHK'] + data_source = ds._generic_check._check_var( + data_source, 'data_source', + default=lok[0], + allowed=lok, + ) + + # lfpe_spect + lpfe_spect = [ + ff for ff in os.listdir(_PATH_INPUTS) + if ff.endswith('_data.npz') + and f"_{data_source}_".lower() in ff.lower() + ] + + # dpfe_spect + dpfe_spect = { + ff.split('_')[-2]: os.path.join(_PATH_INPUTS, ff) + for ff in lpfe_spect + } + + # dmix + dmix = _check_mix(dmix, dpfe_spect=dpfe_spect) + assert all([kk in dpfe_spect.keys() for kk in dmix.keys()]) + + # -------------- + # load files + # -------------- + + dfiles = { + k0: np.load(dpfe_spect[k0], allow_pickle=True)['arr_0'].tolist() + for k0 in dmix.keys() + } + + # --------------- + # Extract common + # --------------- + + # Extract Te + lk = list(dmix.keys()) + lc = ['Te', 'E_photon'] + dcommon = {cc: None for cc in lc} + for cc in lc: + ref = dfiles[lk[0]][cc]['data'] + rshape = ref.shape + for kk, vv in dfiles.items(): + cc_shape = vv[cc]['data'].shape + if vv[cc]['data'].shape != ref.shape: + msg = ( + "Different shapes!\n" + f"\t- dfiles['{kk}']['{cc}']['data'].shape: {cc_shape}\n" + f"\t- dfiles['{lk[0]}']['{cc}']['data'].shape: {rshape}\n" + ) + raise Exception(msg) + if not np.allclose(vv[cc]['data'], ref): + msg = ( + f"Common field '{cc}' not constant accross files!\n" + f"\t- {ref}\n" + f"\t- {vv[cc]['data']}\n" + ) + raise Exception(msg) + dcommon[cc] = dfiles[lk[0]][cc] + + # ne_m3 vs Te + ne_m3, dcommon['Te']['data'] = _check_neTe(ne_m3, dcommon['Te']['data']) + dcommon['ne'] = { + 'data': ne_m3, + 'units': '1/m3', + } + + # -------------- + # shapes + # -------------- + + shape_Te = ne_m3.shape + shape_conc = dmix[lk[0]].shape + try: + shape_plasma = np.broadcast_shapes(shape_Te, shape_conc) + ne_m3 = np.broadcast_to(ne_m3, shape_plasma) + for k0, v0 in dmix.items(): + dmix[k0] = np.broadcast_to(v0, shape_plasma) + except Exception: + msg = "Te and concentrations should be broadcastcable!" + raise Exception(msg) + + shape_spect = dcommon['E_photon']['data'].shape + shape_emiss = shape_plasma + shape_spect + + sli_E = (None,) * len(shape_plasma) + (slice(None),) + sli_ne = (slice(None),) * len(shape_plasma) + (None,) + + # -------------- + # Zeff + # -------------- + + # Zeff + Zeff = np.zeros(ne_m3.shape) + sli_Z = (None,) * len(shape_plasma) + (slice(None),) + sli_cc = (slice(None),) * len(shape_plasma) + (None,) + for k0, v0 in dfiles.items(): + zz = np.arange(0, v0['Xz']['data'].shape[-1]) + + fz = v0['Xz']['data'] + if k0 == 'H': + fz[:, 1] = 1. + + sumfz = np.sum(fz, axis=-1) + if not np.allclose(sumfz, 1, atol=0.01, rtol=0): + lstr = [ + f"\t- Te = {te*1e-3:4.3f} keV: sum(fz) = {sfz}" + for te, sfz in zip(v0['Te']['data'], sumfz) + ] + msg = ( + f"Impurity '{k0}' has wrong fractional abundances:\n" + f"From file: {dpfe_spect[k0]}\n" + + "\n".join(lstr) + ) + raise Exception(msg) + Zeff += np.sum(dmix[k0][sli_cc] * fz * zz[sli_Z]**2, axis=-1) + + dcommon['Zeff'] = {'data': Zeff, 'units': None} + + # -------------- + # dplasma + # -------------- + + lemiss = ['ff', 'fb', 'bb'] + dplasma = { + 'concentration': dmix, + 'emiss': {k0: {cc: None for cc in lemiss} for k0 in dmix.keys()}, + 'common': dcommon, + 'Zmean': {k0: v0[''] for k0, v0 in dfiles.items()}, + } + + # -------------- + # load elements + # -------------- + + units0 = 'J*cm^3/s/eV/atom/electron' + units = "1 / (m3.s.eV.sr)" + E_ph = dcommon['E_photon']['data'] + for k0 in dmix.keys(): + + # find emissivity key + for ke in lemiss: + lk = [kk for kk in dfiles[k0].keys() if f"emis_{ke}" in kk] + if len(lk) == 1: + kemiss = lk[0] + elif len(lk) > 1: + kchianti = [kk for kk in lk if 'chianti' in kk.lower()] + if len(kchianti) > 0: + kemiss = kchianti[0] + else: + msg = "Not sure how to choose between {lk}" + raise Exception(msg) + else: + msg = f"key 'emis_{ke}' not found in {k0}" + raise Exception(msg) + + # check units + uu = str(dfiles[k0][kemiss]['units']).replace('$', '') + assert uu == units0 + + # concentration + cc = np.broadcast_to(dplasma['concentration'][k0], shape_plasma) + + # emiss + emiss = dfiles[k0][kemiss]['data'].squeeze() + + # ph / m3 / s / eV / sr + emiss = ( + np.broadcast_to(emiss, shape_emiss) + * 1e-6 # cm3 => m3 + / (E_ph[sli_E] * scpct.e) # J => ph + * ne_m3[sli_ne]**2 # /electron => /m3 + * cc[sli_ne] # /atom => /m3 + / (4*np.pi) # => /sr + ) + + dplasma['emiss'][k0][ke] = { + 'data': emiss, + 'units': units, + 'ref': 'plasma.shape + (nE,)', + } + + # ----------- + # total + # ----------- + + dplasma['emiss_tot'] = { + ke: { + 'data': np.sum( + [dplasma['emiss'][k0][ke]['data'] for k0 in dmix.keys()], + axis=0, + ), + 'units': units, + } + for ke in lemiss + } + + return dplasma, dmix, data_source + + +# ##################################################### +# ##################################################### +# Check +# ##################################################### + + +def _check_mix( + dmix=None, + dpfe_spect=None, +): + + # ----------- + # dmix + # ----------- + + # None + if dmix is None: + dmix = 'H' + + # str + if isinstance(dmix, str): + dmix = {dmix: 1} + + # --------- + # each dict + + dfail = {} + lok = sorted(dpfe_spect.keys()) + if isinstance(dmix, dict): + + # check each element + for k0, v0 in dmix.items(): + + # key + if not (isinstance(k0, str) and k0 in lok): + dfail[k0] = f"key not allowed {lok}" + continue + + # value + v0 = np.atleast_1d(v0) + iok = np.isfinite(v0) + iok[iok] = (v0[iok] >= 0) & (v0[iok] <= 1) + if not np.all(iok): + dfail[k0] = "non-finite or values not in [0, 1]" + continue + dmix[k0] = v0 + + else: + dfail['dmix'] = "not a dict!" + + # ------------------- + # overall consistency + + if len(dfail) == 0 and len(dmix) > 1: + + # broadcastable + try: + shape = np.broadcast_shapes(*[v0.shape for v0 in dmix.values()]) + + # broadcast + for k0, v0 in dmix.items(): + dmix[k0] = np.broadcast_to(v0, shape) + + except Exception: + dfail['dmisc'] = "values not broadcastable!" + + # ------------- + # sum = 1 + + total = np.sum([v0 for v0 in dmix.values()], axis=0) + for k0, v0 in dmix.items(): + dmix[k0] = v0 / total + + # ----------- + # any error + + if len(dfail) > 0: + lstr = [f"\t- {k0}: {v0}" for k0, v0 in dfail.items()] + msg = ( + "\nArg 'dmix' must be a dict with of the form:\n" + "\t- {'H': cH, 'Ni': cNi, ...}\n" + "where:\n" + "\t- each key is an element\n" + "\t- each value is a concentration float / array in [0, 1]\n" + "All values must be broadcastable together\n" + "All values will be normalized such that their sum = 1\n" + "\nIdentified issues:\n" + + "\n".join(lstr) + + f"\nProvided:\n{dmix}\n" + ) + raise Exception(msg) + + return dmix + + +def _check_neTe( + ne_m3=None, + Te=None, +): + + # ----------- + # ne, Te + # ----------- + + # ne + if ne_m3 is None: + ne_m3 = _NE + ne_m3 = np.atleast_1d(ne_m3) + + # Te + Te = np.atleast_1d(Te) + + # ------------- + # broadcastable + + return np.broadcast_arrays(ne_m3, Te) diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_load_spect_anis.py b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_load_spect_anis.py new file mode 100644 index 0000000..d2c844a --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_load_spect_anis.py @@ -0,0 +1,364 @@ +import os + + +import numpy as np +import scipy.interpolate as scpinterp +import astropy.units as asunits +import tofu as tf + + +from ._load_spect import main as load_spect +from ._fig02_dist import _DDIST + + +tfphysemis = tf.physics_tools.electrons.emission + + +# ##################################################### +# ##################################################### +# DEFAULTS +# ##################################################### + + +# PATHS +_PATH_HERE = os.path.dirname(__file__) +_PATH_PAPER = os.path.dirname(_PATH_HERE) +_PATH_INPUTS = os.path.join(_PATH_PAPER, 'inputs') + + +def _D2CROSS_PHI_FN(nE, ntheta): + return ( + f"d2cross_phi_Ee01eV-100MeV-{nE}log_Eph1eV-100MeV-{nE+1}log" + f"_nthetaph{ntheta+1}_nthetae0{ntheta}_EH.npz" + ) + + +_DPFE_D2CROSS_PHI = { + 'EH0': os.path.join(_PATH_INPUTS, _D2CROSS_PHI_FN(80, 60)), + 'EH1': os.path.join(_PATH_INPUTS, _D2CROSS_PHI_FN(240, 60)), +} +_D2CROSS_PHI = 'EH0' + + +_TE = 1e3 * np.linspace(0.1, 2.5, 25) +_JP_FRAC = np.linspace(0.1, 0.9, 9) +# _JP_FRAC = np.linspace(0.025, 0.975, 39) +_EKIN_MAX_EV = np.r_[100e3, 10e6] +_PNORMW = np.r_[0.1, 5] + + +# ##################################################### +# ##################################################### +# Main +# ##################################################### + + +def main( + dmix=None, + ne_m3=None, + # ddist + dominant=None, + Efield_par_Vm=None, # pcrit + jp_Am2=None, + jp_fraction_re=None, + Ekin_max_eV=None, + Ekin_min_eV=None, + sigmap=None, + pnormW=None, + # intergation method + integration=None, + # d2cross_phi + d2cross_phi=None, + # data_source + data_source=None, +): + """ For a given impurity mix and ne, Te, return a dict of emissivities + + Emissivities are broken down into 2 distributions + - Maxwell: + - ff: anisotro'upper + - fb: isotropic (from SCRAM / FLYCHK / CHIANTI) + - bb: isotropic + - RE: + - ff: anisotropic + + Returns dplasma: + { + 'concentration': {'H': float / array, 'Ni': float / array, ...}, + 'emiss_tot': { + 'maxwell': { + 'ff': {'data': array, 'units': str}}, + 'fb': {'data': array, 'units': str}}, + 'bb': {'data': array, 'units': str}}, + }, + 'RE': { + 'ff': {'data': array, 'units': str}}, + }, + }, + 'common': { + 'Te': {'data': array, 'units': str}, + 'ne': {'data': array, 'units': str}, + 'E_ph': {'data': array, 'units': str}, + 'theta_ph': {'data': array, 'units': str}, + }, + } + + """ + + # -------------- + # isotropic + # -------------- + + # load + dplasma, dmix, data_source = load_spect( + dmix=dmix, + ne_m3=ne_m3, + # data_source + data_source=data_source, + ) + + # extract + + # --------------------- + # extract Te and ddist + # --------------------- + + Te_eV = dplasma['common']['Te']['data'] + Z_eff = dplasma['common']['Zeff']['data'] + + # ------------ + # inputs + # ------------ + + if jp_fraction_re is None: + jp_fraction_re = _JP_FRAC + jp_fraction_re = np.atleast_1d(jp_fraction_re) + + if Ekin_max_eV is None: + Ekin_max_eV = _EKIN_MAX_EV + Ekin_max_eV = np.atleast_1d(Ekin_max_eV) + + if pnormW is None: + pnormW = _PNORMW + pnormW = np.atleast_1d(pnormW) + + # ------------ + # ddist + # ------------ + + ddist = locals() + ddist = { + kk: vv if ddist.get(kk) is None else ddist[kk] + for kk, vv in _DDIST.items() + if kk not in ['E_eV', 'theta'] + } + if ne_m3 is not None: + ddist['ne_m3'] = ne_m3 + + if Ekin_max_eV is not None: + ddist['Ekin_max_eV'] = Ekin_max_eV + + if pnormW is not None: + ddist['pnormW'] = pnormW + + if Efield_par_Vm is not None: + ddist['Efield_par_Vm'] = Efield_par_Vm + + if Ekin_min_eV is not None: + ddist['Ekin_min_eV'] = Ekin_min_eV + + if sigmap is not None: + ddist['sigmap'] = sigmap + + # shape + ddist['Ekin_max_eV'] = ddist['Ekin_max_eV'][:, None, None] + ddist['pnormW'] = ddist['pnormW'][:, None, None] + if ddist['pnormW'].size != ddist['Ekin_max_eV'].size: + raise Exception() + ddist['Te_eV'] = Te_eV[None, :, None] + if jp_Am2 is not None: + ddist['jp_Am2'] = jp_Am2 + ddist['jp_fraction_re'] = jp_fraction_re[None, None, :] + + nEkin = ddist['Ekin_max_eV'].shape[0] + + # ------------ + # d2cross_phi + # ------------ + + if d2cross_phi is None: + d2cross_phi = _D2CROSS_PHI + + if not d2cross_phi.endswith('.npz'): + d2cross_phi = _DPFE_D2CROSS_PHI[d2cross_phi] + + # -------------- + # anisotropic + # -------------- + + # ------------------------ + # integrated cross-section + + demiss, ddist, d2cross_phi = tfphysemis.get_xray_thin_integ_dist( + # ---------------- + # cross-section + # tabulated d2cross_phi + d2cross_phi=d2cross_phi, + # d2cross_phi computation + E_ph_eV=None, + E_e0_eV=None, + E_e0_eV_npts=None, + theta_e0_vsB_npts=None, + phi_e0_vsB_npts=None, + theta_ph_vsB=None, + # integration + integration=integration, + # ----------- + # verb + debug=False, + verb=True, + # ---------------- + # electron distribution + **ddist, + ) + + # add check units + units0 = dplasma['emiss_tot']['ff']['units'] + units1 = demiss['emiss']['maxwell']['emiss']['units'] + assert asunits.Unit(units0) == asunits.Unit(units1) + units = units0 + + # -------------- + # shapes + # -------------- + + shape_Te = dplasma['common']['Te']['data'].shape + shape_plasma = ddist['plasma']['Te_eV']['data'].shape + assert shape_Te[0] == shape_plasma[1] + + shape_emiss_anis = demiss['emiss']['maxwell']['emiss']['data'].shape + nEph1 = demiss['E_ph_eV']['data'].size + ntheta = demiss['theta_ph_vsB']['data'].size + shape_check = (nEkin, shape_Te[0], jp_fraction_re.size, nEph1, ntheta) + assert shape_emiss_anis == shape_check + + # ---------------------- + # add Zeff to emiss_anis + # ---------------------- + + # reshape Zeff + Zeff = dplasma['common']['Zeff']['data'][None, :, None, None, None] + + # add + for k0, v0 in demiss['emiss'].items(): + demiss['emiss'][k0]['emiss']['data'] = v0['emiss']['data'] * Zeff + + # -------------- + # uniformize + # -------------- + + E_ph0 = dplasma['common']['E_photon']['data'] + E_ph1 = demiss['E_ph_eV']['data'] + ilow = E_ph1 < E_ph0.min() + iup = E_ph1 > E_ph0.max() + E_ph = np.r_[E_ph1[ilow], E_ph0, E_ph1[iup]] + + shape_new_anis = shape_emiss_anis[:-2] + (E_ph.size, ntheta) + shape_new_iso = (1, shape_Te[0], 1, E_ph.size, 1) + + iin = np.r_[ + np.zeros(ilow.sum(), dtype=bool), + np.ones(E_ph0.size, dtype=bool), + np.zeros(iup.sum(), dtype=bool), + ] + iout = ~iin + + # --------------------- + # isotropic - 0-padding + + diso = {} + sli_in = (slice(None),) * 3 + (iin, slice(None)) + sli_broad = (None, slice(None), None, slice(None), None) + for ff in ['ff', 'fb', 'bb']: + diso[ff] = np.zeros(shape_new_iso, dtype=float) + diso[ff][sli_in] = dplasma['emiss_tot'][ff]['data'][sli_broad] + + # --------------------- + # anisotropic: interpolate + + danis = {} + sli_out = (slice(None),) * 3 + (~iin, slice(None)) + iedges = (ilow | iup) + sli_edges = (slice(None),) * 3 + (iedges, slice(None)) + for k0, v0 in demiss['emiss'].items(): + + # tabulated + danis[k0] = np.zeros(shape_new_anis, dtype=float) + danis[k0][sli_out] = v0['emiss']['data'][sli_edges] + + # interpolated + iok = v0['emiss']['data'] > 0. + for ind in np.ndindex(v0['emiss']['data'].shape[:-2]): + + sli = ind + (slice(None), slice(None)) + iok = np.all(v0['emiss']['data'][sli] > 0, axis=-1) + if not np.any(iok): + continue + sli = ind + (iok, slice(None)) + + sli_inii = ind + (iin, slice(None)) + + danis[k0][sli_inii] = np.power( + 10, + scpinterp.make_interp_spline( + np.log10(E_ph1[iok]), + np.log10(v0['emiss']['data'][sli]), + k=1, + axis=-2, + check_finite=True, + )(np.log10(E_ph0)), + ) + + # -------------- + # output + # -------------- + + demiss = { + 'emiss': { + 'maxwell': { + 'ff': { + 'data': danis['maxwell'], + 'units': units, + }, + 'ff_iso': { + 'data': diso['ff'], + 'units': units, + }, + 'fb': { + 'data': diso['fb'], + 'units': units, + }, + 'bb': { + 'data': diso['bb'], + 'units': units, + }, + }, + 'RE': { + 'ff': { + 'data': danis['RE'], + 'units': units, + }, + }, + }, + 'E_ph': { + 'data': E_ph, + 'units': 'eV', + }, + 'theta_ph_vsB': demiss['theta_ph_vsB'], + 'Te': ddist['plasma']['Te_eV'], + 'ne': ddist['plasma']['ne_m3'], + 'jp_fraction_re': ddist['plasma']['jp_fraction_re'], + 'Ekin_max_eV': ddist['plasma']['Ekin_max_eV'], + } + + return demiss, ddist, dmix, data_source diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_main.py b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_main.py new file mode 100644 index 0000000..881c7c6 --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_main.py @@ -0,0 +1,50 @@ +import os +import sys + + +_PATH_HERE = os.path.dirname(__file__) +_PATH_PUBLI = os.path.dirname(os.path.dirname(os.path.dirname(_PATH_HERE))) +_PATH_PROJECTS = os.path.dirname(_PATH_PUBLI) +_PATH_TF = os.path.join(_PATH_PROJECTS, 'tofu') +sys.path.insert(0, _PATH_TF) +import tofu as tf +sys.path.pop(0) + + +# ##################################################### +# ##################################################### +# Local imports +# ##################################################### + + +from ._fig01_cross import main as fig01 +from ._fig02_dist import main as fig02 +from ._fig03_validate_ff import main as fig03 +from ._fig04_bremsstrahlung import main as fig04 +from ._fig05_emiss import main as fig05 +from ._fig06_tokamak import main as fig06 +from ._fig07_responsivities import main as fig07 + + +# ##################################################### +# ##################################################### +# main +# ##################################################### + + +def main( + ne_m3=None, + pfe_save=None, +): + + # --------------- + # plot all + # --------------- + + for fig in [fig01, fig02, fig04, fig05, fig06]: + fig( + ne_m3=ne_m3, + pfe_save=pfe_save, + ) + + return diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_perfs.py b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_perfs.py new file mode 100644 index 0000000..c674926 --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_perfs.py @@ -0,0 +1,412 @@ +import copy + + +import numpy as np +import scipy.interpolate as scpinterp +import scipy.integrate as scpinteg +import scipy.constants as scpct +import astropy.units as asunits +import datastock as ds + + +from . import _load_spect_anis +from ._fig02_dist_type import _DDIST +from ._fig05_emiss import _RE, _DDMIX +from ._fig07_responsivities import _PFE_RESPONSIVITIES, _LRESP + + +# ####################################### +# ####################################### +# DEFAULTS +# ####################################### + + +_DANGLES = { + 'in-vessel': { + 'head-on': np.r_[0, 40]*np.pi/180, + 'back': np.r_[140, 180]*np.pi/180, + 'resp': ['cvd_bare', 'cvd_filter', 'bolo'], + 'color': 'tab:orange', + 'alpha': 0.4, + }, + 'ex-cryostat': { + 'head-on': np.r_[60, 70]*np.pi/180, + 'back': np.r_[110, 120]*np.pi/180, + 'resp': [ + 'mesxr_11_keV', 'mesxr_18_keV', + 'mehxr_20_keV', 'mehxr_60_keV', + 'spectro', + ], + 'color': 'tab:green', + 'alpha': 0.4, + }, +} + + +# ####################################### +# ####################################### +# Main +# ####################################### + + +def main( + dmix=None, + # cases + Te_eV=None, + # d2cross + d2cross_phi=None, + # dist + ne_m3=None, + jp_Am2=None, + # RE + re=None, + dominant=None, + jp_fraction_re=None, + Efield_par_Vm=None, + Ekin_min_eV=None, + Ekin_max_eV=None, + sigmap=None, + pnormW=0, + # jp_fraction_re=np.linspace(0.025, 0.975, 39), + # unused + **kwdargs, +): + + # ----------- + # inputs + # ----------- + + if dmix is None: + dmix = _DDMIX[0] + + # Maxwell + kwd_max = {'ne_m3': ne_m3, 'jp_Am2': jp_Am2} + + # RE + re = ds._generic_check._check_var( + re, 're', + types=str, + default=_RE, + allowed=sorted(_DDIST['RE'].keys()), + ) + + lRE = [ + 'dominant', 'jp_fraction_re', 'Efield_par_Vm', + 'Efield_par_Vm', 'Ekin_max_eV', 'Ekin_min_eV', + 'sigmap', 'pnormW' + ] + kwd_RE = { + kk: _DDIST['RE'][re].get(kk) if vv is None else vv + for kk, vv in locals().items() + if kk in lRE + } + + # -------------- + # load elements + # -------------- + + demiss = {} + kwd = dict(kwd_max) + kwd.update(**kwd_RE) + demiss, ddist, dmix, data_source = _load_spect_anis.main( + dmix=dmix, + # d2cross + d2cross_phi=d2cross_phi, + # dist + **kwd, + ) + + # extract + # E_ph = demiss[0]['common']['E_photon']['data'] + Teu = np.unique(ddist['plasma']['Te_eV']['data']) + # ne = np.unique(ddist['plasma']['ne_m3']['data'])[0] + # jp = np.unique(ddist['plasma']['jp_Am2']['data'])[0] + # units = demiss['emiss']['maxwell']['ff']['units'] + + if Te_eV is None: + sli_emiss = (0, slice(None), slice(None), slice(None), slice(None)) + axis_emiss = (0, 1, -1) + sli_resp = (None, None, slice(None), None) + else: + indTe = np.argmin(np.abs(Teu - Te_eV)) + Te_eV = Teu[indTe] + sli_emiss = (0, indTe, 0, slice(None), slice(None)) + axis_emiss = (-1,) + sli_resp = (slice(None), None) + + # ------------------- + # load responsivities + # ------------------- + + # load + dresp = { + k0: v0.tolist() + for k0, v0 in np.load(_PFE_RESPONSIVITIES, allow_pickle=True).items() + if k0 in _LRESP + } + + # ----------------------------- + # Integrate over responsivity + # ----------------------------- + + demiss_integ = {kk: {'maxwell': {}, 'RE': {}} for kk in dresp.keys()} + for kresp, vresp in dresp.items(): + + # Eph_resp + iok_resp = np.isfinite(vresp['responsivity']['data']) + iok_resp[iok_resp] = vresp['responsivity']['data'][iok_resp] > 0. + Eph_resp = vresp['E_eV']['data'][iok_resp] + + # loop on dist + for kdist, vdist in demiss['emiss'].items(): + + # loop on emiss type + for kemiss, vemiss in vdist.items(): + + # skip + if kemiss == 'ff_iso': + continue + + # Eph_emiss + data = demiss['emiss'][kdist][kemiss]['data'][sli_emiss] + iok_emiss = ( + (demiss['E_ph']['data'] > Eph_resp.min()) + & (demiss['E_ph']['data'] < Eph_resp.max()) + & np.all(np.isfinite(data), axis=axis_emiss) + & np.all(data > 0., axis=axis_emiss) + ) + Eph_emiss = demiss['E_ph']['data'][iok_emiss] + + # Eph + Eph = np.unique(np.r_[Eph_resp, Eph_emiss]) + + # interpolate responsivity + resp = scpinterp.make_interp_spline( + Eph_resp, + vresp['responsivity']['data'][iok_resp], + k=1, + axis=0, + check_finite=True, + )(Eph) + + # units + unitsE = asunits.Unit(demiss['E_ph']['units']) + unitsR = asunits.Unit(vresp['responsivity']['units']) + unitsD = asunits.Unit(demiss['emiss'][kdist][kemiss]['units']) + units = unitsE * unitsR * unitsD + + # ph vs energy + if 'W' in str(unitsR): + data = data * demiss['E_ph']['data'][:, None] * scpct.e + units = units * asunits.Unit('W.s') + else: + units = units * asunits.Unit('ph') + + # interpolate emissivity + if np.any(iok_emiss): + sli = (slice(None),) * (data.ndim-2) + (iok_emiss, slice(None)) + + if kemiss == 'bb': + emiss = scpinterp.make_interp_spline( + Eph_emiss, + data[sli], + k=1, + axis=-2, + bc_type=None, + check_finite=True, + )(Eph) + + else: + emiss = np.power( + 10, + scpinterp.make_interp_spline( + np.log10(Eph_emiss), + np.log10(data[sli]), + k=1, + axis=-2, + bc_type=None, + check_finite=True, + )(np.log10(Eph)), + ) + + # adjust out of bounds + iout = (Eph < Eph_emiss[0]) | (Eph > Eph_emiss[-1]) + sli_out = sli[:-2] + (iout, slice(None)) + emiss[sli_out] = 0. + + # sanity check + inan = (~np.isfinite(emiss)) | (emiss < 0.) + if np.any(inan): + msg = ( + "Interpolation prior to energy integration: " + "non-finite values in:\n" + f" \t- {kresp} {kdist} {kemiss}: {inan.sum()}" + ) + raise Exception(msg) + + # integrate + demiss_integ[kresp][kdist][kemiss] = { + 'data': scpinteg.trapezoid( + resp[sli_resp] * emiss, + x=Eph, + axis=-2, + ), + 'units': units, + } + else: + shape = data.shape[:-2] + data.shape[-1:] + zeros = np.zeros(shape, dtype=float) + demiss_integ[kresp][kdist][kemiss] = { + 'data': zeros, + 'units': units, + } + + # ----------------------------- + # Integrate over angles + # ----------------------------- + + dang = {} + dsignal = copy.deepcopy(demiss_integ) + theta = demiss['theta_ph_vsB']['data'] + for kresp, vresp in dresp.items(): + + # get relevant angle + kang = [ + kk for kk, vv in _DANGLES.items() + if kresp in vv['resp'] + ][0] + dang[kresp] = kang + + # integrate each direction + for idir, kdir in enumerate(['head-on', 'back']): + + # iang + iang = ( + (theta >= _DANGLES[kang][kdir][0]) + & (theta <= _DANGLES[kang][kdir][1]) + ) + + # solid angle assuming cone + dcos = np.diff(np.cos(_DANGLES[kang][kdir])[::-1])[0] + sang = 2*np.pi * dcos + + # loop on dist + for kdist, vdist in demiss_integ[kresp].items(): + + # loop on emiss type + for kemiss, vemiss in vdist.items(): + + ndim = vemiss['data'].ndim + sli_emiss = (slice(None),) * (ndim - 1) + (iang,) + sli_ang = (None,) * (ndim - 1) + (iang,) + + # integrate + if vemiss['data'].shape[-1] == 1: + data = vemiss['data'][..., 0] * sang + else: + data = scpinteg.trapezoid( + vemiss['data'][sli_emiss] + * np.sin(theta[sli_ang]), + x=theta[iang], + axis=-1, + ) * 2*np.pi + + # sanity check + inan = (~np.isfinite(data)) | (data < 0.) + if np.any(inan): + msg = ( + "Angle integration: non-finite values in:\n" + f" \t- {kresp} {kdist} {kemiss} {kdir}: " + f"{inan.sum()}" + ) + raise Exception(msg) + + # store + if idir == 0: + dsignal[kresp][kdist][kemiss] = { + kdir: { + 'data': data, + 'units': vemiss['units'] * asunits.Unit('sr'), + } + } + else: + dsignal[kresp][kdist][kemiss][kdir] = { + 'data': data, + 'units': vemiss['units'] * asunits.Unit('sr'), + } + + # -------------- + # total_headon + # -------------- + + lresp = _LRESP + ldist = sorted(dsignal[lresp[0]].keys()) + + # details + shape = (len(lresp),) + data.shape + shape_bb = dsignal[kresp]['maxwell']['bb']['head-on']['data'].shape + if len(shape_bb) == 0: + zeros = 0. + else: + zeros = np.zeros(shape_bb, dtype=float) + total_headon = np.zeros(shape, dtype=float) + total_back = np.zeros(shape, dtype=float) + for kdist in ldist: + for kemiss in sorted(dsignal[lresp[0]][kdist].keys()): + total_headon[...] += np.array([ + zeros if (kresp == 'spectro' and kemiss == 'bb') + else dsignal[kresp][kdist][kemiss]['head-on']['data'] + for kresp in lresp + ]) + total_back[...] += np.array([ + zeros if (kresp == 'spectro' and kemiss == 'bb') + else dsignal[kresp][kdist][kemiss]['back']['data'] + for kresp in lresp + ]) + + # diff_RE + diff_RE = ( + np.array([ + dsignal[kresp]['RE']['ff']['head-on']['data'] + for kresp in lresp + ]) + - np.array([ + dsignal[kresp]['RE']['ff']['back']['data'] + for kresp in lresp + ]) + ) + + # diff_max + diff_max = ( + np.array([ + dsignal[kresp]['maxwell']['ff']['head-on']['data'] + for kresp in lresp + ]) + - np.array([ + dsignal[kresp]['maxwell']['ff']['back']['data'] + for kresp in lresp + ]) + ) + + # sanity check + error = (diff_RE + diff_max) - (total_headon - total_back) + error_percent = np.zeros(diff_RE.shape, dtype=float) + iok = (diff_RE > 0.) & np.isfinite(total_headon) + error_percent[iok] = 100 * error[iok] / diff_RE[iok] + if np.any(error_percent > 0.01): + lstr = [ + f"\t- {ss}: {error_percent[ii]:2.1e} %" + for ii, ss in enumerate(lresp) + ] + msg = ( + "Something wrong with fb or bb:\n" + + "\n".join(lstr) + ) + raise Exception(msg) + + return ( + demiss_integ, dsignal, ddist, dmix, + total_headon, diff_RE, diff_max, + dang, theta, + lresp, ldist, + ) diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_savefig.py b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_savefig.py new file mode 100644 index 0000000..75cef75 --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_savefig.py @@ -0,0 +1,57 @@ + + +import os +import datastock as ds + + +# ####################################### +# ####################################### +# Default +# ####################################### + + +# PATHS +_PATH_HERE = os.path.dirname(__file__) +_PATH_PAPER = os.path.dirname(_PATH_HERE) +_PATH_SAVE = os.path.join(_PATH_PAPER, 'figures') + + +# ####################################### +# ####################################### +# Main +# ####################################### + + +def main( + fig=None, + pfe_save=None, + path_save=None, + file=None, +): + + # ------------- + # inputs + # ------------- + + # pfe_save + pfe_save = ds._generic_check._check_var( + pfe_save, 'pfe_save', + types=(bool, str, type(None)), + default=False, + ) + + # ------------- + # saving + # ------------- + + if pfe_save is not False: + if pfe_save in [None, True]: + name = f"{os.path.split(file)[-1][1:].replace('.py', '')}.png" + if path_save is None: + path_save = _PATH_SAVE + pfe_save = os.path.join(_PATH_SAVE, name) + fig.savefig(pfe_save, format='png', dpi=300) + msg = f"Saved figure in:\n\t{pfe_save}\n" + print(msg) + + return diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_singlept_sensor.py b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_singlept_sensor.py new file mode 100644 index 0000000..5f16305 --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_singlept_sensor.py @@ -0,0 +1,495 @@ +import copy + + +import numpy as np +import scipy.integrate as scpinteg +import astropy.units as asunits +import tofu as tf + + +from ._fig11_hxrva_cad import _PFE_CONFIG, _PFE_COLL, _add_ptcam +from . import _perfs + + +# ####################################### +# ####################################### +# DEFAULTS +# ####################################### + + +_RE = ['avalanche 10 MeV', 'dreicer'] +_JP_FRAC = 0.9 +_TE = 0.5e3 + + +_FONTSIZE = 14 + + +# ####################################### +# ####################################### +# Main +# ####################################### + + +def main( + # coll + coll=None, + key_cam=None, + config=None, + res=None, + # responsivity + key_resp=None, + # equilibrium + t=None, + # dmix + dmix=None, + # d2cross + d2cross_phi=None, + # dist + ne_m3=None, + jp_Am2=None, + Te_eV=None, + jp_fraction_re=None, + # RE + re=None, + # assumptions + cases=None, + # ptcam + angle0=None, + angle1=None, + # unused + **kwdargs, +): + + # ----------- + # inputs + # ----------- + + if Te_eV is None: + Te_eV = _TE + + if jp_fraction_re is None: + jp_fraction_re = _JP_FRAC + + if re is None: + re = _RE + if isinstance(re, str): + re = [re] + + if key_resp is None: + key_resp = 'cvd_filter' + + # -------------- + # compute + # -------------- + + demiss_integ = {} + for rei in re: + ( + demiss_integ[rei], dsignal, ddist, dmix, + total_headon, diff_RE, diff_max, + dang, theta, + lresp, ldist, + ) = _perfs.main( + dmix=dmix, + ne_m3=ne_m3, + jp_Am2=jp_Am2, + d2cross_phi=d2cross_phi, + re=rei, + jp_fraction_re=jp_fraction_re, + ) + + # -------------- + # extract + # -------------- + + Teu = np.unique(ddist['plasma']['Te_eV']['data']) + neu = np.unique(ddist['plasma']['ne_m3']['data']) + jpu = np.unique(ddist['plasma']['jp_Am2']['data']) + jp_fracu = np.unique(ddist['plasma']['jp_fraction_re']['data']) + assert neu.size == jpu.size == jp_fracu.size == 1 + + indTe = np.argmin(np.abs(Teu - Te_eV)) + sli_Te = (indTe, 0, slice(None)) + Te_eV = Teu[indTe] + ne = neu[0] + jp = jpu[0] + jp_frac = jp_fracu[0] + + # -------------- + # load coll + # -------------- + + # config + if config is None: + config = _PFE_CONFIG + if isinstance(config, str): + config = tf.load(config) + + # -------------- + # load cameras + # -------------- + + # coll + if coll is None: + coll = _PFE_COLL + if isinstance(coll, str): + coll = tf.data.load(coll) + + # -------------- + # add pt_cam rays + # -------------- + + dkrays, angle0, angle1 = _add_ptcam( + coll=coll, + key_cam=key_cam, + angle0=angle0, + angle1=angle1, + config=config, + ) + assert len(dkrays) == 2 + + # -------------- + # get angle vs B + # -------------- + + msg = "\nCompute angles vs B..." + print(msg) + + dangles, axis_samp = _get_angles_vs_B( + coll=coll, + cases=cases, + dkrays=dkrays, + res=res, + t=t, + ) + + # ################ + # interpolate ne + # ################ + + dne_coef = {} + for kcam, krays in dkrays.items(): + + R = dangles[krays]['helicity']['R']['data'] + Z = dangles[krays]['helicity']['Z']['data'] + kmR = [kk for kk in coll.ddata.keys() if kk.endswith('_magaxR')][0] + kmZ = [kk for kk in coll.ddata.keys() if kk.endswith('_magaxZ')][0] + magR = coll.ddata[kmR]['data'] + magZ = coll.ddata[kmZ]['data'] + + sli = (slice(None),) + (None,)*R.ndim + ne_coef = np.exp( + - (R[None, ...] - magR[sli])**2/0.5**2 + - (Z[None, ...] - magZ[sli])**2/0.7**2 + ) + + if t is not None: + kt = f"{kmR.split('_')[0]}_t" + indt = np.argmin(np.abs(coll.ddata[kt]['data'] - t)) + ne_coef = ne_coef[indt, ...] + dne_coef[krays] = ne_coef + + # ################ + # interpolate emissivity + # ################ + + msg = "Interpolate emissivity..." + print(msg) + + dsig_los = { + kcase: { + krays: { + rei: {kdist: {} for kdist in ['maxwell', 'RE']} + for rei in re + } + for krays in dkrays.values() + } + for kcase in cases.keys() + } + for kcase, vcase in cases.items(): + for kcam, krays in dkrays.items(): + + # ------------------ + # find closest angle + + khel = 'helicity' if vcase['helicity'] is True else 'nohelicity' + angles = dangles[krays][khel]['angle']['data'] + R = dangles[krays][khel]['angle']['data'] + Z = dangles[krays][khel]['angle']['data'] + iok = np.isfinite(angles) & np.isfinite(R) + ref_ang = dangles[krays][khel]['angle']['ref'] + ref = tuple([ + rr for ii, rr in enumerate(ref_ang) if ii != axis_samp + ]) + + # ------------------ + # find closest angle + + sli = (None,) * angles.ndim + (slice(None),) + delta = angles[..., None] - theta[sli] + iang = np.argmin(np.abs(delta), axis=-1) + shape = tuple([ + ss for ii, ss in enumerate(iang.shape) + if ii != axis_samp + ]) + + length = dangles[krays][khel]['length']['data'] + length[~iok] = np.nan + ndimd = iang.ndim - length.ndim + sli_ndimd = (slice(None),)*ndimd + + # ----------- + # units, ref + + for rei in re: + for kdist, vdist in demiss_integ[rei][key_resp].items(): + for kemiss, vemiss in vdist.items(): + emiss = vemiss['data'][sli_Te] + + # ne profile + if vcase['ne'] is True: + lengthi = length * dne_coef[krays] + else: + lengthi = length + + # units + initialize + units_emiss = asunits.Unit(vemiss['units']) + emiss_los = np.zeros(shape, dtype=float) + if kemiss in ['fb', 'bb']: + emiss_los[...] = ( + emiss * np.nansum(lengthi, axis=0) + )[None, ...] + else: + for ind in np.ndindex(lengthi.shape[1:]): + ioki = iok[(slice(None),) + ind] + if not np.any(ioki): + continue + ll = lengthi[(ioki,) + ind] + sli_emiss = sli_ndimd + (ioki,) + ind + sli_los = sli_ndimd + ind + emiss_los[sli_los] = scpinteg.trapezoid( + emiss[iang[sli_emiss]], + x=ll, + axis=axis_samp, + ) + + dsig_los[kcase][krays][rei][kdist][kemiss] = { + 'data': emiss_los, + 'units': units_emiss * asunits.Unit('m'), + 'ref': ref, + } + + # store + coll.add_data( + key=f"{kcase}_{krays}_{rei}_{kdist}_{kemiss}", + **dsig_los[kcase][krays][rei][kdist][kemiss], + ) + + # ################ + # merits + # ################ + + kcase = [ + kcase for kcase, vcase in cases.items() + if all([vv for vv in vcase.values()]) + ][0] + + lkrays = sorted(dsig_los[kcase].keys()) + krays_mean = [ + np.nanmean(dsig_los[kcase][kk][re[0]]['RE']['ff']['data']) + for kk in lkrays + ] + krays_max = lkrays[np.argmax(krays_mean)] + krays_min = lkrays[np.argmin(krays_mean)] + + dmetrics = {} + for rei in re: + + # RE + RE_headon = dsig_los[kcase][krays_max][rei]['RE']['ff']['data'] + RE_back = dsig_los[kcase][krays_min][rei]['RE']['ff']['data'] + diffRE = RE_headon - RE_back + + # Max + kd = 'maxwell' + Max_headon_ff = dsig_los[kcase][krays_max][rei][kd]['ff']['data'] + Max_back_ff = dsig_los[kcase][krays_min][rei][kd]['ff']['data'] + diffMax = Max_headon_ff - Max_back_ff + + Max_headon_tot = ( + dsig_los[kcase][krays_max][rei][kd]['ff']['data'] + + dsig_los[kcase][krays_max][rei][kd]['fb']['data'] + + dsig_los[kcase][krays_max][rei][kd]['bb']['data'] + ) + total_headon = Max_headon_tot + RE_headon + + # metrics + units = dsig_los[kcase][krays_max][rei]['RE']['ff']['units'] + ref = dsig_los[kcase][krays_max][rei]['RE']['ff']['ref'] + dmetrics[rei] = { + 'xi': { + 'data': diffRE / total_headon, + 'units': '', + 'ref': ref, + }, + 'kappa': { + 'data': diffRE / (diffRE + diffMax), + 'units': '', + 'ref': ref, + }, + 'meas_RE_diff': { + 'data': diffRE, + 'units': units, + 'ref': ref, + }, + 'meas_RE_headon': { + 'data': RE_headon, + 'units': units, + 'ref': ref, + }, + 'meas_headon': { + 'data': total_headon, + 'units': units, + 'ref': ref, + }, + } + + return ( + coll, config, + dangles, dsig_los, dmetrics, + krays_max, krays_min, + angle0, angle1, + ne, jp, jp_frac, Te_eV, + ) + + +# ############################################ +# ############################################ +# angles vs B +# ############################################ + + +def _get_angles_vs_B( + coll=None, + cases=None, + dkrays=None, + res=None, + t=None, +): + + # ---------------- + # keys of interest + # ---------------- + + # B-field keys + kBR, kBZ, kBphi = [ + [kk for kk in coll.ddata.keys() if kk.endswith(f'_B{k0}')][0] + for k0 in ['R', 'Z', 'phi'] + ] + + # Sep + ksepR, ksepZ = [ + [kk for kk in coll.ddata.keys() if kk.endswith(f'_sep{k0}')][0] + for k0 in ["R", "Z"] + ] + + # time + kt = f"{kBR.split('_')[0]}_t" + + # ---------------- + # loop on cases / rays + # ---------------- + + dangles = { + krays: {'helicity': {}, 'nohelicity': {}} + for krays in dkrays.values() + } + + for kcam, krays in dkrays.items(): + + # --------------------------------------- + # compute angles along rays - helicity + + dout = coll.get_rays_angle_vs_vect( + # rays + key_rays=krays, + res=res, + segment=-1, + # vector components + key_XR=kBR, + key_YZ=kBZ, + key_Zphi=kBphi, + geometry='toroidal', + # separatrix + key_sepR=True, + key_sepZ=True, + # verb + verb=None, + )[krays] + + # select time + if t is not None: + indt = np.argmin(np.abs(coll.ddata[kt]['data'] - t)) + t = coll.ddata[kt]['data'][indt] + + dout['angle']['data'] = dout['angle']['data'][indt, ...] + dout['angle']['ref'] = dout['angle']['ref'][1:] + + dangles[krays]['helicity'] = dout + + # --------------------------------------- + # compute angles along rays - helicity + + phi = dangles[krays]['helicity']['phi']['data'] + ephix = -np.sin(phi) + ephiy = np.cos(phi) + if np.nanmean(coll.ddata[kBphi]['data']) < 0: + ephix, ephix = -ephix, -ephix + vx, vy, vz = coll.get_rays_vect(krays) + cos = ephix * (-vx) + ephiy * (-vy) + ang = np.arccos(cos) + if ang.ndim < dout['angle']['data'].ndim: + nt = dout['angle']['data'].shape[0] + ang = np.array([ang for ii in range(nt)]) + ang[~np.isfinite(dout['angle']['data'])] = np.nan + + dangles[krays]['nohelicity'] = copy.deepcopy(dout) + dangles[krays]['nohelicity']['angle']['data'] = ang + + # ------------ + # add missing ref + + axis_samp = dout['angle']['ref'].index(None) + + kref = f'{krays}_nsamp0' + coll.add_ref( + key=kref, + size=ang.shape[axis_samp], + ) + + # ------------ + # store + + ref = dout['angle']['ref'] + ref[ref.index(None)] = kref + + key = f"{krays}_helicity" + coll.add_data( + key=key, + data=dangles[krays]['helicity']['angle']['data'], + ref=tuple(ref), + units=dangles[krays]['helicity']['angle']['units'], + dim=dangles[krays]['helicity']['angle']['dim'], + ) + + key = f"{krays}_nohelicity" + coll.add_data( + key=key, + data=dangles[krays]['nohelicity']['angle']['data'], + ref=tuple(ref), + units=dangles[krays]['nohelicity']['angle']['units'], + dim=dangles[krays]['nohelicity']['angle']['dim'], + ) + + return dangles, axis_samp diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/figures/HXR_MPP200_zoom.png b/articles/RSI_2026_RunawayBremsstrahlungDetection/figures/HXR_MPP200_zoom.png new file mode 100644 index 0000000..cda7a7f --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/figures/HXR_MPP200_zoom.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:d740b14d239753cb278fa11f7318f99b1324c5e88cc929c3111768b7a43d7016 +size 1157098 diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/figures/fig01_cross.png b/articles/RSI_2026_RunawayBremsstrahlungDetection/figures/fig01_cross.png new file mode 100644 index 0000000..7f84726 --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/figures/fig01_cross.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:5dcc4ef9989f783fd6862992637db5366669530b105ffab6ca87ab8e68189159 +size 511464 diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/figures/fig02_RE_ff.png b/articles/RSI_2026_RunawayBremsstrahlungDetection/figures/fig02_RE_ff.png new file mode 100644 index 0000000..599f794 --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/figures/fig02_RE_ff.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:fc1a6ca4c7b9bdfeca44663eb6d09fb33b0946f49100746bfeb26a28cda33e87 +size 350298 diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/figures/fig02_dist.png b/articles/RSI_2026_RunawayBremsstrahlungDetection/figures/fig02_dist.png new file mode 100644 index 0000000..a5672d9 --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/figures/fig02_dist.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:01905bf070436cae533d5054e729555761361ede9bfa226d9549fbf2813171e4 +size 341741 diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/figures/fig02_dist_type.png b/articles/RSI_2026_RunawayBremsstrahlungDetection/figures/fig02_dist_type.png new file mode 100644 index 0000000..e8ba40c --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/figures/fig02_dist_type.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:6cdd302af9a66c56e0bb282e88f0fe96cab14461a1dc169c275d2e977ceae39d +size 309704 diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/figures/fig03_validate_ff.png b/articles/RSI_2026_RunawayBremsstrahlungDetection/figures/fig03_validate_ff.png new file mode 100644 index 0000000..0cb90dd --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/figures/fig03_validate_ff.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:1f12498000348ed6c8fd270024c7129031f9073366e8dce093ba0fb9e22af2a1 +size 324034 diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/figures/fig04_bremsstrahlung.png b/articles/RSI_2026_RunawayBremsstrahlungDetection/figures/fig04_bremsstrahlung.png new file mode 100644 index 0000000..3752584 --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/figures/fig04_bremsstrahlung.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:1da777136108a83ae6a2ad736cad5bcfeec8e7756d791e81def8e87c21f078c6 +size 494583 diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/figures/fig05_tokamak.png b/articles/RSI_2026_RunawayBremsstrahlungDetection/figures/fig05_tokamak.png new file mode 100644 index 0000000..064382f --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/figures/fig05_tokamak.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:cef441d701316cb62d3be568418d8bb64b36ecec53b9748343e918d0fd2184f2 +size 267360 diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/figures/fig06_responsivities.png b/articles/RSI_2026_RunawayBremsstrahlungDetection/figures/fig06_responsivities.png new file mode 100644 index 0000000..963ad42 --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/figures/fig06_responsivities.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:8dbc75675a2ae86e47273887c18a53ea2f49fe5dd39804678c743b4aa86b7de8 +size 144863 diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/figures/fig07_responsivities.png b/articles/RSI_2026_RunawayBremsstrahlungDetection/figures/fig07_responsivities.png new file mode 100644 index 0000000..f2c397a --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/figures/fig07_responsivities.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:c33828cabe2c74d40283e9a30b8afd3dc4d1c3d4e99580d28cc390b148bb5028 +size 136528 diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/figures/fig08_perfs_single.png b/articles/RSI_2026_RunawayBremsstrahlungDetection/figures/fig08_perfs_single.png new file mode 100644 index 0000000..83840a6 --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/figures/fig08_perfs_single.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:2dccc5a7a810e2576ce7b1d9a0145e125a759393833e2bf60d3123a995858dc3 +size 562497 diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/figures/fig11_hxrva_cad.png b/articles/RSI_2026_RunawayBremsstrahlungDetection/figures/fig11_hxrva_cad.png new file mode 100644 index 0000000..750ce87 --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/figures/fig11_hxrva_cad.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:277c885624fdf24611a70abf5846e6e9db35d7929bd739d97d6d00f9e8bbe5c7 +size 450244 diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/figures/fig14_bolo_cad.png b/articles/RSI_2026_RunawayBremsstrahlungDetection/figures/fig14_bolo_cad.png new file mode 100644 index 0000000..37a62c7 --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/figures/fig14_bolo_cad.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:b31ff6962af8ccca83f0a4a9cb2a23ce45779e1e29c245e9595e88cb7f948c35 +size 467193 diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/260616_SCRAM86_Ar_data.npz b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/260616_SCRAM86_Ar_data.npz new file mode 100644 index 0000000..f7b0b27 --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/260616_SCRAM86_Ar_data.npz @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:8ade9d576ca6b6275a37173637a5c3f46f9892b06abe28359c4a11d9b9086ccb +size 16174102 diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/260616_SCRAM86_Fe_data.npz b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/260616_SCRAM86_Fe_data.npz new file mode 100644 index 0000000..ab17ebe --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/260616_SCRAM86_Fe_data.npz @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:fe1399a531f279c4266e91e03f4b953b4a896105b10bf4f2009119faa61d60c9 +size 16178902 diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/260616_SCRAM86_Kr_data.npz b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/260616_SCRAM86_Kr_data.npz new file mode 100644 index 0000000..ce64172 --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/260616_SCRAM86_Kr_data.npz @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:19e41ec071a61319a3093f42d69d43abf6aa4237b45e8c0c0d489f3442b8d85e +size 80824912 diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/260616_SCRAM86_O_data.npz b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/260616_SCRAM86_O_data.npz new file mode 100644 index 0000000..e6a2b9b --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/260616_SCRAM86_O_data.npz @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:fbdae5a5337669df2cb1be49d6094b4040cb09d2608525e0ce36f518ea0f719c +size 16168102 diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/260616_SCRAM86_W_data.npz b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/260616_SCRAM86_W_data.npz new file mode 100644 index 0000000..c1d2360 --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/260616_SCRAM86_W_data.npz @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:16a8d13a73f0ea5bd836e6496fe27213eb37987aa307dcd407dc176f82a11795 +size 1506241 diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/260616_plt_scram.py b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/260616_plt_scram.py new file mode 100644 index 0000000..4108b93 --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/260616_plt_scram.py @@ -0,0 +1,364 @@ +''' + +Look at SCRAM runs for Ca + +cjperks +June 16, 2026 + +''' + + +# Modules +import numpy as np +import sys, os +from matplotlib.cm import get_cmap +from matplotlib.colors import Normalize +import scipy.constants as cnt +from scipy.interpolate import interp1d +import ast + +from atomic_world.run_SCRAM.main import read_v86 as r86 + +colors = plt.rcParams['axes.prop_cycle'].by_key()['color'] + +# Species +elem = 'O_' +Znuc = 8 + +# Plasma +ne_cm3 = 1e14 +Te_eV = 1e3 + +Te_long = np.logspace(np.log10(1000), np.log10(10000), 51) # [eV] +ne_long = 1e14*np.ones_like(Te_long) + +charges = np.arange(0,Znuc+1)[None,:] + +################################################### +# +# SCRAM data +# +################################################### + +# Common file path +fol = os.path.join( + '/home/cjperks', + 'sharing/for_dVEZINET', + '260616_SXR_emis' + ) +# Data file +fil = '260616_SCRAM86_%s_data.npz'%(elem.split('_')[0]) + +# Loads pickle object +dscram = np.load( + os.path.join(fol, fil), + allow_pickle=True + )['arr_0'][()] + + +######################################### +# +# Plotting +# +######################################### + +Znuc = r86._get_nele(sym=elem.split('_')[0]) + +### --- Plots ion/rec rates for a charge state --- ### + +fig1, ax1 = plt.subplots( + 1,2, + figsize = (10,6) + ) +nele = 2 + +# Rates +ax1[0].plot( + dscram['Te']['data']/1e3, + dscram['scd']['data'][:,Znuc-nele], + 'b-', + #label = 'ci+ea; SCRAM' + ) +ax1[0].plot( + dscram['Te']['data']/1e3, + dscram['acd']['data'][:,Znuc-nele], + 'b--', + #label = 'rr+dr; SCRAM' + ) + +ax1[0].plot( + np.nan, np.nan, + 'k-', + label = 'ioniz.' + ) +ax1[0].plot( + np.nan, np.nan, + 'k--', + label = 'recomb.' + ) + +ax1[0].grid('on') +ax1[0].set_xscale('log') +ax1[0].set_xlabel(r'$T_e$ [keV]') +ax1[0].set_yscale('log') +ax1[0].set_ylabel(r'rate [$cm^3/s/electron$]') +ax1[0].set_xlim(0.1, 2.5) +ax1[0].set_ylim(1e-14,1e-9) + +leg = ax1[0].legend( + loc = 'lower right', + labelcolor='linecolor', + borderaxespad=0, + ) + +# Abundance +ax1[1].plot( + dscram['Te']['data']/1e3, + dscram['Xz']['data'][:,Znuc-nele], + 'b-', + label = 'SCRAM' + ) + +ax1[1].grid('on') +ax1[1].set_xscale('log') +ax1[1].set_xlabel(r'$T_e$ [keV]') +ax1[1].set_yscale('log') +ax1[1].set_ylabel(r'abundance [frac]') +ax1[1].set_xlim(0.1, 2.5) +ax1[1].set_ylim(1e-3, 1e0) + +leg = ax1[1].legend( + loc = 'lower left', + labelcolor='linecolor', + borderaxespad=0, + ) + + +fig1.suptitle( + r'%s; nele=%i'%( + elem, nele + ) + ) + +plt.tight_layout() + + + +### --- Plots and cooling curve --- ### + +fig2, ax2 = plt.subplots( + 1,2, + #figsize = (12,6) + ) + +# +ax2[0].plot( + dscram['Te']['data']/1e3, + dscram['']['data'], + 'b-', + label = 'SCRAM' + ) + + +ax2[0].plot( + [0.1, 30], + [Znuc-2, Znuc-2], + 'k--' + ) +ax2[0].plot( + [0.1, 30], + [Znuc-10, Znuc-10], + 'k--' + ) +if Znuc > 28: + ax2[0].plot( + [0.1, 30], + [Znuc-28, Znuc-28], + 'k--' + ) + +ax2[0].grid('on') +ax2[0].set_xscale('log') +ax2[0].set_xlabel(r'$T_e$ [keV]') +ax2[0].set_ylabel('') +ax2[0].set_xlim(0.1, 2.5) +ax2[0].set_ylim(max(Znuc-10,0),Znuc) + +leg = ax2[0].legend( + loc = 'lower right', + labelcolor='linecolor', + borderaxespad=0, + ) + + + +# Prad +ax2[1].plot( + dscram['Te']['data']/1e3, + dscram['Prad']['data'], + 'b-', + label = 'SCRAM' + ) + + +ax2[1].grid('on') +ax2[1].set_xscale('log') +ax2[1].set_xlabel(r'$T_e$ [keV]') +ax2[1].set_yscale('log') +ax2[1].set_ylabel(r'$P_{rad}$ [$W*m^3/atom/electron$]') +ax2[1].set_xlim(0.1, 2.5) +ax2[1].set_ylim(1e-35,1e-32) + + +fig2.suptitle( + r'%s'%( + elem + ) + ) + +plt.tight_layout() + + + +### --- Plots the spectrum at a single Te --- ### + + +fig3, ax3 = plt.subplots() + +ind = np.argmin(abs( + Te_eV - dscram['Te']['data'] + )) + + +ax3.plot( + dscram['E_photon']['data']/1e3, + dscram['emis_tot']['data'][ind,:], + 'b-', + label='tot' + ) + +if True: + ax3.plot( + dscram['E_photon']['data']/1e3, + dscram['emis_ff']['data'][ind,:], + 'b--', + label='ff' + ) + + ax3.plot( + dscram['E_photon']['data']/1e3, + dscram['emis_fb']['data'][ind,:], + 'b:', + label='fb' + ) + + ax3.plot( + dscram['E_photon']['data']/1e3, + dscram['emis_bb']['data'][ind,:], + 'b-.', + label='bb' + ) + +if False: + + # List of ions of interest + ss = dscram['emis_ion']['name_long'] + neles = ast.literal_eval(ss[ss.index('['):]) + for nn, nele in enumerate(neles): + ax3.plot( + dscram['E_photon']['data']/1e3, + ( + dscram['emis_ion']['data'][ind,:,nn] + *dscram['Xz']['data'][ind,Znuc-nele] # [1/ion] --> [1/atom], applies coronal charge balance + ), + '-.', + color = colors[nn], + label='nele=%i'%(nele) + ) + + +if False: + ax3.plot( + dscram['E_photon']['data']/1e3, + dscram['emis_tot']['data'][ind,:]*ff_fil_fine, + 'r-', + label='filtered' + ) + + +ax3.grid('on') +ax3.set_yscale('log') +ax3.set_xscale('log') +ax3.set_ylim(1e-34, 1e-27) + +ax3.set_xlabel('Photon energy [keV]') +ax3.set_ylabel(r'Emissivity [$J*cm^3/s/eV/atom/electron$]') + +leg = ax3.legend( + loc = 'upper right', + labelcolor='linecolor', + borderaxespad=0, + ) + + +ax3.set_title( + r'%s; $T_e$=%0.1f keV'%( + elem, Te_eV/1e3 + ) + ) + +plt.tight_layout() + + + + +### --- Plots the spectrum at v. Te --- ### + + +fig4, ax4 = plt.subplots( + #figsize = (10,8) + ) + + +cmap = get_cmap('rainbow') +norm = Normalize( + vmin=dscram['Te']['data'].min()/1e3, + vmax=dscram['Te']['data'].max()/1e3, + ) + +for ind_te in np.arange(len(dscram['Te']['data'])): + + color = cmap(norm( + dscram['Te']['data'][ind_te]/1e3 + )) + + ax4.plot( + dscram['E_photon']['data']/1e3, + dscram['emis_tot']['data'][ind_te,:], + color = color + ) + +sm = plt.cm.ScalarMappable(norm=norm, cmap=cmap) +sm.set_array([]) # required for Matplotlib < 3.7 +cbar = plt.colorbar(sm, ax=ax4) +cbar.set_label(r'$T_e$ [keV]') + +ax4.grid('on') +ax4.set_yscale('log') +ax4.set_xscale('log') + +ax4.set_ylim(1e-34, 1e-27) + +ax4.set_xlabel('Photon energy [keV]') +ax4.set_ylabel(r'Emissivity [$J*cm^3/s/eV/atom/electron$]') + + +ax4.set_title( + r'%s; SCRAM v8.6'%( + elem + ) + ) + +plt.tight_layout() + + diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/260625_FLYCHK_H_data.npz b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/260625_FLYCHK_H_data.npz new file mode 100644 index 0000000..f4ed6c9 --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/260625_FLYCHK_H_data.npz @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:75f3590ba5d5e77c3bc1ccc75ab3fba898cf0e1690283556dadfae35001f5f0e +size 4402856 diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/260625_FLYCHK_He_data.npz b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/260625_FLYCHK_He_data.npz new file mode 100644 index 0000000..409644a --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/260625_FLYCHK_He_data.npz @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:ff030933fc66befe3f7708fafd8a273705a946efafa4871bd9f1793346f5c625 +size 4402928 diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/260715_SCRAM86_H_data.npz b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/260715_SCRAM86_H_data.npz new file mode 100644 index 0000000..f672ffb --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/260715_SCRAM86_H_data.npz @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:49a127881ca75af4d9345134b92e21db5dfcd96a438a06abf7bd4c427618bd2c +size 16163897 diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/260715_SCRAM86_He_data.npz b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/260715_SCRAM86_He_data.npz new file mode 100644 index 0000000..df280ee --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/260715_SCRAM86_He_data.npz @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:549b950760824c3bdc96124793f31a6ffd3122d58c5c6982b8155a754e612776 +size 16164497 diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/HXR_MPP200.png b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/HXR_MPP200.png new file mode 100644 index 0000000..da5f7ef --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/HXR_MPP200.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:3f838e6f4cb33bb665d0f092810716858e46f6af9383804365ed46a29be16682 +size 1727023 diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/HXR_MPP200_zoom.png b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/HXR_MPP200_zoom.png new file mode 100644 index 0000000..cda7a7f --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/HXR_MPP200_zoom.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:d740b14d239753cb278fa11f7318f99b1324c5e88cc929c3111768b7a43d7016 +size 1157098 diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/Inversion_BOLO_DMS_ptcam_dvezinet_20260709-205535.npz b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/Inversion_BOLO_DMS_ptcam_dvezinet_20260709-205535.npz new file mode 100644 index 0000000..cb020c9 --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/Inversion_BOLO_DMS_ptcam_dvezinet_20260709-205535.npz @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:1b6b9e24650487eabbde34d2edf20aa14c7b799899fbdd05143824fd7524b8f9 +size 46197071 diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/Inversion_HXRVA_dvezinet_20260713-154229.npz b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/Inversion_HXRVA_dvezinet_20260713-154229.npz new file mode 100644 index 0000000..362926e --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/Inversion_HXRVA_dvezinet_20260713-154229.npz @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:1bfb2f34a5ebaad0d06bffe636fdaf4b45e0b14335b0c47c352b2ed7c847ae99 +size 4353458 diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/Inversion_HXRVA_ptcam_dvezinet_20260709-204104.npz b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/Inversion_HXRVA_ptcam_dvezinet_20260709-204104.npz new file mode 100644 index 0000000..11e469f --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/Inversion_HXRVA_ptcam_dvezinet_20260709-204104.npz @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:1817fbae99d413f678f048fa2ecc78ae3446c2f25aa13bf4c3c99fbbd1dfa74b +size 67110740 diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/OMP000U.png b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/OMP000U.png new file mode 100644 index 0000000..0c9cd6a --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/OMP000U.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:4a79be23956c374ac995195d9d34ebd55f9732b7c433d4e52b0c175ff02d37eb +size 295849 diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/OMP000U_Zoom.png b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/OMP000U_Zoom.png new file mode 100644 index 0000000..d50353e --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/OMP000U_Zoom.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:7a21a06cc19be69c68983a103278cd8bc1d96d1dd1c00a80d7b6a3d005a53153 +size 329352 diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/OMP000U_Zoom_inside.png b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/OMP000U_Zoom_inside.png new file mode 100644 index 0000000..c53e3cf --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/OMP000U_Zoom_inside.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:7b2d378ed873bdf221f213f3c839fd97b4e95c54b1dc4a3f9436f673436c9b87 +size 404726 diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/TFG_Config_ExpSPARC_SPARC-V2_sh00000_Vers1.8.18.npz b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/TFG_Config_ExpSPARC_SPARC-V2_sh00000_Vers1.8.18.npz new file mode 100644 index 0000000..3a4e9f3 --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/TFG_Config_ExpSPARC_SPARC-V2_sh00000_Vers1.8.18.npz @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:1657048de43f24a1a062c424b18941b27537f22e24240f64c06fa1c9f466cdde +size 965951 diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/d2cross_Ee0100eV-10MeV-80log_Eph100eV-10MeV-81log_ntheta61_EH.npz b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/d2cross_Ee0100eV-10MeV-80log_Eph100eV-10MeV-81log_ntheta61_EH.npz new file mode 100644 index 0000000..6bee4e3 --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/d2cross_Ee0100eV-10MeV-80log_Eph100eV-10MeV-81log_ntheta61_EH.npz @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:88fa2297b575de69336f358ad34ee7a7f8e1970d2e7cf4c7b09a742941f54ecc +size 3169754 diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/d2cross_Ee01eV-100MeV-80log_Eph1eV-100MeV-81log_ntheta61_EH.npz b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/d2cross_Ee01eV-100MeV-80log_Eph1eV-100MeV-81log_ntheta61_EH.npz new file mode 100644 index 0000000..b71c7ea --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/d2cross_Ee01eV-100MeV-80log_Eph1eV-100MeV-81log_ntheta61_EH.npz @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:ff78b2cd9e8cec94e666e32ccc9b9d48dc061f21758cf7ec39785e754718cac0 +size 3169340 diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/d2cross_phi_Ee01eV-100MeV-240log_Eph1eV-100MeV-241log_nthetaph61_nthetae060_EH.npz b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/d2cross_phi_Ee01eV-100MeV-240log_Eph1eV-100MeV-241log_nthetaph61_nthetae060_EH.npz new file mode 100644 index 0000000..58ba71c --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/d2cross_phi_Ee01eV-100MeV-240log_Eph1eV-100MeV-241log_nthetaph61_nthetae060_EH.npz @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:86f104d3274f6fc9a7d569053c32d281a8e8799594086c34079cf6f61f11733b +size 1693566266 diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/d2cross_phi_Ee01eV-100MeV-80log_Eph1eV-100MeV-81log_nthetaph61_nthetae060_EH.npz b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/d2cross_phi_Ee01eV-100MeV-80log_Eph1eV-100MeV-81log_nthetaph61_nthetae060_EH.npz new file mode 100644 index 0000000..17025fb --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/d2cross_phi_Ee01eV-100MeV-80log_Eph1eV-100MeV-81log_nthetaph61_nthetae060_EH.npz @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:9f59473986eda33f68cef72def3901eebe80c858bceef9dd57af9de0aa143b50 +size 189742906 diff --git a/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/responsivities.npz b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/responsivities.npz new file mode 100644 index 0000000..f96f57c --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/responsivities.npz @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:14fefdfe2f2517470d72db879f9306d5a18020bb706c8b7a73ff2ac02d55840a +size 367243 diff --git a/articles/__init__.py b/articles/__init__.py new file mode 100644 index 0000000..9d06607 --- /dev/null +++ b/articles/__init__.py @@ -0,0 +1 @@ +from . import RSI_2026_RunawayBremsstrahlungDetection