diff --git a/named_arrays/_scalars/scalars.py b/named_arrays/_scalars/scalars.py
index 8625c87..be927d4 100644
--- a/named_arrays/_scalars/scalars.py
+++ b/named_arrays/_scalars/scalars.py
@@ -127,6 +127,12 @@ def __array_matmul__(
         if result is not NotImplemented:
             return result
 
+        if out is None:
+            if np.all(x1 == 0) or np.all(x2 == 0):
+                unit_1 = na.unit(x1, unit_dimensionless=1)
+                unit_2 = na.unit(x2, unit_dimensionless=1)
+                return 0 * unit_1 * unit_2
+
         return np.multiply(
             x1,
             x2,
diff --git a/named_arrays/_scalars/tests/test_scalars.py b/named_arrays/_scalars/tests/test_scalars.py
index 695c65d..cde181c 100644
--- a/named_arrays/_scalars/tests/test_scalars.py
+++ b/named_arrays/_scalars/tests/test_scalars.py
@@ -180,7 +180,7 @@ def test_matmul(
 
             result_expected = np.multiply(array, array_2)
 
-            out = 0 * result
+            out = 0 * array * array_2
             result_out = np.matmul(array, array_2, out=out)
 
             assert np.all(result == result_expected)
diff --git a/named_arrays/_vectors/cartesian/vectors_cartesian_3d.py b/named_arrays/_vectors/cartesian/vectors_cartesian_3d.py
index 3cbf00e..977cce5 100644
--- a/named_arrays/_vectors/cartesian/vectors_cartesian_3d.py
+++ b/named_arrays/_vectors/cartesian/vectors_cartesian_3d.py
@@ -130,7 +130,7 @@ def _sold_angle(
         d3 = (b @ c) * a_
         denomerator = d0 + d1 + d2 + d3
 
-        unit = numerator.unit
+        unit = na.unit(numerator)
 
         if unit is not None:
             numerator = numerator.to(unit).value
diff --git a/named_arrays/_vectors/tests/test_vectors_temporal_spectral.py b/named_arrays/_vectors/tests/test_vectors_temporal_spectral.py
index 76557c4..e101c97 100644
--- a/named_arrays/_vectors/tests/test_vectors_temporal_spectral.py
+++ b/named_arrays/_vectors/tests/test_vectors_temporal_spectral.py
@@ -211,7 +211,7 @@ class TestTemporalSpectralVectorLinearSpace(
                 wavelength=1 * u.nm,
             ),
             pc=na.SpectralMatrixArray(
-                wavelength=na.CartesianNdVectorArray(dict(wavelength=1, x=0, y=0)),
+                wavelength=na.CartesianNdVectorArray(dict(wavelength=1, x=0, y=0.1)),
             ),
             shape_wcs=dict(wavelength=5, x=_num_x, y=_num_y),
         ),