diff --git a/lts_array/classes/lts_classes.py b/lts_array/classes/lts_classes.py
index 45142c7..7acc818 100644
--- a/lts_array/classes/lts_classes.py
+++ b/lts_array/classes/lts_classes.py
@@ -722,10 +722,12 @@ def post_process(dimension_number, co_array_num, alpha, h, nits, tau, xij, coeff
         # Calculate the sigma_tau value (Szuberla et al. 2006).
         residuals = tau[weights, jj, :] - (xij[weights, :] @ z_final)
         m_w, _ = np.shape(xij[weights, :])
-        with np.errstate(invalid='raise'):
+
+        with np.errstate(invalid="raise"):
             try:
-                sigma_tau[jj] = np.sqrt(tau[weights, jj, :].T @ residuals / (
-                    m_w - dimension_number))[0]
+                sigma_tau_value = tau[weights, jj, :].T @ residuals
+                sigma_tau_value = sigma_tau_value / (m_w - dimension_number)
+                sigma_tau[jj] = np.sqrt(float(np.asarray(sigma_tau_value).squeeze()))
             except FloatingPointError:
                 pass
 
@@ -755,8 +757,8 @@ def post_process(dimension_number, co_array_num, alpha, h, nits, tau, xij, coeff
         # Halving here s.t. +/- value expresses uncertainty bounds.
         # Remove the 1/2's to get full values to express
         # coverage ellipse area.
-        conf_int_baz[jj] = 0.5 * sig_theta
-        conf_int_vel[jj] = 0.5 * np.abs(np.diff(1 / eExtrm[:2]))
+        conf_int_baz[jj] = 0.5 * float(np.asarray(sig_theta).squeeze())
+        conf_int_vel[jj] = 0.5 * float(np.asarray(np.abs(np.diff(1 / eExtrm[:2]))).squeeze())
 
         # Cast weights to int for output
         element_weights[:, jj] = weights * 1
@@ -932,8 +934,13 @@ def solve(self, data):
 
             # Calculate the sigma_tau value (Szuberla et al. 2006).
             residuals = self.tau[:, jj, :] - (self.xij @ z_final)
-            self.sigma_tau[jj] = np.sqrt(self.tau[:, jj, :].T @ residuals / (
-                self.co_array_num - self.dimension_number))[0]
+
+            sigma_tau_value = self.tau[:, jj, :].T @ residuals
+            sigma_tau_value = sigma_tau_value / (
+                    self.co_array_num - self.dimension_number
+            )
+
+            self.sigma_tau[jj] = np.sqrt(float(np.asarray(sigma_tau_value).squeeze()))
 
             # Calculate uncertainties from Szuberla & Olson, 2004
             # Equation 16
@@ -960,8 +967,8 @@ def solve(self, data):
                 # Halving here s.t. +/- value expresses uncertainty bounds.
                 # Remove the 1/2's to get full values to express
                 # coverage ellipse area.
-                self.conf_int_baz[jj] = 0.5 * sig_theta
-                self.conf_int_vel[jj] = 0.5 * np.abs(np.diff(1 / eExtrm[:2]))
+                self.conf_int_baz[jj] = 0.5 * float(np.asarray(sig_theta).squeeze())
+                self.conf_int_vel[jj] = 0.5 * float(np.asarray(np.abs(np.diff(1 / eExtrm[:2]))).squeeze())
 
             except ValueError:
                 self.conf_int_baz[jj] = np.nan