moving to eV for nesst

pyproject.toml

@@ -23,7 +23,7 @@ keywords = ["python", "neutron", "fusion", "source", "openmc", "energy", "tokama
 dependencies = [
     "numpy>=1.9",
     "matplotlib>=3.2.2",
-    "NeSST"
+    "NeSST@git+https://github.com/aidancrilly/NeSST",
 ]
 
 [project.urls]

src/openmc_plasma_source/fuel_types.py

@@ -112,8 +112,6 @@ def get_neutron_energy_distribution(
         energy distribution
     """
 
-    ion_temperature_kev = ion_temperature / 1e3  # convert eV to keV
-
     sum_fuel_isotopes = sum(fuel.values())
     if sum_fuel_isotopes > 1.0:
         raise ValueError(
@@ -151,15 +149,15 @@ def get_neutron_energy_distribution(
     if ["D", "T"] == sorted(set(fuel.keys())):
         strength_DT = 1.0
         strength_DD = nst.yield_from_dt_yield_ratio(
-            "dd", strength_DT, ion_temperature_kev, fuel["D"], fuel["T"]
+            "dd", strength_DT, ion_temperature, fuel["D"], fuel["T"]
         )
         strength_TT = nst.yield_from_dt_yield_ratio(
-            "tt", strength_DT, ion_temperature_kev, fuel["D"], fuel["T"]
+            "tt", strength_DT, ion_temperature, fuel["D"], fuel["T"]
         )
 
         total_strength = sum([strength_TT, strength_DD, strength_DT])
 
-        dNdE_TT = strength_TT * nst.dNdE_TT(E_pspec, ion_temperature_kev)
+        dNdE_TT = strength_TT * nst.dNdE_TT(E_pspec, ion_temperature)
         tt_source = openmc.stats.Tabular(E_pspec * 1e6, dNdE_TT)
 
         dd_source = openmc.stats.Normal(mean_value=DDmean, std_dev=DD_std_dev)
@@ -184,6 +182,6 @@ def get_neutron_energy_distribution(
 
     elif ["T"] == sorted(set(fuel.keys())):
         strength_TT = 1.0
-        dNdE_TT = strength_TT * nst.dNdE_TT(E_pspec, ion_temperature_kev)
+        dNdE_TT = strength_TT * nst.dNdE_TT(E_pspec, ion_temperature)
         tt_source = openmc.stats.Tabular(E_pspec * 1e6, dNdE_TT)
         return [tt_source], [strength_TT]

src/openmc_plasma_source/point_source.py

@@ -9,7 +9,7 @@ def fusion_point_source(
     coordinate: Tuple[float, float, float] = (0.0, 0.0, 0.0),
     temperature: float = 20000.0,
     fuel: dict = {"D": 0.5, "T": 0.5},
-):
+) -> list[openmc.IndependentSource]:
     """Creates a list of openmc.IndependentSource objects representing an ICF source.
 
     Resulting ICF (Inertial Confinement Fusion) source will have an energy

src/openmc_plasma_source/tokamak_source.py

@@ -118,8 +118,7 @@ def tokamak_source(
     a = np.random.random(sample_size) * minor_radius
     alpha = np.random.random(sample_size) * 2 * np.pi
 
-    # compute densities, temperatures, neutron source densities and
-    # convert coordinates
+    # compute densities, temperatures
     densities = tokamak_ion_density(
         mode=mode,
         ion_density_centre=ion_density_centre,
@@ -130,6 +129,8 @@ def tokamak_source(
         ion_density_separatrix=ion_density_separatrix,
         r=a,
     )
+
+    # compute temperatures
     temperatures = tokamak_ion_temperature(
         r=a,
         mode=mode,
@@ -142,10 +143,7 @@ def tokamak_source(
         major_radius=major_radius,
     )
 
-    neutron_source_density = tokamak_neutron_source_density(densities, temperatures)
-
-    strengths = neutron_source_density / sum(neutron_source_density)
-
+    # convert coordinates
     RZ = tokamak_convert_a_alpha_to_R_Z(
         a=a,
         alpha=alpha,
@@ -156,6 +154,14 @@ def tokamak_source(
         elongation=elongation,
     )
 
+    #TODO loop through the fuel reactions
+
+    # compute neutron source densities
+    neutron_source_density = tokamak_neutron_source_density(densities, temperatures)
+
+    strengths = neutron_source_density / sum(neutron_source_density)
+
+
     sources = tokamak_make_openmc_sources(
         strengths=strengths,
         angles=angles,
@@ -313,8 +319,12 @@ def tokamak_make_openmc_sources(
     is based on .strengths.
 
     Args:
+        strengths
         angles ((float, float), optional): rotation of the ring source.
-        Defaults to (0, 2*np.pi).
+            Defaults to (0, 2*np.pi).
+        temperatures
+        fuel
+        RZ
 
     Returns:
         list: list of openmc.IndependentSource()
@@ -357,7 +367,7 @@ def tokamak_make_openmc_sources(
     return sources
 
 
-def tokamak_neutron_source_density(ion_density, ion_temperature):
+def tokamak_neutron_source_density(ion_density, ion_temperature, reaction='DT'):
     """Computes the neutron source density given ion density and ion
     temperature.
 
@@ -371,23 +381,26 @@ def tokamak_neutron_source_density(ion_density, ion_temperature):
 
     ion_density = np.asarray(ion_density)
     ion_temperature = np.asarray(ion_temperature)
-
     return ion_density**2 * _DT_xs(ion_temperature)
+    #TODO account for reaction
+    if reaction == 'DT':
+        nst.reac_DT(ion_temperature)
+    if reaction == 'DD'
+        nst.reac_DD(ion_temperature)
+    if reaction == 'TT'
+        nst.reac_TT(ion_temperature)
 
 
+#TODO consider replace with NeSST or getting DD version as well
 def _DT_xs(ion_temperature):
     """Sadler–Van Belle formula
     Ref : https://doi.org/10.1016/j.fusengdes.2012.02.025
-
     Args:
-        ion_temperature (float, ndarray): ion temperature in keV
-
+        ion_temperature (float, ndarray): ion temperature in eV
     Returns:
         float, ndarray: the DT cross section at the given temperature
     """
-
-    ion_temperature = np.asarray(ion_temperature)
-
+    ion_temperature_kev = np.asarray(ion_temperature/1e3)
     c = [
         2.5663271e-18,
         19.983026,
@@ -397,14 +410,12 @@ def _DT_xs(ion_temperature):
         6.6024089e-2,
         8.1215505e-3,
     ]
-
-    U = 1 - ion_temperature * (
-        c[2] + ion_temperature * (c[3] - c[4] * ion_temperature)
-    ) / (1.0 + ion_temperature * (c[5] + c[6] * ion_temperature))
-
+    U = 1 - ion_temperature_kev * (
+        c[2] + ion_temperature_kev * (c[3] - c[4] * ion_temperature_kev)
+    ) / (1.0 + ion_temperature_kev * (c[5] + c[6] * ion_temperature_kev))
     val = (
         c[0]
-        * np.exp(-c[1] * (U / ion_temperature) ** (1 / 3))
-        / (U ** (5 / 6) * ion_temperature ** (2 / 3))
+        * np.exp(-c[1] * (U / ion_temperature_kev) ** (1 / 3))
+        / (U ** (5 / 6) * ion_temperature_kev ** (2 / 3))
     )
     return val