From 783062a286f667f16a524a36c7f84b06349d96fe Mon Sep 17 00:00:00 2001
From: Jonathan Shimwell <mail@jshimwell.com>
Date: Wed, 31 Jan 2024 14:02:17 +0000
Subject: [PATCH 1/9] draft source plotting

---
 openmc/model/model.py | 131 ++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 131 insertions(+)

diff --git a/openmc/model/model.py b/openmc/model/model.py
index 5194fec1daa..92c3963d77b 100644
--- a/openmc/model/model.py
+++ b/openmc/model/model.py
@@ -11,6 +11,7 @@
 
 import h5py
 import lxml.etree as ET
+import numpy as np
 
 import openmc
 import openmc._xml as xml
@@ -802,6 +803,133 @@ def calculate_volumes(self, threads=None, output=True, cwd='.',
                             openmc.lib.materials[domain_id].volume = \
                                 vol_calc.volumes[domain_id].n
 
+    def plot(
+            self,
+            source: dict = {'n_samples': 1, 'plane_tolerance': 1},
+            origin: Optional[Iterable[float]] = None,
+            basis: str = 'xy',
+            *args,
+            **kwargs
+        ):
+        """Display a slice plot of the geometry.
+
+        .. versionadded:: 0.14.1
+
+        Parameters
+        ----------
+        origin : iterable of float
+            Coordinates at the origin of the plot. If left as None then the
+            bounding box center will be used to attempt to ascertain the origin.
+            Defaults to (0, 0, 0) if the bounding box is not finite
+        width : iterable of float
+            Width of the plot in each basis direction. If left as none then the
+            bounding box width will be used to attempt to ascertain the plot
+            width. Defaults to (10, 10) if the bounding box is not finite
+        pixels : Iterable of int or int
+            If iterable of ints provided, then this directly sets the number of
+            pixels to use in each basis direction. If int provided, then this
+            sets the total number of pixels in the plot and the number of pixels
+            in each basis direction is calculated from this total and the image
+            aspect ratio.
+        basis : {'xy', 'xz', 'yz'}
+            The basis directions for the plot
+        color_by : {'cell', 'material'}
+            Indicate whether the plot should be colored by cell or by material
+        colors : dict
+            Assigns colors to specific materials or cells. Keys are instances of
+            :class:`Cell` or :class:`Material` and values are RGB 3-tuples, RGBA
+            4-tuples, or strings indicating SVG color names. Red, green, blue,
+            and alpha should all be floats in the range [0.0, 1.0], for
+            example::
+
+               # Make water blue
+               water = openmc.Cell(fill=h2o)
+               universe.plot(..., colors={water: (0., 0., 1.))
+        seed : int
+            Seed for the random number generator
+        openmc_exec : str
+            Path to OpenMC executable.
+        axes : matplotlib.Axes
+            Axes to draw to
+        legend : bool
+            Whether a legend showing material or cell names should be drawn
+        legend_kwargs : dict
+            Keyword arguments passed to :func:`matplotlib.pyplot.legend`.
+        outline : bool
+            Whether outlines between color boundaries should be drawn
+        axis_units : {'km', 'm', 'cm', 'mm'}
+            Units used on the plot axis
+        source : dict
+            A dictionary containing details for plotting the source locations.
+            Acceptable keys are n_samples and  which is the number of source
+            particles to sample and add to plot. Defaults
+            to None which doesn't plot any particles on the plot. Note particles
+            are drawn
+        **kwargs
+            Keyword arguments passed to :func:`matplotlib.pyplot.imshow`
+        Returns
+        -------
+        matplotlib.axes.Axes
+            Axes containing resulting image
+        """
+
+        plot = self.geometry.plot(*args, **kwargs)
+        if source:
+            particles = self.sample_external_source(source['n_samples'])
+
+            particle_indices = {'xy': (0, 1), 'xz': (0, 2), 'yz': (1, 2)}[basis]
+
+            if np.isinf(self.geometry.bounding_box.extent[basis]).any():
+                if origin is None:
+                    origin = (0, 0, 0)
+            else:
+                if origin is None:
+                    # if nan values in the bb.center they get replaced with 0.0
+                    # this happens when the bounding_box contains inf values
+                    with warnings.catch_warnings():
+                        warnings.simplefilter("ignore", RuntimeWarning)
+                        origin = np.nan_to_num(self.geometry.bounding_box.center)
+
+            slice_axis_index = {'xy': 2, 'xz': 1, 'yz': 0}[basis]
+            slice_axis_value = origin[slice_axis_index]
+
+            x_positions = []
+            y_positions = []
+            for particle in particles:
+                if particle.r[slice_axis_index] < slice_axis_value + source['plane_tolerance']:# and particle.r[slice_axis_index] > slice_axis_value - source['plane_tolerance']:
+                    x_positions.append(particle.r[particle_indices[0]])
+                    y_positions.append(particle.r[particle_indices[1]])
+            plot.plot(x_positions, y_positions, marker='x', color='red', label='neutron source')
+        return plot
+
+    def sample_external_source(self, n_samples: int=1, prn_seed: Optional[int]=None):
+        """Sample external source
+
+        .. versionadded:: 0.14.1
+
+        Parameters
+        ----------
+        n_samples : int
+            Number of samples
+        prn_seed : int
+            Pseudorandom number generator (PRNG) seed; if None, one will be
+            generated randomly.
+
+        Returns
+        -------
+        list of openmc.SourceParticle
+            List of samples source particles
+        """
+        
+        self.export_to_xml()
+        import openmc.lib
+        with openmc.lib.run_in_memory(output=False):
+            particles = openmc.lib.sample_external_source(
+                n_samples=n_samples, prn_seed=prn_seed
+            )
+
+        return particles
+
     def plot_geometry(self, output=True, cwd='.', openmc_exec='openmc'):
         """Creates plot images as specified by the Model.plots attribute
 
@@ -820,6 +948,9 @@ def plot_geometry(self, output=True, cwd='.', openmc_exec='openmc'):
 
         """
 
+        # warnings.warn("plot_geometry method will be deprecated soon, please "
+        #               "make use of the plot method instead.")
+
         if len(self.plots) == 0:
             # Then there is no volume calculation specified
             raise ValueError("The Model.plots attribute must be specified "

From f88e2165885cd1ca22cc25f7f389eb0ae012e26e Mon Sep 17 00:00:00 2001
From: Jonathan Shimwell <mail@jshimwell.com>
Date: Wed, 31 Jan 2024 14:35:29 +0000
Subject: [PATCH 2/9] scatter not plot

---
 openmc/model/model.py | 23 ++++++++++++++---------
 1 file changed, 14 insertions(+), 9 deletions(-)

diff --git a/openmc/model/model.py b/openmc/model/model.py
index 92c3963d77b..8664e8ef97e 100644
--- a/openmc/model/model.py
+++ b/openmc/model/model.py
@@ -804,13 +804,13 @@ def calculate_volumes(self, threads=None, output=True, cwd='.',
                                 vol_calc.volumes[domain_id].n
 
     def plot(
-            self,
-            source: dict = {'n_samples': 1, 'plane_tolerance': 1},
-            origin: Optional[Iterable[float]] = None,
-            basis: str = 'xy',
-            *args,
-            **kwargs
-        ):
+        self,
+        source: dict = None,
+        origin: Optional[Iterable[float]] = None,
+        basis: str = "xy",
+        *args,
+        **kwargs,
+    ):
         """Display a slice plot of the geometry.
 
         .. versionadded:: 0.14.1
@@ -896,10 +896,15 @@ def plot(
             x_positions = []
             y_positions = []
             for particle in particles:
-                if particle.r[slice_axis_index] < slice_axis_value + source['plane_tolerance']:# and particle.r[slice_axis_index] > slice_axis_value - source['plane_tolerance']:
+                if (
+                    particle.r[slice_axis_index]
+                    < slice_axis_value + source["plane_tolerance"]
+                    and particle.r[slice_axis_index]
+                    > slice_axis_value - source["plane_tolerance"]
+                ):
                     x_positions.append(particle.r[particle_indices[0]])
                     y_positions.append(particle.r[particle_indices[1]])
-            plot.plot(x_positions, y_positions, marker='x', color='red', label='neutron source')
+            plot.scatter(x_positions, y_positions, marker="x", label="source")
         return plot
 
     def sample_external_source(self, n_samples: int=1, prn_seed: Optional[int]=None):

From baa3326c5bb8be2311d788b4e5e0f500c7f7ad1e Mon Sep 17 00:00:00 2001
From: Jonathan Shimwell <mail@jshimwell.com>
Date: Wed, 31 Jan 2024 15:26:52 +0000
Subject: [PATCH 3/9] remove warning

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

diff --git a/openmc/model/model.py b/openmc/model/model.py
index 8664e8ef97e..4e9a65e9116 100644
--- a/openmc/model/model.py
+++ b/openmc/model/model.py
@@ -873,7 +873,7 @@ def plot(
             Axes containing resulting image
         """
 
-        plot = self.geometry.plot(*args, **kwargs)
+        plot = self.geometry.plot(basis=basis, origin=origin, *args, **kwargs)
         if source:
             particles = self.sample_external_source(source['n_samples'])
 
@@ -953,9 +953,6 @@ def plot_geometry(self, output=True, cwd='.', openmc_exec='openmc'):
 
         """
 
-        # warnings.warn("plot_geometry method will be deprecated soon, please "
-        #               "make use of the plot method instead.")
-
         if len(self.plots) == 0:
             # Then there is no volume calculation specified
             raise ValueError("The Model.plots attribute must be specified "

From b53a8300031dd7901d0fae655b4e2bb565350c9e Mon Sep 17 00:00:00 2001
From: Jonathan Shimwell <mail@jshimwell.com>
Date: Wed, 31 Jan 2024 15:35:47 +0000
Subject: [PATCH 4/9] added kwarg for source scatter plot

---
 openmc/model/model.py | 8 +++++++-
 1 file changed, 7 insertions(+), 1 deletion(-)

diff --git a/openmc/model/model.py b/openmc/model/model.py
index 4e9a65e9116..6fbac602d36 100644
--- a/openmc/model/model.py
+++ b/openmc/model/model.py
@@ -803,11 +803,15 @@ def calculate_volumes(self, threads=None, output=True, cwd='.',
                             openmc.lib.materials[domain_id].volume = \
                                 vol_calc.volumes[domain_id].n
 
+    # default kwargs that are passed to plt.scatter in the plot method below.
+    _default_source_kwargs = {'marker': 'x'}
+
     def plot(
         self,
         source: dict = None,
         origin: Optional[Iterable[float]] = None,
         basis: str = "xy",
+        source_kwargs: dict = _default_source_kwargs,
         *args,
         **kwargs,
     ):
@@ -855,6 +859,8 @@ def plot(
             Whether a legend showing material or cell names should be drawn
         legend_kwargs : dict
             Keyword arguments passed to :func:`matplotlib.pyplot.legend`.
+        source_kwargs : dict
+            Keyword arguments passed to :func:`matplotlib.pyplot.scatter`.
         outline : bool
             Whether outlines between color boundaries should be drawn
         axis_units : {'km', 'm', 'cm', 'mm'}
@@ -904,7 +910,7 @@ def plot(
                 ):
                     x_positions.append(particle.r[particle_indices[0]])
                     y_positions.append(particle.r[particle_indices[1]])
-            plot.scatter(x_positions, y_positions, marker="x", label="source")
+            plot.scatter(x_positions, y_positions, **source_kwargs)
         return plot
 
     def sample_external_source(self, n_samples: int=1, prn_seed: Optional[int]=None):

From ae0f49223726d597474735929d66b47b94cd46a3 Mon Sep 17 00:00:00 2001
From: shimwell <mail@jshimwell.com>
Date: Wed, 31 Jan 2024 21:06:50 +0000
Subject: [PATCH 5/9] sep args for source

---
 openmc/model/model.py | 29 ++++++++++++++++-------------
 1 file changed, 16 insertions(+), 13 deletions(-)

diff --git a/openmc/model/model.py b/openmc/model/model.py
index 6fbac602d36..98dad33c076 100644
--- a/openmc/model/model.py
+++ b/openmc/model/model.py
@@ -808,9 +808,10 @@ def calculate_volumes(self, threads=None, output=True, cwd='.',
 
     def plot(
         self,
-        source: dict = None,
         origin: Optional[Iterable[float]] = None,
         basis: str = "xy",
+        n_samples: Optional[int] = None,
+        plane_tolerance: float = 1.,
         source_kwargs: dict = _default_source_kwargs,
         *args,
         **kwargs,
@@ -859,20 +860,22 @@ def plot(
             Whether a legend showing material or cell names should be drawn
         legend_kwargs : dict
             Keyword arguments passed to :func:`matplotlib.pyplot.legend`.
-        source_kwargs : dict
-            Keyword arguments passed to :func:`matplotlib.pyplot.scatter`.
         outline : bool
             Whether outlines between color boundaries should be drawn
         axis_units : {'km', 'm', 'cm', 'mm'}
             Units used on the plot axis
-        source : dict
-            A dictionary containing details for plotting the source locations.
-            Acceptable keys are n_samples and  which is the number of source
-            particles to sample and add to plot. Defaults
-            to None which doesn't plot any particles on the plot. Note particles
-            are drawn
+        n_samples : dict
+            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
+            Keyword arguments passed to :func:`matplotlib.pyplot.scatter`.
         **kwargs
             Keyword arguments passed to :func:`matplotlib.pyplot.imshow`
+
         Returns
         -------
         matplotlib.axes.Axes
@@ -880,8 +883,8 @@ def plot(
         """
 
         plot = self.geometry.plot(basis=basis, origin=origin, *args, **kwargs)
-        if source:
-            particles = self.sample_external_source(source['n_samples'])
+        if n_samples:
+            particles = self.sample_external_source(n_samples)
 
             particle_indices = {'xy': (0, 1), 'xz': (0, 2), 'yz': (1, 2)}[basis]
 
@@ -904,9 +907,9 @@ def plot(
             for particle in particles:
                 if (
                     particle.r[slice_axis_index]
-                    < slice_axis_value + source["plane_tolerance"]
+                    < slice_axis_value + plane_tolerance
                     and particle.r[slice_axis_index]
-                    > slice_axis_value - source["plane_tolerance"]
+                    > slice_axis_value - plane_tolerance
                 ):
                     x_positions.append(particle.r[particle_indices[0]])
                     y_positions.append(particle.r[particle_indices[1]])

From 34b605eec33c445e22bc039fda74324a846e563d Mon Sep 17 00:00:00 2001
From: shimwell <mail@jshimwell.com>
Date: Wed, 31 Jan 2024 21:10:03 +0000
Subject: [PATCH 6/9] added check type to args

---
 openmc/model/model.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/openmc/model/model.py b/openmc/model/model.py
index 98dad33c076..43fade6053f 100644
--- a/openmc/model/model.py
+++ b/openmc/model/model.py
@@ -882,6 +882,9 @@ def plot(
             Axes containing resulting image
         """
 
+        check_type('n_samples', n_samples, int)
+        check_type('plane_tolerance', plane_tolerance, float)
+
         plot = self.geometry.plot(basis=basis, origin=origin, *args, **kwargs)
         if n_samples:
             particles = self.sample_external_source(n_samples)

From 57d9a11c550314d7ba60de923ee2071b33195e5d Mon Sep 17 00:00:00 2001
From: Jonathan Shimwell <mail@jshimwell.com>
Date: Thu, 1 Feb 2024 09:07:04 +0000
Subject: [PATCH 7/9] allowing none type

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

diff --git a/openmc/model/model.py b/openmc/model/model.py
index 43fade6053f..b491ab178e1 100644
--- a/openmc/model/model.py
+++ b/openmc/model/model.py
@@ -882,7 +882,7 @@ def plot(
             Axes containing resulting image
         """
 
-        check_type('n_samples', n_samples, int)
+        check_type('n_samples', n_samples, (int, type(None)))
         check_type('plane_tolerance', plane_tolerance, float)
 
         plot = self.geometry.plot(basis=basis, origin=origin, *args, **kwargs)

From 658a3813b55049526eb97e0a6d93bf7c7b67e6c7 Mon Sep 17 00:00:00 2001
From: Jonathan Shimwell <mail@jshimwell.com>
Date: Thu, 1 Feb 2024 09:34:19 +0000
Subject: [PATCH 8/9] added unit test for coord locations

---
 tests/unit_tests/test_model.py | 31 +++++++++++++++++++++++++++++++
 1 file changed, 31 insertions(+)

diff --git a/tests/unit_tests/test_model.py b/tests/unit_tests/test_model.py
index 24f136c8c4d..0336850714c 100644
--- a/tests/unit_tests/test_model.py
+++ b/tests/unit_tests/test_model.py
@@ -585,3 +585,34 @@ def test_single_xml_exec(run_in_tmpdir):
 
     with pytest.raises(RuntimeError, match='input_dir'):
         openmc.run(path_input='input_dir/pincell.xml')
+
+def test_model_plot():
+    # plots the geometry with source location and checks the resulting matplotlib
+    # includes the correct coordinates for the scatter plot for all basis.
+
+    surface = openmc.Sphere(r=600, boundary_type="vacuum")
+    cell = openmc.Cell(region=-surface)
+    geometry = openmc.Geometry([cell])
+    settings = openmc.Settings()
+    settings.particles = 1
+    settings.batches = 1
+    source = openmc.IndependentSource()
+    source.space = openmc.stats.Point((1, 2, 3))
+    settings.source = source
+    model = openmc.Model(geometry, None, settings)
+
+    plot = model.plot(n_samples=1, plane_tolerance=4.0, basis="xy")
+    coords = plot.axes.collections[0].get_offsets().data.flatten()
+    assert (coords == np.array([1.0, 2.0])).all()
+
+    plot = model.plot(n_samples=1, plane_tolerance=4.0, basis="xz")
+    coords = plot.axes.collections[0].get_offsets().data.flatten()
+    assert (coords == np.array([1.0, 3.0])).all()
+
+    plot = model.plot(n_samples=1, plane_tolerance=4.0, basis="yz")
+    coords = plot.axes.collections[0].get_offsets().data.flatten()
+    assert (coords == np.array([2.0, 3.0])).all()
+
+    plot = model.plot(n_samples=1, plane_tolerance=0.1, basis="xy")
+    coords = plot.axes.collections[0].get_offsets().data.flatten()
+    assert (coords == np.array([])).all()

From d9c046c77565a9f7d14561f6600ae720a7a988e2 Mon Sep 17 00:00:00 2001
From: Paul Romano <paul.k.romano@gmail.com>
Date: Thu, 26 Sep 2024 09:12:05 -0500
Subject: [PATCH 9/9] Run Model.sample_external_source in temporary directory

---
 openmc/geometry.py             |   5 +-
 openmc/lib/core.py             |   9 ++-
 openmc/model/model.py          | 139 +++++++++++++--------------------
 openmc/universe.py             |   3 +-
 tests/unit_tests/test_model.py |  15 ++--
 5 files changed, 71 insertions(+), 100 deletions(-)

diff --git a/openmc/geometry.py b/openmc/geometry.py
index 6cce4c18c70..a52dc4418d2 100644
--- a/openmc/geometry.py
+++ b/openmc/geometry.py
@@ -7,8 +7,6 @@
 import warnings
 import lxml.etree as ET
 
-import numpy as np
-
 import openmc
 import openmc._xml as xml
 from .checkvalue import check_type, check_less_than, check_greater_than, PathLike
@@ -67,7 +65,7 @@ def root_universe(self, root_universe):
         self._root_universe = root_universe
 
     @property
-    def bounding_box(self) -> np.ndarray:
+    def bounding_box(self) -> openmc.BoundingBox:
         return self.root_universe.bounding_box
 
     @property
@@ -800,6 +798,7 @@ def plot(self, *args, **kwargs):
             Units used on the plot axis
         **kwargs
             Keyword arguments passed to :func:`matplotlib.pyplot.imshow`
+
         Returns
         -------
         matplotlib.axes.Axes
diff --git a/openmc/lib/core.py b/openmc/lib/core.py
index a9a549fa05a..e646a9ae1df 100644
--- a/openmc/lib/core.py
+++ b/openmc/lib/core.py
@@ -474,8 +474,11 @@ def run(output=True):
         _dll.openmc_run()
 
 
-def sample_external_source(n_samples=1, prn_seed=None):
-    """Sample external source
+def sample_external_source(
+        n_samples: int = 1000,
+        prn_seed: int | None = None
+) -> list[openmc.SourceParticle]:
+    """Sample external source and return source particles.
 
     .. versionadded:: 0.13.1
 
@@ -490,7 +493,7 @@ def sample_external_source(n_samples=1, prn_seed=None):
     Returns
     -------
     list of openmc.SourceParticle
-        List of samples source particles
+        List of sampled source particles
 
     """
     if n_samples <= 0:
diff --git a/openmc/model/model.py b/openmc/model/model.py
index 729b42ead99..78f743a03ed 100644
--- a/openmc/model/model.py
+++ b/openmc/model/model.py
@@ -794,67 +794,19 @@ def calculate_volumes(self, threads=None, output=True, cwd='.',
                             openmc.lib.materials[domain_id].volume = \
                                 vol_calc.volumes[domain_id].n
 
-    # default kwargs that are passed to plt.scatter in the plot method below.
-    _default_source_kwargs = {'marker': 'x'}
-
     def plot(
         self,
-        origin: Optional[Iterable[float]] = None,
-        basis: str = "xy",
-        n_samples: Optional[int] = None,
+        n_samples: int | None = None,
         plane_tolerance: float = 1.,
-        source_kwargs: dict = _default_source_kwargs,
-        *args,
+        source_kwargs: dict | None = None,
         **kwargs,
     ):
         """Display a slice plot of the geometry.
 
-        .. versionadded:: 0.14.1
+        .. versionadded:: 0.15.1
 
         Parameters
         ----------
-        origin : iterable of float
-            Coordinates at the origin of the plot. If left as None then the
-            bounding box center will be used to attempt to ascertain the origin.
-            Defaults to (0, 0, 0) if the bounding box is not finite
-        width : iterable of float
-            Width of the plot in each basis direction. If left as none then the
-            bounding box width will be used to attempt to ascertain the plot
-            width. Defaults to (10, 10) if the bounding box is not finite
-        pixels : Iterable of int or int
-            If iterable of ints provided, then this directly sets the number of
-            pixels to use in each basis direction. If int provided, then this
-            sets the total number of pixels in the plot and the number of pixels
-            in each basis direction is calculated from this total and the image
-            aspect ratio.
-        basis : {'xy', 'xz', 'yz'}
-            The basis directions for the plot
-        color_by : {'cell', 'material'}
-            Indicate whether the plot should be colored by cell or by material
-        colors : dict
-            Assigns colors to specific materials or cells. Keys are instances of
-            :class:`Cell` or :class:`Material` and values are RGB 3-tuples, RGBA
-            4-tuples, or strings indicating SVG color names. Red, green, blue,
-            and alpha should all be floats in the range [0.0, 1.0], for
-            example::
-
-               # Make water blue
-               water = openmc.Cell(fill=h2o)
-               universe.plot(..., colors={water: (0., 0., 1.))
-        seed : int
-            Seed for the random number generator
-        openmc_exec : str
-            Path to OpenMC executable.
-        axes : matplotlib.Axes
-            Axes to draw to
-        legend : bool
-            Whether a legend showing material or cell names should be drawn
-        legend_kwargs : dict
-            Keyword arguments passed to :func:`matplotlib.pyplot.legend`.
-        outline : bool
-            Whether outlines between color boundaries should be drawn
-        axis_units : {'km', 'm', 'cm', 'mm'}
-            Units used on the plot axis
         n_samples : dict
             The number of source particles to sample and add to plot. Defaults
             to None which doesn't plot any particles on the plot.
@@ -865,7 +817,7 @@ def plot(
         source_kwargs : dict
             Keyword arguments passed to :func:`matplotlib.pyplot.scatter`.
         **kwargs
-            Keyword arguments passed to :func:`matplotlib.pyplot.imshow`
+            Keyword arguments passed to :func:`openmc.Universe.plot`
 
         Returns
         -------
@@ -873,47 +825,57 @@ def plot(
             Axes containing resulting image
         """
 
-        check_type('n_samples', n_samples, (int, type(None)))
+        check_type('n_samples', n_samples, int | None)
         check_type('plane_tolerance', plane_tolerance, float)
+        if source_kwargs is None:
+            source_kwargs = {}
+        source_kwargs.setdefault('marker', 'x')
 
-        plot = self.geometry.plot(basis=basis, origin=origin, *args, **kwargs)
+        ax = self.geometry.plot(**kwargs)
         if n_samples:
+            # Sample external source particles
             particles = self.sample_external_source(n_samples)
 
-            particle_indices = {'xy': (0, 1), 'xz': (0, 2), 'yz': (1, 2)}[basis]
+            # Determine plotting parameters and bounding box of geometry
+            bbox = self.geometry.bounding_box
+            origin = kwargs.get('origin', None)
+            basis = kwargs.get('basis', 'xy')
+            indices = {'xy': (0, 1, 2), 'xz': (0, 2, 1), 'yz': (1, 2, 0)}[basis]
 
-            if np.isinf(self.geometry.bounding_box.extent[basis]).any():
+            # Infer origin if not provided
+            if np.isinf(bbox.extent[basis]).any():
                 if origin is None:
                     origin = (0, 0, 0)
             else:
                 if origin is None:
-                    # if nan values in the bb.center they get replaced with 0.0
+                    # if nan values in the bbox.center they get replaced with 0.0
                     # this happens when the bounding_box contains inf values
                     with warnings.catch_warnings():
                         warnings.simplefilter("ignore", RuntimeWarning)
-                        origin = np.nan_to_num(self.geometry.bounding_box.center)
+                        origin = np.nan_to_num(bbox.center)
 
-            slice_axis_index = {'xy': 2, 'xz': 1, 'yz': 0}[basis]
-            slice_axis_value = origin[slice_axis_index]
+            slice_index = indices[2]
+            slice_value = origin[slice_index]
 
-            x_positions = []
-            y_positions = []
+            xs = []
+            ys = []
+            tol = plane_tolerance
             for particle in particles:
-                if (
-                    particle.r[slice_axis_index]
-                    < slice_axis_value + plane_tolerance
-                    and particle.r[slice_axis_index]
-                    > slice_axis_value - plane_tolerance
-                ):
-                    x_positions.append(particle.r[particle_indices[0]])
-                    y_positions.append(particle.r[particle_indices[1]])
-            plot.scatter(x_positions, y_positions, **source_kwargs)
-        return plot
-
-    def sample_external_source(self, n_samples: int=1, prn_seed: Optional[int]=None):
-        """Sample external source
-
-        .. versionadded:: 0.14.1
+                if (slice_value - tol < particle.r[slice_index] < slice_value + tol):
+                    xs.append(particle.r[indices[0]])
+                    ys.append(particle.r[indices[1]])
+            ax.scatter(xs, ys, **source_kwargs)
+        return ax
+
+    def sample_external_source(
+            self,
+            n_samples: int = 1000,
+            prn_seed: int | None = None,
+            **init_kwargs
+    ) -> list[openmc.SourceParticle]:
+        """Sample external source and return source particles.
+
+        .. versionadded:: 0.15.1
 
         Parameters
         ----------
@@ -922,21 +884,30 @@ def sample_external_source(self, n_samples: int=1, prn_seed: Optional[int]=None)
         prn_seed : int
             Pseudorandom number generator (PRNG) seed; if None, one will be
             generated randomly.
+        **init_kwargs
+            Keyword arguments passed to :func:`openmc.lib.init`
 
         Returns
         -------
         list of openmc.SourceParticle
             List of samples source particles
         """
-        
-        self.export_to_xml()
         import openmc.lib
-        with openmc.lib.run_in_memory(output=False):
-            particles = openmc.lib.sample_external_source(
-                n_samples=n_samples, prn_seed=prn_seed
-            )
 
-        return particles
+        # Silence output by default. Also set arguments to start in volume
+        # calculation mode to avoid loading cross sections
+        init_kwargs.setdefault('output', False)
+        init_kwargs.setdefault('args', ['-c'])
+
+        with change_directory(tmpdir=True):
+            # Export model within temporary directory
+            self.export_to_model_xml()
+
+            # Sample external source sites
+            with openmc.lib.run_in_memory(**init_kwargs):
+                return openmc.lib.sample_external_source(
+                    n_samples=n_samples, prn_seed=prn_seed
+                )
 
     def plot_geometry(self, output=True, cwd='.', openmc_exec='openmc'):
         """Creates plot images as specified by the Model.plots attribute
diff --git a/openmc/universe.py b/openmc/universe.py
index 9fab9ae51b6..d424c243bd4 100644
--- a/openmc/universe.py
+++ b/openmc/universe.py
@@ -1,3 +1,4 @@
+from __future__ import annotations
 import math
 from abc import ABC, abstractmethod
 from collections.abc import Iterable
@@ -230,7 +231,7 @@ def cells(self):
         return self._cells
 
     @property
-    def bounding_box(self):
+    def bounding_box(self) -> openmc.BoundingBox:
         regions = [c.region for c in self.cells.values()
                    if c.region is not None]
         if regions:
diff --git a/tests/unit_tests/test_model.py b/tests/unit_tests/test_model.py
index ee0fe40a670..4b567c56d62 100644
--- a/tests/unit_tests/test_model.py
+++ b/tests/unit_tests/test_model.py
@@ -594,19 +594,16 @@ def test_single_xml_exec(run_in_tmpdir):
 
 
 def test_model_plot():
-    # plots the geometry with source location and checks the resulting matplotlib
-    # includes the correct coordinates for the scatter plot for all basis.
+    # plots the geometry with source location and checks the resulting
+    # matplotlib includes the correct coordinates for the scatter plot for all
+    # basis.
 
     surface = openmc.Sphere(r=600, boundary_type="vacuum")
     cell = openmc.Cell(region=-surface)
     geometry = openmc.Geometry([cell])
-    settings = openmc.Settings()
-    settings.particles = 1
-    settings.batches = 1
-    source = openmc.IndependentSource()
-    source.space = openmc.stats.Point((1, 2, 3))
-    settings.source = source
-    model = openmc.Model(geometry, None, settings)
+    source = openmc.IndependentSource(space=openmc.stats.Point((1, 2, 3)))
+    settings = openmc.Settings(particles=1, batches=1, source=source)
+    model = openmc.Model(geometry, settings=settings)
 
     plot = model.plot(n_samples=1, plane_tolerance=4.0, basis="xy")
     coords = plot.axes.collections[0].get_offsets().data.flatten()