diff --git a/openmc/lib/math.py b/openmc/lib/math.py
index 781402c46..64225868d 100644
--- a/openmc/lib/math.py
+++ b/openmc/lib/math.py
@@ -31,7 +31,7 @@ _dll.maxwell_spectrum.restype = c_double
 _dll.maxwell_spectrum.argtypes = [c_double, POINTER(c_uint64), c_int]
 
 _dll.watt_spectrum.restype = c_double
-_dll.watt_spectrum.argtypes = [c_double, c_double]
+_dll.watt_spectrum.argtypes = [c_double, c_double, POINTER(c_uint64), c_int]
 
 _dll.broaden_wmp_polynomials.restype = None
 _dll.broaden_wmp_polynomials.argtypes = [c_double, c_double, c_int,
@@ -144,7 +144,7 @@ def calc_zn(n, rho, phi):
     rho : float
         Radial location in the unit disk
     phi : float
-        Theta (radians) location in the unit disk
+		Theta (radians) location in the unit disk
 
     Returns
     -------
@@ -185,7 +185,7 @@ def calc_zn_rad(n, rho):
     return zn_rad
 
 
-def rotate_angle(uvw0, mu, phi=None):
+def rotate_angle(uvw0, mu, phi, prn_seeds, stream ):
     """ Rotates direction cosines through a polar angle whose cosine is
     mu and through an azimuthal angle sampled uniformly.
 
@@ -195,8 +195,12 @@ def rotate_angle(uvw0, mu, phi=None):
         Original direction cosine
     mu : float
         Polar angle cosine to rotate
-    phi : float, optional
-        Azimuthal angle; if None, one will be sampled uniformly
+    phi : iterable of float
+        Azimuthal angle; if NULL, one will be sampled uniformly
+    prn_seeds : iterable of int
+        PRNG seed array
+    stream : int
+        PRNG stream index
 
     Returns
     -------
@@ -206,17 +210,18 @@ def rotate_angle(uvw0, mu, phi=None):
     """
 
     uvw0_arr = np.array(uvw0, dtype=np.float64)
+    phi_arr  = np.array(phi, dtype=np.float64)
+    if phi_arr.size == 0:
+        phi_arr = None
+    prn_seeds_arr = np.array(prn_seeds, dtype=np.uint64)
 
-    if phi is None:
-        _dll.rotate_angle_c(uvw0_arr, mu, None)
-    else:
-        _dll.rotate_angle_c(uvw0_arr, mu, c_double(phi))
+    dll.rotate_angle_c(uvw0_arr, mu, phi_arr, prn_seeds_arr, stream)
     uvw = uvw0_arr
 
     return uvw
 
 
-def maxwell_spectrum(T):
+def maxwell_spectrum(T, prn_seeds, stream):
     """ Samples an energy from the Maxwell fission distribution based
     on a direct sampling scheme.
 
@@ -224,6 +229,10 @@ def maxwell_spectrum(T):
     ----------
     T : float
         Spectrum parameter
+    prn_seeds : iterable of int
+        PRNG seed array
+    stream : int
+        PRNG stream index
 
     Returns
     -------
@@ -231,11 +240,13 @@ def maxwell_spectrum(T):
         Sampled outgoing energy
 
     """
+	
+    prn_seeds_arr = np.array(prn_seeds, dtype=np.uint64)
 
-    return _dll.maxwell_spectrum(T)
+    return _dll.maxwell_spectrum(T, prn_seeds_arr, stream)
 
 
-def watt_spectrum(a, b):
+def watt_spectrum(a, b, prn_seeds, stream):
     """ Samples an energy from the Watt energy-dependent fission spectrum.
 
     Parameters
@@ -244,6 +255,10 @@ def watt_spectrum(a, b):
         Spectrum parameter a
     b : float
         Spectrum parameter b
+    prn_seeds : iterable of int
+        PRNG seed array
+    stream : int
+        PRNG stream index
 
     Returns
     -------
@@ -251,11 +266,13 @@ def watt_spectrum(a, b):
         Sampled outgoing energy
 
     """
+	
+    prn_seeds_arr = np.array(prn_seeds, dtype=np.uint64)
 
-    return _dll.watt_spectrum(a, b)
+    return _dll.watt_spectrum(a, b, prn_seeds_arr, stream)
 
 
-def normal_variate(mean_value, std_dev):
+def normal_variate(mean_value, std_dev, prn_seeds, stream):
     """ Samples an energy from the Normal distribution.
 
     Parameters
@@ -264,6 +281,10 @@ def normal_variate(mean_value, std_dev):
         Mean of the Normal distribution
     std_dev : float
         Standard deviation of the normal distribution
+    prn_seeds : iterable of int
+        PRNG seed array
+    stream : int
+        PRNG stream index
 
     Returns
     -------
@@ -271,6 +292,8 @@ def normal_variate(mean_value, std_dev):
         Sampled outgoing normally distributed value
 
     """
+	
+    prn_seeds_arr = np.array(prn_seeds, dtype=np.uint64)
 
     return _dll.normal_variate(mean_value, std_dev)
 
diff --git a/tests/unit_tests/test_math.py b/tests/unit_tests/test_math.py
index 359181d23..60c872431 100644
--- a/tests/unit_tests/test_math.py
+++ b/tests/unit_tests/test_math.py
@@ -154,9 +154,9 @@ def test_calc_zn_rad():
 
 def test_rotate_angle():
     uvw0 = np.array([1., 0., 0.])
-    phi = 0.
+    phi = [0.]
     mu = 0.
-    prn_seeds = np.array([1, 2, 3, 4, 5, 6], dtype=np.uint64)
+	prn_seeds = [1, 2, 3, 4, 5, 6]
     stream = 0
 
     # reference: mu of 0 pulls the vector the bottom, so:
@@ -176,6 +176,7 @@ def test_rotate_angle():
 
     # Now to test phi is None
     mu = 0.9
+	phi = []
 
     # When seed = 1, phi will be sampled as 1.9116495709698769
     # The resultant reference is from hand-calculations given the above
@@ -186,7 +187,7 @@ def test_rotate_angle():
 
 
 def test_maxwell_spectrum():
-    prn_seeds = np.array([1, 2, 3, 4, 5, 6], dtype=np.uint64)
+	prn_seeds = [1, 2, 3, 4, 5, 6]
     stream = 0
     T = 0.5
     ref_val = 0.6129982175261098
@@ -196,7 +197,7 @@ def test_maxwell_spectrum():
 
 
 def test_watt_spectrum():
-    prn_seeds = np.array([1, 2, 3, 4, 5, 6], dtype=np.uint64)
+	prn_seeds = [1, 2, 3, 4, 5, 6]
     stream = 0
     a = 0.5
     b = 0.75
@@ -207,7 +208,7 @@ def test_watt_spectrum():
 
 
 def test_normal_dist():
-    prn_seeds = np.array([1, 2, 3, 4, 5, 6], dtype=np.uint64)
+	prn_seeds = [1, 2, 3, 4, 5, 6]
     stream = 0
     a = 14.08
     b = 0.0
@@ -216,7 +217,7 @@ def test_normal_dist():
 
     assert ref_val == pytest.approx(test_val)
 
-    prn_seeds = np.array([1, 2, 3, 4, 5, 6], dtype=np.uint64)
+	prn_seeds = [1, 2, 3, 4, 5, 6]
     stream = 0
     a = 14.08
     b = 1.0