diff --git a/tests/test_nuclide.py b/tests/test_nuclide.py index b3b761d7..8f3124c5 100644 --- a/tests/test_nuclide.py +++ b/tests/test_nuclide.py @@ -1,6 +1,7 @@ # Copyright 2024, Battelle Energy Alliance, LLC All Rights Reserved. import pytest -from hypothesis import assume, given, note, strategies as st, settings +from hypothesis import assume, example, given, note, strategies as st, settings +import numpy as np import montepy @@ -383,10 +384,56 @@ def test_particle_str(_): class TestNucleus: + @staticmethod + def _semi_emper_binding_energy_p_nucleon(Z, A): + """ + Calculates the base nuclide's binding energy using the semi-empirical mass-formula. + + Based on Wikipedia: + """ + CONSTANTS = { + "A_V": 15.75, + "A_S": -17.8, + "A_C": -0.711, + "A_A": -23.7, + "A_P": 11.18, + } + """ + Units MeV + based on Rohlf + """ + FUNCTIONS = { + "A_V": lambda z, a: a, + "A_S": lambda z, a: np.pow(a, 2.0 / 3), + "A_C": lambda z, a: (z * (z - 1)) / (np.pow(a, 1.0 / 3)), + "A_A": lambda z, a: (a - 2 * z) ** 2 / a, + } + energy = 0.0 + for (name, constant), func in zip(CONSTANTS.items(), FUNCTIONS.values()): + energy += constant * func(Z, A) + # handle even odd stuff + N = A - Z + delta = CONSTANTS["A_P"] / (np.pow(A, 1 / 2)) + even_n = N % 2 == 0 + even_z = Z % 2 == 0 + parity_match = even_z == even_z + # 0 for parity mismatch + if not parity_match: + return energy / A + energy += delta * 1 if even_n else -1 + return energy / A + + @example(Z=97, A=185, meta=2) @given(Z=st.integers(1, 99), A=st.integers(0, 300), meta=st.integers(0, 4)) def test_nucleus_init_eq_hash(_, Z, A, meta): # avoid metastable elemental assume((A == 0) == (meta == 0)) + # use SEMF to rule out non-physical nuclides + # start caring after iron + if Z > 26 and A > 0: + assume( + _._semi_emper_binding_energy_p_nucleon(Z, A) > 7.5 + ) # 7.5 MeV per nucleon cutoff nucleus = Nucleus(Element(Z), A, meta) assert nucleus.Z == Z assert nucleus.A == A