diff --git a/pyproject.toml b/pyproject.toml index 4a12828..9977621 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -40,6 +40,9 @@ docs = [ "myst-parser>=2.0", "nbsphinx>=0.9", ] +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 e439089..111c5ac 100644 --- a/src/NeSST/__init__.py +++ b/src/NeSST/__init__.py @@ -1,6 +1,10 @@ ################## from .constants import * from .core import * -from .dress_interface import * + +try: + from .dress_interface import * +except ImportError: + pass from .fitting import * from .time_of_flight import * diff --git a/src/NeSST/core.py b/src/NeSST/core.py index 4c94559..c99e394 100644 --- a/src/NeSST/core.py +++ b/src/NeSST/core.py @@ -99,6 +99,41 @@ def QBallabio(Ein: npt.NDArray, mean: float, variance: float) -> npt.NDArray: return spec +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: + 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..ba7eb67 100644 --- a/src/NeSST/dress_interface.py +++ b/src/NeSST/dress_interface.py @@ -1,22 +1,95 @@ -import dress +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 diff --git a/tests/test_dress.py b/tests/test_dress.py new file mode 100644 index 0000000..fd9bd89 --- /dev/null +++ b/tests/test_dress.py @@ -0,0 +1,94 @@ +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 % +# 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)