diff --git a/openmc/lib/core.py b/openmc/lib/core.py
index de5f4adf21..3f4a301cd8 100644
--- a/openmc/lib/core.py
+++ b/openmc/lib/core.py
@@ -1,8 +1,10 @@
 from contextlib import contextmanager
 from ctypes import (c_bool, c_int, c_int32, c_int64, c_double, c_char_p,
-                    c_char, POINTER, Structure, c_void_p, create_string_buffer)
+                    c_char, POINTER, Structure, c_void_p, create_string_buffer,
+                    c_uint64)
 import sys
 import os
+from random import getrandbits
 
 import numpy as np
 from numpy.ctypeslib import as_array
@@ -10,6 +12,7 @@ from numpy.ctypeslib import as_array
 from . import _dll
 from .error import _error_handler
 import openmc.lib
+import openmc
 
 
 class _SourceSite(Structure):
@@ -95,6 +98,9 @@ _dll.openmc_global_bounding_box.argtypes = [POINTER(c_double),
                                             POINTER(c_double)]
 _dll.openmc_global_bounding_box.restype = c_int
 _dll.openmc_global_bounding_box.errcheck = _error_handler
+_dll.openmc_sample_external_source.argtypes = [POINTER(c_uint64), POINTER(_SourceSite)]
+_dll.openmc_sample_external_source.restype = c_int
+_dll.openmc_sample_external_source.errcheck = _error_handler
 
 
 def global_bounding_box():
@@ -415,6 +421,37 @@ def run(output=True):
         _dll.openmc_run()
 
 
+def sample_external_source(prn_seed=None):
+    """Sample external source
+
+    .. versionadded:: 0.13.1
+
+    Parameters
+    ----------
+    prn_seed : int
+        PRNG seed; if None, one will be generated randomly
+
+    Returns
+    -------
+    openmc.SourceParticle
+        Sampled source particle
+
+    """
+    if prn_seed is None:
+        prn_seed = getrandbits(63)
+
+    # Call into C API to sample source
+    site = _SourceSite()
+    _dll.openmc_sample_external_source(c_uint64(prn_seed), site)
+
+    # Convert to SourceParticle and return
+    return openmc.SourceParticle(
+        r=site.r, u=site.u, E=site.E, time=site.time, wgt=site.wgt,
+        delayed_group=site.delayed_group, surf_id=site.surf_id,
+        particle=openmc.ParticleType(site.particle)
+    )
+
+
 def simulation_init():
     """Initialize simulation"""
     _dll.openmc_simulation_init()
diff --git a/openmc/lib/math.py b/openmc/lib/math.py
index bdb09b2d1f..d70c20a9f6 100644
--- a/openmc/lib/math.py
+++ b/openmc/lib/math.py
@@ -1,12 +1,11 @@
 from ctypes import c_int, c_double, POINTER, c_uint64
+from random import getrandbits
 
 import numpy as np
 from numpy.ctypeslib import ndpointer
 
 from . import _dll
 
-from random import getrandbits
-
 
 _dll.t_percentile.restype = c_double
 _dll.t_percentile.argtypes = [c_double, c_int]
@@ -240,10 +239,10 @@ def maxwell_spectrum(T, prn_seed=None):
         Sampled outgoing energy
 
     """
-    
+
     if prn_seed is None:
         prn_seed = getrandbits(63)
-	
+
     return _dll.maxwell_spectrum(T, c_uint64(prn_seed))
 
 
@@ -265,7 +264,7 @@ def watt_spectrum(a, b, prn_seed=None):
         Sampled outgoing energy
 
     """
-    
+
     if prn_seed is None:
         prn_seed = getrandbits(63)
 
@@ -290,7 +289,7 @@ def normal_variate(mean_value, std_dev, prn_seed=None):
         Sampled outgoing normally distributed value
 
     """
-    
+
     if prn_seed is None:
         prn_seed = getrandbits(63)
 
diff --git a/src/source.cpp b/src/source.cpp
index 11ede5fc4e..8b634ec459 100644
--- a/src/source.cpp
+++ b/src/source.cpp
@@ -399,4 +399,14 @@ void free_memory_source()
   model::external_sources.clear();
 }
 
+//==============================================================================
+// C API
+//==============================================================================
+
+extern "C" int openmc_sample_external_source(uint64_t* seed, SourceSite* site)
+{
+  *site = sample_external_source(seed);
+  return 0;
+}
+
 } // namespace openmc