Simplify angle_between + refactor quadrilateral_area

regional_mom6/regional_mom6.py

@@ -261,7 +261,6 @@ def dz(npoints, ratio, target_depth, min_dz=0.0001, tolerance=1):
 
     Returns:
         numpy.array: An array containing the thickness profile.
-
     """
 
     profile = min_dz + 0.5 * (np.abs(ratio) * min_dz - min_dz) * (
@@ -279,24 +278,16 @@ def dz(npoints, ratio, target_depth, min_dz=0.0001, tolerance=1):
     return dz(npoints, ratio, target_depth, min_dz * err_ratio)
 
 
-# Borrowed from grid tools (GFDL)
 def angle_between(v1, v2, v3):
-    """Returns angle v2-v1-v3 i.e betweeen v1-v2 and v1-v3."""
-
-    # vector product between v1 and v2
-    px = v1[1] * v2[2] - v1[2] * v2[1]
-    py = v1[2] * v2[0] - v1[0] * v2[2]
-    pz = v1[0] * v2[1] - v1[1] * v2[0]
-    # vector product between v1 and v3
-    qx = v1[1] * v3[2] - v1[2] * v3[1]
-    qy = v1[2] * v3[0] - v1[0] * v3[2]
-    qz = v1[0] * v3[1] - v1[1] * v3[0]
-
-    ddd = (px * px + py * py + pz * pz) * (qx * qx + qy * qy + qz * qz)
-    ddd = (px * qx + py * qy + pz * qz) / np.sqrt(ddd)
-    angle = np.arccos(ddd)
-
-    return angle
+    """Returns angle v2-v1-v3 that is between the vectors v1-v2 and v1-v3."""
+
+    v1xv2 = np.cross(v1, v2)
+    v1xv3 = np.cross(v1, v3)
+    cosangle = np.dot(v1xv2, v1xv3) / np.sqrt(
+        np.dot(v1xv2, v1xv2) * np.dot(v1xv3, v1xv3)
+    )
+
+    return np.arccos(cosangle)
 
 
 # Borrowed from grid tools (GFDL)

tests/test_grid_generation.py

@@ -6,6 +6,14 @@
 import xarray as xr
 
 
+def random_unit_vector():
+    """Return a random unit vector on the unit sphere."""
+    z = 2 * np.random.rand(1)[0] - 1
+    z
+    θ = 2 * np.pi * np.random.rand(1)[0]
+    return [np.sqrt(1 - z**2) * np.cos(θ), np.sqrt(1 - z**2) * np.sin(θ), z]
+
+
 @pytest.mark.parametrize(
     ("v1", "v2", "v3", "true_angle"),
     [