diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst
index d2c15c56fd..a5cf36ece7 100644
--- a/docs/source/io_formats/settings.rst
+++ b/docs/source/io_formats/settings.rst
@@ -490,9 +490,10 @@ attributes/sub-elements:
       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-, theta-, and phi-coordinates where theta
-      is the angle with respect to the z-axis, phi is the azimuthal angle, and
-      the sphere is centered on the coordinate (x0,y0,z0).
+      independent distributions of r-, cos_theta-, and phi-coordinates where
+      cos_theta is the cosinus 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).
 
       *Default*: None
 
diff --git a/openmc/stats/multivariate.py b/openmc/stats/multivariate.py
index 8635bcd71c..c1de22bb6f 100644
--- a/openmc/stats/multivariate.py
+++ b/openmc/stats/multivariate.py
@@ -399,9 +399,9 @@ class SphericalIndependent(Spatial):
     ----------
     r : openmc.stats.Univariate
         Distribution of r-coordinates in the local reference frame
-    theta : openmc.stats.Univariate
-        Distribution of theta-coordinates (angle relative to the z-axis) in the
-        local reference frame
+    cos_theta : openmc.stats.Univariate
+        Distribution of the cosinus of the theta-coordinates (angle relative to
+        the z-axis) in the local reference frame
     phi : openmc.stats.Univariate
         Distribution of phi-coordinates (azimuthal angle) in the local
         reference frame
@@ -411,9 +411,9 @@ class SphericalIndependent(Spatial):
 
     """
 
-    def __init__(self, r, theta, phi, origin=(0.0, 0.0, 0.0)):
+    def __init__(self, r, cos_theta, phi, origin=(0.0, 0.0, 0.0)):
         self.r = r
-        self.theta = theta
+        self.cos_theta = cos_theta
         self.phi = phi
         self.origin = origin
 
@@ -422,8 +422,8 @@ class SphericalIndependent(Spatial):
         return self._r
 
     @property
-    def theta(self):
-        return self._theta
+    def cos_theta(self):
+        return self._cos_theta
 
     @property
     def phi(self):
@@ -438,10 +438,10 @@ class SphericalIndependent(Spatial):
         cv.check_type('r coordinate', r, Univariate)
         self._r = r
 
-    @theta.setter
-    def theta(self, theta):
-        cv.check_type('theta coordinate', theta, Univariate)
-        self._theta = theta
+    @cos_theta.setter
+    def cos_theta(self, cos_theta):
+        cv.check_type('cos_theta coordinate', cos_theta, Univariate)
+        self._cos_theta = cos_theta
 
     @phi.setter
     def phi(self, phi):
@@ -466,7 +466,7 @@ class SphericalIndependent(Spatial):
         element = ET.Element('space')
         element.set('type', 'spherical')
         element.append(self.r.to_xml_element('r'))
-        element.append(self.theta.to_xml_element('theta'))
+        element.append(self.cos_theta.to_xml_element('cos_theta'))
         element.append(self.phi.to_xml_element('phi'))
         element.set("origin", ' '.join(map(str, self.origin)))
         return element
diff --git a/src/distribution_spatial.cpp b/src/distribution_spatial.cpp
index 6a15992e10..2542842c32 100644
--- a/src/distribution_spatial.cpp
+++ b/src/distribution_spatial.cpp
@@ -132,15 +132,15 @@ SphericalIndependent::SphericalIndependent(pugi::xml_node node)
     r_ = make_unique<Discrete>(x, p, 1);
   }
 
-  // Read distribution for theta-coordinate
-  if (check_for_node(node, "theta")) {
-    pugi::xml_node node_dist = node.child("theta");
-    theta_ = distribution_from_xml(node_dist);
+  // Read distribution for cos_theta-coordinate
+  if (check_for_node(node, "cos_theta")) {
+    pugi::xml_node node_dist = node.child("cos_theta");
+    cos_theta_ = distribution_from_xml(node_dist);
   } else {
-    // If no distribution was specified, default to a single point at theta=0
+    // If no distribution was specified, default to a single point at cos_theta=0
     double x[] {0.0};
     double p[] {1.0};
-    theta_ = make_unique<Discrete>(x, p, 1);
+    cos_theta_ = make_unique<Discrete>(x, p, 1);
   }
 
   // Read distribution for phi-coordinate
@@ -171,15 +171,12 @@ SphericalIndependent::SphericalIndependent(pugi::xml_node node)
 Position SphericalIndependent::sample(uint64_t* seed) const
 {
   double r = r_->sample(seed);
-  double theta = theta_->sample(seed);
-  // To distribute the points uniformly in the sphere, we must ensure that
-  // cos(theta) and not theta is uniformly distributed.
-  // This is done with theta_cos.
-  double theta_cos = std::acos((theta / (0.5 * PI)) - 1);
+  double cos_theta = cos_theta_->sample(seed);
   double phi = phi_->sample(seed);
-  double x = r * sin(theta_cos) * cos(phi) + origin_.x;
-  double y = r * sin(theta_cos) * sin(phi) + origin_.y;
-  double z = r * cos(theta_cos) + origin_.z;
+  // sin(theta) by sin**2 + cos**2 = 1
+  double x = r * std::pow(1 - std::pow(cos_theta, 2.0), 0.5) * cos(phi) + origin_.x;
+  double y = r * std::pow(1 - std::pow(cos_theta, 2.0), 0.5) * sin(phi) + origin_.y;
+  double z = r * cos_theta + origin_.z;
   return {x, y, z};
 }