From 87e99ac4ceab6d6daecb600d1c25d208934379e7 Mon Sep 17 00:00:00 2001
From: Leni Aniva <v@leni.sh>
Date: Fri, 28 Jun 2024 21:36:33 -0400
Subject: [PATCH] fix: Collision problem with Hirth joints

---
 nhf/joints.py | 65 ++++++++++++++++++++++++++++++++++++---------------
 nhf/test.py   |  6 ++++-
 2 files changed, 51 insertions(+), 20 deletions(-)

diff --git a/nhf/joints.py b/nhf/joints.py
index 7c227b4..43c5169 100644
--- a/nhf/joints.py
+++ b/nhf/joints.py
@@ -23,50 +23,69 @@ def hirth_joint(radius=60,
     Creates a cylindrical Hirth Joint
 
     is_mated: If set to true, rotate the teeth so they line up at 0 degrees.
+
+    FIXME: The curves don't mate perfectly. See if non-planar lofts can solve
+    this issue.
     """
-    # ensures secant doesn't blow up
+    # ensures tangent doesn't blow up
     assert n_tooth >= 5
     assert radius > radius_inner
+    assert tooth_height >= tooth_height_inner
 
     # angle of half of a single tooth
     theta = math.pi / n_tooth
 
-    # Generate a tooth by lofting between two curves
-
+    c, s, t = math.cos(theta), math.sin(theta), math.tan(theta)
+    span = radius * t
+    radius_proj = radius / c
+    span_inner = radius_inner * s
+    # 2 * raise + (inner tooth height) = (tooth height)
     inner_raise = (tooth_height - tooth_height_inner) / 2
-    # Outer tooth triangle spans a curve of length `2 pi r / n_tooth`. This
-    # creates the side profile (looking radially inwards) of each of the
-    # triangles.
+    # Outer tooth triangle spans 2*theta radians. This profile is the radial
+    # profile projected onto a plane `radius` away from the centre of the
+    # cylinder. The y coordinates on the edge must drop to compensate.
+
+    # The drop is equal to, via similar triangles
+    drop = inner_raise * (radius_proj - radius) / (radius - radius_inner)
     outer = [
-        (radius * math.tan(theta), 0),
-        (radius * math.tan(theta) - tol, 0),
+        (span, -tol - drop),
+        (span, -drop),
         (0, tooth_height),
-        (-radius * math.tan(theta) + tol, 0),
-        (-radius * math.tan(theta), 0),
+        (-span, -drop),
+        (-span, -tol - drop),
     ]
+    adj = radius_inner * c
+    # In the case of the inner triangle, it is projected onto a plane `adj` away
+    # from the centre. The apex must extrapolate
+
+    # Via similar triangles
+    #
+    #  (inner_raise + tooth_height_inner) -
+    #    (tooth_height - inner_raise - tooth_height_inner) * ((radius_inner - adj) / (radius - radius_inner))
+    apex = (inner_raise + tooth_height_inner)  - \
+        inner_raise * (radius_inner - adj) / (radius - radius_inner)
     inner = [
-        (radius_inner * math.sin(theta), 0),
-        (radius_inner * math.sin(theta), inner_raise),
-        (0, inner_raise + tooth_height_inner),
-        (-radius_inner * math.sin(theta), inner_raise),
-        (-radius_inner * math.sin(theta), 0),
+        (span_inner, -tol),
+        (span_inner, inner_raise),
+        (0, apex),
+        (-span_inner, inner_raise),
+        (-span_inner, -tol),
     ]
     tooth = (
         Cq.Workplane('YZ')
         .polyline(inner)
         .close()
-        .workplane(offset=radius - radius_inner)
+        .workplane(offset=radius - adj)
         .polyline(outer)
         .close()
-        .loft(ruled=True, combine=True)
+        .loft(ruled=False, combine=True)
         .val()
     )
-    tooth_centre_radius = radius_inner * math.cos(theta)
     angle_offset = hirth_tooth_angle(n_tooth) / 2 if is_mated else 0
     teeth = (
         Cq.Workplane('XY')
         .polarArray(
-            radius=tooth_centre_radius,
+            radius=adj,
             startAngle=angle_offset,
             angle=360,
             count=n_tooth)
@@ -75,6 +94,14 @@ def hirth_joint(radius=60,
             height=base_height + tooth_height,
             radius=radius,
         ))
+        .intersect(Cq.Solid.makeCylinder(
+            height=base_height + tooth_height,
+            radius=radius,
+        ))
+        .cut(Cq.Solid.makeCylinder(
+            height=base_height + tooth_height,
+            radius=radius_inner,
+        ))
     )
     base = (
         Cq.Workplane('XY')
diff --git a/nhf/test.py b/nhf/test.py
index e29b781..b0ab29a 100644
--- a/nhf/test.py
+++ b/nhf/test.py
@@ -3,6 +3,7 @@ import cadquery as Cq
 import nhf.joints
 import nhf.handle
 import nhf.metric_threads as NMt
+from nhf.checks import binary_intersection
 
 class TestJoints(unittest.TestCase):
 
@@ -12,7 +13,10 @@ class TestJoints(unittest.TestCase):
             j.val().solids(), Cq.Solid,
             msg="Hirth joint must be in one piece")
     def test_joints_hirth_assembly(self):
-        nhf.joints.hirth_assembly()
+        assembly = nhf.joints.hirth_assembly()
+        isect = binary_intersection(assembly)
+        self.assertLess(isect.Volume(), 1e-6,
+                        "Hirth joint assembly must not have intersection")
     def test_joints_comma_assembly(self):
         nhf.joints.comma_assembly()
     def test_torsion_joint(self):