Skip to content
Open
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
49 changes: 48 additions & 1 deletion tests/test_nuclide.py
Original file line number Diff line number Diff line change
@@ -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

Expand Down Expand Up @@ -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: <https://en.wikipedia.org/wiki/Semi-empirical_mass_formula>
"""
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
Expand Down