added neutronics tools

libra_toolbox/__init__.py

@@ -1,2 +1,3 @@
 from . import tritium
 from . import neutron_detection
+from . import neutronics

libra_toolbox/_version.py

@@ -0,0 +1,16 @@
+# file generated by setuptools_scm
+# don't change, don't track in version control
+TYPE_CHECKING = False
+if TYPE_CHECKING:
+    from typing import Tuple, Union
+    VERSION_TUPLE = Tuple[Union[int, str], ...]
+else:
+    VERSION_TUPLE = object
+
+version: str
+__version__: str
+__version_tuple__: VERSION_TUPLE
+version_tuple: VERSION_TUPLE
+
+__version__ = version = '0.1.dev87+geef7af6.d20241209'
+__version_tuple__ = version_tuple = (0, 1, 'dev87', 'geef7af6.d20241209')

libra_toolbox/neutronics/__init__.py

@@ -0,0 +1 @@
+from .neutron_source import *

libra_toolbox/neutronics/neutron_source.py

@@ -0,0 +1,61 @@
+# building the MIT-VaultLab neutron generator
+# angular and energy distribution
+
+from pathlib import Path
+import numpy as np
+import h5py
+import openmc
+
+
+def neutron_source_diamond(center=(0, 0, 0), reference_uvw=(0, 0, 1)):
+    '''method for building the MIT-VaultLab neutron generator in OpenMC
+    with data tabulated from John Ball and Shon Mackie characterization
+    via diamond detectors
+
+    Parameters
+    ----------
+    center : coordinate position of the source (it is a point source)
+
+    reference_uvw : direction for the polar angle (tuple or list of versors)
+    it is the same for the openmc.PolarAzimuthal class
+    more specifically, polar angle = 0 is the direction of the D accelerator
+    towards the Zr-T target
+    '''
+
+    angles = ["0", "15", "30", "45", "60", "75", "90", "105", "120", "135", "150"]
+
+    filename = 'neutron_source_diamond.h5'
+    filename = str(Path(__file__).parent) / Path(filename)
+
+    with h5py.File(filename, 'r') as mvng_source:
+        # energy values
+        energies = mvng_source['values/table']['Energy (MeV)'] * 1e6
+        # angular bins in [0, pi)
+        pbins = np.cos([np.deg2rad(float(a)) for a in angles] + [np.pi])
+        spectra = [mvng_source['values/table'][col] for col in angles]
+
+    # yield values for strengths
+    yields = np.sum(spectra, axis=-1) * np.diff(pbins)
+    yields /= np.sum(yields)
+
+    # azimuthal values
+    phi = openmc.stats.Uniform(a=0, b=2*np.pi)
+
+    all_sources = []
+    for i, angle in enumerate(pbins[:-1]):
+
+        mu = openmc.stats.Uniform(a=pbins[i+1], b=pbins[i])
+
+        space = openmc.stats.Point(center)
+        angle = openmc.stats.PolarAzimuthal(
+            mu=mu, phi=phi, reference_uvw=reference_uvw)
+        energy = openmc.stats.Tabular(
+            energies, spectra[i], interpolation='linear-linear')
+        strength = yields[i]
+
+        my_source = openmc.Source(
+            space=space, angle=angle, energy=energy, strength=strength, particle='neutron')
+
+        all_sources.append(my_source)
+
+    return all_sources

pyproject.toml

@@ -21,4 +21,7 @@ dependencies = ["numpy", "pint", "scipy", "matplotlib", "sympy", "pandas"]
 tests = ["pytest>=5.4.3", "pytest-cov", "nbconvert", "ipykernel"]
 
 [tool.setuptools_scm]
-write_to = "libra_toolbox/_version.py"
+write_to = "libra_toolbox/_version.py"
+
+[tool.setuptools.package-data]
+mvng_source = ["libra_toolbox/neutronics/neutron_source_diamond.h5"]