diff --git a/conda_package/tests/test_conversion.py b/conda_package/tests/test_conversion.py
index 3a02eb98e..23377d622 100755
--- a/conda_package/tests/test_conversion.py
+++ b/conda_package/tests/test_conversion.py
@@ -47,6 +47,41 @@ def test_conversion_angle_edge():
     assert np.max(np.abs(angle_diff)) < 1.0e-10
 
 
+def test_conversion_triangle_angle_quality():
+    ds_mesh = xarray.open_dataset(
+        get_test_data_file('mesh.QU.1920km.151026.nc')
+    )
+    ds_mesh = convert(dsIn=ds_mesh)
+
+    edges_on_vertex = ds_mesh.edgesOnVertex.values - 1
+    assert np.all(edges_on_vertex >= 0)
+    dc_edge = ds_mesh.dcEdge.values
+    a_len = dc_edge[edges_on_vertex[:, 0]]
+    b_len = dc_edge[edges_on_vertex[:, 1]]
+    c_len = dc_edge[edges_on_vertex[:, 2]]
+
+    # law of cosines for the angle opposite each side of the dual triangle
+    angle1 = np.arccos(
+        np.clip((b_len**2 + c_len**2 - a_len**2) / (2 * b_len * c_len), -1, 1)
+    )
+    angle2 = np.arccos(
+        np.clip((a_len**2 + c_len**2 - b_len**2) / (2 * a_len * c_len), -1, 1)
+    )
+    angle3 = np.arccos(
+        np.clip((a_len**2 + b_len**2 - c_len**2) / (2 * a_len * b_len), -1, 1)
+    )
+    angles = np.stack([angle1, angle2, angle3], axis=1)
+    assert np.max(np.abs(angles.sum(axis=1) - np.pi)) < 1.0e-10
+
+    min_angle = angles.min(axis=1)
+    max_angle = angles.max(axis=1)
+    quality = ds_mesh.triangleAngleQuality.values
+    assert np.max(np.abs(quality - min_angle / max_angle)) < 1.0e-10
+
+    obtuse = (max_angle > 0.5 * np.pi).astype(int)
+    assert np.array_equal(ds_mesh.obtuseTriangle.values, obtuse)
+
+
 def test_masks_to_int_dataset_copy():
     ds_in = xarray.Dataset(
         data_vars={
@@ -71,3 +106,4 @@ def test_masks_to_int_dataset_copy():
 if __name__ == '__main__':
     test_conversion()
     test_conversion_angle_edge()
+    test_conversion_triangle_angle_quality()
diff --git a/mesh_tools/mesh_conversion_tools/mpas_mesh_converter.cpp b/mesh_tools/mesh_conversion_tools/mpas_mesh_converter.cpp
index 83d052aea..a578eb0b0 100644
--- a/mesh_tools/mesh_conversion_tools/mpas_mesh_converter.cpp
+++ b/mesh_tools/mesh_conversion_tools/mpas_mesh_converter.cpp
@@ -2314,7 +2314,7 @@ int buildMeshQualities(){/*{{{*/
 #endif
 
 					angle1 = acos( max(-1.0, min(1.0, (b_len * b_len + c_len * c_len - a_len * a_len) / (2 * b_len * c_len))));
-					angle2 = acos( max(-1.0, min(1.0, (a_len * a_len + c_len * c_len - b_len * c_len) / (2 * a_len * c_len))));
+					angle2 = acos( max(-1.0, min(1.0, (a_len * a_len + c_len * c_len - b_len * b_len) / (2 * a_len * c_len))));
 					angle3 = acos( max(-1.0, min(1.0, (a_len * a_len + b_len * b_len - c_len * c_len) / (2 * a_len * b_len))));
 
 					minAngle = min(angle1, min(angle2, angle3));
diff --git a/mesh_tools/mesh_conversion_tools_netcdf_c/mpas_mesh_converter.cpp b/mesh_tools/mesh_conversion_tools_netcdf_c/mpas_mesh_converter.cpp
index 6f38e7424..aafeda8c4 100755
--- a/mesh_tools/mesh_conversion_tools_netcdf_c/mpas_mesh_converter.cpp
+++ b/mesh_tools/mesh_conversion_tools_netcdf_c/mpas_mesh_converter.cpp
@@ -2395,7 +2395,7 @@ int buildMeshQualities(){/*{{{*/
                     angle1 = acos( max(-1.0, min(1.0,
                         (b_len * b_len + c_len * c_len - a_len * a_len) / (2 * b_len * c_len))));
                     angle2 = acos( max(-1.0, min(1.0,
-                        (a_len * a_len + c_len * c_len - b_len * c_len) / (2 * a_len * c_len))));
+                        (a_len * a_len + c_len * c_len - b_len * b_len) / (2 * a_len * c_len))));
                     angle3 = acos( max(-1.0, min(1.0,
                         (a_len * a_len + b_len * b_len - c_len * c_len) / (2 * a_len * b_len))));