diff --git a/.gitattributes b/.gitattributes new file mode 100644 index 0000000..540cecf --- /dev/null +++ b/.gitattributes @@ -0,0 +1,12 @@ +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 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..66591aa --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/__init__.py @@ -0,0 +1,9 @@ +from ._main import main +from ._load_spect import main as load_spect +from ._fig01_cross import main as fig01_cross +from ._fig02_dist import main as fig02_dist +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 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..e639099 --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig01_cross.py @@ -0,0 +1,295 @@ + + +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=(15, 7), + 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 = { + 'left': 0.06, 'right': 0.99, + 'bottom': 0.08, 'top': 0.93, + 'wspace': 0.25, 'hspace': 0.10, + } + + fig = plt.figure(figsize=figsize) + + gs = gridspec.GridSpec(ncols=4, nrows=2, **dmargin) + dax = {} + + # -------------- + # prepare axes + # -------------- + + # -------------- + # ax - isolines + + ax = fig.add_subplot( + gs[:, -2:], + 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[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', + ) + + # store + dax['theta_norm0'] = {'handle': ax, 'type': 'isolines'} + + # theta_norm1 + ax = fig.add_subplot( + gs[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', + ) + + # store + dax['theta_norm1'] = {'handle': ax, 'type': 'isolines'} + + # -------------- + # ax - theta_abs + + # theta_abs0 + ax = fig.add_subplot( + gs[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', + ) + + # store + dax['theta_abs0'] = {'handle': ax, 'type': 'isolines'} + + # theta_abs1 + ax = fig.add_subplot( + gs[1, 1], + sharex=dax['theta_norm0']['handle'], + sharey=dax['theta_abs0']['handle'], + ) + ax.set_xlabel( + r"$\theta_{ph}$ (deg)", + size=fontsize, + fontweight='bold', + ) + + # store + dax['theta_abs1'] = {'handle': ax, 'type': 'isolines'} + + # ------------------ + # 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_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/_fig03_validate_ff.py b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig03_validate_ff.py new file mode 100644 index 0000000..e47e109 --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig03_validate_ff.py @@ -0,0 +1,342 @@ + + +import os + + +import numpy as np +import scipy.constants as scpct +import matplotlib.pyplot as plt +import matplotlib.gridspec as gridspec +import astropy.units as asunits +import datastock as ds +import tofu as tf + + +from ._load_spect import main as load_spect +from ._savefig import main as savefig + + +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') + + +# SPECTRAL MODELLING FILES +_LPFE_SPECT = [ + ff for ff in os.listdir(_PATH_INPUTS) + if ff.endswith('_data.npz') + and any([ss in ff for ss in ['_SCRAM86_', '_FLYCHK_']]) +] +_DPFE_SPECT = { + ff.split('_')[-2]: os.path.join(_PATH_INPUTS, ff) + for ff in _LPFE_SPECT +} + + +# CROSS-SECTION FILES +_DPFE_DCROSS = { + 'EH0': os.path.join( + _PATH_INPUTS, + 'd2cross_phi_Ee01eV-100MeV-240log_Eph1eV-100MeV-241log_nthetaph61_nthetae060_EH.npz', + ), + 'EH1': os.path.join( + _PATH_INPUTS, + 'd2cross_phi_Ee01eV-100MeV-80log_Eph1eV-100MeV-81log_nthetaph61_nthetae060_EH.npz' + ), +} + + +# ##################################################### +# ##################################################### +# Main +# ##################################################### + + +def main( + d2cross_phi='EH1', + elements='H', + ne_m3=1e19, + Te_plot=None, + # plot + figsize=(6, 8), + fontsize=14, + # save + path_save=None, + pfe_save=None, + # unused + **kwdargs, +): + """ Validate Bremsstrahlung vs SCRAM and FLYCHK + + """ + + # -------------- + # inputs + # -------------- + + # d2cross_phi + d2cross_phi = ds._generic_check._check_var( + d2cross_phi, 'd2cross_phi', + types=str, + allowed=list(_DPFE_DCROSS.keys()), + default='EH0', + ) + d2cross_phi = _DPFE_DCROSS[d2cross_phi] + + # -------------- + # load elements + # -------------- + + dplasma = load_spect( + dmix='H', + ne_m3=ne_m3, + ) + + # extract + E_ph = dplasma['common']['E_photon']['data'] + Te = dplasma['common']['Te']['data'] + units = dplasma['emiss_tot']['ff']['units'] + + # -------------- + # load cross + # -------------- + + demiss, ddist, d2cross_phi = tfphysemis.get_xray_thin_integ_dist( + # dist + Te_eV=Te, + ne_m3=ne_m3, + jp_Am2=0, + Zeff=1, + jp_fraction_re=0, + # ---------------- + # cross-section + # tabulated d2cross_phi + d2cross_phi=d2cross_phi, + # d2cross_phi computation + E_ph_eV=dplasma['common']['E_photon']['data'], + # ----------- + # verb + debug=False, + verb=True, + ) + + # check units + units2 = demiss['emiss']['maxwell']['emiss']['units'] + assert asunits.Unit(units) == asunits.Unit(units2) + + # check isotropy + emiss = demiss['emiss']['maxwell']['emiss']['data'] + mean = np.nanmean(emiss, axis=-1) + diff = emiss - mean[:, :, None] + iok = emiss > 0. + c0 = np.abs(diff[iok] / emiss[iok]) + assert np.all(c0 < 1e-2) + + # Interpolate + interp = np.full(dplasma['emiss_tot']['ff']['data'].shape, np.nan) + for ind in np.ndindex(mean.shape[:-1]): + sli0 = ind + (slice(None),) + iok = mean[sli0] > 0. + sli = ind + (iok,) + + interp[sli0] = 10**(np.interp( + np.log10(E_ph), + np.log10(demiss['E_ph_eV']['data'][iok]), + np.log10(mean[sli]), + )) + + # diff + emiss_min = np.minimum(interp, dplasma['emiss_tot']['ff']['data']) + iok = emiss_min > 0. + diff = np.full(interp.shape, np.nan) + diff[iok] = 100 * ( + np.abs(interp - dplasma['emiss_tot']['ff']['data'])[iok] + / emiss_min[iok] + ) + + # -------------- + # prepare axes + # -------------- + + dmargin = { + 'left': 0.13, 'right': 0.99, + 'bottom': 0.10, 'top': 0.92, + 'wspace': 0.25, 'hspace': 0.10, + } + + fig = plt.figure(figsize=figsize) + + gs = gridspec.GridSpec(ncols=1, nrows=2, **dmargin) + dax = {} + + # -------------- + # prepare axes + # -------------- + + # -------------- + # ax - abs + + ax = fig.add_subplot( + gs[0, 0], + xscale='log', + yscale='log', + aspect='auto', + ) + ax.set_xlabel( + r"$E_{ph}$ (keV)", + size=fontsize, + fontweight='bold', + ) + ax.set_ylabel( + r"$\epsilon^{Max}_{ff}$" + f" ({asunits.Unit(units)})", + size=fontsize, + fontweight='bold', + ) + tit = ( + "Validation of Bremstrahlung implemented from EH cross-section\n" + "Maxwellian distribution with " + + r"$j_{P}$" + " = 0 A/m2, " + + r"$n_e$" + f" = {ne_m3:1.0e}" + ) + ax.set_title( + tit, + size=fontsize, + fontweight='bold', + ) + + # store + dax['abs'] = {'handle': ax} + + # -------------- + # ax - diff + + ax0 = ax + ax = fig.add_subplot( + gs[1, 0], + sharex=ax0, + yscale='log', + aspect='auto', + ) + ax.set_xlabel( + r"$E_{ph}$ (keV)", + size=fontsize, + fontweight='bold', + ) + ax.set_ylabel( + "error (\%)", + size=fontsize, + fontweight='bold', + ) + + # store + dax['diff'] = {'handle': ax} + + # standardize + dax = ds._generic_check._check_dax(dax) + + # --------------------- + # Absolute values + # --------------------- + + kax = 'abs' + dcolor = {} + if dax.get(kax) is not None: + ax = dax[kax]['handle'] + + # -------------- + # loop on Te + + for ii, te in enumerate(dplasma['common']['Te']['data']): + + # label + lab = r"$T_e$" + f" = {te*1e-3:3.2f} keV" + + # slice + sli = (ii, slice(None)) + + # loop on elements + for kk, vv in dplasma['emiss'].items(): + + # SCRAM / FLYCHK + l0, = ax.loglog( + E_ph*1e-3, + vv['ff']['data'][sli], + ls='-', + lw=1, + label=lab, + ) + dcolor[ii] = l0.get_color() + + # ff + l0, = ax.loglog( + E_ph*1e-3, + interp[sli], + ls='--', + lw=1, + color=dcolor[ii], + # label=lab, + ) + + vmax = np.nanmax(np.maximum(vv['ff']['data'], interp)) + vmax_plot = 10**np.ceil(np.log10(vmax)) + + ax.legend() + ax.set_xlim(E_ph.min()*1e-3, E_ph.max()*1e-3) + ax.set_ylim(vmax_plot/1e15, vmax_plot) + + # --------------- + # plot diff + # --------------- + + kax = 'diff' + if dax.get(kax) is not None: + ax = dax[kax]['handle'] + + # -------------- + # loop on Te + + for ii, te in enumerate(Te): + + # label + lab = r"$T_e$" + f" = {te*1e-3:3.2f} keV" + + # slice + sli = (ii, slice(None)) + + # loop on elements + for kk, vv in dplasma['emiss'].items(): + + # SCRAM / FLYCHK + l0, = ax.loglog( + E_ph*1e-3, + diff[sli], + ls='-', + lw=1, + color=dcolor[ii], + label=lab, + ) + + ax.set_ylim(1e-4, 1e4) + + # -------------- + # save + # -------------- + + savefig( + fig=fig, + pfe_save=pfe_save, + path_save=path_save, + file=__file__, + ) + + return dax, demiss, ddist, dplasma 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..8b6c981 --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig04_bremsstrahlung.py @@ -0,0 +1,536 @@ +import os + + +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 ._fig02_dist import _DDIST +from ._savefig import main as savefig + + +tfphysemis = tf.physics_tools.electrons.emission + + +# ##################################################### +# ##################################################### +# DEFAULTS +# ##################################################### + + +_PATH_HERE = os.path.dirname(__file__) +_PATH_INPUTS = os.path.join(os.path.dirname(_PATH_HERE), 'inputs') +_PATH_SAVE = os.path.join(os.path.dirname(_PATH_HERE), 'figures') + + +_PFE_D2CROSS_PHI = os.path.join( + _PATH_INPUTS, + 'd2cross_phi_Ee01eV-100MeV-80log_Eph1eV-100MeV-81log_nthetaph61_nthetae060_EH.npz', + # 'd2cross_phi_Ee01eV-100MeV-240log_Eph1eV-100MeV-241log_nthetaph61_nthetae060_EH.npz', +) + + +_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=np.r_[0.1, 5], + Ekin_max_eV=np.r_[100e3, 10e6], + # Te_eV=1e3 * np.linspace(0.1, 2.5, 25), + Te_eV=1e3 * np.linspace(0.1, 2.5, 11), + # jp_fraction_re=np.linspace(0.025, 0.975, 39), + jp_fraction_re=np.linspace(0.1, 0.9, 9), + # 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 + + """ + + # ------------ + # d2cross_phi + # ------------ + + if d2cross_phi is None: + d2cross_phi = _PFE_D2CROSS_PHI + + # ------------ + # 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 pnormW is not None: + ddist['pnormW'] = pnormW + + if Ekin_max_eV is not None: + ddist['Ekin_max_eV'] = Ekin_max_eV + + if pnormW is not None: + ddist['pnormW'] = pnormW + + # shape + ddist['Ekin_max_eV'] = np.atleast_1d(ddist['Ekin_max_eV'])[:, None, None] + ddist['pnormW'] = np.atleast_1d(ddist['pnormW'])[:, None, None] + if ddist['pnormW'].size != ddist['Ekin_max_eV'].size: + raise Exception() + ddist['Te_eV'] = Te_eV[None, :, None] + ddist['jp_fraction_re'] = jp_fraction_re[None, None, :] + + nEkin = ddist['Ekin_max_eV'].shape[0] + + # -------------- + # 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=E_ph_eV, + E_e0_eV=None, + E_e0_eV_npts=None, + theta_e0_vsB_npts=None, + phi_e0_vsB_npts=None, + theta_ph_vsB=None, + # inputs + Z=None, + # hypergeometric parameter + ninf=None, + source=None, + # integration parameters + nthetae=None, + ndphi=None, + # output customization + version_cross=None, + # save / load + save_d2cross_phi=False, + # --------------------- + # optional responsivity + dresponsivity=None, + plot_responsivity_integration=None, + # ----------- + # verb + debug=False, + verb=True, + # ---------------- + # electron distribution + **ddist, + ) + + units = demiss['emiss']['RE']['emiss']['units'] + ne = np.unique(ddist['plasma']['ne_m3']['data'])[0] + jp = np.unique(ddist['plasma']['jp_Am2']['data'])[0] + + tit = ( + r"$n_e$" + f" = {ne:1.0e} /m3\n" + + r"$j_{P,RE}$" + f" = {jp*1e-6:2.1f} MA/m2" + ) + + # -------------- + # cases + # -------------- + + if cases is None: + cases = _CASES + + # -------------- + # Elim + # -------------- + + shape = ddist['plasma']['Te_eV']['data'].shape + Elim = np.full(shape, np.nan) + theta_Elim = 0 + for ii, ind in enumerate(np.ndindex(shape)): + + sli_emiss = ind + (slice(None), 0) + emiss_RE = demiss['emiss']['RE']['emiss']['data'][sli_emiss] + emiss_max = demiss['emiss']['maxwell']['emiss']['data'][sli_emiss] + + ilim = (emiss_RE > emiss_max) + if np.any(ilim): + Elim[ind] = np.min(demiss['E_ph_eV']['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 + + 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]['emiss']['data'][sli] + + # plot + ax.fill_between( + demiss['E_ph_eV']['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_eV']['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_eV']['data'] - cc)) + sli = ic + (iE, slice(None)) + + for kdist in demiss['emiss'].keys(): + emiss_theta = demiss['emiss'][kdist]['emiss']['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_eV']['data'] - cc)) + sli = ic + (iE, slice(None)) + + for kdist in demiss['emiss'].keys(): + emiss_theta = demiss['emiss'][kdist]['emiss']['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..1097486 --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig05_emiss.py @@ -0,0 +1,316 @@ +import os + + +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 main as load_spect +from ._fig04_bremsstrahlung import _PFE_D2CROSS_PHI +from ._savefig import main as savefig + + +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') + + +# ##################################################### +# ##################################################### +# Main +# ##################################################### + + +def main( + d2cross_phi='EH1', + dmix={ + 0: 'H', + 1: {'H': 0.95, 'O': 0.04, 'Fe': 0.01}, + }, + ne_m3=1e19, + # plot + figsize=(15, 8), + fontsize=14, + # save + path_save=None, + pfe_save=None, + # unused + **kwdargs, +): + """ Validate Bremsstrahlung vs SCRAM and FLYCHK + + """ + + # -------------- + # inputs + # -------------- + + nmix = len(dmix) + + # -------------- + # load elements + # -------------- + + dplasma = {} + for ii in dmix.keys(): + dplasma[ii] = load_spect( + dmix=dmix[ii], + ne_m3=ne_m3, + ) + + # extract + E_ph = dplasma[0]['common']['E_photon']['data'] + Te = dplasma[0]['common']['Te']['data'] + units = dplasma[0]['emiss_tot']['ff']['units'] + + # -------------- + # 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=E_ph_eV, + E_e0_eV=None, + E_e0_eV_npts=None, + theta_e0_vsB_npts=None, + phi_e0_vsB_npts=None, + theta_ph_vsB=None, + # inputs + Z=None, + # hypergeometric parameter + ninf=None, + source=None, + # integration parameters + nthetae=None, + ndphi=None, + # output customization + version_cross=None, + # save / load + save_d2cross_phi=False, + # --------------------- + # optional responsivity + dresponsivity=None, + plot_responsivity_integration=None, + # ----------- + # verb + debug=False, + verb=True, + # ---------------- + # electron distribution + **ddist, + ) + + # -------------- + # prepare axes + # -------------- + + dmargin = { + 'left': 0.06, 'right': 0.98, + 'bottom': 0.06, 'top': 0.93, + 'wspace': 0.25, 'hspace': 0.30, + } + + fig = plt.figure(figsize=figsize) + + gs = gridspec.GridSpec(ncols=3, nrows=nmix, **dmargin) + dax = {} + + # ---------------- + # ax - spectra + # ---------------- + + ax0_spect = None + ax0_map = None + for ii in sorted(dmix.keys()): + + ax = fig.add_subplot( + gs[ii, :2], + sharex=ax0_spect, + sharey=ax0_spect, + aspect='auto', + ) + ax.set_ylabel( + r"$\epsilon$" + f' ({units})', + fontsize=fontsize, + fontweight='bold', + ) + if ii == 0: + ax0_spect = ax + elif ii == nmix - 1: + ax.set_xlabel( + r"$E_{ph}$" + ' (keV)', + fontsize=fontsize, + fontweight='bold', + ) + dax[f'spect_{ii}'] = ax + + # ---------------- + # ax - Elim + # ---------------- + + ax = fig.add_subplot( + gs_map[ii, 2:], + aspect='auto', + sharex=ax0, + sharey=ax0, + ) + if ii == 0: + ax0_map = ax + ax.set_title( + tit, + fontsize=fontsize, + fontweight='bold', + ) + elif ii == nmix - 1: + ax.set_xlabel('Te (keV)', fontsize=fontsize, fontweight='bold') + ax.set_ylabel('jp_frac', fontsize=fontsize, fontweight='bold') + + dax[f'Elim_{ii}'] = ax + + # ---------------- + # check dax format + + dax = ds._generic_check._check_dax(dax) + + # -------------- + # plot - spectra + # -------------- + + for ii, vplasma in sorted(dplasma.keys()): + + kax = f"spect_{ii}" + if dax.get(kax) is not None: + ax = dax[kax]['handle'] + + # ---------------------- + # total emiss Maxwellian + + emiss_tot = np.sum( + [vv['data'] for vv in plasma['emiss_tot'].values()], + axis=0, + ) + + ax.plot( + vplasma['common']['E_photon']['data'], + emiss_tot, + color=None, + ls='-', + lw=1, + ) + + # --------- + # ff - RE + + # plot + ax.fill_between( + demiss['E_ph_eV']['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) + + # Elim + iE = np.argmin(np.abs(demiss['E_ph_eV']['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) + + # -------------- + # 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 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..0259354 --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_fig07_responsivities.py @@ -0,0 +1,148 @@ +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', +) + + +# ##################################################### +# ##################################################### +# 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, v0 in dresp.items(): + + # 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-4, 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/_load_spect.py b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_load_spect.py new file mode 100644 index 0000000..3f71dd1 --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_load_spect.py @@ -0,0 +1,331 @@ +import os + + +import numpy as np +import scipy.constants as scpct +import astropy.units as asunits +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 +_LPFE_SPECT = [ + ff for ff in os.listdir(_PATH_INPUTS) + if ff.endswith('_data.npz') + and any([ss in ff for ss in ['_SCRAM86_', '_FLYCHK_']]) +] +_DPFE_SPECT = { + ff.split('_')[-2]: os.path.join(_PATH_INPUTS, ff) + for ff in _LPFE_SPECT +} + + +# PLASMA +_NE = 1e19 # /m3 + + +# ##################################################### +# ##################################################### +# Main +# ##################################################### + + +def main( + dmix=None, + ne_m3=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 + # -------------- + + dmix = _check_mix(dmix) + + # -------------- + # 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'] + for kk, vv in dfiles.items(): + assert np.allclose(vv[cc]['data'], ref) + 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) + 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,) + + # -------------- + # dplasma + # -------------- + + lemiss = ['ff', 'fb', 'bb'] + dplasma = { + 'concentration': dmix, + 'emiss': {k0: {cc: None for cc in lemiss} for k0 in dmix.keys()}, + 'common': dcommon, + } + + # -------------- + # 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 + + +# ##################################################### +# ##################################################### +# Check +# ##################################################### + + +def _check_mix( + dmix=None, +): + + # ----------- + # dmix + # ----------- + + # None + if dmix is None: + dmix = 'H' + + # str + if isinstance(dmix, str): + dmix = {dmix: 1} + + # --------- + # each dict + + lok = sorted(_DPFE_SPECT.keys()) + if isinstance(dmix, dict): + + # check each element + dfail = {} + for k0, v0 in dmix.items(): + + # key + if not (isinstance(k0, str) and k0 in lok): + dfail[k0] = "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 for v0 in dmix.values()]) + + # broadcast + for k0, v0 in dmix.items(): + dmix[k0] = np.broadcast_to(v0, shape) + + # sum = 1 + total = np.sum([v0 for v0 in dmix.values()], axis=0) + for k0, v0 in dmix.items(): + dmix[k0] = v0 / total + + except Exception: + dfail['dmisc'] = "values not broadcastable!" + + # ----------- + # 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) + ) + 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 + if Te is None: + Te = _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..e728e08 --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/code/_load_spect_anis.py @@ -0,0 +1,98 @@ +import os + + +import numpy as np +import astropy.units as asunits + + +from ._load_spect import main as load_spect + + +# ##################################################### +# ##################################################### +# DEFAULTS +# ##################################################### + + +# PATHS +_PATH_HERE = os.path.dirname(__file__) +_PATH_PAPER = os.path.dirname(_PATH_HERE) +_PATH_INPUTS = os.path.join(_PATH_PAPER, 'inputs') + + +# ##################################################### +# ##################################################### +# Main +# ##################################################### + + +def main( + dmix=None, + ne_m3=None, +): + """ For a given impurity mix and ne, Te, return a dict of emissivities + + Emissivities are broken down into 2 distributions + - Maxwell: + - ff: anisotropic + - 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 + # -------------- + + dplasma = load_spect( + dmix=dmix, + ne_m3=ne_m3, + ) + + # -------------- + # anisotropic + # -------------- + + + + # -------------- + # output + # -------------- + + demiss = { + 'emiss': { + 'maxwell': { + 'ff': {}, + 'fb': dplasma['emiss_tot']['fb'], + 'bb': dplasma['emiss_tot']['bb'], + }, + 'RE': { + 'ff': {}, + }, + }, + } + + return demiss 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/_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/figures/fig01_cross.png b/articles/RSI_2026_RunawayBremsstrahlungDetection/figures/fig01_cross.png new file mode 100644 index 0000000..2915d61 --- /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:4167f592a00d053bae903aa0f5886dc5f52bbcd57a9d97598b952c52fb0f39f6 +size 630505 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/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/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/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..811ce40 --- /dev/null +++ b/articles/RSI_2026_RunawayBremsstrahlungDetection/inputs/responsivities.npz @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:4c275a8bc8cb2109b18a6595f136964653d4a4ecfecfdd0f8a28604850a1ebc3 +size 365835 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