From 54ec402320a8dd8fd044d6688155784704d0a7ca Mon Sep 17 00:00:00 2001
From: magnoxemo <walid.ahammed01@gmail.com>
Date: Thu, 7 Nov 2024 15:34:33 -0600
Subject: [PATCH 01/24] trying to understand what it does!

---
 src/source.cpp | 11 ++++++++++-
 1 file changed, 10 insertions(+), 1 deletion(-)

diff --git a/src/source.cpp b/src/source.cpp
index 15fe8433ba5..ae77cc1e06c 100644
--- a/src/source.cpp
+++ b/src/source.cpp
@@ -616,8 +616,17 @@ SourceSite sample_external_source(uint64_t* seed)
   // Sample from among multiple source distributions
   int i = 0;
   if (model::external_sources.size() > 1) {
-    double xi = prn(seed) * total_strength;
+    double xi = prn(seed) * total_strength; //random strength from the total source strength 
     double c = 0.0;
+
+    /* First, it checks if there are multiple sources in model::external_sources. 
+    If so, it proceeds to sample one. xi is a random value scaled by total_strength.
+    This random value (xi) is used to decide which source distribution to sample from.
+    c accumulates the strengths of each source sequentially. The loop increments c with 
+    the strength of each source until xi < c, indicating that the randomly chosen point
+    falls within the range associated with a specific source. The index i at this point
+    identifies the selected source.*/
+    
     for (; i < model::external_sources.size(); ++i) {
       c += model::external_sources[i]->strength();
       if (xi < c)

From 693302bd2c636de55917de91f2b401ddc87f3fd2 Mon Sep 17 00:00:00 2001
From: Patrick Shriwise <pshriwise@gmail.com>
Date: Fri, 8 Nov 2024 15:18:15 -0600
Subject: [PATCH 02/24] Reset values of lattice offset tables when allocated
 (#3188)

---
 include/openmc/lattice.h | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/include/openmc/lattice.h b/include/openmc/lattice.h
index 429d81ddda8..fb374ab75c3 100644
--- a/include/openmc/lattice.h
+++ b/include/openmc/lattice.h
@@ -71,7 +71,8 @@ class Lattice {
   //! Allocate offset table for distribcell.
   void allocate_offset_table(int n_maps)
   {
-    offsets_.resize(n_maps * universes_.size(), C_NONE);
+    offsets_.resize(n_maps * universes_.size());
+    std::fill(offsets_.begin(), offsets_.end(), C_NONE);
   }
 
   //! Populate the distribcell offset tables.

From 26f08cbae17e80939d40e4c7701402d0271f81b1 Mon Sep 17 00:00:00 2001
From: Jonathan Shimwell <mail@jshimwell.com>
Date: Sat, 9 Nov 2024 00:52:40 +0100
Subject: [PATCH 03/24] allowing varible offsets for polygon.offset (#3120)

Co-authored-by: Paul Romano <paul.k.romano@gmail.com>
---
 openmc/model/surface_composite.py          | 27 ++++++++++++++++++----
 tests/unit_tests/test_surface_composite.py | 11 +++++----
 2 files changed, 30 insertions(+), 8 deletions(-)

diff --git a/openmc/model/surface_composite.py b/openmc/model/surface_composite.py
index df290329647..7a5ab8f1002 100644
--- a/openmc/model/surface_composite.py
+++ b/openmc/model/surface_composite.py
@@ -1,3 +1,4 @@
+from __future__ import annotations
 from abc import ABC, abstractmethod
 from collections.abc import Iterable, Sequence
 from copy import copy
@@ -1316,25 +1317,43 @@ def _decompose_polygon_into_convex_sets(self):
             surfsets.append(surf_ops)
         return surfsets
 
-    def offset(self, distance):
+    def offset(self, distance: float | Sequence[float] | np.ndarray) -> Polygon:
         """Offset this polygon by a set distance
 
         Parameters
         ----------
-        distance : float
+        distance : float or sequence of float or np.ndarray
             The distance to offset the polygon by. Positive is outward
-            (expanding) and negative is inward (shrinking).
+            (expanding) and negative is inward (shrinking). If a float is
+            provided, the same offset is applied to all vertices. If a list or
+            tuple is provided, each vertex gets a different offset. If an
+            iterable or numpy array is provided, each vertex gets a different
+            offset.
 
         Returns
         -------
         offset_polygon : openmc.model.Polygon
         """
+
+        if isinstance(distance, float):
+            distance = np.full(len(self.points), distance)
+        elif isinstance(distance, Sequence):
+            distance = np.array(distance)
+        elif not isinstance(distance, np.ndarray):
+            raise TypeError("Distance must be a float or sequence of float.")
+
+        if len(distance) != len(self.points):
+            raise ValueError(
+                f"Length of distance {len(distance)} array must "
+                f"match number of polygon points {len(self.points)}"
+            )
+
         normals = np.insert(self._normals, 0, self._normals[-1, :], axis=0)
         cos2theta = np.sum(normals[1:, :]*normals[:-1, :], axis=-1, keepdims=True)
         costheta = np.cos(np.arccos(cos2theta) / 2)
         nvec = (normals[1:, :] + normals[:-1, :])
         unit_nvec = nvec / np.linalg.norm(nvec, axis=-1, keepdims=True)
-        disp_vec = distance / costheta * unit_nvec
+        disp_vec = distance[:, np.newaxis] / costheta * unit_nvec
 
         return type(self)(self.points + disp_vec, basis=self.basis)
 
diff --git a/tests/unit_tests/test_surface_composite.py b/tests/unit_tests/test_surface_composite.py
index 963bbe00d19..015ac667e45 100644
--- a/tests/unit_tests/test_surface_composite.py
+++ b/tests/unit_tests/test_surface_composite.py
@@ -347,10 +347,13 @@ def test_polygon():
         assert any([points_in[i] in reg for reg in star_poly.regions])
         assert points_in[i] not in +star_poly
         assert (0, 0, 0) not in -star_poly
-        if basis != 'rz':
-            offset_star = star_poly.offset(.6)
-            assert (0, 0, 0) in -offset_star
-            assert any([(0, 0, 0) in reg for reg in offset_star.regions])
+        if basis != "rz":
+            for offsets in [0.6, np.array([0.6] * 10), [0.6] * 10]:
+                offset_star = star_poly.offset(offsets)
+                assert (0, 0, 0) in -offset_star
+                assert any([(0, 0, 0) in reg for reg in offset_star.regions])
+            with pytest.raises(ValueError):
+                star_poly.offset([0.6, 0.6])
 
     # check invalid Polygon input points
     # duplicate points not just at start and end

From 5fe31912911db43b48c859ed8d8e8a4aa29f7ae1 Mon Sep 17 00:00:00 2001
From: azimG <azimgivron@yahoo.com>
Date: Sat, 9 Nov 2024 17:56:37 +0100
Subject: [PATCH 04/24] Update surface_composite.py (#3189)

---
 openmc/model/surface_composite.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/openmc/model/surface_composite.py b/openmc/model/surface_composite.py
index 7a5ab8f1002..cc570cc56e5 100644
--- a/openmc/model/surface_composite.py
+++ b/openmc/model/surface_composite.py
@@ -27,7 +27,7 @@ def translate(self, vector, inplace=False):
         return surf
 
     def rotate(self, rotation, pivot=(0., 0., 0.), order='xyz', inplace=False):
-        surf = copy(self)
+        surf = self if inplace else copy(self)
         for name in self._surface_names:
             s = getattr(surf, name)
             setattr(surf, name, s.rotate(rotation, pivot, order, inplace))

From c87f749fc421e31d4692027c980a06cc66432649 Mon Sep 17 00:00:00 2001
From: Nicolas Linden <nicolas@nlinden.com>
Date: Sat, 9 Nov 2024 20:33:14 +0100
Subject: [PATCH 05/24] add export_model_xml arguments to Model.plot_geometry
 and Model.calculate_volumes (#3190)

Co-authored-by: Nicolas Linden <n.linde@naarea.fr>
Co-authored-by: Paul Romano <paul.k.romano@gmail.com>
---
 openmc/model/model.py | 41 +++++++++++++++++++++++++++++++++--------
 1 file changed, 33 insertions(+), 8 deletions(-)

diff --git a/openmc/model/model.py b/openmc/model/model.py
index 78f743a03ed..0ea6db8db88 100644
--- a/openmc/model/model.py
+++ b/openmc/model/model.py
@@ -719,7 +719,8 @@ def run(self, particles=None, threads=None, geometry_debug=False,
 
     def calculate_volumes(self, threads=None, output=True, cwd='.',
                           openmc_exec='openmc', mpi_args=None,
-                          apply_volumes=True):
+                          apply_volumes=True, export_model_xml=True,
+                          **export_kwargs):
         """Runs an OpenMC stochastic volume calculation and, if requested,
         applies volumes to the model
 
@@ -748,6 +749,13 @@ def calculate_volumes(self, threads=None, output=True, cwd='.',
         apply_volumes : bool, optional
             Whether apply the volume calculation results from this calculation
             to the model. Defaults to applying the volumes.
+        export_model_xml : bool, optional
+            Exports a single model.xml file rather than separate files. Defaults
+            to True.
+        **export_kwargs
+            Keyword arguments passed to either :meth:`Model.export_to_model_xml`
+            or :meth:`Model.export_to_xml`.
+
         """
 
         if len(self.settings.volume_calculations) == 0:
@@ -769,10 +777,15 @@ def calculate_volumes(self, threads=None, output=True, cwd='.',
                 openmc.lib.calculate_volumes(output)
 
             else:
-                self.export_to_xml()
-                openmc.calculate_volumes(threads=threads, output=output,
-                                         openmc_exec=openmc_exec,
-                                         mpi_args=mpi_args)
+                if export_model_xml:
+                    self.export_to_model_xml(**export_kwargs)
+                else:
+                    self.export_to_xml(**export_kwargs)
+                path_input = export_kwargs.get("path", None)
+                openmc.calculate_volumes(
+                    threads=threads, output=output, openmc_exec=openmc_exec,
+                    mpi_args=mpi_args, path_input=path_input
+                )
 
             # Now we apply the volumes
             if apply_volumes:
@@ -909,7 +922,8 @@ def sample_external_source(
                     n_samples=n_samples, prn_seed=prn_seed
                 )
 
-    def plot_geometry(self, output=True, cwd='.', openmc_exec='openmc'):
+    def plot_geometry(self, output=True, cwd='.', openmc_exec='openmc',
+                      export_model_xml=True, **export_kwargs):
         """Creates plot images as specified by the Model.plots attribute
 
         .. versionadded:: 0.13.0
@@ -924,6 +938,12 @@ def plot_geometry(self, output=True, cwd='.', openmc_exec='openmc'):
         openmc_exec : str, optional
             Path to OpenMC executable. Defaults to 'openmc'.
             This only applies to the case when not using the C API.
+        export_model_xml : bool, optional
+            Exports a single model.xml file rather than separate files. Defaults
+            to True.
+        **export_kwargs
+            Keyword arguments passed to either :meth:`Model.export_to_model_xml`
+            or :meth:`Model.export_to_xml`.
 
         """
 
@@ -937,8 +957,13 @@ def plot_geometry(self, output=True, cwd='.', openmc_exec='openmc'):
                 # Compute the volumes
                 openmc.lib.plot_geometry(output)
             else:
-                self.export_to_xml()
-                openmc.plot_geometry(output=output, openmc_exec=openmc_exec)
+                if export_model_xml:
+                    self.export_to_model_xml(**export_kwargs)
+                else:
+                    self.export_to_xml(**export_kwargs)
+                path_input = export_kwargs.get("path", None)
+                openmc.plot_geometry(output=output, openmc_exec=openmc_exec,
+                                     path_input=path_input)
 
     def _change_py_lib_attribs(self, names_or_ids, value, obj_type,
                                attrib_name, density_units='atom/b-cm'):

From 454df35979d9efa21800b9ed1542859e79b40454 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Mon, 11 Nov 2024 16:03:31 -0600
Subject: [PATCH 06/24] Fixes in MicroXS.from_multigroup_flux (#3192)

---
 openmc/deplete/independent_operator.py                |  5 +++--
 openmc/deplete/microxs.py                             | 10 ++++++----
 tests/unit_tests/test_deplete_decay.py                |  4 ++--
 tests/unit_tests/test_deplete_independent_operator.py |  6 +++---
 tests/unit_tests/test_deplete_microxs.py              |  8 +++++---
 5 files changed, 19 insertions(+), 14 deletions(-)

diff --git a/openmc/deplete/independent_operator.py b/openmc/deplete/independent_operator.py
index cc1a05bd99f..226682ffa53 100644
--- a/openmc/deplete/independent_operator.py
+++ b/openmc/deplete/independent_operator.py
@@ -44,7 +44,7 @@ class IndependentOperator(OpenMCOperator):
 
     Parameters
     ----------
-    materials : openmc.Materials
+    materials : iterable of openmc.Material
         Materials to deplete.
     fluxes : list of numpy.ndarray
         Flux in each group in [n-cm/src] for each domain
@@ -127,8 +127,9 @@ def __init__(self,
                  reduce_chain_level=None,
                  fission_yield_opts=None):
         # Validate micro-xs parameters
-        check_type('materials', materials, openmc.Materials)
+        check_type('materials', materials, Iterable, openmc.Material)
         check_type('micros', micros, Iterable, MicroXS)
+        materials = openmc.Materials(materials)
 
         if not (len(fluxes) == len(micros) == len(materials)):
             msg = (f'The length of fluxes ({len(fluxes)}) should be equal to '
diff --git a/openmc/deplete/microxs.py b/openmc/deplete/microxs.py
index 5be9875f304..bbe9b2a4322 100644
--- a/openmc/deplete/microxs.py
+++ b/openmc/deplete/microxs.py
@@ -198,6 +198,8 @@ class MicroXS:
     """
     def __init__(self, data: np.ndarray, nuclides: list[str], reactions: list[str]):
         # Validate inputs
+        if len(data.shape) != 3:
+            raise ValueError('Data array must be 3D.')
         if data.shape[:2] != (len(nuclides), len(reactions)):
             raise ValueError(
                 f'Nuclides list of length {len(nuclides)} and '
@@ -291,11 +293,11 @@ def from_multigroup_flux(
         mts = [REACTION_MT[name] for name in reactions]
 
         # Normalize multigroup flux
-        multigroup_flux = np.asarray(multigroup_flux)
+        multigroup_flux = np.array(multigroup_flux)
         multigroup_flux /= multigroup_flux.sum()
 
-        # Create 2D array for microscopic cross sections
-        microxs_arr = np.zeros((len(nuclides), len(mts)))
+        # Create 3D array for microscopic cross sections
+        microxs_arr = np.zeros((len(nuclides), len(mts), 1))
 
         # Create a material with all nuclides
         mat_all_nucs = openmc.Material()
@@ -327,7 +329,7 @@ def from_multigroup_flux(
                         xs = lib_nuc.collapse_rate(
                             mt, temperature, energies, multigroup_flux
                         )
-                        microxs_arr[nuc_index, mt_index] = xs
+                        microxs_arr[nuc_index, mt_index, 0] = xs
 
         return cls(microxs_arr, nuclides, reactions)
 
diff --git a/tests/unit_tests/test_deplete_decay.py b/tests/unit_tests/test_deplete_decay.py
index cebbc16fa79..6e7b0b101ec 100644
--- a/tests/unit_tests/test_deplete_decay.py
+++ b/tests/unit_tests/test_deplete_decay.py
@@ -20,7 +20,7 @@ def test_deplete_decay_products(run_in_tmpdir):
         """)
 
     # Create MicroXS object with no cross sections
-    micro_xs = openmc.deplete.MicroXS(np.empty((0, 0)), [], [])
+    micro_xs = openmc.deplete.MicroXS(np.empty((0, 0, 0)), [], [])
 
     # Create depletion operator with no reactions
     op = openmc.deplete.IndependentOperator.from_nuclides(
@@ -59,7 +59,7 @@ def test_deplete_decay_step_fissionable(run_in_tmpdir):
     """
 
     # Set up a pure decay operator
-    micro_xs = openmc.deplete.MicroXS(np.empty((0, 0)), [], [])
+    micro_xs = openmc.deplete.MicroXS(np.empty((0, 0, 0)), [], [])
     mat = openmc.Material()
     mat.name = 'I do not decay.'
     mat.add_nuclide('U238', 1.0, 'ao')
diff --git a/tests/unit_tests/test_deplete_independent_operator.py b/tests/unit_tests/test_deplete_independent_operator.py
index 9129cf0642f..c765d065009 100644
--- a/tests/unit_tests/test_deplete_independent_operator.py
+++ b/tests/unit_tests/test_deplete_independent_operator.py
@@ -6,7 +6,7 @@
 
 import pytest
 
-from openmc import Material, Materials
+from openmc import Material
 from openmc.deplete import IndependentOperator, MicroXS
 
 CHAIN_PATH = Path(__file__).parents[1] / "chain_simple.xml"
@@ -34,7 +34,7 @@ def test_operator_init():
     fuel.set_density("g/cc", 10.4)
     fuel.depletable = True
     fuel.volume = 1
-    materials = Materials([fuel])
+    materials = [fuel]
     fluxes = [1.0]
     micros = [micro_xs]
     IndependentOperator(materials, fluxes, micros, CHAIN_PATH)
@@ -47,7 +47,7 @@ def test_error_handling():
     fuel.set_density("g/cc", 1)
     fuel.depletable = True
     fuel.volume = 1
-    materials = Materials([fuel])
+    materials = [fuel]
     fluxes = [1.0, 2.0]
     micros = [micro_xs]
     with pytest.raises(ValueError, match=r"The length of fluxes \(2\)"):
diff --git a/tests/unit_tests/test_deplete_microxs.py b/tests/unit_tests/test_deplete_microxs.py
index ad54026f014..073b3f162d1 100644
--- a/tests/unit_tests/test_deplete_microxs.py
+++ b/tests/unit_tests/test_deplete_microxs.py
@@ -46,9 +46,7 @@ def test_from_array():
     data.shape = (12, 2, 1)
 
     MicroXS(data, nuclides, reactions)
-    with pytest.raises(ValueError, match=r'Nuclides list of length \d* and '
-                       r'reactions array of length \d* do not '
-                       r'match dimensions of data array of shape \(\d*\, \d*\)'):
+    with pytest.raises(ValueError, match='Data array must be 3D'):
         MicroXS(data[:, 0], nuclides, reactions)
 
 
@@ -96,9 +94,13 @@ def test_multigroup_flux_same():
     energies = [0., 6.25e-1, 5.53e3, 8.21e5, 2.e7]
     flux_per_ev = [0.3, 0.3, 1.0, 1.0]
     flux = flux_per_ev * np.diff(energies)
+    flux_sum = flux.sum()
     microxs_4g = MicroXS.from_multigroup_flux(
         energies=energies, multigroup_flux=flux, chain_file=chain_file)
 
+    # from_multigroup_flux should not modify the flux
+    assert flux.sum() == flux_sum
+
     # Generate micro XS based on 2-group flux, where the boundaries line up with
     # the 4 group flux and have the same flux per eV across the full energy
     # range

From 3ad81c5e59ae02e005ec6acedb23358c52960c0d Mon Sep 17 00:00:00 2001
From: Patrick Shriwise <pshriwise@gmail.com>
Date: Tue, 12 Nov 2024 14:58:15 -0600
Subject: [PATCH 07/24] Fix documentation typo in `boundary_type` (#3196)

---
 openmc/surface.py | 66 +++++++++++++++++++++++------------------------
 1 file changed, 33 insertions(+), 33 deletions(-)

diff --git a/openmc/surface.py b/openmc/surface.py
index 537823ba16a..c95025f949b 100644
--- a/openmc/surface.py
+++ b/openmc/surface.py
@@ -119,7 +119,7 @@ class Surface(IDManagerMixin, ABC):
     surface_id : int, optional
         Unique identifier for the surface. If not specified, an identifier will
         automatically be assigned.
-    boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic', 'white'}, optional
+    boundary_type : {'transmission', 'vacuum', 'reflective', 'periodic', 'white'}, optional
         Boundary condition that defines the behavior for particles hitting the
         surface. Defaults to transmissive boundary condition where particles
         freely pass through the surface. Note that periodic boundary conditions
@@ -135,7 +135,7 @@ class Surface(IDManagerMixin, ABC):
 
     Attributes
     ----------
-    boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic', 'white'}
+    boundary_type : {'transmission', 'vacuum', 'reflective', 'periodic', 'white'}
         Boundary condition that defines the behavior for particles hitting the
         surface.
     albedo : float
@@ -683,7 +683,7 @@ class Plane(PlaneMixin, Surface):
         The 'C' parameter for the plane. Defaults to 0.
     d : float, optional
         The 'D' parameter for the plane. Defaults to 0.
-    boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic', 'white'}, optional
+    boundary_type : {'transmission', 'vacuum', 'reflective', 'periodic', 'white'}, optional
         Boundary condition that defines the behavior for particles hitting the
         surface. Defaults to transmissive boundary condition where particles
         freely pass through the surface.
@@ -707,7 +707,7 @@ class Plane(PlaneMixin, Surface):
         The 'C' parameter for the plane
     d : float
         The 'D' parameter for the plane
-    boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic', 'white'}
+    boundary_type : {'transmission', 'vacuum', 'reflective', 'periodic', 'white'}
         Boundary condition that defines the behavior for particles hitting the
         surface.
     albedo : float
@@ -818,7 +818,7 @@ class XPlane(PlaneMixin, Surface):
     ----------
     x0 : float, optional
         Location of the plane in [cm]. Defaults to 0.
-    boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic', 'white'}, optional
+    boundary_type : {'transmission', 'vacuum', 'reflective', 'periodic', 'white'}, optional
         Boundary condition that defines the behavior for particles hitting the
         surface. Defaults to transmissive boundary condition where particles
         freely pass through the surface. Only axis-aligned periodicity is
@@ -837,7 +837,7 @@ class XPlane(PlaneMixin, Surface):
     ----------
     x0 : float
         Location of the plane in [cm]
-    boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic', 'white'}
+    boundary_type : {'transmission', 'vacuum', 'reflective', 'periodic', 'white'}
         Boundary condition that defines the behavior for particles hitting the
         surface.
     albedo : float
@@ -883,7 +883,7 @@ class YPlane(PlaneMixin, Surface):
     ----------
     y0 : float, optional
         Location of the plane in [cm]
-    boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic', 'white'}, optional
+    boundary_type : {'transmission', 'vacuum', 'reflective', 'periodic', 'white'}, optional
         Boundary condition that defines the behavior for particles hitting the
         surface. Defaults to transmissive boundary condition where particles
         freely pass through the surface. Only axis-aligned periodicity is
@@ -902,7 +902,7 @@ class YPlane(PlaneMixin, Surface):
     ----------
     y0 : float
         Location of the plane in [cm]
-    boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic', 'white'}
+    boundary_type : {'transmission', 'vacuum', 'reflective', 'periodic', 'white'}
         Boundary condition that defines the behavior for particles hitting the
         surface.
     albedo : float
@@ -948,7 +948,7 @@ class ZPlane(PlaneMixin, Surface):
     ----------
     z0 : float, optional
         Location of the plane in [cm]. Defaults to 0.
-    boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic', 'white'}, optional
+    boundary_type : {'transmission', 'vacuum', 'reflective', 'periodic', 'white'}, optional
         Boundary condition that defines the behavior for particles hitting the
         surface. Defaults to transmissive boundary condition where particles
         freely pass through the surface. Only axis-aligned periodicity is
@@ -967,7 +967,7 @@ class ZPlane(PlaneMixin, Surface):
     ----------
     z0 : float
         Location of the plane in [cm]
-    boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic', 'white'}
+    boundary_type : {'transmission', 'vacuum', 'reflective', 'periodic', 'white'}
         Boundary condition that defines the behavior for particles hitting the
         surface.
     albedo : float
@@ -1197,7 +1197,7 @@ class Cylinder(QuadricMixin, Surface):
     dz : float, optional
         z-component of the vector representing the axis of the cylinder.
         Defaults to 1.
-    boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}, optional
+    boundary_type : {'transmission', 'vacuum', 'reflective', 'white'}, optional
         Boundary condition that defines the behavior for particles hitting the
         surface. Defaults to transmissive boundary condition where particles
         freely pass through the surface.
@@ -1228,7 +1228,7 @@ class Cylinder(QuadricMixin, Surface):
         y-component of the vector representing the axis of the cylinder
     dz : float
         z-component of the vector representing the axis of the cylinder
-    boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}
+    boundary_type : {'transmission', 'vacuum', 'reflective', 'white'}
         Boundary condition that defines the behavior for particles hitting the
         surface.
     albedo : float
@@ -1372,7 +1372,7 @@ class XCylinder(QuadricMixin, Surface):
         z-coordinate for the origin of the Cylinder in [cm]. Defaults to 0
     r : float, optional
         Radius of the cylinder in [cm]. Defaults to 1.
-    boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}, optional
+    boundary_type : {'transmission', 'vacuum', 'reflective', 'white'}, optional
         Boundary condition that defines the behavior for particles hitting the
         surface. Defaults to transmissive boundary condition where particles
         freely pass through the surface.
@@ -1395,7 +1395,7 @@ class XCylinder(QuadricMixin, Surface):
         z-coordinate for the origin of the Cylinder in [cm]
     r : float
         Radius of the cylinder in [cm]
-    boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}
+    boundary_type : {'transmission', 'vacuum', 'reflective', 'white'}
         Boundary condition that defines the behavior for particles hitting the
         surface.
     albedo : float
@@ -1470,7 +1470,7 @@ class YCylinder(QuadricMixin, Surface):
         z-coordinate for the origin of the Cylinder in [cm]. Defaults to 0
     r : float, optional
         Radius of the cylinder in [cm]. Defaults to 1.
-    boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}, optional
+    boundary_type : {'transmission', 'vacuum', 'reflective', 'white'}, optional
         Boundary condition that defines the behavior for particles hitting the
         surface. Defaults to transmissive boundary condition where particles
         freely pass through the surface.
@@ -1493,7 +1493,7 @@ class YCylinder(QuadricMixin, Surface):
         z-coordinate for the origin of the Cylinder in [cm]
     r : float
         Radius of the cylinder in [cm]
-    boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}
+    boundary_type : {'transmission', 'vacuum', 'reflective', 'white'}
         Boundary condition that defines the behavior for particles hitting the
         surface.
     albedo : float
@@ -1568,7 +1568,7 @@ class ZCylinder(QuadricMixin, Surface):
         y-coordinate for the origin of the Cylinder in [cm]. Defaults to 0
     r : float, optional
         Radius of the cylinder in [cm]. Defaults to 1.
-    boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}, optional
+    boundary_type : {'transmission', 'vacuum', 'reflective', 'white'}, optional
         Boundary condition that defines the behavior for particles hitting the
         surface. Defaults to transmissive boundary condition where particles
         freely pass through the surface.
@@ -1591,7 +1591,7 @@ class ZCylinder(QuadricMixin, Surface):
         y-coordinate for the origin of the Cylinder in [cm]
     r : float
         Radius of the cylinder in [cm]
-    boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}
+    boundary_type : {'transmission', 'vacuum', 'reflective', 'white'}
         Boundary condition that defines the behavior for particles hitting the
         surface.
     albedo : float
@@ -1667,7 +1667,7 @@ class Sphere(QuadricMixin, Surface):
         z-coordinate of the center of the sphere in [cm]. Defaults to 0.
     r : float, optional
         Radius of the sphere in [cm]. Defaults to 1.
-    boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}, optional
+    boundary_type : {'transmission', 'vacuum', 'reflective', 'white'}, optional
         Boundary condition that defines the behavior for particles hitting the
         surface. Defaults to transmissive boundary condition where particles
         freely pass through the surface.
@@ -1691,7 +1691,7 @@ class Sphere(QuadricMixin, Surface):
         z-coordinate of the center of the sphere in [cm]
     r : float
         Radius of the sphere in [cm]
-    boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}
+    boundary_type : {'transmission', 'vacuum', 'reflective', 'white'}
         Boundary condition that defines the behavior for particles hitting the
         surface.
     albedo : float
@@ -1784,7 +1784,7 @@ class Cone(QuadricMixin, Surface):
     surface_id : int, optional
         Unique identifier for the surface. If not specified, an identifier will
         automatically be assigned.
-    boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}, optional
+    boundary_type : {'transmission', 'vacuum', 'reflective', 'white'}, optional
         Boundary condition that defines the behavior for particles hitting the
         surface. Defaults to transmissive boundary condition where particles
         freely pass through the surface.
@@ -1812,7 +1812,7 @@ class Cone(QuadricMixin, Surface):
         y-component of the vector representing the axis of the cone.
     dz : float
         z-component of the vector representing the axis of the cone.
-    boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}
+    boundary_type : {'transmission', 'vacuum', 'reflective', 'white'}
         Boundary condition that defines the behavior for particles hitting the
         surface.
     albedo : float
@@ -1930,7 +1930,7 @@ class XCone(QuadricMixin, Surface):
         Parameter related to the aperture [:math:`\\rm cm^2`].
         It can be interpreted as the increase in the radius squared per cm along
         the cone's axis of revolution.
-    boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}, optional
+    boundary_type : {'transmission', 'vacuum', 'reflective', 'white'}, optional
         Boundary condition that defines the behavior for particles hitting the
         surface. Defaults to transmissive boundary condition where particles
         freely pass through the surface.
@@ -1954,7 +1954,7 @@ class XCone(QuadricMixin, Surface):
         z-coordinate of the apex in [cm]
     r2 : float
         Parameter related to the aperature
-    boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}
+    boundary_type : {'transmission', 'vacuum', 'reflective', 'white'}
         Boundary condition that defines the behavior for particles hitting the
         surface.
     albedo : float
@@ -2031,7 +2031,7 @@ class YCone(QuadricMixin, Surface):
         Parameter related to the aperture [:math:`\\rm cm^2`].
         It can be interpreted as the increase in the radius squared per cm along
         the cone's axis of revolution.
-    boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}, optional
+    boundary_type : {'transmission', 'vacuum', 'reflective', 'white'}, optional
         Boundary condition that defines the behavior for particles hitting the
         surface. Defaults to transmissive boundary condition where particles
         freely pass through the surface.
@@ -2055,7 +2055,7 @@ class YCone(QuadricMixin, Surface):
         z-coordinate of the apex in [cm]
     r2 : float
         Parameter related to the aperature
-    boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}
+    boundary_type : {'transmission', 'vacuum', 'reflective', 'white'}
         Boundary condition that defines the behavior for particles hitting the
         surface.
     albedo : float
@@ -2132,7 +2132,7 @@ class ZCone(QuadricMixin, Surface):
         Parameter related to the aperature [cm^2].
         This is the square of the radius of the cone 1 cm from.
         This can also be treated as the square of the slope of the cone relative to its axis.
-    boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}, optional
+    boundary_type : {'transmission', 'vacuum', 'reflective', 'white'}, optional
         Boundary condition that defines the behavior for particles hitting the
         surface. Defaults to transmissive boundary condition where particles
         freely pass through the surface.
@@ -2156,7 +2156,7 @@ class ZCone(QuadricMixin, Surface):
         z-coordinate of the apex in [cm]
     r2 : float
         Parameter related to the aperature
-    boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}
+    boundary_type : {'transmission', 'vacuum', 'reflective', 'white'}
         Boundary condition that defines the behavior for particles hitting the
         surface.
     albedo : float
@@ -2221,7 +2221,7 @@ class Quadric(QuadricMixin, Surface):
     ----------
     a, b, c, d, e, f, g, h, j, k : float, optional
         coefficients for the surface. All default to 0.
-    boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}, optional
+    boundary_type : {'transmission', 'vacuum', 'reflective', 'white'}, optional
         Boundary condition that defines the behavior for particles hitting the
         surface. Defaults to transmissive boundary condition where particles
         freely pass through the surface.
@@ -2239,7 +2239,7 @@ class Quadric(QuadricMixin, Surface):
     ----------
     a, b, c, d, e, f, g, h, j, k : float
         coefficients for the surface
-    boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}
+    boundary_type : {'transmission', 'vacuum', 'reflective', 'white'}
         Boundary condition that defines the behavior for particles hitting the
         surface.
     albedo : float
@@ -2391,7 +2391,7 @@ class XTorus(TorusMixin, Surface):
         Minor radius of the torus in [cm] (parallel to axis of revolution)
     c : float
         Minor radius of the torus in [cm] (perpendicular to axis of revolution)
-    boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}
+    boundary_type : {'transmission', 'vacuum', 'reflective', 'white'}
         Boundary condition that defines the behavior for particles hitting the
         surface.
     albedo : float
@@ -2466,7 +2466,7 @@ class YTorus(TorusMixin, Surface):
         Minor radius of the torus in [cm] (parallel to axis of revolution)
     c : float
         Minor radius of the torus (perpendicular to axis of revolution)
-    boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}
+    boundary_type : {'transmission', 'vacuum', 'reflective', 'white'}
         Boundary condition that defines the behavior for particles hitting the
         surface.
     albedo : float
@@ -2541,7 +2541,7 @@ class ZTorus(TorusMixin, Surface):
         Minor radius of the torus in [cm] (parallel to axis of revolution)
     c : float
         Minor radius of the torus in [cm] (perpendicular to axis of revolution)
-    boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}
+    boundary_type : {'transmission', 'vacuum', 'reflective', 'white'}
         Boundary condition that defines the behavior for particles hitting the
         surface.
     albedo : float

From 70347ff6a90549b558e4f7c0f0c00d990bfbdc9f Mon Sep 17 00:00:00 2001
From: Seda Yilmaz <133794379+nsedayilmaz@users.noreply.github.com>
Date: Wed, 13 Nov 2024 02:11:59 +0300
Subject: [PATCH 08/24] Add a vessel composite surface with ellipsoids on top
 and bottom. (#3168)

Co-authored-by: Patrick Shriwise <pshriwise@gmail.com>
Co-authored-by: Paul Romano <paul.k.romano@gmail.com>
---
 docs/source/pythonapi/model.rst            |   1 +
 openmc/model/surface_composite.py          | 105 ++++++++++++++++++++-
 tests/unit_tests/test_surface_composite.py |  49 ++++++++++
 3 files changed, 150 insertions(+), 5 deletions(-)

diff --git a/docs/source/pythonapi/model.rst b/docs/source/pythonapi/model.rst
index e7d6d320f1e..3034826bddd 100644
--- a/docs/source/pythonapi/model.rst
+++ b/docs/source/pythonapi/model.rst
@@ -32,6 +32,7 @@ Composite Surfaces
    openmc.model.RectangularParallelepiped
    openmc.model.RectangularPrism
    openmc.model.RightCircularCylinder
+   openmc.model.Vessel
    openmc.model.XConeOneSided
    openmc.model.YConeOneSided
    openmc.model.ZConeOneSided
diff --git a/openmc/model/surface_composite.py b/openmc/model/surface_composite.py
index cc570cc56e5..1f25b9ca636 100644
--- a/openmc/model/surface_composite.py
+++ b/openmc/model/surface_composite.py
@@ -41,9 +41,11 @@ def boundary_type(self):
     def boundary_type(self, boundary_type):
         # Set boundary type on underlying surfaces, but not for ambiguity plane
         # on one-sided cones
+        classes = (XConeOneSided, YConeOneSided, ZConeOneSided, Vessel)
         for name in self._surface_names:
-            if name != 'plane':
-                getattr(self, name).boundary_type = boundary_type
+            if isinstance(self, classes) and name.startswith('plane'):
+                continue
+            getattr(self, name).boundary_type = boundary_type
 
     def __repr__(self):
         return f"<{type(self).__name__} at 0x{id(self):x}>"
@@ -90,8 +92,8 @@ class CylinderSector(CompositeSurface):
         counterclockwise direction with respect to the first basis axis
         (+y, +z, or +x). Must be greater than :attr:`theta1`.
     center : iterable of float
-       Coordinate for central axes of cylinders in the (y, z), (x, z), or (x, y)
-       basis. Defaults to (0,0).
+        Coordinate for central axes of cylinders in the (y, z), (x, z), or (x, y)
+        basis. Defaults to (0,0).
     axis : {'x', 'y', 'z'}
         Central axis of the cylinders defining the inner and outer surfaces of
         the sector. Defaults to 'z'.
@@ -119,7 +121,7 @@ def __init__(self,
                  r2,
                  theta1,
                  theta2,
-                 center=(0.,0.),
+                 center=(0., 0.),
                  axis='z',
                  **kwargs):
 
@@ -1856,3 +1858,96 @@ def __init__(self, center_base: Sequence[float], axis: Sequence[float],
 
     def __neg__(self) -> openmc.Region:
         return +self.plane_bottom & -self.plane_top & -self.cone
+
+
+class Vessel(CompositeSurface):
+    """Vessel composed of cylinder with semi-ellipsoid top and bottom.
+
+    This composite surface is represented by a finite cylinder with ellipsoidal
+    top and bottom surfaces. This surface is equivalent to the 'vesesl' surface
+    in Serpent.
+
+    .. versionadded:: 0.15.1
+
+    Parameters
+    ----------
+    r : float
+        Radius of vessel.
+    p1 : float
+        Minimum coordinate for cylindrical part of vessel.
+    p2 : float
+        Maximum coordinate for cylindrical part of vessel.
+    h1 : float
+        Height of bottom ellipsoidal part of vessel.
+    h2 : float
+        Height of top ellipsoidal part of vessel.
+    center : 2-tuple of float
+        Coordinate for central axis of the cylinder in the (y, z), (x, z), or
+        (x, y) basis. Defaults to (0,0).
+    axis : {'x', 'y', 'z'}
+        Central axis of the cylinder.
+
+    """
+
+    _surface_names = ('cyl', 'plane_bottom', 'plane_top', 'bottom', 'top')
+
+    def __init__(self, r: float, p1: float, p2: float, h1: float, h2: float,
+                 center: Sequence[float] = (0., 0.), axis: str = 'z', **kwargs):
+        if p1 >= p2:
+            raise ValueError('p1 must be less than p2')
+        check_value('axis', axis, {'x', 'y', 'z'})
+
+        c1, c2 = center
+        cyl_class = getattr(openmc, f'{axis.upper()}Cylinder')
+        plane_class = getattr(openmc, f'{axis.upper()}Plane')
+        self.cyl = cyl_class(c1, c2, r, **kwargs)
+        self.plane_bottom = plane_class(p1)
+        self.plane_top = plane_class(p2)
+
+        # General equation for an ellipsoid:
+        #   (x-x₀)²/r² + (y-y₀)²/r² + (z-z₀)²/h² = 1
+        #   (x-x₀)² + (y-y₀)² + (z-z₀)²s² = r²
+        # Let s = r/h:
+        #   (x² - 2x₀x + x₀²) + (y² - 2y₀y + y₀²) + (z² - 2z₀z + z₀²)s² = r²
+        #   x² + y² + s²z² - 2x₀x - 2y₀y - 2s²z₀z + (x₀² + y₀² + z₀²s² - r²) = 0
+
+        sb = (r/h1)
+        st = (r/h2)
+        kwargs['a'] = kwargs['b'] = kwargs['c'] = 1.0
+        kwargs_bottom = kwargs
+        kwargs_top = kwargs.copy()
+
+        sb2 = sb*sb
+        st2 = st*st
+        kwargs_bottom['k'] = c1*c1 + c2*c2 + p1*p1*sb2 - r*r
+        kwargs_top['k'] = c1*c1 + c2*c2 + p2*p2*st2 - r*r
+
+        if axis == 'x':
+            kwargs_bottom['a'] *= sb2
+            kwargs_top['a'] *= st2
+            kwargs_bottom['g'] = -2*p1*sb2
+            kwargs_top['g'] = -2*p2*st2
+            kwargs_top['h'] = kwargs_bottom['h'] = -2*c1
+            kwargs_top['j'] = kwargs_bottom['j'] = -2*c2
+        elif axis == 'y':
+            kwargs_bottom['b'] *= sb2
+            kwargs_top['b'] *= st2
+            kwargs_top['g'] = kwargs_bottom['g'] = -2*c1
+            kwargs_bottom['h'] = -2*p1*sb2
+            kwargs_top['h'] = -2*p2*st2
+            kwargs_top['j'] = kwargs_bottom['j'] = -2*c2
+        elif axis == 'z':
+            kwargs_bottom['c'] *= sb2
+            kwargs_top['c'] *= st2
+            kwargs_top['g'] = kwargs_bottom['g'] = -2*c1
+            kwargs_top['h'] = kwargs_bottom['h'] = -2*c2
+            kwargs_bottom['j'] = -2*p1*sb2
+            kwargs_top['j'] = -2*p2*st2
+
+        self.bottom = openmc.Quadric(**kwargs_bottom)
+        self.top = openmc.Quadric(**kwargs_top)
+
+    def __neg__(self):
+        return ((-self.cyl & +self.plane_bottom & -self.plane_top) |
+                (-self.bottom & -self.plane_bottom) |
+                (-self.top & +self.plane_top))
diff --git a/tests/unit_tests/test_surface_composite.py b/tests/unit_tests/test_surface_composite.py
index 015ac667e45..da7ffbcc470 100644
--- a/tests/unit_tests/test_surface_composite.py
+++ b/tests/unit_tests/test_surface_composite.py
@@ -599,3 +599,52 @@ def test_conical_frustum():
     # Denegenerate case with r1 = r2
     s = openmc.model.ConicalFrustum(center_base, axis, r1, r1)
     assert (1., 1., -0.01) in -s
+
+
+def test_vessel():
+    center = (3.0, 2.0)
+    r = 1.0
+    p1, p2 = -5.0, 5.0
+    h1 = h2 = 1.0
+    s = openmc.model.Vessel(r, p1, p2, h1, h2, center)
+    assert isinstance(s.cyl, openmc.Cylinder)
+    assert isinstance(s.plane_bottom, openmc.Plane)
+    assert isinstance(s.plane_top, openmc.Plane)
+    assert isinstance(s.bottom, openmc.Quadric)
+    assert isinstance(s.top, openmc.Quadric)
+
+    # Make sure boundary condition propagates (but not for planes)
+    s.boundary_type = 'reflective'
+    assert s.boundary_type == 'reflective'
+    assert s.cyl.boundary_type == 'reflective'
+    assert s.bottom.boundary_type == 'reflective'
+    assert s.top.boundary_type == 'reflective'
+    assert s.plane_bottom.boundary_type == 'transmission'
+    assert s.plane_top.boundary_type == 'transmission'
+
+    # Check bounding box
+    ll, ur = (+s).bounding_box
+    assert np.all(np.isinf(ll))
+    assert np.all(np.isinf(ur))
+    ll, ur = (-s).bounding_box
+    assert np.all(np.isinf(ll))
+    assert np.all(np.isinf(ur))
+
+    # __contains__ on associated half-spaces
+    assert (3., 2., 0.) in -s
+    assert (3., 2., -5.0) in -s
+    assert (3., 2., 5.0) in -s
+    assert (3., 2., -5.9) in -s
+    assert (3., 2., 5.9) in -s
+    assert (3., 2., -6.1) not in -s
+    assert (3., 2., 6.1) not in -s
+    assert (4.5, 2., 0.) in +s
+    assert (3., 3.2, 0.) in +s
+    assert (3., 2., 7.) in +s
+
+    # translate method
+    s_t = s.translate((0., 0., 1.))
+    assert (3., 2., 6.1) in -s_t
+
+    # Make sure repr works
+    repr(s)

From 030b3d80bab11a4a66b016c774df28c079fd8d67 Mon Sep 17 00:00:00 2001
From: Paul Wilson <paul.wilson@wisc.edu>
Date: Wed, 13 Nov 2024 15:39:41 -0600
Subject: [PATCH 09/24] Add PointCloud spatial distribution (#3161)

Co-authored-by: Patrick Shriwise <pshriwise@gmail.com>
Co-authored-by: Paul Romano <paul.k.romano@gmail.com>
---
 docs/source/io_formats/settings.rst   |  56 ++++++++++---
 docs/source/pythonapi/stats.rst       |   1 +
 docs/source/usersguide/settings.rst   |   4 +-
 include/openmc/distribution_spatial.h |  20 +++++
 include/openmc/xml_interface.h        |   3 +
 openmc/stats/multivariate.py          | 114 ++++++++++++++++++++++++++
 src/distribution_spatial.cpp          |  45 ++++++++++
 src/xml_interface.cpp                 |  19 +++++
 tests/unit_tests/test_source.py       |  56 +++++++++++++
 9 files changed, 306 insertions(+), 12 deletions(-)

diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst
index 2c6c3265649..34b73fc55c5 100644
--- a/docs/source/io_formats/settings.rst
+++ b/docs/source/io_formats/settings.rst
@@ -579,24 +579,38 @@ attributes/sub-elements:
 
     :type:
       The type of spatial distribution. Valid options are "box", "fission",
-      "point", "cartesian", "cylindrical", and "spherical". A "box" spatial
-      distribution has coordinates sampled uniformly in a parallelepiped. A
-      "fission" spatial distribution samples locations from a "box"
+      "point", "cartesian", "cylindrical", "spherical", "mesh", and "cloud".
+
+      A "box" spatial distribution has coordinates sampled uniformly in a
+      parallelepiped.
+
+      A "fission" spatial distribution samples locations from a "box"
       distribution but only locations in fissionable materials are accepted.
+
       A "point" spatial distribution has coordinates specified by a triplet.
+
       A "cartesian" spatial distribution specifies independent distributions of
-      x-, y-, and z-coordinates. A "cylindrical" spatial distribution specifies
-      independent distributions of r-, phi-, and z-coordinates where phi is the
-      azimuthal angle and the origin for the cylindrical coordinate system is
-      specified by origin. A "spherical" spatial distribution specifies
-      independent distributions of r-, cos_theta-, and phi-coordinates where
-      cos_theta is the cosine of the angle with respect to the z-axis, phi is
-      the azimuthal angle, and the sphere is centered on the coordinate
-      (x0,y0,z0). A "mesh" spatial distribution samples source sites from a mesh element
+      x-, y-, and z-coordinates.
+
+      A "cylindrical" spatial distribution specifies independent distributions
+      of r-, phi-, and z-coordinates where phi is the azimuthal angle and the
+      origin for the cylindrical coordinate system is specified by origin.
+
+      A "spherical" spatial distribution specifies independent distributions of
+      r-, cos_theta-, and phi-coordinates where cos_theta is the cosine of the
+      angle with respect to the z-axis, phi is the azimuthal angle, and the
+      sphere is centered on the coordinate (x0,y0,z0).
+
+      A "mesh" spatial distribution samples source sites from a mesh element
       based on the relative strengths provided in the node. Source locations
       within an element are sampled isotropically. If no strengths are provided,
       the space within the mesh is uniformly sampled.
 
+      A "cloud" spatial distribution samples source sites from a list of spatial
+      positions provided in the node, based on the relative strengths provided
+      in the node. If no strengths are provided, the positions are uniformly
+      sampled.
+
       *Default*: None
 
     :parameters:
@@ -662,6 +676,26 @@ attributes/sub-elements:
       For "cylindrical and "spherical" distributions, this element specifies
       the coordinates for the origin of the coordinate system.
 
+    :mesh_id:
+      For "mesh" spatial distributions, this element specifies which mesh ID to
+      use for the geometric description of the mesh.
+
+    :coords:
+      For "cloud" distributions, this element specifies a list of coordinates
+      for each of the points in the cloud.
+
+    :strengths:
+      For "mesh" and "cloud" spatial distributions, this element specifies the
+      relative source strength of each mesh element or each point in the cloud.
+
+    :volume_normalized:
+      For "mesh" spatial distrubtions, this optional boolean element specifies
+      whether the vector of relative strengths should be multiplied by the mesh
+      element volume. This is most common if the strengths represent a source
+      per unit volume.
+
+      *Default*: false
+
   :angle:
     An element specifying the angular distribution of source sites. This element
     has the following attributes:
diff --git a/docs/source/pythonapi/stats.rst b/docs/source/pythonapi/stats.rst
index b72896c1860..c8318ba8620 100644
--- a/docs/source/pythonapi/stats.rst
+++ b/docs/source/pythonapi/stats.rst
@@ -59,6 +59,7 @@ Spatial Distributions
    openmc.stats.Box
    openmc.stats.Point
    openmc.stats.MeshSpatial
+   openmc.stats.PointCloud
 
 .. autosummary::
    :toctree: generated
diff --git a/docs/source/usersguide/settings.rst b/docs/source/usersguide/settings.rst
index 4111481a9ea..349aa34d07a 100644
--- a/docs/source/usersguide/settings.rst
+++ b/docs/source/usersguide/settings.rst
@@ -192,7 +192,9 @@ distributions using spherical or cylindrical coordinates, you can use
 :class:`openmc.stats.SphericalIndependent` or
 :class:`openmc.stats.CylindricalIndependent`, respectively. Meshes can also be
 used to represent spatial distributions with :class:`openmc.stats.MeshSpatial`
-by specifying a mesh and source strengths for each mesh element.
+by specifying a mesh and source strengths for each mesh element. It is also
+possible to define a "cloud" of source points, each with a different relative
+probability, using :class:`openmc.stats.PointCloud`.
 
 The angular distribution can be set equal to a sub-class of
 :class:`openmc.stats.UnitSphere` such as :class:`openmc.stats.Isotropic`,
diff --git a/include/openmc/distribution_spatial.h b/include/openmc/distribution_spatial.h
index 28568c890d8..8ff766d1333 100644
--- a/include/openmc/distribution_spatial.h
+++ b/include/openmc/distribution_spatial.h
@@ -136,6 +136,26 @@ class MeshSpatial : public SpatialDistribution {
                                 //!< mesh element indices
 };
 
+//==============================================================================
+//! Distribution of points
+//==============================================================================
+
+class PointCloud : public SpatialDistribution {
+public:
+  explicit PointCloud(pugi::xml_node node);
+  explicit PointCloud(
+    std::vector<Position> point_cloud, gsl::span<const double> strengths);
+
+  //! Sample a position from the distribution
+  //! \param seed Pseudorandom number seed pointer
+  //! \return Sampled position
+  Position sample(uint64_t* seed) const override;
+
+private:
+  std::vector<Position> point_cloud_;
+  DiscreteIndex point_idx_dist_; //!< Distribution of Position indices
+};
+
 //==============================================================================
 //! Uniform distribution of points over a box
 //==============================================================================
diff --git a/include/openmc/xml_interface.h b/include/openmc/xml_interface.h
index bd6554c134a..f49613ecde1 100644
--- a/include/openmc/xml_interface.h
+++ b/include/openmc/xml_interface.h
@@ -50,6 +50,9 @@ xt::xarray<T> get_node_xarray(
   return xt::adapt(v, shape);
 }
 
+std::vector<Position> get_node_position_array(
+  pugi::xml_node node, const char* name, bool lowercase = false);
+
 Position get_node_position(
   pugi::xml_node node, const char* name, bool lowercase = false);
 
diff --git a/openmc/stats/multivariate.py b/openmc/stats/multivariate.py
index 06c65896f49..cd474fa9b28 100644
--- a/openmc/stats/multivariate.py
+++ b/openmc/stats/multivariate.py
@@ -279,6 +279,8 @@ def from_xml_element(cls, elem, meshes=None):
             return Point.from_xml_element(elem)
         elif distribution == 'mesh':
             return MeshSpatial.from_xml_element(elem, meshes)
+        elif distribution == 'cloud':
+            return PointCloud.from_xml_element(elem)
 
 
 class CartesianIndependent(Spatial):
@@ -756,6 +758,118 @@ def from_xml_element(cls, elem, meshes):
         return cls(meshes[mesh_id], strengths, volume_normalized)
 
 
+class PointCloud(Spatial):
+    """Spatial distribution from a point cloud.
+
+    This distribution specifies a discrete list of points, with corresponding
+    relative probabilities.
+
+    .. versionadded:: 0.15.1
+
+    Parameters
+    ----------
+    positions : iterable of 3-tuples
+        The points in space to be sampled
+    strengths : iterable of float, optional
+        An iterable of values that represents the relative probabilty of each
+        point.
+
+    Attributes
+    ----------
+    positions : numpy.ndarray
+        The points in space to be sampled with shape (N, 3)
+    strengths : numpy.ndarray or None
+        An array of relative probabilities for each mesh point
+    """
+
+    def __init__(
+        self,
+        positions: Sequence[Sequence[float]],
+        strengths: Sequence[float] | None = None
+    ):
+        self.positions = positions
+        self.strengths = strengths
+
+    @property
+    def positions(self) -> np.ndarray:
+        return self._positions
+
+    @positions.setter
+    def positions(self, positions):
+        positions = np.array(positions, dtype=float)
+        if positions.ndim != 2:
+            raise ValueError('positions must be a 2D array')
+        elif positions.shape[1] != 3:
+            raise ValueError('Each position must have 3 values')
+        self._positions = positions
+
+    @property
+    def strengths(self) -> np.ndarray:
+        return self._strengths
+
+    @strengths.setter
+    def strengths(self, strengths):
+        if strengths is not None:
+            strengths = np.array(strengths, dtype=float)
+            if strengths.ndim != 1:
+                raise ValueError('strengths must be a 1D array')
+            elif strengths.size != self.positions.shape[0]:
+                raise ValueError('strengths must have the same length as positions')
+        self._strengths = strengths
+
+    @property
+    def num_strength_bins(self) -> int:
+        if self.strengths is None:
+            raise ValueError('Strengths are not set')
+        return self.strengths.size
+
+    def to_xml_element(self) -> ET.Element:
+        """Return XML representation of the spatial distribution
+
+        Returns
+        -------
+        element : lxml.etree._Element
+            XML element containing spatial distribution data
+
+        """
+        element = ET.Element('space')
+        element.set('type', 'cloud')
+
+        subelement = ET.SubElement(element, 'coords')
+        subelement.text = ' '.join(str(e) for e in self.positions.flatten())
+
+        if self.strengths is not None:
+            subelement = ET.SubElement(element, 'strengths')
+            subelement.text = ' '.join(str(e) for e in self.strengths)
+
+        return element
+
+    @classmethod
+    def from_xml_element(cls, elem: ET.Element) -> PointCloud:
+        """Generate spatial distribution from an XML element
+
+        Parameters
+        ----------
+        elem : lxml.etree._Element
+            XML element
+
+        Returns
+        -------
+        openmc.stats.PointCloud
+            Spatial distribution generated from XML element
+
+
+        """
+        coord_data = get_text(elem, 'coords')
+        positions = np.array([float(b) for b in coord_data.split()]).reshape((-1, 3))
+
+        strengths = get_text(elem, 'strengths')
+        if strengths is not None:
+            strengths = [float(b) for b in strengths.split()]
+
+        return cls(positions, strengths)
+
+
 class Box(Spatial):
     """Uniform distribution of coordinates in a rectangular cuboid.
 
diff --git a/src/distribution_spatial.cpp b/src/distribution_spatial.cpp
index 8dbd7d88706..5c193a95d29 100644
--- a/src/distribution_spatial.cpp
+++ b/src/distribution_spatial.cpp
@@ -26,6 +26,8 @@ unique_ptr<SpatialDistribution> SpatialDistribution::create(pugi::xml_node node)
     return UPtrSpace {new SphericalIndependent(node)};
   } else if (type == "mesh") {
     return UPtrSpace {new MeshSpatial(node)};
+  } else if (type == "cloud") {
+    return UPtrSpace {new PointCloud(node)};
   } else if (type == "box") {
     return UPtrSpace {new SpatialBox(node)};
   } else if (type == "fission") {
@@ -298,6 +300,49 @@ Position MeshSpatial::sample(uint64_t* seed) const
   return this->sample_mesh(seed).second;
 }
 
+//==============================================================================
+// PointCloud implementation
+//==============================================================================
+
+PointCloud::PointCloud(pugi::xml_node node)
+{
+  if (check_for_node(node, "coords")) {
+    point_cloud_ = get_node_position_array(node, "coords");
+  } else {
+    fatal_error("No coordinates were provided for the PointCloud "
+                "spatial distribution");
+  }
+
+  std::vector<double> strengths;
+
+  if (check_for_node(node, "strengths"))
+    strengths = get_node_array<double>(node, "strengths");
+  else
+    strengths.resize(point_cloud_.size(), 1.0);
+
+  if (strengths.size() != point_cloud_.size()) {
+    fatal_error(
+      fmt::format("Number of entries for the strengths array {} does "
+                  "not match the number of spatial points provided {}.",
+        strengths.size(), point_cloud_.size()));
+  }
+
+  point_idx_dist_.assign(strengths);
+}
+
+PointCloud::PointCloud(
+  std::vector<Position> point_cloud, gsl::span<const double> strengths)
+{
+  point_cloud_.assign(point_cloud.begin(), point_cloud.end());
+  point_idx_dist_.assign(strengths);
+}
+
+Position PointCloud::sample(uint64_t* seed) const
+{
+  int32_t index = point_idx_dist_.sample(seed);
+  return point_cloud_[index];
+}
+
 //==============================================================================
 // SpatialBox implementation
 //==============================================================================
diff --git a/src/xml_interface.cpp b/src/xml_interface.cpp
index 6ce465bafb9..840d3f5b871 100644
--- a/src/xml_interface.cpp
+++ b/src/xml_interface.cpp
@@ -4,6 +4,7 @@
 
 #include "openmc/error.h"
 #include "openmc/string_utils.h"
+#include "openmc/vector.h"
 
 namespace openmc {
 
@@ -48,6 +49,24 @@ bool get_node_value_bool(pugi::xml_node node, const char* name)
   return false;
 }
 
+vector<Position> get_node_position_array(
+  pugi::xml_node node, const char* name, bool lowercase)
+{
+  vector<double> coords = get_node_array<double>(node, name, lowercase);
+  if (coords.size() % 3 != 0) {
+    fatal_error(fmt::format(
+      "Incorect number of coordinates in Position array ({}) for \"{}\"",
+      coords.size(), name));
+  }
+  vector<Position> positions;
+  positions.reserve(coords.size() / 3);
+  auto it = coords.begin();
+  for (size_t i = 0; i < coords.size(); i += 3) {
+    positions.push_back({coords[i], coords[i + 1], coords[i + 2]});
+  }
+  return positions;
+}
+
 Position get_node_position(
   pugi::xml_node node, const char* name, bool lowercase)
 {
diff --git a/tests/unit_tests/test_source.py b/tests/unit_tests/test_source.py
index 32650d54936..c88fbcbe62d 100644
--- a/tests/unit_tests/test_source.py
+++ b/tests/unit_tests/test_source.py
@@ -1,3 +1,4 @@
+from collections import Counter
 from math import pi
 
 import openmc
@@ -49,6 +50,61 @@ def test_spherical_uniform():
 
     assert isinstance(sph_indep_function, openmc.stats.SphericalIndependent)
 
+def test_point_cloud():
+    positions = [(1, 0, 2), (0, 1, 0), (0, 0, 3), (4, 9, 2)]
+    strengths = [1, 2, 3, 4]
+
+    space = openmc.stats.PointCloud(positions, strengths)
+    np.testing.assert_equal(space.positions, positions)
+    np.testing.assert_equal(space.strengths, strengths)
+
+    src = openmc.IndependentSource(space=space)
+    assert src.space == space
+    np.testing.assert_equal(src.space.positions, positions)
+    np.testing.assert_equal(src.space.strengths, strengths)
+
+    elem = src.to_xml_element()
+    src = openmc.IndependentSource.from_xml_element(elem)
+    np.testing.assert_equal(src.space.positions, positions)
+    np.testing.assert_equal(src.space.strengths, strengths)
+
+
+def test_point_cloud_invalid():
+    with pytest.raises(ValueError, match='2D'):
+        openmc.stats.PointCloud([1, 0, 2, 0, 1, 0])
+
+    with pytest.raises(ValueError, match='3 values'):
+        openmc.stats.PointCloud([(1, 0, 2, 3), (4, 5, 2, 3)])
+
+    with pytest.raises(ValueError, match='1D'):
+        openmc.stats.PointCloud([(1, 0, 2), (4, 5, 2)], [(1, 2), (3, 4)])
+
+    with pytest.raises(ValueError, match='same length'):
+        openmc.stats.PointCloud([(1, 0, 2), (4, 5, 2)], [1, 2, 4])
+
+
+def test_point_cloud_strengths(run_in_tmpdir, sphere_box_model):
+    positions = [(1., 0., 2.), (0., 1., 0.), (0., 0., 3.), (-1., -1., 2.)]
+    strengths = [1, 2, 3, 4]
+    space = openmc.stats.PointCloud(positions, strengths)
+
+    model = sphere_box_model[0]
+    model.settings.run_mode = 'fixed source'
+    model.settings.source = openmc.IndependentSource(space=space)
+
+    try:
+        model.init_lib()
+        n_samples = 50_000
+        sites = openmc.lib.sample_external_source(n_samples)
+    finally:
+        model.finalize_lib()
+
+    count = Counter(s.r for s in sites)
+    for i, (strength, position) in enumerate(zip(strengths, positions)):
+        sampled_strength = count[position] / n_samples
+        expected_strength = pytest.approx(strength/sum(strengths), abs=0.02)
+        assert sampled_strength == expected_strength, f'Strength incorrect for {positions[i]}'
+
 
 def test_source_file():
     filename = 'source.h5'

From 6fd184a5e70fe8bac5a509a04f150167d2819859 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Wed, 13 Nov 2024 15:58:22 -0600
Subject: [PATCH 10/24] Fix docstring for Model.plot (#3198)

---
 openmc/model/model.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/openmc/model/model.py b/openmc/model/model.py
index 0ea6db8db88..ffad8f503c7 100644
--- a/openmc/model/model.py
+++ b/openmc/model/model.py
@@ -820,14 +820,14 @@ def plot(
 
         Parameters
         ----------
-        n_samples : dict
+        n_samples : int, optional
             The number of source particles to sample and add to plot. Defaults
             to None which doesn't plot any particles on the plot.
         plane_tolerance: float
             When plotting a plane the source locations within the plane +/-
             the plane_tolerance will be included and those outside of the
             plane_tolerance will not be shown
-        source_kwargs : dict
+        source_kwargs : dict, optional
             Keyword arguments passed to :func:`matplotlib.pyplot.scatter`.
         **kwargs
             Keyword arguments passed to :func:`openmc.Universe.plot`

From 3cd9ca98c8b39c3728bdfec21b13d69fdefc15fc Mon Sep 17 00:00:00 2001
From: magnoxemo <walid.ahammed01@gmail.com>
Date: Fri, 8 Nov 2024 15:56:33 -0600
Subject: [PATCH 11/24] init_alias is private. I have to figure how to acces
 that in the settings.cpp

---
 src/settings.cpp |  6 ++++++
 src/source.cpp   | 16 ++++++++++++----
 2 files changed, 18 insertions(+), 4 deletions(-)

diff --git a/src/settings.cpp b/src/settings.cpp
index 3a905bdf961..d0af5e17911 100644
--- a/src/settings.cpp
+++ b/src/settings.cpp
@@ -32,11 +32,14 @@
 #include "openmc/volume_calc.h"
 #include "openmc/weight_windows.h"
 #include "openmc/xml_interface.h"
+#include "openmc/distribution.h"
+#include "openmc/source.h"
 
 namespace openmc {
 
 //==============================================================================
 // Global variables
+external_sources_alias_sampler = DiscreteIndex();
 //==============================================================================
 
 namespace settings {
@@ -558,6 +561,9 @@ void read_settings_xml(pugi::xml_node root)
     model::external_sources.push_back(make_unique<FileSource>(path));
   }
 
+  external_sources_alias_sampler = DiscreteIndex(model::externa_sources);
+  external_sources_alias_sampler.init_alias();
+
   // If no source specified, default to isotropic point source at origin with
   // Watt spectrum. No default source is needed in random ray mode.
   if (model::external_sources.empty() &&
diff --git a/src/source.cpp b/src/source.cpp
index ae77cc1e06c..3e6f8d385cd 100644
--- a/src/source.cpp
+++ b/src/source.cpp
@@ -609,12 +609,19 @@ void initialize_source()
 SourceSite sample_external_source(uint64_t* seed)
 {
   // Determine total source strength
+  /*
   double total_strength = 0.0;
   for (auto& s : model::external_sources)
-    total_strength += s->strength();
+  {
+    total_strength += s->strength(); //is that strength normalized? if yes then it can be used as a pdf 
+  }
+  */
+  //external_sources_alias_sampler = DiscreteIndex(model::external_sources);
+  //external_sources_alias_sampler.init_alias();
+  int i=external_sources_alias_sampler.sample();
 
   // Sample from among multiple source distributions
-  int i = 0;
+  /*int i = 0;
   if (model::external_sources.size() > 1) {
     double xi = prn(seed) * total_strength; //random strength from the total source strength 
     double c = 0.0;
@@ -625,13 +632,14 @@ SourceSite sample_external_source(uint64_t* seed)
     c accumulates the strengths of each source sequentially. The loop increments c with 
     the strength of each source until xi < c, indicating that the randomly chosen point
     falls within the range associated with a specific source. The index i at this point
-    identifies the selected source.*/
-    
+    identifies the selected source.
+
     for (; i < model::external_sources.size(); ++i) {
       c += model::external_sources[i]->strength();
       if (xi < c)
         break;
     }
+    */
   }
 
   // Sample source site from i-th source distribution

From ba55c037d35e2e305a0cf5f9226f71f4b3144a40 Mon Sep 17 00:00:00 2001
From: magnoxemo <walid.ahammed01@gmail.com>
Date: Sat, 9 Nov 2024 10:07:07 -0600
Subject: [PATCH 12/24] ig it works?

---
 include/openmc/distribution.h | 1 -
 src/settings.cpp              | 3 +--
 src/source.cpp                | 7 ++++---
 3 files changed, 5 insertions(+), 6 deletions(-)

diff --git a/include/openmc/distribution.h b/include/openmc/distribution.h
index e70c803de2f..0026ba23276 100644
--- a/include/openmc/distribution.h
+++ b/include/openmc/distribution.h
@@ -8,7 +8,6 @@
 
 #include "pugixml.hpp"
 #include <gsl/gsl-lite.hpp>
-
 #include "openmc/constants.h"
 #include "openmc/memory.h" // for unique_ptr
 #include "openmc/vector.h" // for vector
diff --git a/src/settings.cpp b/src/settings.cpp
index d0af5e17911..f7689bfa64d 100644
--- a/src/settings.cpp
+++ b/src/settings.cpp
@@ -561,8 +561,7 @@ void read_settings_xml(pugi::xml_node root)
     model::external_sources.push_back(make_unique<FileSource>(path));
   }
 
-  external_sources_alias_sampler = DiscreteIndex(model::externa_sources);
-  external_sources_alias_sampler.init_alias();
+  external_sources_alias_sampler = DiscreteIndex();
 
   // If no source specified, default to isotropic point source at origin with
   // Watt spectrum. No default source is needed in random ray mode.
diff --git a/src/source.cpp b/src/source.cpp
index 3e6f8d385cd..bd8480a516f 100644
--- a/src/source.cpp
+++ b/src/source.cpp
@@ -609,15 +609,16 @@ void initialize_source()
 SourceSite sample_external_source(uint64_t* seed)
 {
   // Determine total source strength
-  /*
-  double total_strength = 0.0;
+  
+  std::vector <double> source_strength ;
   for (auto& s : model::external_sources)
   {
-    total_strength += s->strength(); //is that strength normalized? if yes then it can be used as a pdf 
+    source_strength.push_back(s->strength()); //is that strength normalized? if yes then it can be used as a pdf 
   }
   */
   //external_sources_alias_sampler = DiscreteIndex(model::external_sources);
   //external_sources_alias_sampler.init_alias();
+  external_sources_alias_sampler=DiscreteIndex(source_strength);
   int i=external_sources_alias_sampler.sample();
 
   // Sample from among multiple source distributions

From 2647dc9492ae1cf886617a87766d207e087a5595 Mon Sep 17 00:00:00 2001
From: magnoxemo <walid.ahammed01@gmail.com>
Date: Sat, 9 Nov 2024 10:45:21 -0600
Subject: [PATCH 13/24] not sure

---
 src/distribution.cpp |  2 +-
 src/settings.cpp     | 17 ++++++++++++++---
 src/source.cpp       | 26 ++++++++++++--------------
 3 files changed, 27 insertions(+), 18 deletions(-)

diff --git a/src/distribution.cpp b/src/distribution.cpp
index a6b4acd58b1..3cc54ee0e6a 100644
--- a/src/distribution.cpp
+++ b/src/distribution.cpp
@@ -101,7 +101,7 @@ size_t DiscreteIndex::sample(uint64_t* seed) const
     return 0;
   }
 }
-
+assign
 void DiscreteIndex::normalize()
 {
   // Renormalize density function so that it sums to unity. Note that we save
diff --git a/src/settings.cpp b/src/settings.cpp
index f7689bfa64d..4e138f55ba2 100644
--- a/src/settings.cpp
+++ b/src/settings.cpp
@@ -39,7 +39,7 @@ namespace openmc {
 
 //==============================================================================
 // Global variables
-external_sources_alias_sampler = DiscreteIndex();
+DiscreteIndex external_sources_alias_sampler;
 //==============================================================================
 
 namespace settings {
@@ -560,8 +560,19 @@ void read_settings_xml(pugi::xml_node root)
     }
     model::external_sources.push_back(make_unique<FileSource>(path));
   }
-
-  external_sources_alias_sampler = DiscreteIndex();
+  std::vector <double> source_strength ;
+  for (auto& s : model::external_sources)
+  {
+    source_strength.push_back(s->strength()); //is that strength normalized? if yes then it can be used as a pdf 
+  }
+  */
+  //external_sources_alias_sampler = DiscreteIndex(model::external_sources);
+  //external_sources_alias_sampler.init_alias();
+  vector<double> strengths;
+  for (auto& s : model::external_sources) {
+    strengths.push_back(s->strength());
+  }
+  external_source_alias_sampler.assign(strengths);
 
   // If no source specified, default to isotropic point source at origin with
   // Watt spectrum. No default source is needed in random ray mode.
diff --git a/src/source.cpp b/src/source.cpp
index bd8480a516f..8737ad68892 100644
--- a/src/source.cpp
+++ b/src/source.cpp
@@ -608,24 +608,20 @@ void initialize_source()
 
 SourceSite sample_external_source(uint64_t* seed)
 {
+  /*
   // Determine total source strength
-  
-  std::vector <double> source_strength ;
+  double total_strength = 0.0;
   for (auto& s : model::external_sources)
-  {
-    source_strength.push_back(s->strength()); //is that strength normalized? if yes then it can be used as a pdf 
-  }
+    total_strength += s->strength();
+
+  "I will be only needing the strength and normalize it in the DiscreteIndex class"
   */
-  //external_sources_alias_sampler = DiscreteIndex(model::external_sources);
-  //external_sources_alias_sampler.init_alias();
-  external_sources_alias_sampler=DiscreteIndex(source_strength);
-  int i=external_sources_alias_sampler.sample();
 
   // Sample from among multiple source distributions
-  /*int i = 0;
+  int i = 0;
   if (model::external_sources.size() > 1) {
-    double xi = prn(seed) * total_strength; //random strength from the total source strength 
-    double c = 0.0;
+    //double xi = prn(seed) * total_strength; //random strength from the total source strength 
+    //double c = 0.0;
 
     /* First, it checks if there are multiple sources in model::external_sources. 
     If so, it proceeds to sample one. xi is a random value scaled by total_strength.
@@ -634,13 +630,15 @@ SourceSite sample_external_source(uint64_t* seed)
     the strength of each source until xi < c, indicating that the randomly chosen point
     falls within the range associated with a specific source. The index i at this point
     identifies the selected source.
-
+    
     for (; i < model::external_sources.size(); ++i) {
       c += model::external_sources[i]->strength();
       if (xi < c)
         break;
     }
-    */
+  */
+  int i = source_alias_sampler.sample(seed);
+
   }
 
   // Sample source site from i-th source distribution

From 648510ce8fbaf7709c3342835f2ad0d6cf8e7487 Mon Sep 17 00:00:00 2001
From: magnoxemo <walid.ahammed01@gmail.com>
Date: Sat, 9 Nov 2024 20:03:45 -0600
Subject: [PATCH 14/24] implementing the suggested changes

---
 include/openmc/distribution.h |  4 ++--
 src/distribution.cpp          |  2 +-
 src/settings.cpp              | 19 ++++++-------------
 src/source.cpp                | 29 +----------------------------
 4 files changed, 10 insertions(+), 44 deletions(-)

diff --git a/include/openmc/distribution.h b/include/openmc/distribution.h
index 0026ba23276..3b3acf0a0b2 100644
--- a/include/openmc/distribution.h
+++ b/include/openmc/distribution.h
@@ -6,11 +6,11 @@
 
 #include <cstddef> // for size_t
 
-#include "pugixml.hpp"
-#include <gsl/gsl-lite.hpp>
 #include "openmc/constants.h"
 #include "openmc/memory.h" // for unique_ptr
 #include "openmc/vector.h" // for vector
+#include "pugixml.hpp"
+#include <gsl/gsl-lite.hpp>
 
 namespace openmc {
 
diff --git a/src/distribution.cpp b/src/distribution.cpp
index 3cc54ee0e6a..a6b4acd58b1 100644
--- a/src/distribution.cpp
+++ b/src/distribution.cpp
@@ -101,7 +101,7 @@ size_t DiscreteIndex::sample(uint64_t* seed) const
     return 0;
   }
 }
-assign
+
 void DiscreteIndex::normalize()
 {
   // Renormalize density function so that it sums to unity. Note that we save
diff --git a/src/settings.cpp b/src/settings.cpp
index 4e138f55ba2..4b660b09214 100644
--- a/src/settings.cpp
+++ b/src/settings.cpp
@@ -32,18 +32,18 @@
 #include "openmc/volume_calc.h"
 #include "openmc/weight_windows.h"
 #include "openmc/xml_interface.h"
-#include "openmc/distribution.h"
-#include "openmc/source.h"
 
 namespace openmc {
 
 //==============================================================================
 // Global variables
-DiscreteIndex external_sources_alias_sampler;
 //==============================================================================
 
 namespace settings {
 
+// DiscreteIndex class for the alias sampling
+
+DiscreteIndex external_sources_alias_sampler;
 // Default values for boolean flags
 bool assume_separate {false};
 bool check_overlaps {false};
@@ -560,19 +560,12 @@ void read_settings_xml(pugi::xml_node root)
     }
     model::external_sources.push_back(make_unique<FileSource>(path));
   }
-  std::vector <double> source_strength ;
-  for (auto& s : model::external_sources)
-  {
-    source_strength.push_back(s->strength()); //is that strength normalized? if yes then it can be used as a pdf 
-  }
-  */
-  //external_sources_alias_sampler = DiscreteIndex(model::external_sources);
-  //external_sources_alias_sampler.init_alias();
-  vector<double> strengths;
+
+  vector<double> source_strengths;
   for (auto& s : model::external_sources) {
     strengths.push_back(s->strength());
   }
-  external_source_alias_sampler.assign(strengths);
+  external_source_alias_sampler.assign(source_strengths);
 
   // If no source specified, default to isotropic point source at origin with
   // Watt spectrum. No default source is needed in random ray mode.
diff --git a/src/source.cpp b/src/source.cpp
index 8737ad68892..0a58c67952c 100644
--- a/src/source.cpp
+++ b/src/source.cpp
@@ -608,37 +608,10 @@ void initialize_source()
 
 SourceSite sample_external_source(uint64_t* seed)
 {
-  /*
-  // Determine total source strength
-  double total_strength = 0.0;
-  for (auto& s : model::external_sources)
-    total_strength += s->strength();
-
-  "I will be only needing the strength and normalize it in the DiscreteIndex class"
-  */
-
   // Sample from among multiple source distributions
   int i = 0;
   if (model::external_sources.size() > 1) {
-    //double xi = prn(seed) * total_strength; //random strength from the total source strength 
-    //double c = 0.0;
-
-    /* First, it checks if there are multiple sources in model::external_sources. 
-    If so, it proceeds to sample one. xi is a random value scaled by total_strength.
-    This random value (xi) is used to decide which source distribution to sample from.
-    c accumulates the strengths of each source sequentially. The loop increments c with 
-    the strength of each source until xi < c, indicating that the randomly chosen point
-    falls within the range associated with a specific source. The index i at this point
-    identifies the selected source.
-    
-    for (; i < model::external_sources.size(); ++i) {
-      c += model::external_sources[i]->strength();
-      if (xi < c)
-        break;
-    }
-  */
-  int i = source_alias_sampler.sample(seed);
-
+    int i = externaal_source_alias_sampler.sample(seed);
   }
 
   // Sample source site from i-th source distribution

From 7a1b4ad3f8f0076cde15ef5110deabd93d25a4df Mon Sep 17 00:00:00 2001
From: magnoxemo <walid.ahammed01@gmail.com>
Date: Sat, 9 Nov 2024 20:14:23 -0600
Subject: [PATCH 15/24] don't need the int i again cause already decleared the
 data type

---
 src/distribution.cpp | 1 +
 src/source.cpp       | 2 +-
 2 files changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/distribution.cpp b/src/distribution.cpp
index a6b4acd58b1..fc20e17a2ee 100644
--- a/src/distribution.cpp
+++ b/src/distribution.cpp
@@ -16,6 +16,7 @@
 #include "openmc/random_lcg.h"
 #include "openmc/xml_interface.h"
 
+
 namespace openmc {
 
 //==============================================================================
diff --git a/src/source.cpp b/src/source.cpp
index 0a58c67952c..94d2e60f249 100644
--- a/src/source.cpp
+++ b/src/source.cpp
@@ -611,7 +611,7 @@ SourceSite sample_external_source(uint64_t* seed)
   // Sample from among multiple source distributions
   int i = 0;
   if (model::external_sources.size() > 1) {
-    int i = externaal_source_alias_sampler.sample(seed);
+    i = externaal_source_alias_sampler.sample(seed);
   }
 
   // Sample source site from i-th source distribution

From 93afa218436ac3668e189a1674342d05981daba4 Mon Sep 17 00:00:00 2001
From: magnoxemo <walid.ahammed01@gmail.com>
Date: Sat, 9 Nov 2024 22:13:49 -0600
Subject: [PATCH 16/24] fixing the typos

---
 include/openmc/distribution.h | 6 ++++--
 src/settings.cpp              | 5 +----
 src/source.cpp                | 5 ++++-
 3 files changed, 9 insertions(+), 7 deletions(-)

diff --git a/include/openmc/distribution.h b/include/openmc/distribution.h
index 3b3acf0a0b2..8b360c8d664 100644
--- a/include/openmc/distribution.h
+++ b/include/openmc/distribution.h
@@ -6,11 +6,13 @@
 
 #include <cstddef> // for size_t
 
+#include "pugixml.hpp"
+#include <gsl/gsl-lite.hpp>
+
 #include "openmc/constants.h"
 #include "openmc/memory.h" // for unique_ptr
 #include "openmc/vector.h" // for vector
-#include "pugixml.hpp"
-#include <gsl/gsl-lite.hpp>
+
 
 namespace openmc {
 
diff --git a/src/settings.cpp b/src/settings.cpp
index 4b660b09214..828ee1b3160 100644
--- a/src/settings.cpp
+++ b/src/settings.cpp
@@ -41,9 +41,6 @@ namespace openmc {
 
 namespace settings {
 
-// DiscreteIndex class for the alias sampling
-
-DiscreteIndex external_sources_alias_sampler;
 // Default values for boolean flags
 bool assume_separate {false};
 bool check_overlaps {false};
@@ -563,7 +560,7 @@ void read_settings_xml(pugi::xml_node root)
 
   vector<double> source_strengths;
   for (auto& s : model::external_sources) {
-    strengths.push_back(s->strength());
+    source_strengths.push_back(s->strength());
   }
   external_source_alias_sampler.assign(source_strengths);
 
diff --git a/src/source.cpp b/src/source.cpp
index 94d2e60f249..eaf8c9beebd 100644
--- a/src/source.cpp
+++ b/src/source.cpp
@@ -44,6 +44,9 @@ namespace openmc {
 namespace model {
 
 vector<unique_ptr<Source>> external_sources;
+
+// DiscreteIndex class for the alias sampling
+DiscreteIndex external_sources_alias_sampler;
 }
 
 //==============================================================================
@@ -611,7 +614,7 @@ SourceSite sample_external_source(uint64_t* seed)
   // Sample from among multiple source distributions
   int i = 0;
   if (model::external_sources.size() > 1) {
-    i = externaal_source_alias_sampler.sample(seed);
+    i = external_source_alias_sampler.sample(seed);
   }
 
   // Sample source site from i-th source distribution

From 7dccde7fc037e764f7c35a27c3f522e716971e1b Mon Sep 17 00:00:00 2001
From: magnoxemo <walid.ahammed01@gmail.com>
Date: Mon, 11 Nov 2024 09:06:52 -0600
Subject: [PATCH 17/24] removing unnecessary spaces

---
 include/openmc/distribution.h | 1 -
 src/distribution.cpp          | 1 -
 2 files changed, 2 deletions(-)

diff --git a/include/openmc/distribution.h b/include/openmc/distribution.h
index 8b360c8d664..e70c803de2f 100644
--- a/include/openmc/distribution.h
+++ b/include/openmc/distribution.h
@@ -13,7 +13,6 @@
 #include "openmc/memory.h" // for unique_ptr
 #include "openmc/vector.h" // for vector
 
-
 namespace openmc {
 
 //==============================================================================
diff --git a/src/distribution.cpp b/src/distribution.cpp
index fc20e17a2ee..a6b4acd58b1 100644
--- a/src/distribution.cpp
+++ b/src/distribution.cpp
@@ -16,7 +16,6 @@
 #include "openmc/random_lcg.h"
 #include "openmc/xml_interface.h"
 
-
 namespace openmc {
 
 //==============================================================================

From 16431501f0f38104b626ab0c7bd503b9b768a425 Mon Sep 17 00:00:00 2001
From: magnoxemo <walid.ahammed01@gmail.com>
Date: Fri, 15 Nov 2024 15:49:05 -0600
Subject: [PATCH 18/24] defining the external_sources_alias_sample in souce.cpp
 cause it's failing the tests

---
 include/openmc/source.h | 2 ++
 src/source.cpp          | 2 --
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/include/openmc/source.h b/include/openmc/source.h
index 5cc0356e3e9..253f5a457f3 100644
--- a/include/openmc/source.h
+++ b/include/openmc/source.h
@@ -35,6 +35,8 @@ class Source;
 namespace model {
 
 extern vector<unique_ptr<Source>> external_sources;
+// DiscreteIndex class for the alias sampling
+DiscreteIndex external_sources_alias_sampler;
 
 } // namespace model
 
diff --git a/src/source.cpp b/src/source.cpp
index eaf8c9beebd..a70f2bd38bd 100644
--- a/src/source.cpp
+++ b/src/source.cpp
@@ -45,8 +45,6 @@ namespace model {
 
 vector<unique_ptr<Source>> external_sources;
 
-// DiscreteIndex class for the alias sampling
-DiscreteIndex external_sources_alias_sampler;
 }
 
 //==============================================================================

From f0656bd6b6870a7f755ef7af100f79a1157c8ef5 Mon Sep 17 00:00:00 2001
From: magnoxemo <walid.ahammed01@gmail.com>
Date: Fri, 15 Nov 2024 23:08:35 -0600
Subject: [PATCH 19/24] needed to use the extern in source.h and redefine in
 the source.cpp in openmc model namespace. ig it passes all the tests now

---
 include/openmc/source.h | 2 +-
 src/settings.cpp        | 2 +-
 src/source.cpp          | 5 ++++-
 3 files changed, 6 insertions(+), 3 deletions(-)

diff --git a/include/openmc/source.h b/include/openmc/source.h
index 253f5a457f3..2b8839abc62 100644
--- a/include/openmc/source.h
+++ b/include/openmc/source.h
@@ -36,7 +36,7 @@ namespace model {
 
 extern vector<unique_ptr<Source>> external_sources;
 // DiscreteIndex class for the alias sampling
-DiscreteIndex external_sources_alias_sampler;
+extern DiscreteIndex external_sources_alias_sampler;
 
 } // namespace model
 
diff --git a/src/settings.cpp b/src/settings.cpp
index 828ee1b3160..c198f3d65db 100644
--- a/src/settings.cpp
+++ b/src/settings.cpp
@@ -562,7 +562,7 @@ void read_settings_xml(pugi::xml_node root)
   for (auto& s : model::external_sources) {
     source_strengths.push_back(s->strength());
   }
-  external_source_alias_sampler.assign(source_strengths);
+  openmc::model::external_sources_alias_sampler.assign(source_strengths);
 
   // If no source specified, default to isotropic point source at origin with
   // Watt spectrum. No default source is needed in random ray mode.
diff --git a/src/source.cpp b/src/source.cpp
index a70f2bd38bd..a5dbb35bc37 100644
--- a/src/source.cpp
+++ b/src/source.cpp
@@ -45,6 +45,9 @@ namespace model {
 
 vector<unique_ptr<Source>> external_sources;
 
+// DiscreteIndex class for the alias sampling
+DiscreteIndex external_sources_alias_sampler;
+
 }
 
 //==============================================================================
@@ -612,7 +615,7 @@ SourceSite sample_external_source(uint64_t* seed)
   // Sample from among multiple source distributions
   int i = 0;
   if (model::external_sources.size() > 1) {
-    i = external_source_alias_sampler.sample(seed);
+    i = model::external_sources_alias_sampler.sample(seed);
   }
 
   // Sample source site from i-th source distribution

From b8e0ca6ec5b130682ab83caef93c6b0331c0d303 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Fri, 22 Nov 2024 19:31:45 -0600
Subject: [PATCH 20/24] Fix uniform source sampling

---
 src/settings.cpp | 2 +-
 src/source.cpp   | 8 ++++++--
 2 files changed, 7 insertions(+), 3 deletions(-)

diff --git a/src/settings.cpp b/src/settings.cpp
index 3fdbb3cd512..cee5b830518 100644
--- a/src/settings.cpp
+++ b/src/settings.cpp
@@ -567,7 +567,7 @@ void read_settings_xml(pugi::xml_node root)
   for (auto& s : model::external_sources) {
     source_strengths.push_back(s->strength());
   }
-  openmc::model::external_sources_alias_sampler.assign(source_strengths);
+  model::external_sources_alias_sampler.assign(source_strengths);
 
   // If no source specified, default to isotropic point source at origin with
   // Watt spectrum. No default source is needed in random ray mode.
diff --git a/src/source.cpp b/src/source.cpp
index 1d9a0af736c..e4303c1fa5a 100644
--- a/src/source.cpp
+++ b/src/source.cpp
@@ -48,7 +48,7 @@ vector<unique_ptr<Source>> external_sources;
 // DiscreteIndex class for the alias sampling
 DiscreteIndex external_sources_alias_sampler;
 
-}
+} // namespace model
 
 //==============================================================================
 // Source implementation
@@ -615,7 +615,11 @@ SourceSite sample_external_source(uint64_t* seed)
   // Sample from among multiple source distributions
   int i = 0;
   if (model::external_sources.size() > 1) {
-    i = model::external_sources_alias_sampler.sample(seed);
+    if (settings::uniform_source_sampling) {
+      i = prn(seed) * model::external_sources.size();
+    } else {
+      i = model::external_sources_alias_sampler.sample(seed);
+    }
   }
 
   // Sample source site from i-th source distribution

From 8ef337969e86b995f91ee84a1961140ba8a377a8 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Fri, 22 Nov 2024 19:44:41 -0600
Subject: [PATCH 21/24] Add test for multiple sources

---
 openmc/model/model.py                            |  4 ++--
 tests/unit_tests/test_uniform_source_sampling.py | 12 ++++++++++++
 2 files changed, 14 insertions(+), 2 deletions(-)

diff --git a/openmc/model/model.py b/openmc/model/model.py
index ffad8f503c7..49cb0a83ef5 100644
--- a/openmc/model/model.py
+++ b/openmc/model/model.py
@@ -885,7 +885,7 @@ def sample_external_source(
             n_samples: int = 1000,
             prn_seed: int | None = None,
             **init_kwargs
-    ) -> list[openmc.SourceParticle]:
+    ) -> openmc.ParticleList:
         """Sample external source and return source particles.
 
         .. versionadded:: 0.15.1
@@ -902,7 +902,7 @@ def sample_external_source(
 
         Returns
         -------
-        list of openmc.SourceParticle
+        openmc.ParticleList
             List of samples source particles
         """
         import openmc.lib
diff --git a/tests/unit_tests/test_uniform_source_sampling.py b/tests/unit_tests/test_uniform_source_sampling.py
index ece8afe8946..7f805e37d27 100644
--- a/tests/unit_tests/test_uniform_source_sampling.py
+++ b/tests/unit_tests/test_uniform_source_sampling.py
@@ -61,3 +61,15 @@ def test_tally_mean(run_in_tmpdir, sphere_model):
 
     # Check that tally means match
     assert mean == pytest.approx(reference_mean)
+
+
+def test_multiple_sources(sphere_model):
+    low_strength_src = openmc.IndependentSource(
+        energy=openmc.stats.delta_function(1.0e6), strength=1e-7)
+    sphere_model.settings.source.append(low_strength_src)
+    sphere_model.settings.uniform_source_sampling = True
+
+    # Sample particles from source and make sure 1 MeV shows up despite
+    # negligible strength
+    particles = sphere_model.sample_external_source(100)
+    assert {p.E for p in particles} == {1.0e3, 1.0e6}

From 47db5a77bd248992e4dc8c544c85b540405c0758 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Fri, 22 Nov 2024 19:46:25 -0600
Subject: [PATCH 22/24] Change name to external_sources_probability

---
 include/openmc/source.h | 5 +++--
 src/settings.cpp        | 2 +-
 src/source.cpp          | 4 ++--
 3 files changed, 6 insertions(+), 5 deletions(-)

diff --git a/include/openmc/source.h b/include/openmc/source.h
index 2b8839abc62..6733eaeffca 100644
--- a/include/openmc/source.h
+++ b/include/openmc/source.h
@@ -35,8 +35,9 @@ class Source;
 namespace model {
 
 extern vector<unique_ptr<Source>> external_sources;
-// DiscreteIndex class for the alias sampling
-extern DiscreteIndex external_sources_alias_sampler;
+
+// Probability distribution for selecting external sources
+extern DiscreteIndex external_sources_probability;
 
 } // namespace model
 
diff --git a/src/settings.cpp b/src/settings.cpp
index cee5b830518..eda0f8ee71c 100644
--- a/src/settings.cpp
+++ b/src/settings.cpp
@@ -567,7 +567,7 @@ void read_settings_xml(pugi::xml_node root)
   for (auto& s : model::external_sources) {
     source_strengths.push_back(s->strength());
   }
-  model::external_sources_alias_sampler.assign(source_strengths);
+  model::external_sources_probability.assign(source_strengths);
 
   // If no source specified, default to isotropic point source at origin with
   // Watt spectrum. No default source is needed in random ray mode.
diff --git a/src/source.cpp b/src/source.cpp
index e4303c1fa5a..d7bd021da20 100644
--- a/src/source.cpp
+++ b/src/source.cpp
@@ -46,7 +46,7 @@ namespace model {
 vector<unique_ptr<Source>> external_sources;
 
 // DiscreteIndex class for the alias sampling
-DiscreteIndex external_sources_alias_sampler;
+DiscreteIndex external_sources_probability;
 
 } // namespace model
 
@@ -618,7 +618,7 @@ SourceSite sample_external_source(uint64_t* seed)
     if (settings::uniform_source_sampling) {
       i = prn(seed) * model::external_sources.size();
     } else {
-      i = model::external_sources_alias_sampler.sample(seed);
+      i = model::external_sources_probability.sample(seed);
     }
   }
 

From 8dba43e74b3d2144e82984212ab0f28ddc29224d Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Fri, 22 Nov 2024 19:48:43 -0600
Subject: [PATCH 23/24] Two comments

---
 src/settings.cpp | 1 +
 src/source.cpp   | 1 -
 2 files changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/settings.cpp b/src/settings.cpp
index eda0f8ee71c..61eda79967a 100644
--- a/src/settings.cpp
+++ b/src/settings.cpp
@@ -563,6 +563,7 @@ void read_settings_xml(pugi::xml_node root)
     model::external_sources.push_back(make_unique<FileSource>(path));
   }
 
+  // Build probability mass function for sampling external sources
   vector<double> source_strengths;
   for (auto& s : model::external_sources) {
     source_strengths.push_back(s->strength());
diff --git a/src/source.cpp b/src/source.cpp
index d7bd021da20..2116809f2dd 100644
--- a/src/source.cpp
+++ b/src/source.cpp
@@ -45,7 +45,6 @@ namespace model {
 
 vector<unique_ptr<Source>> external_sources;
 
-// DiscreteIndex class for the alias sampling
 DiscreteIndex external_sources_probability;
 
 } // namespace model

From 2fdf150712e3aaed8e1a760e74fe2e7330a2ba94 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Sat, 23 Nov 2024 11:01:58 -0600
Subject: [PATCH 24/24] Update result for source regression test

---
 tests/regression_tests/source/results_true.dat | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/regression_tests/source/results_true.dat b/tests/regression_tests/source/results_true.dat
index 94d9ced1551..673d27c8af8 100644
--- a/tests/regression_tests/source/results_true.dat
+++ b/tests/regression_tests/source/results_true.dat
@@ -1,2 +1,2 @@
 k-combined:
-3.054101E-01 1.167865E-03
+2.959436E-01 2.782384E-03