Added test to compare with NeSST results

.github/workflows/python-publish.yml

@@ -60,7 +60,7 @@ jobs:
           set -e
           python3 -m venv .venv
           source .venv/bin/activate
-          pip install fusion_neutron_utils --find-links dist --force-reinstall
+          pip install fusion_neutron_utils nesst --find-links dist --force-reinstall
           pip install pytest
           pytest
       - name: pytest
@@ -76,7 +76,7 @@ jobs:
             pip3 install -U pip pytest
           run: |
             set -e
-            pip3 install fusion_neutron_utils --find-links dist --force-reinstall
+            pip3 install fusion_neutron_utils nesst --find-links dist --force-reinstall
             pytest
 
   windows:
@@ -112,7 +112,7 @@ jobs:
           set -e
           python3 -m venv .venv
           source .venv/Scripts/activate
-          pip install fusion_neutron_utils --find-links dist --force-reinstall
+          pip install fusion_neutron_utils nesst --find-links dist --force-reinstall
           pip install pytest
           pytest
 

src/lib.rs

@@ -11,18 +11,17 @@ fn neutron_energy_mean_and_std_dev(
 ) -> 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),
+        Some("D+D=n+He3") => (4.69515, -0.040729, 0.47, 0.81844, 2.4486858678216934e6),
+        Some("D+T=n+a") => (5.30509, 0.0024736, 1.84, 1.3818, 14028394.744466662),
         _ => return Err(PyErr::new::<pyo3::exceptions::PyValueError, _>("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;
+    // units of mean_delta are in put into ev with the 1000 multiplication
+    let mean_delta = 1000.0 *( a_1 * ion_temperature_kev.powf(0.66666666) / (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);
 

tests/test_average_neutron_energy.py

@@ -1,6 +1,32 @@
+import NeSST as nst
 from fusion_neutron_utils import neutron_energy_mean_and_std_dev
 from pytest import approx
 
+def test_mean_energy_with_nesst():
+
+    ion_temperature = 19e3
+
+    nesst_dt_mean, nest_dt_std_dev, nest_dt_var = nst.DTprimspecmoments(ion_temperature)
+    nesst_dd_mean, nest_dd_std_dev, nest_dd_var = nst.DDprimspecmoments(ion_temperature)
+
+    fnu_dd_mean, fnu_dd_std_dev = neutron_energy_mean_and_std_dev(
+        reaction='D+D=n+He3',
+        ion_temperature=ion_temperature,
+        temperature_units='eV',
+        neutron_energy_units='eV'
+    )
+
+    fnu_dt_mean, fnu_dt_std_dev = neutron_energy_mean_and_std_dev(
+        reaction='D+T=n+a',
+        ion_temperature=ion_temperature,
+        temperature_units='eV',
+        neutron_energy_units='eV'
+    )
+    assert nesst_dt_mean == approx(fnu_dt_mean, rel=1e-6)
+    assert nesst_dd_mean == approx(fnu_dd_mean, rel=1e-6)
+    assert nest_dt_std_dev == approx(fnu_dt_std_dev, rel=1e-6)
+    assert nest_dd_std_dev == approx(fnu_dd_std_dev, rel=1e-6)
+
 def test_mean_energy():
     mean, std_dev = neutron_energy_mean_and_std_dev(
         reaction='D+D=n+He3',