fix: Hirth joint mating

nhf/joints.py

@@ -4,126 +4,86 @@ import cadquery as Cq
 import nhf.springs as NS
 from nhf import Role
 
-def hirth_tooth_angle(n_tooth):
-    """
-    Angle of one whole tooth
-    """
-    return 360 / n_tooth
-
 @dataclass(frozen=True)
 class HirthJoint:
     """
     A Hirth joint attached to a cylindrical base
     """
+
+    # r
     radius: float = 60
+    # r_i
     radius_inner: float = 40
     base_height: float = 20
     n_tooth: float = 16
+    # h_o
     tooth_height: float = 16
-    tooth_height_inner: float = 2
 
     def __post_init__(self):
         # Ensures tangent doesn't blow up
         assert self.n_tooth >= 5
         assert self.radius > self.radius_inner
-        assert self.tooth_height >= self.tooth_height_inner
-
-    @property
-    def _theta(self):
-        return math.pi / self.n_tooth
 
     @property
     def tooth_angle(self):
-        return hirth_tooth_angle(self.n_tooth)
+        return 360 / self.n_tooth
 
 
     def generate(self, tag_prefix="", is_mated=False, tol=0.01):
         """
         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
+        FIXME: Mate is not exact when number of tooth is low
-        this issue.
         """
-        c, s, t = math.cos(self._theta), math.sin(self._theta), math.tan(self._theta)
+        phi = math.radians(self.tooth_angle)
-        span = self.radius * t
+        alpha = 2 * math.atan(self.radius / self.tooth_height * math.tan(phi/2))
-        radius_proj = self.radius / c
+        #alpha = math.atan(self.radius * math.radians(180 / self.n_tooth) / self.tooth_height)
-        span_inner = self.radius_inner * s
+        gamma = math.radians(90 / self.n_tooth)
-        # 2 * raise + (inner tooth height) = (tooth height)
+        # Tooth half height
-        inner_raise = (self.tooth_height - self.tooth_height_inner) / 2
+        l = self.radius * math.cos(gamma)
-        # Outer tooth triangle spans 2*theta radians. This profile is the radial
+        a = self.radius * math.sin(gamma)
-        # profile projected onto a plane `radius` away from the centre of the
+        t = a / math.tan(alpha / 2)
-        # cylinder. The y coordinates on the edge must drop to compensate.
+        beta = math.asin(t / l)
-
+        dx = self.tooth_height * math.tan(alpha / 2)
-        # The drop is equal to, via similar triangles
+        profile = (
-        drop = inner_raise * (radius_proj - self.radius) / (self.radius - self.radius_inner)
-        outer = [
-            (span, -tol - drop),
-            (span, -drop),
-            (0, self.tooth_height),
-            (-span, -drop),
-            (-span, -tol - drop),
-        ]
-        adj = self.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 + self.tooth_height_inner)  - \
-            inner_raise * (self.radius_inner - adj) / (self.radius - self.radius_inner)
-        inner = [
-            (span_inner, -tol),
-            (span_inner, inner_raise),
-            (0, apex),
-            (-span_inner, inner_raise),
-            (-span_inner, -tol),
-        ]
-        tooth = (
             Cq.Workplane('YZ')
-            .polyline(inner)
+            .polyline([
+                (0, 0),
+                (dx, self.tooth_height),
+                (-dx, self.tooth_height),
+            ])
             .close()
-            .workplane(offset=self.radius - adj)
+            .extrude(-self.radius)
-            .polyline(outer)
-            .close()
-            .loft(ruled=False, combine=True)
             .val()
+            .rotate((0, 0, 0), (0, 1, 0), math.degrees(beta))
+            .moved(Cq.Location((0, 0, self.base_height)))
         )
-        angle_offset = hirth_tooth_angle(self.n_tooth) / 2 if is_mated else 0
+        core = Cq.Solid.makeCylinder(
-        h = self.base_height + self.tooth_height
+            radius=self.radius_inner,
-        teeth = (
+            height=self.tooth_height,
+            pnt=(0, 0, self.base_height),
+        )
+        angle_offset = self.tooth_angle / 2 if is_mated else 0
+        result = (
             Cq.Workplane('XY')
+            .cylinder(
+                radius=self.radius,
+                height=self.base_height + self.tooth_height,
+                centered=(True, True, False))
+            .faces(">Z")
+            .tag(f"{tag_prefix}bore")
+            .cut(core)
             .polarArray(
-                radius=adj,
+                radius=self.radius,
                 startAngle=angle_offset,
                 angle=360,
                 count=self.n_tooth)
-            .eachpoint(lambda loc: tooth.located(loc))
+            .cutEach(
-            .intersect(Cq.Solid.makeCylinder(
+                lambda loc: profile.moved(loc),
-                height=h,
+            )
-                radius=self.radius,
-            ))
-            .cut(Cq.Solid.makeCylinder(
-                height=h,
-                radius=self.radius_inner,
-            ))
         )
-        base = (
-            Cq.Workplane('XY')
-            .cylinder(
-                height=self.base_height,
-                radius=self.radius,
-                centered=(True, True, False))
-            .faces(">Z").tag(f"{tag_prefix}bore")
-            .union(
-                teeth.val().move(Cq.Location((0,0,self.base_height))),
-                tol=tol)
-            .clean()
-        )
-        #base.workplane(offset=tooth_height/2).circle(radius=radius,forConstruction=True).tag("mate")
         (
-            base
+            result
             .polyline([
                 (0, 0, self.base_height),
                 (0, 0, self.base_height + self.tooth_height)
@@ -131,11 +91,11 @@ class HirthJoint:
             .tag(f"{tag_prefix}mate")
         )
         (
-            base
+            result
             .polyline([(0, 0, 0), (1, 0, 0)], forConstruction=True)
             .tag(f"{tag_prefix}directrix")
         )
-        return base
+        return result
 
     def assembly(self):
         """

nhf/test.py

@@ -13,12 +13,15 @@ class TestJoints(unittest.TestCase):
         self.assertIsInstance(
             obj.val().solids(), Cq.Solid,
             msg="Hirth joint must be in one piece")
+
     def test_joints_hirth_assembly(self):
-        j = nhf.joints.HirthJoint()
+        for n_tooth in [16, 20, 24]:
-        assembly = j.assembly()
+            with self.subTest(n_tooth=n_tooth):
-        isect = binary_intersection(assembly)
+                j = nhf.joints.HirthJoint()
-        self.assertLess(isect.Volume(), 1e-6,
+                assembly = j.assembly()
-                        "Hirth joint assembly must not have intersection")
+                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):