From c914e08a74025af7fd1d406a512c4e9c0edec39b Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sat, 28 Feb 2026 18:10:50 +0000 Subject: [PATCH 1/4] Initial plan From 0acc05ab8b1e9ad4ad1804592f25fac13bad6902 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Sat, 28 Feb 2026 18:18:25 +0000 Subject: [PATCH 2/4] Add DRESS integration for DT and DD thermal spectra Co-authored-by: aidancrilly <20111148+aidancrilly@users.noreply.github.com> --- pyproject.toml | 3 ++ src/NeSST/__init__.py | 5 +- src/NeSST/core.py | 31 ++++++++++++ src/NeSST/dress_interface.py | 88 +++++++++++++++++++++++++++++---- tests/test_dress.py | 95 ++++++++++++++++++++++++++++++++++++ 5 files changed, 212 insertions(+), 10 deletions(-) create mode 100644 tests/test_dress.py diff --git a/pyproject.toml b/pyproject.toml index 7d2352d..55b0dae 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -34,6 +34,9 @@ dependencies = [ test = [ "pytest", "jupyter", "matplotlib" ] +dress = [ + "pydress" +] [tool.setuptools.package-data] "NeSST.data" = ["*.dat", "*.txt","*.csv","*.json"] diff --git a/src/NeSST/__init__.py b/src/NeSST/__init__.py index d545839..b37fa73 100644 --- a/src/NeSST/__init__.py +++ b/src/NeSST/__init__.py @@ -1,6 +1,9 @@ ################## from .constants import * from .core import * -from .dress_interface import * +try: + from .dress_interface import * +except ImportError: + pass import NeSST.fitting import NeSST.time_of_flight \ No newline at end of file diff --git a/src/NeSST/core.py b/src/NeSST/core.py index 6b6a9a5..b2a427b 100644 --- a/src/NeSST/core.py +++ b/src/NeSST/core.py @@ -63,6 +63,37 @@ def QBallabio(Ein : npt.NDArray, mean : float, variance : float) -> npt.NDArray: spec = np.exp(-2.0*Ebar*(np.sqrt(Ein)-np.sqrt(Ebar))**2/sig2)/norm return spec +def QDress_DT(Ein : npt.NDArray, T_D : float, T_T : float = None, n_samples : int = int(1e6)) -> npt.NDArray: + """Calculates the DT primary neutron spectrum using the DRESS Monte Carlo code. + + Args: + Ein (numpy.array): array of energy bin edges on which to compute the spectrum (eV) + T_D (float): temperature of deuterons in eV + T_T (float, optional): temperature of tritons in eV. If None, uses T_D for both. + n_samples (int): number of Monte Carlo samples (default 1e6) + + Returns: + numpy.array: normalised DT spectrum (1/eV) evaluated at the bin centres of Ein + """ + from NeSST.dress_interface import DRESS_DT_spec + if T_T is None: + T_T = T_D + return DRESS_DT_spec(T_D, T_T, n_samples, Ein) + +def QDress_DD(Ein : npt.NDArray, Tion : float, n_samples : int = int(1e6)) -> npt.NDArray: + """Calculates the DD primary neutron spectrum using the DRESS Monte Carlo code. + + Args: + Ein (numpy.array): array of energy bin edges on which to compute the spectrum (eV) + Tion (float): temperature of the ions in eV + n_samples (int): number of Monte Carlo samples (default 1e6) + + Returns: + numpy.array: normalised DD spectrum (1/eV) evaluated at the bin centres of Ein + """ + from NeSST.dress_interface import DRESS_DD_spec + return DRESS_DD_spec(Tion, n_samples, Ein) + # TT spectral shape def dNdE_TT(E : npt.NDArray, Tion : float) -> npt.NDArray: """Calculates the TT primary spectrum with Doppler broadening effect as diff --git a/src/NeSST/dress_interface.py b/src/NeSST/dress_interface.py index 2bf37dc..7d15e6e 100644 --- a/src/NeSST/dress_interface.py +++ b/src/NeSST/dress_interface.py @@ -1,22 +1,92 @@ -import dress +import numpy as np +import numpy.typing as npt -def DRESS_DT_spec_single_T(T,n_samples,bins): - # Conversion to keV units of DRESS - T_keV = T/1e3 - bins /= 1e3 +try: + import dress + dress_available = True +except ImportError: + dress_available = False + +def _check_dress(): + if not dress_available: + raise ImportError("pydress is required for DRESS interface functions. " + "Install it with: pip install pydress") + +def DRESS_DT_spec(T_D : float, T_T : float, n_samples : int, bins : npt.NDArray) -> npt.NDArray: + """Computes DT primary neutron spectrum using DRESS with separate D and T temperatures. + + Args: + T_D (float): temperature of deuterons in eV + T_T (float): temperature of tritons in eV + n_samples (int): number of Monte Carlo samples + bins (numpy.array): energy bin edges in eV + + Returns: + numpy.array: normalised DT spectrum (1/eV) on the bin centres + """ + _check_dress() + T_D_keV = T_D / 1e3 + T_T_keV = T_T / 1e3 + bins_keV = bins / 1e3 reaction = dress.reactions.DTNHe4Reaction() spec_calc = dress.SpectrumCalculator(reaction, n_samples=n_samples) + dist_a = dress.dists.MaxwellianDistribution(T_D_keV, spec_calc.reactant_a.particle) + dist_b = dress.dists.MaxwellianDistribution(T_T_keV, spec_calc.reactant_b.particle) + + spec_calc.reactant_a.v = dist_a.sample(spec_calc.n_samples) + spec_calc.reactant_b.v = dist_b.sample(spec_calc.n_samples) + + spec = spec_calc(bins=bins_keV, normalize=True) + + # Convert from unit-sum to unit-area (1/eV) + bins_w = bins[1:] - bins[:-1] + + return spec / bins_w + +def DRESS_DT_spec_single_T(T : float, n_samples : int, bins : npt.NDArray) -> npt.NDArray: + """Computes DT primary neutron spectrum using DRESS with a single ion temperature. + + Wrapper around DRESS_DT_spec for the case where D and T share the same temperature. + + Args: + T (float): ion temperature in eV + n_samples (int): number of Monte Carlo samples + bins (numpy.array): energy bin edges in eV + + Returns: + numpy.array: normalised DT spectrum (1/eV) on the bin centres + """ + return DRESS_DT_spec(T, T, n_samples, bins) + +def DRESS_DD_spec(T : float, n_samples : int, bins : npt.NDArray) -> npt.NDArray: + """Computes DD primary neutron spectrum using DRESS. + + Args: + T (float): temperature of deuterons in eV + n_samples (int): number of Monte Carlo samples + bins (numpy.array): energy bin edges in eV + + Returns: + numpy.array: normalised DD spectrum (1/eV) on the bin centres + """ + _check_dress() + T_keV = T / 1e3 + bins_keV = bins / 1e3 + + reaction = dress.reactions.DDNHe3Reaction() + spec_calc = dress.SpectrumCalculator(reaction, n_samples=n_samples) + dist_a = dress.dists.MaxwellianDistribution(T_keV, spec_calc.reactant_a.particle) dist_b = dress.dists.MaxwellianDistribution(T_keV, spec_calc.reactant_b.particle) spec_calc.reactant_a.v = dist_a.sample(spec_calc.n_samples) spec_calc.reactant_b.v = dist_b.sample(spec_calc.n_samples) - spec = spec_calc(bins=bins,normalize=True) + spec = spec_calc(bins=bins_keV, normalize=True) - # Ensure unit area for eV - bins_w = 1e3*(bins[1:]-bins[:-1]) + # Convert from unit-sum to unit-area (1/eV) + bins_w = bins[1:] - bins[:-1] - return spec/bins_w \ No newline at end of file + return spec / bins_w \ No newline at end of file diff --git a/tests/test_dress.py b/tests/test_dress.py new file mode 100644 index 0000000..70cde00 --- /dev/null +++ b/tests/test_dress.py @@ -0,0 +1,95 @@ +import pytest +import numpy as np + +pytest.importorskip("dress", reason="pydress is required for DRESS tests") + +import NeSST as nst + + +# Number of Monte Carlo samples used in tests - large enough for accurate +# statistics but small enough to keep test runtime reasonable. +_N_SAMPLES = int(1e5) + +# Tolerance for relative comparison of mean energy against Ballabio (%) +_MEAN_RTOL = 0.01 # 1 % +# Tolerance for relative comparison of standard deviation against Ballabio (%) +_STDDEV_RTOL = 0.05 # 5 % + + +def _spectrum_moments(Ein, spec): + """Return the mean and standard deviation of a normalised spectrum.""" + Emid = 0.5 * (Ein[:-1] + Ein[1:]) + dE = Ein[1:] - Ein[:-1] + total = np.sum(spec * dE) + mean = np.sum(Emid * spec * dE) / total + stddev = np.sqrt(np.sum((Emid - mean) ** 2 * spec * dE) / total) + return mean, stddev + + +def test_QDress_DT_mean_matches_Ballabio(): + """DRESS DT mean neutron energy should match Ballabio within 1 %.""" + Tion = 5e3 # 5 keV in eV + mean_ball, stddev_ball, _ = nst.DTprimspecmoments(Tion) + + Ein = np.linspace(mean_ball - 10 * stddev_ball, mean_ball + 10 * stddev_ball, 500) + spec = nst.QDress_DT(Ein, Tion, n_samples=_N_SAMPLES) + + mean_dress, _ = _spectrum_moments(Ein, spec) + + assert mean_dress == pytest.approx(mean_ball, rel=_MEAN_RTOL) + + +def test_QDress_DT_stddev_matches_Ballabio(): + """DRESS DT spectral width should match Ballabio within 5 %.""" + Tion = 5e3 # 5 keV in eV + mean_ball, stddev_ball, _ = nst.DTprimspecmoments(Tion) + + Ein = np.linspace(mean_ball - 10 * stddev_ball, mean_ball + 10 * stddev_ball, 500) + spec = nst.QDress_DT(Ein, Tion, n_samples=_N_SAMPLES) + + _, stddev_dress = _spectrum_moments(Ein, spec) + + assert stddev_dress == pytest.approx(stddev_ball, rel=_STDDEV_RTOL) + + +def test_QDress_DT_separate_temperatures(): + """Using different D and T temperatures should produce a different spectrum + from a single-temperature run.""" + Tion = 5e3 # 5 keV in eV + T_D = 5e3 + T_T = 10e3 # tritons hotter than deuterons + + mean_ball, stddev_ball, _ = nst.DTprimspecmoments(Tion) + Ein = np.linspace(mean_ball - 10 * stddev_ball, mean_ball + 10 * stddev_ball, 500) + + spec_single = nst.QDress_DT(Ein, Tion, n_samples=_N_SAMPLES) + spec_diff = nst.QDress_DT(Ein, T_D, T_T, n_samples=_N_SAMPLES) + + # The spectra should not be identical when temperatures differ + assert not np.allclose(spec_single, spec_diff) + + +def test_QDress_DD_mean_matches_Ballabio(): + """DRESS DD mean neutron energy should match Ballabio within 1 %.""" + Tion = 5e3 # 5 keV in eV + mean_ball, stddev_ball, _ = nst.DDprimspecmoments(Tion) + + Ein = np.linspace(mean_ball - 10 * stddev_ball, mean_ball + 10 * stddev_ball, 500) + spec = nst.QDress_DD(Ein, Tion, n_samples=_N_SAMPLES) + + mean_dress, _ = _spectrum_moments(Ein, spec) + + assert mean_dress == pytest.approx(mean_ball, rel=_MEAN_RTOL) + + +def test_QDress_DD_stddev_matches_Ballabio(): + """DRESS DD spectral width should match Ballabio within 5 %.""" + Tion = 5e3 # 5 keV in eV + mean_ball, stddev_ball, _ = nst.DDprimspecmoments(Tion) + + Ein = np.linspace(mean_ball - 10 * stddev_ball, mean_ball + 10 * stddev_ball, 500) + spec = nst.QDress_DD(Ein, Tion, n_samples=_N_SAMPLES) + + _, stddev_dress = _spectrum_moments(Ein, spec) + + assert stddev_dress == pytest.approx(stddev_ball, rel=_STDDEV_RTOL) From 010c07aae28c948068a628551a0cb51b5a94bad4 Mon Sep 17 00:00:00 2001 From: aidancrilly Date: Sat, 28 Feb 2026 19:19:26 +0000 Subject: [PATCH 3/4] Lint --- src/NeSST/__init__.py | 1 + src/NeSST/core.py | 9 +++++++-- src/NeSST/dress_interface.py | 17 ++++++++++------- tests/test_dress.py | 13 ++++++------- 4 files changed, 24 insertions(+), 16 deletions(-) diff --git a/src/NeSST/__init__.py b/src/NeSST/__init__.py index 37f535e..111c5ac 100644 --- a/src/NeSST/__init__.py +++ b/src/NeSST/__init__.py @@ -1,6 +1,7 @@ ################## from .constants import * from .core import * + try: from .dress_interface import * except ImportError: diff --git a/src/NeSST/core.py b/src/NeSST/core.py index 54c6649..c99e394 100644 --- a/src/NeSST/core.py +++ b/src/NeSST/core.py @@ -98,7 +98,8 @@ def QBallabio(Ein: npt.NDArray, mean: float, variance: float) -> npt.NDArray: spec = np.exp(-2.0 * Ebar * (np.sqrt(Ein) - np.sqrt(Ebar)) ** 2 / sig2) / norm return spec -def QDress_DT(Ein : npt.NDArray, T_D : float, T_T : float = None, n_samples : int = int(1e6)) -> npt.NDArray: + +def QDress_DT(Ein: npt.NDArray, T_D: float, T_T: float | None = None, n_samples: int = int(1e6)) -> npt.NDArray: """Calculates the DT primary neutron spectrum using the DRESS Monte Carlo code. Args: @@ -111,11 +112,13 @@ def QDress_DT(Ein : npt.NDArray, T_D : float, T_T : float = None, n_samples : in numpy.array: normalised DT spectrum (1/eV) evaluated at the bin centres of Ein """ from NeSST.dress_interface import DRESS_DT_spec + if T_T is None: T_T = T_D return DRESS_DT_spec(T_D, T_T, n_samples, Ein) -def QDress_DD(Ein : npt.NDArray, Tion : float, n_samples : int = int(1e6)) -> npt.NDArray: + +def QDress_DD(Ein: npt.NDArray, Tion: float, n_samples: int = int(1e6)) -> npt.NDArray: """Calculates the DD primary neutron spectrum using the DRESS Monte Carlo code. Args: @@ -127,8 +130,10 @@ def QDress_DD(Ein : npt.NDArray, Tion : float, n_samples : int = int(1e6)) -> np numpy.array: normalised DD spectrum (1/eV) evaluated at the bin centres of Ein """ from NeSST.dress_interface import DRESS_DD_spec + return DRESS_DD_spec(Tion, n_samples, Ein) + # TT spectral shape def dNdE_TT(E: npt.NDArray, Tion: float) -> npt.NDArray: """Calculates the TT primary spectrum with Doppler broadening effect as diff --git a/src/NeSST/dress_interface.py b/src/NeSST/dress_interface.py index 7d15e6e..ba7eb67 100644 --- a/src/NeSST/dress_interface.py +++ b/src/NeSST/dress_interface.py @@ -1,18 +1,19 @@ -import numpy as np import numpy.typing as npt try: import dress + dress_available = True except ImportError: dress_available = False + def _check_dress(): if not dress_available: - raise ImportError("pydress is required for DRESS interface functions. " - "Install it with: pip install pydress") + raise ImportError("pydress is required for DRESS interface functions. Install it with: pip install pydress") + -def DRESS_DT_spec(T_D : float, T_T : float, n_samples : int, bins : npt.NDArray) -> npt.NDArray: +def DRESS_DT_spec(T_D: float, T_T: float, n_samples: int, bins: npt.NDArray) -> npt.NDArray: """Computes DT primary neutron spectrum using DRESS with separate D and T temperatures. Args: @@ -45,7 +46,8 @@ def DRESS_DT_spec(T_D : float, T_T : float, n_samples : int, bins : npt.NDArray) return spec / bins_w -def DRESS_DT_spec_single_T(T : float, n_samples : int, bins : npt.NDArray) -> npt.NDArray: + +def DRESS_DT_spec_single_T(T: float, n_samples: int, bins: npt.NDArray) -> npt.NDArray: """Computes DT primary neutron spectrum using DRESS with a single ion temperature. Wrapper around DRESS_DT_spec for the case where D and T share the same temperature. @@ -60,7 +62,8 @@ def DRESS_DT_spec_single_T(T : float, n_samples : int, bins : npt.NDArray) -> np """ return DRESS_DT_spec(T, T, n_samples, bins) -def DRESS_DD_spec(T : float, n_samples : int, bins : npt.NDArray) -> npt.NDArray: + +def DRESS_DD_spec(T: float, n_samples: int, bins: npt.NDArray) -> npt.NDArray: """Computes DD primary neutron spectrum using DRESS. Args: @@ -89,4 +92,4 @@ def DRESS_DD_spec(T : float, n_samples : int, bins : npt.NDArray) -> npt.NDArray # Convert from unit-sum to unit-area (1/eV) bins_w = bins[1:] - bins[:-1] - return spec / bins_w \ No newline at end of file + return spec / bins_w diff --git a/tests/test_dress.py b/tests/test_dress.py index 70cde00..fd9bd89 100644 --- a/tests/test_dress.py +++ b/tests/test_dress.py @@ -1,17 +1,16 @@ -import pytest import numpy as np +import pytest pytest.importorskip("dress", reason="pydress is required for DRESS tests") import NeSST as nst - # Number of Monte Carlo samples used in tests - large enough for accurate # statistics but small enough to keep test runtime reasonable. _N_SAMPLES = int(1e5) # Tolerance for relative comparison of mean energy against Ballabio (%) -_MEAN_RTOL = 0.01 # 1 % +_MEAN_RTOL = 0.01 # 1 % # Tolerance for relative comparison of standard deviation against Ballabio (%) _STDDEV_RTOL = 0.05 # 5 % @@ -55,15 +54,15 @@ def test_QDress_DT_stddev_matches_Ballabio(): def test_QDress_DT_separate_temperatures(): """Using different D and T temperatures should produce a different spectrum from a single-temperature run.""" - Tion = 5e3 # 5 keV in eV - T_D = 5e3 - T_T = 10e3 # tritons hotter than deuterons + Tion = 5e3 # 5 keV in eV + T_D = 5e3 + T_T = 10e3 # tritons hotter than deuterons mean_ball, stddev_ball, _ = nst.DTprimspecmoments(Tion) Ein = np.linspace(mean_ball - 10 * stddev_ball, mean_ball + 10 * stddev_ball, 500) spec_single = nst.QDress_DT(Ein, Tion, n_samples=_N_SAMPLES) - spec_diff = nst.QDress_DT(Ein, T_D, T_T, n_samples=_N_SAMPLES) + spec_diff = nst.QDress_DT(Ein, T_D, T_T, n_samples=_N_SAMPLES) # The spectra should not be identical when temperatures differ assert not np.allclose(spec_single, spec_diff) From 64c8aeda5eb9b59a6514fe3c37642a1c31eb6b5b Mon Sep 17 00:00:00 2001 From: aidancrilly Date: Sun, 1 Mar 2026 10:08:12 +0000 Subject: [PATCH 4/4] Lint --- src/NeSST/__init__.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/NeSST/__init__.py b/src/NeSST/__init__.py index 37f535e..111c5ac 100644 --- a/src/NeSST/__init__.py +++ b/src/NeSST/__init__.py @@ -1,6 +1,7 @@ ################## from .constants import * from .core import * + try: from .dress_interface import * except ImportError: