diff --git a/src/sphstat/descriptives.py b/src/sphstat/descriptives.py
index 4db30eb..1b8fa2b 100644
--- a/src/sphstat/descriptives.py
+++ b/src/sphstat/descriptives.py
@@ -128,7 +128,7 @@ def rotationmatrix_withaxis(ua: np.array, ub: np.array) -> np.ndarray:
     ca /= np.linalg.norm(ca)
     B = ub.reshape((3, 1)) @ ca.reshape((1, 3))
     B -= B.T
-    th = np.arccos(ab)
+    th = np.arccos(np.clip(ab, -1.0, 1.0))
     A = np.eye(3) + np.sin(th) * B + (np.cos(th) - 1) * ((ub.reshape((3, 1)) @ ub.reshape((1, 3))) + ca.reshape((3, 1))
                                                          @ ca.reshape((1, 3)))
     return A
@@ -381,7 +381,7 @@ def errmin(thph, sample):
         pts = sample['points']
         err = 0
         for pt in pts:
-            err += np.arccos(np.dot(pt, x))
+            err += np.arccos(np.clip(np.dot(pt, x), -1.0, 1.0))
         return err
     # We will use Nelder-Mead to minimise the arc lengths
     res = minimize(errmin, np.array([ms[0], ms[1]]), args = (sample), method='Nelder-Mead', tol=1e-6)
diff --git a/src/sphstat/distributions.py b/src/sphstat/distributions.py
index 904e01b..4dd8f37 100644
--- a/src/sphstat/distributions.py
+++ b/src/sphstat/distributions.py
@@ -289,7 +289,7 @@ def watson(numsamp: int, lamb: float, mu: float, nu: float, kappa: float) -> dic
                 U = rng.uniform(0, 1, 1)
                 V = rng.uniform(0, 1, 1)
                 S = 1 / kappa * np.log(U / C + 1)
-            the = np.arccos(S)
+            the = np.arccos(np.clip(S, -1.0, 1.0))
             phi = 2 * np.pi * rng.uniform(0, 1, 1)
             sample['tetas'].append(the)
             sample['phis'].append(phi)
@@ -310,7 +310,7 @@ def watson(numsamp: int, lamb: float, mu: float, nu: float, kappa: float) -> dic
                 U = rng.uniform(0, 1, 1)
                 V = rng.uniform(0, 1, 1)
                 S = (1 / c1) * np.tan(c2 * U)
-            the = np.arccos(S)
+            the = np.arccos(np.clip(S, -1.0, 1.0))
             phi = 2 * np.pi * rng.uniform(0, 1, 1)
             sample['tetas'].append(the)
             sample['phis'].append(phi)
diff --git a/src/sphstat/twosample.py b/src/sphstat/twosample.py
index 2e92ea6..80e0542 100644
--- a/src/sphstat/twosample.py
+++ b/src/sphstat/twosample.py
@@ -898,5 +898,5 @@ def errors(samplecart: dict, srcpos: tuple) -> list:
     errs = []
     srcvec = sph2cart(srcpos[0], srcpos[1])
     for pt in samplecart['points']:
-        errs.append(np.arccos(np.dot(pt, srcvec)))
+        errs.append(np.arccos(np.clip(np.dot(pt, srcvec), -1.0, 1.0)))
     return errs
diff --git a/src/sphstat/utils.py b/src/sphstat/utils.py
index 98ffa89..2d25529 100644
--- a/src/sphstat/utils.py
+++ b/src/sphstat/utils.py
@@ -258,7 +258,7 @@ def angle(pt1: np.ndarray, pt2: np.ndarray) -> float:
     assert type(pt2) == np.ndarray
     assert np.isclose(np.linalg.norm(pt1), 1)
     assert np.isclose(np.linalg.norm(pt2), 1)
-    return np.arccos(np.dot(pt1, pt2))
+    return np.arccos(np.clip(np.dot(pt1, pt2), -1.0, 1.0))
 
 
 def cart2sph(pt: np.array) -> tuple:
@@ -270,7 +270,7 @@ def cart2sph(pt: np.array) -> tuple:
     :rtype: tuple
     """
     assert np.isclose(np.linalg.norm(pt), 1)
-    th = np.arccos(pt[2])
+    th = np.arccos(np.clip(pt[2], -1.0, 1.0))
     if pt[1] == 0:
         ph = 0
         return th, ph