From 8936f651cd570e05f6506c823944aa3608086cf8 Mon Sep 17 00:00:00 2001 From: shimwell Date: Mon, 21 Oct 2024 09:13:27 +0100 Subject: [PATCH] end of train journey --- .gitignore | 162 +++++++++++ Cargo.lock | 171 +++++++++++ Cargo.toml | 12 + README.md | 63 +++++ pyproject.toml | 18 ++ src/lib.rs | 393 ++++++++++++++++++++++++++ tests/test_average_neutron_energy.py | 58 ++++ tests/test_relative_reaction_rates.py | 41 +++ 8 files changed, 918 insertions(+) create mode 100644 .gitignore create mode 100644 Cargo.lock create mode 100644 Cargo.toml create mode 100644 README.md create mode 100644 pyproject.toml create mode 100644 src/lib.rs create mode 100644 tests/test_average_neutron_energy.py create mode 100644 tests/test_relative_reaction_rates.py diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..82f9275 --- /dev/null +++ b/.gitignore @@ -0,0 +1,162 @@ +# Byte-compiled / optimized / DLL files +__pycache__/ +*.py[cod] +*$py.class + +# C extensions +*.so + +# Distribution / packaging +.Python +build/ +develop-eggs/ +dist/ +downloads/ +eggs/ +.eggs/ +lib/ +lib64/ +parts/ +sdist/ +var/ +wheels/ +share/python-wheels/ +*.egg-info/ +.installed.cfg +*.egg +MANIFEST + +# PyInstaller +# Usually these files are written by a python script from a template +# before PyInstaller builds the exe, so as to inject date/other infos into it. +*.manifest +*.spec + +# Installer logs +pip-log.txt +pip-delete-this-directory.txt + +# Unit test / coverage reports +htmlcov/ +.tox/ +.nox/ +.coverage +.coverage.* +.cache +nosetests.xml +coverage.xml +*.cover +*.py,cover +.hypothesis/ +.pytest_cache/ +cover/ + +# Translations +*.mo +*.pot + +# Django stuff: +*.log +local_settings.py +db.sqlite3 +db.sqlite3-journal + +# Flask stuff: +instance/ +.webassets-cache + +# Scrapy stuff: +.scrapy + +# Sphinx documentation +docs/_build/ + +# PyBuilder +.pybuilder/ +target/ + +# Jupyter Notebook +.ipynb_checkpoints + +# IPython +profile_default/ +ipython_config.py + +# pyenv +# For a library or package, you might want to ignore these files since the code is +# intended to run in multiple environments; otherwise, check them in: +# .python-version + +# pipenv +# According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control. +# However, in case of collaboration, if having platform-specific dependencies or dependencies +# having no cross-platform support, pipenv may install dependencies that don't work, or not +# install all needed dependencies. +#Pipfile.lock + +# poetry +# Similar to Pipfile.lock, it is generally recommended to include poetry.lock in version control. +# This is especially recommended for binary packages to ensure reproducibility, and is more +# commonly ignored for libraries. +# https://python-poetry.org/docs/basic-usage/#commit-your-poetrylock-file-to-version-control +#poetry.lock + +# pdm +# Similar to Pipfile.lock, it is generally recommended to include pdm.lock in version control. +#pdm.lock +# pdm stores project-wide configurations in .pdm.toml, but it is recommended to not include it +# in version control. +# https://pdm.fming.dev/latest/usage/project/#working-with-version-control +.pdm.toml +.pdm-python +.pdm-build/ + +# PEP 582; used by e.g. github.com/David-OConnor/pyflow and github.com/pdm-project/pdm +__pypackages__/ + +# Celery stuff +celerybeat-schedule +celerybeat.pid + +# SageMath parsed files +*.sage.py + +# Environments +.env +.venv +env/ +venv/ +ENV/ +env.bak/ +venv.bak/ + +# Spyder project settings +.spyderproject +.spyproject + +# Rope project settings +.ropeproject + +# mkdocs documentation +/site + +# mypy +.mypy_cache/ +.dmypy.json +dmypy.json + +# Pyre type checker +.pyre/ + +# pytype static type analyzer +.pytype/ + +# Cython debug symbols +cython_debug/ + +# PyCharm +# JetBrains specific template is maintained in a separate JetBrains.gitignore that can +# be found at https://github.com/github/gitignore/blob/main/Global/JetBrains.gitignore +# and can be added to the global gitignore or merged into this file. For a more nuclear +# option (not recommended) you can uncomment the following to ignore the entire idea folder. +#.idea/ diff --git a/Cargo.lock b/Cargo.lock new file mode 100644 index 0000000..347eed1 --- /dev/null +++ b/Cargo.lock @@ -0,0 +1,171 @@ +# This file is automatically @generated by Cargo. +# It is not intended for manual editing. +version = 3 + +[[package]] +name = "autocfg" +version = "1.3.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "0c4b4d0bd25bd0b74681c0ad21497610ce1b7c91b1022cd21c80c6fbdd9476b0" + +[[package]] +name = "cfg-if" +version = "1.0.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "baf1de4339761588bc0619e3cbc0120ee582ebb74b53b4efbf79117bd2da40fd" + +[[package]] +name = "fusion-reaction-rate" +version = "0.1.0" +dependencies = [ + "pyo3", +] + +[[package]] +name = "heck" +version = "0.5.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "2304e00983f87ffb38b55b444b5e3b60a884b5d30c0fca7d82fe33449bbe55ea" + +[[package]] +name = "indoc" +version = "2.0.5" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b248f5224d1d606005e02c97f5aa4e88eeb230488bcc03bc9ca4d7991399f2b5" + +[[package]] +name = "libc" +version = "0.2.158" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "d8adc4bb1803a324070e64a98ae98f38934d91957a99cfb3a43dcbc01bc56439" + +[[package]] +name = "memoffset" +version = "0.9.1" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "488016bfae457b036d996092f6cb448677611ce4449e970ceaf42695203f218a" +dependencies = [ + "autocfg", +] + +[[package]] +name = "once_cell" +version = "1.19.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3fdb12b2476b595f9358c5161aa467c2438859caa136dec86c26fdd2efe17b92" + +[[package]] +name = "portable-atomic" +version = "1.7.0" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "da544ee218f0d287a911e9c99a39a8c9bc8fcad3cb8db5959940044ecfc67265" + +[[package]] +name = "proc-macro2" +version = "1.0.86" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5e719e8df665df0d1c8fbfd238015744736151d4445ec0836b8e628aae103b77" +dependencies = [ + "unicode-ident", +] + +[[package]] +name = "pyo3" +version = "0.22.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "831e8e819a138c36e212f3af3fd9eeffed6bf1510a805af35b0edee5ffa59433" +dependencies = [ + "cfg-if", + "indoc", + "libc", + "memoffset", + "once_cell", + "portable-atomic", + "pyo3-build-config", + "pyo3-ffi", + "pyo3-macros", + "unindent", +] + +[[package]] +name = "pyo3-build-config" +version = "0.22.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "1e8730e591b14492a8945cdff32f089250b05f5accecf74aeddf9e8272ce1fa8" +dependencies = [ + "once_cell", + "target-lexicon", +] + +[[package]] +name = "pyo3-ffi" +version = "0.22.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "5e97e919d2df92eb88ca80a037969f44e5e70356559654962cbb3316d00300c6" +dependencies = [ + "libc", + "pyo3-build-config", +] + +[[package]] +name = "pyo3-macros" +version = "0.22.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "eb57983022ad41f9e683a599f2fd13c3664d7063a3ac5714cae4b7bee7d3f206" +dependencies = [ + "proc-macro2", + "pyo3-macros-backend", + "quote", + "syn", +] + +[[package]] +name = "pyo3-macros-backend" +version = "0.22.2" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "ec480c0c51ddec81019531705acac51bcdbeae563557c982aa8263bb96880372" +dependencies = [ + "heck", + "proc-macro2", + "pyo3-build-config", + "quote", + "syn", +] + +[[package]] +name = "quote" +version = "1.0.37" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "b5b9d34b8991d19d98081b46eacdd8eb58c6f2b201139f7c5f643cc155a633af" +dependencies = [ + "proc-macro2", +] + +[[package]] +name = "syn" +version = "2.0.77" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "9f35bcdf61fd8e7be6caf75f429fdca8beb3ed76584befb503b1569faee373ed" +dependencies = [ + "proc-macro2", + "quote", + "unicode-ident", +] + +[[package]] +name = "target-lexicon" +version = "0.12.16" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "61c41af27dd6d1e27b1b16b489db798443478cef1f06a660c96db617ba5de3b1" + +[[package]] +name = "unicode-ident" +version = "1.0.12" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "3354b9ac3fae1ff6755cb6db53683adb661634f67557942dea4facebec0fee4b" + +[[package]] +name = "unindent" +version = "0.2.3" +source = "registry+https://github.com/rust-lang/crates.io-index" +checksum = "c7de7d73e1754487cb58364ee906a499937a0dfabd86bcb980fa99ec8c8fa2ce" diff --git a/Cargo.toml b/Cargo.toml new file mode 100644 index 0000000..a2b4fd0 --- /dev/null +++ b/Cargo.toml @@ -0,0 +1,12 @@ +[package] +name = "fusion-reaction-rate" +version = "0.1.0" +edition = "2021" + +# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html +[lib] +name = "fusion_neutron_utils" +crate-type = ["cdylib"] + +[dependencies] +pyo3 = "0.22.2" \ No newline at end of file diff --git a/README.md b/README.md new file mode 100644 index 0000000..a1f0500 --- /dev/null +++ b/README.md @@ -0,0 +1,63 @@ +A package for calculating neutron properties from DT and DD fusion reactions. + +This package accurately calculates the neutron energies and distributions by accounting for the plasma temperature. + +- mean neutron energy +- neutron energy standard deviation () +- thermal reactivity + +A Python :snake: package with a Rust :crab: backend. + +The package makes use of the [Sadler–Van Belle formula](https://doi.org/10.1016/j.fusengdes.2012.02.025) and [Bosch-Hale](https://doi.org/10.1088/0029-5515%2F32%2F4%2FI07) parametrization for reactivity. For energy distributions [Ballabio](https://doi.org/10.1088/0029-5515/38/11/310) is used. + +This package has been inspired by the [NeSST package](https://github.com/aidancrilly/NeSST). + +## Install + +```bash +pip install fusion_neutron_utils +``` + +## Usage + +To find the average neutron energy and the standard deviation of that energy. + +These results could then be input into [openmc.stats.Normal](https://docs.openmc.org/en/v0.12.1/pythonapi/generated/openmc.stats.Normal.html) distribution for a neutron source term in OpenMC + +```python +from fusion_neutron_utils import neutron_energy_mean_and_std_dev +neutron_energy_mean_and_std_dev( + reaction='D+T=n+a', + ion_temperature=30e3, + temperature_units='eV', + neutron_energy_units='eV' +) +>>>>(14021063.802198425, 413861.375751198) +``` + + +The relative reaction rates can be found for the different reactions, this can be useful for setting the relative ```Source.strength``` in OpenMC. Relative reaction rates returns the DT, DD (n+He3), and DD (p+T) reaction rates in that order. Note the DD (p+T) does not emit a neutron but is there for completeness. +```python +from fusion_neutron_utils import relative_reaction_rates +relative_reaction_rates( + ion_temperature=30e3, + temperature_units='eV', + dt_fraction=0.1, + dd_fraction=0.9, +) +>>>[0.9989900631769298, 0.0005595866500573114, 0.00045035017301275736] +``` + +The reactivity can also be found, this can be useful for finding the relative reaction rate for D+D or D+T at a specific temperature + +```python +from fusion_neutron_utils import reactivity +reactivity( + ion_temperature=30e3, + temperature_units='eV', + reactivity_units='m^3/s', + reaction='D+D=p+T', + equation='Bosch-Hale' +) +>>>2.2356373732343755e-53 +``` diff --git a/pyproject.toml b/pyproject.toml new file mode 100644 index 0000000..66b86eb --- /dev/null +++ b/pyproject.toml @@ -0,0 +1,18 @@ +[build-system] +requires = ["maturin>=1.7,<2.0"] +build-backend = "maturin" + +[project] +name = "fusion_neutron_utils" +requires-python = ">=3.8" +classifiers = [ + "Programming Language :: Rust", + "Programming Language :: Python :: Implementation :: CPython", + "Programming Language :: Python :: Implementation :: PyPy", +] +dynamic = ["version"] +[tool.maturin] +features = ["pyo3/extension-module"] + +[project.optional-dependencies] +test = ["pytest", "scipy==1.13.1", "NeSST"] diff --git a/src/lib.rs b/src/lib.rs new file mode 100644 index 0000000..aca2991 --- /dev/null +++ b/src/lib.rs @@ -0,0 +1,393 @@ +use pyo3::prelude::*; + +#[pyfunction(signature = (ion_temperature, temperature_units=None, neutron_energy_units=None, reaction=None))] +/// Calculate the average neutron energy for a given ion temperature and reaction. +#[pyo3(text_signature = "(ion_temperature, temperature_units='eV', neutron_energy_units='eV', reaction='D+T=n+a')")] +fn neutron_energy_mean_and_std_dev( + ion_temperature: f64, + temperature_units: Option<&str>, + neutron_energy_units: Option<&str>, + reaction: Option<&str>, +) -> PyResult<(f64, f64)> { + // values from Ballabio paper + let (a_1, a_2, a_3, a_4, mean) = match reaction { + Some("D+D=n+He3") => (4.69515, -0.040729, 0.47, 0.81844, 2.4495e6), + Some("D+T=n+a") => (5.30509, 2.4736e-3, 1.84, 1.3818, 14.021e6), + _ => return Err(PyErr::new::("reaction must be either 'D+D=n+He3' or 'D+T=n+a'")), + }; + + let ion_temperature_kev: f64 = scale_temperature_units_to_kev(ion_temperature, temperature_units); // Ballabio equation accepts KeV units + + // units of mean_delta are in ev + let mean_delta = a_1 * ion_temperature_kev.powf(2.0 / 3.0) / (1.0 + a_2 * ion_temperature_kev.powf(a_3)) + a_4 * ion_temperature_kev; + + let mean_adjusted = mean + mean_delta; + print!("{}", mean); + + let mean_scaled = scale_energy_in_kev_to_requested_units(mean_adjusted/1e3, neutron_energy_units); + + let (w_0, a_1, a_2, a_3, a_4) = match reaction { + Some("D+D=n+He3") => (82.542, 1.7013e-3, 0.16888, 0.49, 7.9460e-4), + Some("D+T=n+a") => (177.259, 5.1068e-4, 7.6223e-3, 1.78, 8.7691e-5), + _ => unreachable!(), // This case is already handled above + }; + + let delta = a_1 * ion_temperature_kev.powf(2.0 / 3.0) / (1.0 + a_2 * ion_temperature_kev.powf(a_3)) + a_4 * ion_temperature_kev; + + // 2.3548200450309493 on the line below comes from equation 2* math.sqrt(math.log(2)*2) + let variance = ((w_0 * (1.0 + delta)).powi(2) * ion_temperature_kev) / 2.3548200450309493_f64.powi(2); + // let variance = variance * 1e6; // converting keV^2 back to eV^2 + let std_dev = variance.sqrt(); + let std_dev = scale_energy_in_kev_to_requested_units(std_dev, neutron_energy_units); + + Ok((mean_scaled, std_dev)) +} + + + +#[pyfunction(signature = (ion_temperature, temperature_units=None, reactivity_units=None, reaction=None, equation=None))] +/// Bosch-Hale parametrization of D+T thermal reactivity, assuming a Maxwellian ion +/// temperature distribution. If ion_temperature_kev is given in keV, the returned +/// is in m^3/s. This is valid for 0.2 keV <= ion_temperature_kev <= 100 keV. +/// +/// D + T -> T(3.56 MeV) + n(14.03 MeV) +/// +/// Args: +/// ion_temperature_kev (float): Ion temperature. +/// temperature_units +/// +/// Returns: +/// float: The thermal reactivity in m^3/s. +/// +/// Source: +/// Sec. 5.2, Eqn. (12) - (14), Table VII in +/// H.-s. Bosch and G. M. Dale, +/// "Improved Formulas for Fusion Cross-Sections and Thermal Reactivities", +/// Nucl. Fusion 32, 611 (1992) +#[pyo3(text_signature = "(ion_temperature, temperature_units='eV', reactivity_units='m^3/s', reaction='D+T=n+a', equation='Bosch-Hale')")] +fn reactivity( + ion_temperature: f64, + temperature_units: Option<&str>, + reactivity_units: Option<&str>, + reaction: Option<&str>, + equation: Option<&str>, +) -> PyResult { + + if ion_temperature <= 0.0 { + return Err(PyErr::new::("Ion temperature must be positive and non-zero")); + } + + let reaction = reaction.unwrap_or("D+T=n+a"); + let equation = equation.unwrap_or("Bosch-Hale"); + + let ion_temperature_kev: f64 = scale_temperature_units_to_kev(ion_temperature, temperature_units); + + let sigma_thermal_reactivity_scaled = if equation == "Bosch-Hale"{ + let (C1, C2, C3, C4, C5, C6, C7, GAMOV, MRC2) = if reaction == "D+T=n+a" { + ( + 1.17302e-9, 1.51361e-2, 7.51886e-2, 4.60643e-3, 1.35000e-2, -1.06750e-4, 1.36600e-5, + 34.3827, 1_124_656.0, + ) + } else if reaction == "D+D=p+T" { + ( + 5.65718e-12, 3.41267e-3, 1.99167e-3, 0.0, 1.05060e-5, 0.0, 0.0, + 31.3970, 937_814.0, + ) + } else if reaction == "D+D=n+He3" { + ( + 5.43360e-12, 5.85778e-3, 7.68222e-3, 0.0, -2.96400e-6, 0.0, 0.0, + 31.3970, 937_814.0, + ) + } else { + return Err(PyErr::new::("Only 'D+T=n+a', 'D+D=p+T', and 'D+D=n+He3' reactions are supported")); + }; + + let sigma_thermal_reactivity = bosch_and_hale_equations(C1, C2, C3, C4, C5, C6, C7, GAMOV, MRC2, ion_temperature_kev)?; + let scaling_factor = scale_reactivity_units(sigma_thermal_reactivity, reactivity_units); + Ok(sigma_thermal_reactivity*scaling_factor) + }else if equation == "Sadler-Van Belle"{ + if reaction != "D+T=n+a" { + return Err(PyErr::new::("Only 'D+T=n+a' reaction is supported for 'Sadler-Van Belle' equation")); + } + let sigma_thermal_reactivity = sadler_van_belle(ion_temperature_kev)?; + let scaling_factor = scale_reactivity_units(sigma_thermal_reactivity, reactivity_units); + Ok(sigma_thermal_reactivity*scaling_factor) + }else{ + panic!("Only 'Bosch-Hale' and 'Sadler-Van Belle' equations are supported"); + }; + sigma_thermal_reactivity_scaled + +} + + +#[pyfunction(signature = (ion_temperature, temperature_units=None, dt_fraction=None, dd_fraction=None, equation=None))] +#[pyo3(text_signature = "(ion_temperature, temperature_units='eV', dt_fraction=0.5, dd_fraction=0.5 equation='Bosch-Hale')")] +/// Calculate the relative reaction rates for given ion temperature and isotope fractions. +/// +/// Parameters +/// ---------- +/// ion_temperature : float +/// The ion temperature. +/// temperature_units : str, optional +/// The units of the ion temperature. Default is 'eV'. +/// dt_fraction : float, optional +/// The fraction of DT reactions. Default is 0.5. +/// dd_fraction : float, optional +/// The fraction of DD reactions. Default is 0.5. +/// equation : str, optional +/// The equation to use for reactivity calculations. Default is 'Bosch-Hale'. +/// +/// Returns +/// ------- +/// List[float] +/// A list containing the relative reaction rates for DT, DD (n+He3), and DD (p+T) reactions. +/// +/// Examples +/// -------- +/// >>> relative_reaction_rates(10.0) +/// [dt_reactivity, dd_reactivity_1, dd_reactivity_2] +/// +/// >>> relative_reaction_rates(10.0, temperature_units='K', dt_fraction=0.3, dd_fraction=0.7, equation='Custom-Equation') +/// [dt_reactivity, dd_reactivity_1, dd_reactivity_2] +fn relative_reaction_rates( + ion_temperature: f64, + temperature_units: Option<&str>, + dt_fraction: Option, + dd_fraction: Option, + equation: Option<&str>, +) -> Result, PyErr> { + + let ion_temperature_kev: f64 = scale_temperature_units_to_kev(ion_temperature, temperature_units); + + let dt_fraction = dt_fraction.unwrap_or(0.5); + let dd_fraction = dd_fraction.unwrap_or(0.5); + + let equation = equation.unwrap_or("Bosch-Hale"); + + let total_fraction = dt_fraction + dd_fraction; + let tol : f64 = 0.000001; + if !(total_fraction > 1. - 0.000001 && total_fraction < 1. + 0.000001) { + return Err(PyErr::new::("The dt_fraction + dd_fraction do not sum to 1.0 and are not within a small tolerance (+/-0.000001)")); + } + + let dt_reactivity = reactivity( + ion_temperature_kev, + Some("keV"), + Some("m^3/s"), + Some("D+T=n+a"), + Some("Bosch-Hale") + )?; + + let dd_reactivity_1 = reactivity( + ion_temperature_kev, + Some("keV"), + Some("m^3/s"), + Some("D+D=n+He3"), + Some("Bosch-Hale") + )?; + + let dd_reactivity_2 = reactivity( + ion_temperature_kev, + Some("keV"), + Some("m^3/s"), + Some("D+D=p+T"), + Some("Bosch-Hale") + )?; + + let total_reactivity = dt_reactivity * dt_fraction + dd_reactivity_1 * dd_fraction + dd_reactivity_2 * dd_fraction; + let dt_normalized = dt_reactivity * dt_fraction / total_reactivity; + let dd1_normalized = dd_reactivity_1 * dd_fraction / total_reactivity; + let dd2_normalized = dd_reactivity_2 * dd_fraction / total_reactivity; + + Ok(vec![dt_normalized, dd1_normalized, dd2_normalized]) + +} + + +fn sadler_van_belle(ion_temperature: f64) -> Result { + let c = [ + 2.5663271e-18, + 19.983026, + 2.5077133e-2, + 2.5773408e-3, + 6.1880463e-5, + 6.6024089e-2, + 8.1215505e-3, + ]; + + let u = 1.0 - ion_temperature * (c[2] + ion_temperature * (c[3] - c[4] * ion_temperature)) + / (1.0 + ion_temperature * (c[5] + c[6] * ion_temperature)); + + let val = c[0] + * ((-c[1] * (u / ion_temperature).powf(1.0 / 3.0)).exp()) + / (u.powf(5.0 / 6.0) * ion_temperature.powf(2.0 / 3.0)); + + Ok(val) +} + +fn bosch_and_hale_equations(C1: f64, C2: f64, C3: f64, C4: f64, C5: f64, C6: f64, C7: f64, GAMOV: f64, MRC2: f64, ion_temperature_kev: f64) -> Result { + // Equation 13 + let theta: f64 = ion_temperature_kev * (1.0 - (ion_temperature_kev * (C2 + ion_temperature_kev * (C4 + ion_temperature_kev * C6))) / (1.0 + ion_temperature_kev * (C3 + ion_temperature_kev * (C5 + ion_temperature_kev * C7)))) .powi(-1); + + // Equation 14 + let xi: f64 = (GAMOV.powi(2) / (4.0 * theta)).powf(1.0 / 3.0); + + // Equation 12 + let sigma_thermal_reactivity: f64 = C1 * theta * (xi / (MRC2 * ion_temperature_kev.powi(3))).sqrt() * (-3.0 * xi).exp(); + + Ok(sigma_thermal_reactivity * 1.0e-6) +} + +fn scale_reactivity_units(sigma_thermal_reactivity: f64, reactivity_units: Option<&str>) -> f64 { + match reactivity_units.unwrap_or("m^3/s") { + "m^3/s" => sigma_thermal_reactivity * 1.0e-6, + "cm^3/s" => sigma_thermal_reactivity, + "mm^3/s" => sigma_thermal_reactivity * 1.0e3, + _ => panic!("Invalid reaction rate units, accepted values are 'm^3/s', 'cm^3/s', 'mm^3/s'") + } +} + +fn scale_temperature_units_to_kev(ion_temperature: f64, temperature_units: Option<&str>) -> f64 { + match temperature_units.unwrap_or("eV") { + "keV" => ion_temperature, + "eV" => ion_temperature * 1e-3, + "MeV" => ion_temperature * 1e3, + "GeV" => ion_temperature * 1e6, + _ => panic!("Invalid temperature units, accepted values are 'eV', 'keV', 'MeV' or 'GeV'") + } +} + +fn scale_energy_in_kev_to_requested_units(energy_in_kev: f64, temperature_units: Option<&str>) -> f64 { + match temperature_units.unwrap_or("eV") { + "eV" => energy_in_kev * 1e3, // converting keV to eV + "keV" => energy_in_kev, // converting keV to MeV + "MeV" => energy_in_kev / 1e3, // converting keV to MeV + "GeV" => energy_in_kev / 1e6, // converting keV to GeV + _ => panic!("Unsupported temperature units, accepted values are 'eV', 'MeV', 'GeV'"), + } +} + +/// A Python module implemented in Rust. +#[pymodule] +fn fusion_neutron_utils(m: &Bound<'_, PyModule>) -> PyResult<()> { + m.add_function(wrap_pyfunction!(reactivity, m)?)?; + m.add_function(wrap_pyfunction!(relative_reaction_rates, m)?)?; + m.add_function(wrap_pyfunction!(neutron_energy_mean_and_std_dev, m)?)?; + Ok(()) +} + + + + +#[cfg(test)] +mod tests { + use super::*; + + #[test] + fn test_scale_temperature_units_to_kev_eV() { + let ion_temperature = 1000.0; // 1000 eV + let temperature_units = Some("eV"); + let result = scale_temperature_units_to_kev(ion_temperature, temperature_units); + assert_eq!(result, 1.0); // 1 keV + } + + #[test] + fn test_scale_temperature_units_to_kev_keV() { + let ion_temperature = 1.0; // 1 keV + let temperature_units = Some("keV"); + let result = scale_temperature_units_to_kev(ion_temperature, temperature_units); + assert_eq!(result, 1.0); // 1 keV + } + + #[test] + fn test_scale_temperature_units_to_kev_MeV() { + let ion_temperature = 0.001; // 1 MeV + let temperature_units = Some("MeV"); + let result = scale_temperature_units_to_kev(ion_temperature, temperature_units); + assert_eq!(result, 1.0); // 1 keV + } + + #[test] + fn test_scale_temperature_units_to_kev_GeV() { + let ion_temperature = 0.000001; // 1 GeV + let temperature_units = Some("GeV"); + let result = scale_temperature_units_to_kev(ion_temperature, temperature_units); + assert_eq!(result, 1.0); // 1 keV + } + + #[test] + #[should_panic(expected = "Invalid temperature units")] + fn test_scale_temperature_units_to_kev_invalid_units() { + let ion_temperature = 1000.0; + let temperature_units = Some("K"); + scale_temperature_units_to_kev(ion_temperature, temperature_units); + } +} + + +#[cfg(test)] +mod tests2 { + use super::*; + + #[test] + fn test_scale_to_ev() { + let energy_in_kev = 1.0; + let target_unit = "eV"; + let result = scale_energy_in_kev_to_requested_units(energy_in_kev, Some(target_unit)); + assert_eq!(result, 1000.0); // 1 keV = 1000 eV + } + + #[test] + fn test_scale_to_mev() { + let energy_in_kev = 1000.0; + let target_unit = "MeV"; + let result = scale_energy_in_kev_to_requested_units(energy_in_kev, Some(target_unit)); + assert_eq!(result, 1.0); // 1000 keV = 1 MeV + } + + #[test] + fn test_scale_to_gev() { + let energy_in_kev = 1_000_000.0; + let target_unit = "GeV"; + let result = scale_energy_in_kev_to_requested_units(energy_in_kev, Some(target_unit)); + assert_eq!(result, 1.0); // 1,000,000 keV = 1 GeV + } + + #[test] + #[should_panic(expected = "Unsupported temperature units, accepted values are 'eV', 'MeV', 'GeV'")] + fn test_invalid_unit() { + let energy_in_kev = 1.0; + let target_unit = "invalid"; + scale_energy_in_kev_to_requested_units(energy_in_kev, Some(target_unit)); + } +} + +#[cfg(test)] +mod tests3 { + use super::*; + + #[test] + fn test_scale_to_m3_per_s() { + let sigma_thermal_reactivity = 1.0; + let reactivity_units = Some("m^3/s"); + let result = scale_reactivity_units(sigma_thermal_reactivity, reactivity_units); + assert_eq!(result, 1.0 * 1.0e-6); // 1.0 * 1.0e-6 = 1.0e-6 + } + + #[test] + fn test_scale_to_cm3_per_s() { + let sigma_thermal_reactivity = 1.0; + let reactivity_units = Some("cm^3/s"); + let result = scale_reactivity_units(sigma_thermal_reactivity, reactivity_units); + assert_eq!(result, 1.0); // 1.0 cm^3/s = 1.0 cm^3/s + } + + + #[test] + #[should_panic(expected = "Invalid reaction rate units, accepted values are 'm^3/s', 'cm^3/s', 'mm^3/s'")] + fn test_invalid_unit() { + let sigma_thermal_reactivity = 1.0; + let reactivity_units = Some("invalid"); + scale_reactivity_units(sigma_thermal_reactivity, reactivity_units); + } +} \ No newline at end of file diff --git a/tests/test_average_neutron_energy.py b/tests/test_average_neutron_energy.py new file mode 100644 index 0000000..3f7ee3c --- /dev/null +++ b/tests/test_average_neutron_energy.py @@ -0,0 +1,58 @@ +from fusion_neutron_utils import neutron_energy_mean_and_std_dev +from pytest import approx + +def test_mean_energy(): + mean, std_dev = neutron_energy_mean_and_std_dev( + reaction='D+D=n+He3', + ion_temperature=20e3, + temperature_units='eV', + neutron_energy_units='eV' + ) + assert mean == approx(2.5e6, abs = 0.2e6) + + mean, std_dev = neutron_energy_mean_and_std_dev( + reaction='D+D=n+He3', + ion_temperature=20, + temperature_units='keV', + neutron_energy_units='MeV' + ) + assert mean == approx(2.5, abs = 0.2) + + mean, std_dev = neutron_energy_mean_and_std_dev( + reaction='D+T=n+a', + ion_temperature=20e3, + temperature_units='eV', + ) + assert mean == approx(14.06e6, abs = 0.04e6) + +def test_mean_energy_increases_with_ion_temperature(): + mean_cold, std_dev = neutron_energy_mean_and_std_dev( + reaction='D+D=n+He3', + ion_temperature=20e3, + temperature_units='eV', + ) + + mean_hot, std_dev = neutron_energy_mean_and_std_dev( + reaction='D+D=n+He3', + ion_temperature=40e3, + temperature_units='eV', + ) + assert mean_cold < mean_hot + +def test_mean_energy_units(): + mean_kev, std_dev = neutron_energy_mean_and_std_dev( + reaction='D+T=n+a', + ion_temperature=30, + temperature_units='keV', + neutron_energy_units='keV' + ) + + mean_ev, std_dev = neutron_energy_mean_and_std_dev( + reaction='D+T=n+a', + ion_temperature=30e3, + temperature_units='eV', + neutron_energy_units='eV' + ) + assert mean_kev == mean_ev/1e3 + +# def test_ \ No newline at end of file diff --git a/tests/test_relative_reaction_rates.py b/tests/test_relative_reaction_rates.py new file mode 100644 index 0000000..1a97dc0 --- /dev/null +++ b/tests/test_relative_reaction_rates.py @@ -0,0 +1,41 @@ +import pytest +from fusion_neutron_utils import relative_reaction_rates + +def test_relative_reaction_rates_default_fractions(): + result = relative_reaction_rates(ion_temperature=10.0) + assert isinstance(result, list) + assert len(result) == 3 + assert all(isinstance(x, float) for x in result) + +def test_relative_reaction_rates_custom_fractions(): + result = relative_reaction_rates(ion_temperature=10.0, dt_fraction=0.3, dd_fraction=0.7) + assert isinstance(result, list) + assert len(result) == 3 + assert all(isinstance(x, float) for x in result) + +def test_relative_reaction_rates_temperature_units(): + result = relative_reaction_rates(ion_temperature=10.0, temperature_units="keV") + assert isinstance(result, list) + assert len(result) == 3 + assert all(isinstance(x, float) for x in result) + +def test_relative_reaction_rates_custom_equation(): + result = relative_reaction_rates(ion_temperature=10.0, equation="Bosch-Hale") + assert isinstance(result, list) + assert len(result) == 3 + assert all(isinstance(x, float) for x in result) + +def test_relative_reaction_rates_invalid_fractions(): + with pytest.raises(TypeError): + relative_reaction_rates(ion_temperature=10.0, dt_fraction="invalid", dd_fraction="invalid") + +def test_relative_reaction_rates_zero_fractions(): + with pytest.raises(ValueError): + relative_reaction_rates(ion_temperature=10.0, dt_fraction=0.0, dd_fraction=0.0) + +def test_relative_reaction_rates_non_sum_1_fractions(): + with pytest.raises(ValueError): + relative_reaction_rates(ion_temperature=10.0, dt_fraction=0.5, dd_fraction=0.4) + +if __name__ == "__main__": + pytest.main() \ No newline at end of file